function [imp,sim] = iscmemsw_ms(param,shocks);
%
%ISCMEMSW_MS  computes equilibria for the International-Sticky business
%             cycle model with Complete Markets, Exogenous Money growth,
%             and Sticky Wages in Chari, Kehoe, and McGrattan, ``Can 
%             Sticky Price Models Generate Volatile and Persistent Real 
%             Exchange Rates?''  This version includes government spending 
%             and technology shocks.
%
%             The optimization problems being solved are as follows.
%             Home CONSUMER j solves:
%
%               max       sum  beta^t pi(s^t) U[C(j,s^t),L(j,s^t),M(j,s^t)]
%              C(j,s^t)  t,s^t                                    --------
%              M(j,s^t)                                            P(s^t)
%              B(j,s^t)
%              W(j,s^{t-1})
%
%             subject  to
%                 P(s^t)*C(j,s^t) + M(j,s^t) 
%                     + sum_{s_{t+1}} Q(s^{t+1}|s^t)*B(j,s^{t+1})
%                   <=  P(s^t)*W(j,s^{t-1})*L(j,s^t)) + M(j,s^{t-1}) 
%                          + B(j,s^t) + Profits(s^t) - Tax(s^t)
%
%                 L(j,s^t)   = int l(i,j,s^t) di
%                 l(i,j,s^t) = (Wbar(s^t)/W(j,s^{t-1}))^(1/(1-v)) L^d(i,s^t)
%
%
%             Home FINAL GOODS PRODUCERS solve:
%
%               max     P*[ome1*(int Yh(i)^theta di)^(rho/theta)
%               Yh(i)     +ome2*(int Yf(i)^theta di)^(rho/theta)]^(1/rho)
%               Yf(i)         -int Ph(i) Yh(i) di -int Pf(i) Yf(i) di
%
%
%             Home INTERMEDIATE GOODS PRODUCERS solve:
%
%               max     sum  Q(s^t|.)*[Ph(i,s^t)*Yh(i,s^t)
%               Ph     t,s^t         + Ph^*(i,s^t)*Yh^*(i,s^t)*e(s^t)
%               Ph^*                 - P(s^t) int W(j,s^t)l(i,j,s^t) dj
%                                    - P(s^t)X(i,s^t)]
%
%                       s.t. Yh(i,s^t)+Yh^*(i,s^t) 
%                                    = F(K(i,s^{t-1}),A(s^t) L^d(i,s^t))
%
%                            L^d(i,s^t) 
%                                    = [ int l(i,j,s^t)^v dj ] ^(1/v)
%
%                            Yh(i)   = Y [ome1 P]^(1/(1-rho)) 
%                                      Phb(i)^[(rho-theta)/(1-rho)/(theta-1)]
%                                     /Ph(i)^(1/(1-theta))
%
%                            Yh^*(i) = Y^* [ome2 P^*]^(1/(1-rho)) 
%                                      Phb^*(i)^[(rho-theta)/(1-rho)/(theta-1)]
%                                     /Ph^*(i)^(1/(1-theta))
%
%                            K(i,s^t)  = (1-delta) K(i,s^{t-1}) + X(i,s^t)
%                                         - phi(X(i,s^t)/K(i,s^{t-1}))
%                                             * K(i,s^{t-1})
%
%                            Ph(i,s^tau) = Ph(i,s^{t-1}), tau=t,...,t+N-1  
%                            Ph(i,s^tau) = Ph(i,s^{t+N-1}), tau=t+N,...,t+2N-1  
%                                        :
%
%             Foreign consumers and firms solve problems analogous to
%             those given above.
%
%

%             Ellen McGrattan, 3-21-00
%             Revised, ERM, 7-25-01

%______________________________________________________________________________
%
%PARAMETERS:
%
%  Staggered contracts:
%
%    T = fraction of year prices are fixed (T=1 ==> fixed for 1 year)
%    N = number of cohorts
%
%  Preferences:
%
%    Utype=1:
%      U(C,M,L) = [CES(C,M) * (1-L)^psi]^(1-sigma)/(1-sigma)
%      CES(C,M) = [omega*C^nu+(1-omega)*M^nu]^(1/nu), nu=(eta-1)/eta
%
%    Utype=2:
%      U(C,M,L) = CES(C,M)^(1-sigma)/(1-sigma) +psi *(1-L)^(1-xi)/(1-xi)
%
%    Utype=3:
%      U(C,M,L) = [C^(1-kappa)/(1-kappa) -psi * L^(1+xi)/(1+xi)
%                      +omega * M^nu/nu]^(1-sigma)/(1-sigma), nu=(eta-1)/eta
%    Utype=4:
%      U(C,M,L) = [C^(1-kappa)/(1-kappa) +psi * (1-L)^(1-xi)/(1-xi)
%                      +omega * M^nu/nu]^(1-sigma)/(1-sigma), nu=(eta-1)/eta
%
%  Discounting:
%
%    beta     = the discount factor
%
%  Technology:
%
%    F(K,AL)  = K^alpha * (AL)^(1-alpha)
%    Y        = [ome1*(int Yh(i)^theta di)^(rho/theta)
%               +ome2*(int Yf(i)^theta di)^(rho/theta)]^(1/rho)
%
%  Capital accumulation
%
%    K[i,t]   = (1-delta) K[i,t-1] + X[i,t] - phi(X[i,t]/K[i,t-1]) K[i,t-1]
%    phi(X/K) = b/2 * (X/K-delta)^2
%
%  Process for money:
%
%    log mu[t+1] = (1-rhom)*mu + rhom*log mu[t] + eps[t+1],   eps ~ N(0,sigm)
%
%  Timing of shocks:
%
%    if timing = 1, monopolists observe only lagged shocks when choosing prices
%              = 2, monopolists observe lagged shocks and a[t]
%              = 3, monopolists observe lagged shocks and g[t],a[t]
%  
%______________________________________________________________________________
%

T        = param(1);
N        = param(2);
Nw       = param(3);
Utype    = param(4);
ome1     = param(12);
rho      = param(13);
theta    = param(14);
alpha    = param(15);
delta    = param(16);
b        = param(17);
mu       = param(19);
rhom     = param(20);
sigm     = param(21);
corm     = param(22);
gs       = param(36);
rhog     = param(37);
sigg     = param(38);
corg     = param(39);
rhoa     = param(40);
siga     = param(41);
cora     = param(42);
timing   = param(43);
v        = param(46);
nti      = param(47);
nts      = param(48);
nsim     = param(49);
mimp     = 1;
gimp     = 0;
aimp     = 0;

ome2     = 1-ome1;
alph1    = 1-alpha;
thet1    = 1-theta;
v1       = 1-v;
sup1     = 1/(1-rho);
sup2     = rho/(rho-1);
sup3     = (rho-theta)/(1-rho)/(theta-1);

%______________________________________________________________________________



%______________________________________________________________________________
%
% STEADY STATE FOR THE ECONOMY
%

ss       = steadyw_ms(param);
param    = [param;ss];

%______________________________________________________________________________



%______________________________________________________________________________
%
% LINEARIZED FIRST ORDER CONDITIONS
%

% Let mN = max(N,Nw).
% The ordering of variables to the residuals is
%
%     ZZ  = [z_{t+mN-1},z_{t+mN-2},...z_t,
%            ph_{t-2}-m_{t-2},...,ph_{t-N}-m_{t-N},
%            pf_{t-2}-m_{t-2},...,pf_{t-N}-m_{t-N},
%            pf*_{t-2}-m*_{t-2},...,pf*_{t-N}-m*_{t-N},
%            ph*_{t-2}-m*_{t-2},...,ph*_{t-N}-m*_{t-N},
%            k_{1t-1},...,k_{Nt-1},k*_{1t-1},...,k*_{Nt-1},
%            w_{t-2}-m_{t-2},...,w_{t-Nw}-m_{t-Nw},
%            w*_{t-2}-m_{t-2},...,w*_{t-Nw}-m_{t-Nw},
%            mu_{t+mN-1},...,mu_{t-mN+1},mu*_{t+mN-1},...,mu*_{t-mN+1},
%            g_{t+mN-1},...,g_t,g*_{t+mN-1},...,g*_t,
%            a_{t+mN-1},...,a_t,a*_{t+mN-1},...,a*_t]'
%
%     z_t = [ph_{t-1}-m_{t-1}                NOTE: variables are logged 
%            pf_{t-1}-m_{t-1}                      and demeaned 
%            pf*_{t-1}-m*_{t-1}
%            ph*_{t-1}-m*_{t-1}
%              k_{1t}
%                :
%              k_{Nt}
%              k*_{1t}
%                :
%              k*_{Nt}
%              y_t   
%              y*_t 
%              w_{t-1}-m_{t-1}
%              w*_{t-1}-m*_{t-1}]      (2*N+8 x 1)
%


phs      = ss(1);
ys       = ss(2);
ws       = ss(3);
ks       = ss(12*Nw+4:12*Nw+3+N);
z        = [phs;phs;phs;phs;ks;ks;ys;ys;ws;ws];
nz       = 2*N+8;

mN       = max(N,Nw);
tem      = z*ones(1,mN);
ZZ       = log([tem(:);phs*ones(4*N-4,1);ks;ks;ws*ones(2*Nw-2,1);
                mu*ones(4*mN-2,1);gs*ones(2*mN,1);ones(2*mN,1)]);
nZZ      = 2*N*mN+6*N+2*Nw+16*mN-8;
del      = max(abs(ZZ)*1e-4,1e-8);
for i=1:nZZ;
  ZZp     = ZZ;
  ZZm     = ZZ;
  ZZp(i)  = ZZp(i)+del(i);
  ZZm(i)  = ZZm(i)-del(i);
  dR(:,i) = (res_iscmemsw_ms(ZZp,param)-res_iscmemsw_ms(ZZm,param))/(2*del(i));
end;
%______________________________________________________________________________


%______________________________________________________________________________
%
% SOLUTION OF DETERMINISTIC PART OF PROBLEM (A)
%

%
% The residuals are given by
%
%     A1*[X_{t+1},Z_{t+mN-1}]' + A2*[X_t,Z_{t+mN-2}]' + shock terms = 0
%
% where
%
%     Z_t = [z_t,z_{t-1}, ....z_{t-N+2}]'    (nz*(mN-1) x 1)
%
%     z_t = [ph_{t-1}-m_{t-1}
%            pf_{t-1}-m_{t-1}
%            pf*_{t-1}-m*_{t-1}
%            ph*_{t-1}-m*_{t-1}
%              k_{1t}
%                :
%              k_{Nt}
%              k*_{1t}
%                :
%              k*_{Nt}
%              y_t   
%              y*_t 
%              w_{t-1}-m_{t-1}
%              w*_{t-1}-m*_{t-1}]      (2*N+8 x 1)
%
%     X_t = [ph_{t-2}-m_{t-2}
%                :
%            ph_{t-N}-m_{t-N}
%            pf_{t-2}-m_{t-2}
%                :
%            pf_{t-N}-m_{t-N}
%            pf*_{t-2}-m*_{t-2}
%                :
%            pf*_{t-N}-m*_{t-N}
%            ph*_{t-2}-m*_{t-2}
%                :
%            ph*_{t-N}-m*_{t-N}
%              k_{1,t-1}
%                :
%              k_{N,t-1}
%              k*_{1,t-1}
%                :
%              k*_{N,t-1}
%              w_{t-2}-m_{t-2}
%                :
%              w_{t-Nw}-m_{t-Nw}
%              w*_{t-2}-m*_{t-2}
%                :
%              w*_{t-Nw}-m*_{t-Nw}]      (6*N+2*Nw-6 x 1)
%
% The solution to the deterministic part of the problem is given by
%
%     Z_t = A Z_{t-1}
%
%

nX       = 6*N+2*Nw-6;
nZ       = (mN-1)*nz;
nA1      = nX+nZ;
A1       = [eye(nX,nX),zeros(nX,(mN-1)*nz);
            zeros(nz,nX),dR(:,1:(mN-1)*nz);
            zeros((mN-2)*nz,nX+nz),eye((mN-2)*nz)];
I1       = [kron(eye(4),[zeros(1,N-1);eye(N-2),zeros(N-2,1)]), ...
              zeros(4*N-4,2*N+2*Nw-2);
            zeros(2*N+2*Nw-2,nX)];
I2       = [kron(eye(4),[1;zeros(N-2,1)]),zeros(4*N-4,nz-4);
            zeros(2*N,4),eye(2*N),zeros(2*N,4);
            zeros(2*Nw-2,2*N+6),kron(eye(2),[1;zeros(Nw-2,1)])];

A2       = [-I1,zeros(nX,(mN-2)*nz),-I2;
            dR(:,mN*nz+1:mN*nz+nX),zeros(nz,(mN-2)*nz), ...
            dR(:,(mN-1)*nz+1:mN*nz);
            zeros((mN-2)*nz,nX),-eye((mN-2)*nz),zeros((mN-2)*nz,nz)];

[V,L]    = eig(A2,-A1);
L        = diag(L);
[D,I]    = sort(abs(L));
if (~( (D(nX) <= 1) & (D(nX+1) > 1) ) )
  disp('Can''t order eigenvalues')
end;
V11      = V(1:nX,I(1:nX));
V21      = V(nX+1:nA1,I(1:nX));
W        = real([V11*diag(L(I(1:nX)))/V11;
                 V21/V11]);

A                         = zeros((mN-1)*nz);
A(1:nz,1:nz-4)            = W(nA1-nz+1:nA1,[1,N,2*N-1,3*N-2,4*N-3:6*N-4]);
A(1:nz,nz+1:nz:(N-1)*nz)  = W(nA1-nz+1:nA1,2:N-1);
A(1:nz,nz+2:nz:(N-1)*nz)  = W(nA1-nz+1:nA1,N+1:2*N-2);
A(1:nz,nz+3:nz:(N-1)*nz)  = W(nA1-nz+1:nA1,2*N:3*N-3);
A(1:nz,nz+4:nz:(N-1)*nz)  = W(nA1-nz+1:nA1,3*N-1:4*N-4);
A(1:nz,nz-1:nz:(Nw-1)*nz) = W(nA1-nz+1:nA1,6*N-3:6*N+Nw-5);
A(1:nz,nz:nz:(Nw-1)*nz)   = W(nA1-nz+1:nA1,6*N+Nw-4:6*N+2*Nw-6);
A(nz+1:nZ,1:nZ-nz)        = eye(nZ-nz);

%______________________________________________________________________________


%______________________________________________________________________________
%
% SOLUTION OF COEFFICIENTS ON MONETARY GROWTH (B)
%

% The residuals can be written as:
%
%   a_0*Z_{t+mN-1} + a_1*Z_{t+mN-2} .... + a_{mN-1}*Z_t + a_mN*Z_{t-1}
% + b_0*S_{t+mN-1} + b_1*S_{t+mN-2} .... + b_{mN-1}*S_t               = 0
%
% where
%
%     S_t = [mu_t
%             :
%            mu_{t-mN+1}
%            mu*_t
%             :
%            mu*_{t-mN+1}  
%            g_t
%            g_{t-1}
%            g*_t
%            g*_{t-1}
%            a_t
%            a_{t-1}
%            a*_t
%            a*_{t-1}]     (2*mN+8 x 1)
%

nS             = 2*mN+8;
a0             = zeros(nz*nZ,mN+1);
b0             = zeros(nz*nS,mN);

a0(:,1)        = vec(dR(:,1:nZ));
a0(:,mN)       = vec([dR(:,nZ+1:mN*nz),zeros(nz,nZ-nz)]);
b0(:,1)        = vec([dR(:,mN*nz+nX+[1:mN,2*mN:3*mN-1]), ...
                      dR(:,mN*nz+nX+4*mN-1),zeros(nz,1), ...
                      dR(:,mN*nz+nX+5*mN-1),zeros(nz,1), ...
                      dR(:,mN*nz+nX+6*mN-1),zeros(nz,1), ...
                      dR(:,mN*nz+nX+7*mN-1),zeros(nz,1)]);

for i=2:mN-1;
  b0(:,i)      = vec([zeros(nz,2*mN), ...
                      dR(:,mN*nz+nX+4*mN-2+i),zeros(nz,1), ...
                      dR(:,mN*nz+nX+5*mN-2+i),zeros(nz,1), ...
                      dR(:,mN*nz+nX+6*mN-2+i),zeros(nz,1), ...
                      dR(:,mN*nz+nX+7*mN-2+i),zeros(nz,1)]);
end;
b0(:,mN)       = vec([zeros(nz,1),dR(:,mN*nz+nX+[mN+1:2*mN-1]), ...
                      zeros(nz,1),dR(:,mN*nz+nX+[3*mN:4*mN-2]), ...
                      dR(:,mN*nz+nX+5*mN-2),zeros(nz,1), ...
                      dR(:,mN*nz+nX+6*mN-2),zeros(nz,1), ...
                      dR(:,mN*nz+nX+7*mN-2),zeros(nz,1), ...
                      dR(:,mN*nz+nX+8*mN-2),zeros(nz,1)]);


sb             = [eye(nz);zeros(nZ-nz,nz)];
sa             = dR(:,1:nZ);
P              = [kron(eye(2),[rhom,zeros(1,mN-1);
                               eye(mN-1),zeros(mN-1,1)]), zeros(2*mN,8);
                  zeros(8,2*mN),[rhog,0,0,0,0,0,0,0;
                                 1,0,0,0,0,0,0,0;
                                 0,0,rhog,0,0,0,0,0;
                                 0,0,1,0,0,0,0,0;
                                 0,0,0,0,rhoa,0,0,0;
                                 0,0,0,0,1,0,0,0;
                                 0,0,0,0,0,0,rhoa,0;
                                 0,0,0,0,0,0,1,0]];

sum1           = 0;
sum2           = 0;
aa             = zeros(nz,nZ);
bb             = zeros(nz,nS);
for i=1:mN;
  sum1         = sum1 + kron(sa*sb,(P^(mN-i))');
  aa(:)        = a0(:,i+1);
  sa           = sa*A + aa;
  bb(:)        = b0(:,i);
  sum2         = sum2 + (P^(mN-i))'*bb';
end;

Psi1           = [[rhom;zeros(mN-2,1)],eye(mN-1)];
nzero          = ones(nz,nS);
if timing == 2;
  Psi2         = [kron(eye(2),Psi1),zeros(2*mN-2,8);
                  zeros(6,2*mN),[rhog,1,0,0,0,0,0,0;
                                 0,0,rhog,1,0,0,0,0;
                                 0,0,0,0,1,0,0,0;
                                 0,0,0,0,0,1,0,0;
                                 0,0,0,0,0,0,1,0;
                                 0,0,0,0,0,0,0,1]];
  dd           = [kron(eye(4),Psi2), zeros(4*(nS-4),nz*nS-4*nS);
                  zeros(nz*nS-6*nS,4*nS),eye(nz*nS-6*nS),zeros(nz*nS-6*nS,2*nS);
                  zeros(2*(nS-4),nz*nS-2*nS),kron(eye(2),Psi2)];

  nzero([1:4,2*N+7,2*N+8],[1,mN+1,2*mN+1,2*mN+3]) ...
               = zeros(6,4);

elseif timing == 3;
  Psi2         = [kron(eye(2),Psi1),zeros(2*mN-2,8);
                  zeros(8,2*mN),eye(8)];
  dd           = [kron(eye(4),Psi2), zeros(4*(nS-2),nz*nS-4*nS);
                  zeros(nz*nS-6*nS,4*nS),eye(nz*nS-6*nS),zeros(nz*nS-6*nS,2*nS);
                  zeros(2*(nS-2),nz*nS-2*nS),kron(eye(2),Psi2)];

  nzero([1:4,2*N+7,2*N+8],[1,mN+1]) ...
               = zeros(6,2);
else;
  Psi2         = [kron(eye(2),Psi1),zeros(2*mN-2,8);
                  zeros(4,2*mN),[rhog,1,0,0,0,0,0,0;
                                 0,0,rhog,1,0,0,0,0;
                                 0,0,0,0,rhoa,1,0,0;
                                 0,0,0,0,0,0,rhoa,1]];
  dd           = [kron(eye(4),Psi2), zeros(4*(nS-6),nz*nS-4*nS);
                  zeros(nz*nS-6*nS,4*nS),eye(nz*nS-6*nS),zeros(nz*nS-6*nS,2*nS);
                  zeros(2*(nS-6),nz*nS-2*nS),kron(eye(2),Psi2)];

  nzero([1:4,2*N+7,2*N+8],[1,mN+1,2*mN+1,2*mN+3,2*mN+5,2*mN+7]) ...
               = zeros(6);

%
%NOTE: using the following assumes that wages are sticky but not preset
%
% dd           = [kron(eye(4),Psi2), zeros(4*(nS-6),nz*nS-4*nS)
%                 zeros(nz*nS-4*nS,4*nS),eye(nz*nS-4*nS)];
% nzero([1:4],[1,mN+1,2*mN+1,2*mN+3,2*mN+5,2*mN+7]) ...
%              = zeros(4,6);

end;
nzero          = find(nzero')';
qq             = dd*sum1;
rr             = dd*vec(sum2);
B1             = zeros(nS,nz);
B1(nzero)      = -qq(:,nzero)\rr;
B              = sb*B1';
%
% FINAL SOLUTION:  Z[t] = A* Z[t-1] + B* S[t]
%
%______________________________________________________________________________


%______________________________________________________________________________
%
% IMPULSE RESPONSES FOR SYSTEM
%
%       [ Z[t]   ] = [ A  B Means]  [ Z[t-1] ]      [0    ]
%       [ S[t+1] ]   [ 0  P Means]  [ S[t]   ]   +  [Sigma] * e[t+1]
%       [   1    ]   [ 0  0   1  ]  [   1    ]      [0    ]
%
%                      (G)                           (H)

tem       = z*ones(1,mN-1);
Xbar      = log([tem(:);mu*ones(2*mN,1);gs;gs;gs;gs;1;1;1;1]);
G         = [A, B;
             zeros(nS,nZ), P];
G         = [G, (eye(nZ+nS)-G)*Xbar;
             zeros(1,nZ+nS),1];
H         = [zeros(nZ,6);
             sigm,corm*sigm/sqrt(1-corm*corm),0,0,0,0;
             zeros(mN-1,6);
             0,sigm/sqrt(1-corm*corm),0,0,0,0;
             zeros(mN-1,6);
             0,0,sigg,corg*sigg/sqrt(1-corg*corg),0,0;
             0,0,0,0,0,0;
             0,0,0,sigg/sqrt(1-corg*corg),0,0;
             0,0,0,0,0,0;
             0,0,0,0,siga,cora*siga/sqrt(1-cora*cora);
             0,0,0,0,0,0;
             0,0,0,0,0,siga/sqrt(1-cora*cora);
             0,0,0,0,0,0;
             0,0,0,0,0,0];
cjs       = ss(4:3+Nw);
mdjs      = ss(4+Nw:3+2*Nw);
Lsjs      = ss(4+2*Nw:3+3*Nw);
UCC       = ss(4+6*Nw:3+7*Nw);
UMM       = ss(4+7*Nw:3+8*Nw);
UCM       = ss(4+9*Nw:3+10*Nw);
UCL       = ss(4+10*Nw:3+11*Nw);
ULM       = ss(4+11*Nw:3+12*Nw);
xis       = ss(4+12*Nw+N:3+12*Nw+2*N);
Aa        = [cjs'/Nw,zeros(1,Nw);
             -UCC(1)*cjs(1)*ones(Nw-1,1),diag(UCC(2:Nw).*cjs(2:Nw)), ...
             -UCM(1)*mdjs(1)*ones(Nw-1,1),diag(UCM(2:Nw).*mdjs(2:Nw));
             zeros(1,Nw),ones(1,Nw)/Nw;
             -UCM(1)*cjs(1)*ones(Nw-1,1),diag(UCM(2:Nw).*cjs(2:Nw)), ...
             -UMM(1)*mdjs(1)*ones(Nw-1,1),diag(UMM(2:Nw).*mdjs(2:Nw))];

if nti > 0;
  
  nimp     = mN-1;
  nt       = nti+nimp;
  inputs   = [zeros(nimp-1,6);mimp,0,gimp,0,aimp,0;zeros(nti,6)];
  [Xt,Xtl] = dlsim(G, H, G, H, inputs, [Xbar;1]);
  pht      = exp(Xt(:,1));
  pft      = exp(Xt(:,2));
  pfat     = exp(Xt(:,3));
  phat     = exp(Xt(:,4));
  Kit      = exp(Xt(:,5:N+4));
  Kait     = exp(Xt(:,N+5:2*N+4));
  Kitl     = exp(Xtl(:,5:N+4));
  Kaitl    = exp(Xtl(:,N+5:2*N+4));
  Yt       = exp(Xt(:,2*N+5));
  Yat      = exp(Xt(:,2*N+6));
  wt       = exp(Xt(:,2*N+7));
  wat      = exp(Xt(:,2*N+8));
  Mut      = exp(Xtl(:,nZ+1)); 
  Muat     = exp(Xtl(:,nZ+mN+1));
  Gt       = exp(Xtl(:,nZ+2*mN+1));
  Gat      = exp(Xtl(:,nZ+2*mN+3));
  At       = exp(Xtl(:,nZ+2*mN+5));
  Aat      = exp(Xtl(:,nZ+2*mN+7));
  Mt       = cumprod(Mut);
  Mat      = cumprod(Muat);
  Pht      = pht.*Mt./Mut;
  Pft      = pft.*Mt./Mut;
  Pfat     = pfat.*Mat./Muat;
  Phat     = phat.*Mat./Muat;
  Wt       = wt.*Mt./Mut;
  Wat      = wat.*Mat./Mut;
  phvec    = [phs*ones(1,N-1)./mu.^[N-1:-1:1],Pht']';
  pfvec    = [phs*ones(1,N-1)./mu.^[N-1:-1:1],Pft']';
  pfavec   = [phs*ones(1,N-1)./mu.^[N-1:-1:1],Pfat']';
  phavec   = [phs*ones(1,N-1)./mu.^[N-1:-1:1],Phat']';
  wvec     = [ws*ones(1,Nw-1)./mu.^[Nw-1:-1:1],Wt']';
  wavec    = [ws*ones(1,Nw-1)./mu.^[Nw-1:-1:1],Wat']';
  Phit     = zeros(nt,N);
  Pfit     = zeros(nt,N);
  Pfait    = zeros(nt,N);
  Phait    = zeros(nt,N);
  Wit      = zeros(nt,Nw);
  Wait     = zeros(nt,Nw);
  for i=1:N;
    Phit(:,i)  = phvec(N-i+1:nt-i+N);
    Pfit(:,i)  = pfvec(N-i+1:nt-i+N);
    Pfait(:,i) = pfavec(N-i+1:nt-i+N);
    Phait(:,i) = phavec(N-i+1:nt-i+N);
  end;
  for i=1:Nw;
    Wit(:,i)   = wvec(Nw-i+1:nt-i+Nw);
    Wait(:,i)  = wavec(Nw-i+1:nt-i+Nw);
  end; 
  Phbt     = (mean((Phit').^(-theta/thet1))').^(-thet1/theta);    
  Pfbt     = (mean((Pfit').^(-theta/thet1))').^(-thet1/theta);    
  Pfabt    = (mean((Pfait').^(-theta/thet1))').^(-thet1/theta);    
  Phabt    = (mean((Phait').^(-theta/thet1))').^(-thet1/theta);    
  Wbt      = (mean((Wit').^(-v/v1))').^(-v1/v);
  Wabt     = (mean((Wait').^(-v/v1))').^(-v1/v);
  Pt       = (ome1^sup1*Phbt.^sup2+ome2^sup1*Pfbt.^sup2).^(1/sup2);
  Pat      = (ome1^sup1*Pfabt.^sup2+ome2^sup1*Phabt.^sup2).^(1/sup2);
  Yhit     = (((ome1*Pt).^sup1.*Phbt.^sup3.*Yt)*ones(1,N))./Phit.^(1/thet1);
  Yfit     = (((ome2*Pt).^sup1.*Pfbt.^sup3.*Yt)*ones(1,N))./Pfit.^(1/thet1);
  Yfait    = (((ome1*Pat).^sup1.*Pfabt.^sup3.*Yat)*ones(1,N))./Pfait.^(1/thet1);
  Yhait    = (((ome2*Pat).^sup1.*Phabt.^sup3.*Yat)*ones(1,N))./Phait.^(1/thet1);
  Ldit     = ((Yhit+Yhait)./(Kitl(:,[N,1:N-1]).^alpha)).^(1/alph1);
  Ldait    = ((Yfait+Yfit)./(Kaitl(:,[N,1:N-1]).^alpha)).^(1/alph1);
  Ldit     = Ldit./(At*ones(1,N));
  Ldait    = Ldait./(Aat*ones(1,N));
  Ldt      = mean(Ldit')';            
  Ldat     = mean(Ldait')';          
  KK       = Kit./Kitl(:,[N,1:N-1]);
  KKa      = Kait./Kaitl(:,[N,1:N-1]);
  Iit      = (1+b*delta-sqrt(1+2*b-2*b*KK))/b .* Kitl(:,[N,1:N-1]);     
  Iait     = (1+b*delta-sqrt(1+2*b-2*b*KKa))/b.* Kaitl(:,[N,1:N-1]);     
  It       = mean(Iit')';          
  Iat      = mean(Iait')';         
  Lsjt     = ( (Ldt.*Wbt.^(1/v1))*ones(1,Nw) )./ Wit.^(1/v1);
  Lsajt    = ( (Ldat.*Wabt.^(1/v1))*ones(1,Nw) )./ Wait.^(1/v1);
  Lt       = mean(Lsjt')';
  Lat      = mean(Lsajt')';

  for i=1:nt;
    tem       = Aa\[ys*log(Yt(i)/ys)-gs*log(Gt(i)/gs)-log(Iit(i,:)./xis')*xis/N;
                   UCL(1)*Lsjs(1)*log(Lsjt(i,1)/Lsjs(1))- ...
                     UCL(2:Nw).*Lsjs(2:Nw).*log(Lsjt(i,2:Nw)'./Lsjs(2:Nw));
                   log(Mt(i)/Pt(i) /mean(mdjs));
                   ULM(1)*Lsjs(1)*log(Lsjt(i,1)/Lsjs(1))- ...
                     ULM(2:Nw).*Lsjs(2:Nw).*log(Lsjt(i,2:Nw)'./Lsjs(2:Nw))];
    cjt(i,:)  = tem(1:Nw)';
    mdjt(i,:) = tem(Nw+1:2*Nw)';
    tem       = Aa\[ys*log(Yat(i)/ys)-gs*log(Gat(i)/gs)- ...
                       log(Iait(i,:)./xis')*xis/N;
                   UCL(1)*Lsjs(1)*log(Lsajt(i,1)/Lsjs(1))- ...
                     UCL(2:Nw).*Lsjs(2:Nw).*log(Lsajt(i,2:Nw)'./Lsjs(2:Nw));
                   log(Mat(i)/Pat(i)/mean(mdjs));
                   ULM(1)*Lsjs(1)*log(Lsajt(i,1)/Lsjs(1))- ...
                     ULM(2:Nw).*Lsjs(2:Nw).*log(Lsajt(i,2:Nw)'./Lsjs(2:Nw))];
    cajt(i,:) = tem(1:Nw)';
    mdajt(i,:)= tem(Nw+1:2*Nw)';
  end;
  Cjt      = exp(cjt).*( ones(nt,1)*cjs');
  Cajt     = exp(cajt).*( ones(nt,1)*cjs');
  Ct       = mean(Cjt')';
  Cat      = mean(Cajt')';
  Mdjt     = exp(mdjt).*( ones(nt,1)*mdjs');
  Mdajt    = exp(mdajt).*( ones(nt,1)*mdjs');
  RGDP     = Yt + mean(Yhait')' - mean(Yfit')';
  RGDPa    = Yat + mean(Yfit')' - mean(Yhait')';
  [U,UC,UM]= util(Cjt,Mdjt,Lsjt,param);
  [Ua,UCa] = util(Cajt,Mdajt,Lsajt,param);
  qt       = UCa(:,1)./UC(:,1);
  et       = qt.*Pt./Pat;
  pit      = Pt(2:nt) ./Pt(1:nt-1);  
  pit      = [pit(1);pit];
  rt       = UM(:,1) ./(UC(:,1) -UM(:,1))*100;
  rrt      = ((1+rt(1:nt-1)/100) ./pit(2:nt)-1)*100;
  rrt(nt)  = rrt(nt-1);

  imp      = [Ct,Lt,RGDP,Cat,Lat,RGDPa,Mt./(mu.^[1:nt]'), ...
             (Pt./mu.^[0:nt-1]'),et,qt,rt,rrt];                 % in levels

  imp(:,1:10)  ...
           = (imp(:,1:10)./(ones(nt,1)*imp(1,1:10)) - 1)*100;   % put in %
  imp      = imp(nimp+1:nimp+nti,:);

else
  
  imp      = 0;

end;


if (nts > 0 & nsim > 0);
  
  mshocks  = shocks(:,1:nsim);
  mashocks = shocks(:,nsim+1:2*nsim);
  gshocks  = shocks(:,2*nsim+1:3*nsim);
  gashocks = shocks(:,3*nsim+1:4*nsim);
  ashocks  = shocks(:,4*nsim+1:5*nsim);
  aashocks = shocks(:,5*nsim+1:6*nsim);
  Phit     = zeros(nts,N);
  Pfit     = zeros(nts,N);
  Pfait    = zeros(nts,N);
  Phait    = zeros(nts,N);
  Wit      = zeros(nts,Nw);
  Wait     = zeros(nts,Nw);

  for n = 1:nsim;

    inputs   = [mshocks(:,n),mashocks(:,n),gshocks(:,n),gashocks(:,n), ...
                ashocks(:,n),aashocks(:,n)];
    [Xt,Xtl] = dlsim(G, H, G, H, inputs, [Xbar;1]);
    pht      = exp(Xt(:,1));
    pft      = exp(Xt(:,2));
    pfat     = exp(Xt(:,3));
    phat     = exp(Xt(:,4));
    Kit      = exp(Xt(:,5:N+4));
    Kait     = exp(Xt(:,N+5:2*N+4));
    Kitl     = exp(Xtl(:,5:N+4));
    Kaitl    = exp(Xtl(:,N+5:2*N+4));
    Yt       = exp(Xt(:,2*N+5));
    Yat      = exp(Xt(:,2*N+6));
    wt       = exp(Xt(:,2*N+7));
    wat      = exp(Xt(:,2*N+8));
    Mut      = exp(Xtl(:,nZ+1)); 
    Muat     = exp(Xtl(:,nZ+N+1));
    Gt       = exp(Xtl(:,nZ+2*N+1));
    Gat      = exp(Xtl(:,nZ+2*N+3));
    At       = exp(Xtl(:,nZ+2*N+5));
    Aat      = exp(Xtl(:,nZ+2*N+7));
    Mt       = cumprod(Mut);
    Mat      = cumprod(Muat);
    Pht      = pht.*Mt./Mut;
    Pft      = pft.*Mt./Mut;
    Pfat     = pfat.*Mat./Muat;
    Phat     = phat.*Mat./Muat;
    Wt       = wt.*Mt./Mut;
    Wat      = wat.*Mat./Mut;
    phvec    = [phs*ones(1,N-1)./mu.^[N-1:-1:1],Pht']';
    pfvec    = [phs*ones(1,N-1)./mu.^[N-1:-1:1],Pft']';
    pfavec   = [phs*ones(1,N-1)./mu.^[N-1:-1:1],Pfat']';
    phavec   = [phs*ones(1,N-1)./mu.^[N-1:-1:1],Phat']';
    wvec     = [ws*ones(1,Nw-1)./mu.^[Nw-1:-1:1],Wt']';
    wavec    = [ws*ones(1,Nw-1)./mu.^[Nw-1:-1:1],Wat']';
    for i=1:N;
      Phit(:,i)  = phvec(N-i+1:nts-i+N);
      Pfit(:,i)  = pfvec(N-i+1:nts-i+N);
      Pfait(:,i) = pfavec(N-i+1:nts-i+N);
      Phait(:,i) = phavec(N-i+1:nts-i+N);
    end;
    for i=1:Nw;
      Wit(:,i)   = wvec(Nw-i+1:nts-i+Nw);
      Wait(:,i)  = wavec(Nw-i+1:nts-i+Nw);
    end;
    Phbt     = (mean((Phit').^(-theta/thet1))').^(-thet1/theta);    
    Pfbt     = (mean((Pfit').^(-theta/thet1))').^(-thet1/theta);    
    Pfabt    = (mean((Pfait').^(-theta/thet1))').^(-thet1/theta);    
    Phabt    = (mean((Phait').^(-theta/thet1))').^(-thet1/theta);    
    Wbt      = (mean((Wit').^(-v/v1))').^(-v1/v);
    Wabt     = (mean((Wait').^(-v/v1))').^(-v1/v);
    Pt       = (ome1^sup1*Phbt.^sup2+ome2^sup1*Pfbt.^sup2).^(1/sup2);
    Pat      = (ome1^sup1*Pfabt.^sup2+ome2^sup1*Phabt.^sup2).^(1/sup2);
    Yhit     = (((ome1*Pt).^sup1.*Phbt.^sup3.*Yt)*ones(1,N))./Phit.^(1/thet1);
    Yfit     = (((ome2*Pt).^sup1.*Pfbt.^sup3.*Yt)*ones(1,N))./Pfit.^(1/thet1);
    Yfait    = (((ome1*Pat).^sup1.*Pfabt.^sup3.*Yat)*ones(1,N))./Pfait.^ ...
               (1/thet1);
    Yhait    = (((ome2*Pat).^sup1.*Phabt.^sup3.*Yat)*ones(1,N))./Phait.^ ...
               (1/thet1);
    Ldit     = ((Yhit+Yhait)./(Kitl(:,[N,1:N-1]).^alpha)).^(1/alph1);
    Ldait    = ((Yfait+Yfit)./(Kaitl(:,[N,1:N-1]).^alpha)).^(1/alph1);
    Ldit     = Ldit./(At*ones(1,N));
    Ldait    = Ldait./(Aat*ones(1,N));
    Ldt      = mean(Ldit')';            
    Ldat     = mean(Ldait')';          
    KK       = Kit./Kitl(:,[N,1:N-1]);
    KKa      = Kait./Kaitl(:,[N,1:N-1]);
    KK       = min(KK,(1+2*b)/(2*b)-1e-8);
    KKa      = min(KKa,(1+2*b)/(2*b)-1e-8);
    Iit      = (1+b*delta-sqrt(1+2*b-2*b*KK))/b .* Kitl(:,[N,1:N-1]);     
    Iait     = (1+b*delta-sqrt(1+2*b-2*b*KKa))/b.* Kaitl(:,[N,1:N-1]);     
    It       = mean(Iit')';          
    Iat      = mean(Iait')';         
    Lsjt     = ( (Ldt.*Wbt.^(1/v1))*ones(1,Nw) )./ Wit.^(1/v1);
    Lsajt    = ( (Ldat.*Wabt.^(1/v1))*ones(1,Nw) )./ Wait.^(1/v1);
    Lt       = mean(Lsjt')';
    Lat      = mean(Lsajt')';


    for i=1:nts;
      tem       = Aa\[ys*log(Yt(i))-gs*log(Gt(i))-log(Iit(i,:))*xis/N;
                     UCL(1)*Lsjs(1)*log(Lsjt(i,1))- ...
                       UCL(2:Nw).*Lsjs(2:Nw).*log(Lsjt(i,2:Nw)');
                     log(Mt(i)/Pt(i));
                     ULM(1)*Lsjs(1)*log(Lsjt(i,1))- ...
                       ULM(2:Nw).*Lsjs(2:Nw).*log(Lsjt(i,2:Nw)')];
      cjt(i,:)  = tem(1:Nw)';
      mdjt(i,:) = tem(Nw+1:2*Nw)';
      tem       = Aa\[ys*log(Yat(i))-gs*log(Gat(i))-log(Iait(i,:))*xis/N;
                     UCL(1)*Lsjs(1)*log(Lsajt(i,1))- ...
                       UCL(2:Nw).*Lsjs(2:Nw).*log(Lsajt(i,2:Nw)');
                     log(Mat(i)/Pat(i));
                     ULM(1)*Lsjs(1)*log(Lsajt(i,1))- ...
                       ULM(2:Nw).*Lsjs(2:Nw).*log(Lsajt(i,2:Nw)')];
      cajt(i,:) = tem(1:Nw)';
      mdajt(i,:)= tem(Nw+1:2*Nw)';
    end;
    Cjt      = exp(cjt);
    Cajt     = exp(cajt);
    Ct       = mean(Cjt')'; 
    Cat      = mean(Cajt')'; 
    Mdjt     = exp(mdjt);
    Mdajt    = exp(mdajt);
    RGDP     = Yt + mean(Yhait')' - mean(Yfit')';
    RGDPa    = Yat + mean(Yfit')' - mean(Yhait')';
    [U,UC]   = util(Cjt,Mdjt,Lsjt,param);
    [Ua,UCa] = util(Cajt,Mdajt,Lsajt,param);
    qt       = UCa(:,1)./UC(:,1);
    et       = qt.*Pt./Pat;
    nex1     = et.*(mean(Phait'.*Yhait')');
    nex2     = mean(Phit'.*Yhait')';
    nim1     = mean(Pfit'.*Yfit')';
    nim2     = et.*mean(Pfait'.*Yfit')';
    rex1     = mean(Yhait')';
    rim1     = mean(Yfit')';
    tot1     = Pfbt./(et.*Phabt);
    tot2     = et.*Pfabt./Phbt;
    tot3     = Pfbt./Phbt;
    nnxt     = nex1-nim1;
    rnxt     = rex1-rim1;
    NGDP     = Pt .*Yt +nnxt;
    nnxy     = nnxt ./NGDP;
    rnxy     = rnxt ./RGDP;
    phi      = b/2*(Iit./Kitl(:,[N,1:N-1])-delta).^2 .* Kitl(:,[N,1:N-1]);
    adj      = phi./Iit * 100;
    adj      = mean(adj')';

    sim((n-1)*nts+1:n*nts,:)  ...
             = [Ct,It,Lt,RGDP,NGDP,Cat,Iat,Lat,RGDPa,Pat./Pt,et,qt,rnxy, ...
                adj,nex1,nex2,nim1,nim2,rex1,rim1,tot1,tot2,tot3];

  end;

else
    
  sim      = 0;

end;

return

% A check on solution of A:
%   (1) add the following lines after a0(:,1) =  ...
%   (2) check that sa is roughly 0.

tem                  = zeros(nz,nZ);
tem(:,1:nz:(N-1)*nz) = dR(:,mN*nz+1:mN*nz+N-1);
tem(:,2:nz:(N-1)*nz) = dR(:,mN*nz+N:mN*nz+2*N-2);
tem(:,3:nz:(N-1)*nz) = dR(:,mN*nz+2*N-1:mN*nz+3*N-3);
tem(:,4:nz:(N-1)*nz) = dR(:,mN*nz+3*N-2:mN*nz+4*N-4);
tem(:,5:nz-4)        = dR(:,mN*nz+4*N-3:mN*nz+6*N-4);
tem(:,nz-1:nz:(Nw-1)*nz) ....
                     = dR(:,mN*nz+6*N-3:mN*nz+6*N+Nw-5);
tem(:,nz:nz:(Nw-1)*nz) ....
                     = dR(:,mN*nz+6*N+Nw-4:mN*nz+nX);
a0(:,mN+1)           = vec(tem); 


