function [imp,sim] = isimem(param,shocks);
%
%ISIMEM     computes equilibria for the International-Sticky business
%           cycle model with Incomplete Markets and Exogenous Money growth 
%           in Chari, Kehoe, and McGrattan, ``Can Sticky Price Models 
%           Generate Volatile and Persistent Real Exchange Rates?''
%
%           The optimization problems being solved are as follows.
%           Home CONSUMERS solve:
%
%               max           sum sum beta^t pi(s^t) U[ C(s^t,L(s^t),M(s^t)]
%             C(s^t),L(s^t)    t  s^t                                ------
%             M(s^t),B(s^t)                                          P(s^t)
%
%             subject  to
%                 P(s^t) C(s^t) + M(s^t) 
%                   +  sum_{s_{t+1}} Q(s^{t+1}|s^t) B(s^{t+1}) + V(s^t) D(s^t) 
%                  <=  P(s^t) W(s^t) L(s^t) + M(s^{t-1}) + B(s^t) 
%                           + D(s^{t-1}) + Profits(s^t) - Tax(s^t)
%
%           Foreign CONSUMERS budget:
%
%                 P*(s^t) C*(s^t) + M*(s^t) 
%                   +  sum_{s_{t+1}} Q*(s^{t+1}|s^t) B*(s^{t+1})
%                   +  V(s^t) D*(s^t)/e(s^t) 
%                  <=  P*(s^t) W*(s^t) L*(s^t) + M*(s^{t-1}) + B*(s^t) 
%                          + D*(s^{t-1})/e(s^t) + Profits*(s^t) - Tax*(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)W(s^t)L(i,s^t)
%                                    - P(s^t)X(i,s^t)]
%
%                       s.t. Yh(i,s^t)+Yh^*(i,s^t) = F(K(i,s^{t-1}),L(i,s^t))
%
%                            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, 11-2-99
%           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,L)   = K^alpha * L^(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)
%
%  
%______________________________________________________________________________
%

T        = param(1);
N        = param(2);
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);
nti      = param(47);
nts      = param(48);
nsim     = param(49);

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

%______________________________________________________________________________



%______________________________________________________________________________
%
% STEADY STATE FOR THE ECONOMY
%

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

%______________________________________________________________________________



%______________________________________________________________________________
%
% LINEARIZED FIRST ORDER CONDITIONS
%

% The ordering of variables to the residuals is
%
%     ZZ  = [z_{t+N-1},z_{t+N-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},d*_{t-1},
%            mu_{t+N-1},...,mu_{t-N+1},mu*_{t+N-1},...,mu*_{t-N+1}]'
%
%     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}
%              d*_t
%              y_t   
%              y*_t   ]      (2*N+7 x 1)
%

phs      = ss(1);
ys       = ss(2);
ks       = ss(14:13+N);
z        = [phs*ones(4,1);ks;ks;1;ys;ys];
nz       = 2*N+7;
tem      = z*ones(1,N);
ZZ       = log([tem(:);phs*ones(4*N-4,1);ks;ks;1;mu*ones(4*N-2,1)]);
nZZ      = 2*N*N+17*N-5;
del      = max(abs(ZZ)*1e-5,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_isimem(ZZp,param)-res_isimem(ZZm,param))/(2*del(i));
end;
%______________________________________________________________________________


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

%
% The residuals are given by
%
%     A1*[X_{t+1},Z_{t+N-1}]' + A2*[X_t,Z_{t+N-2}]' + shock terms = 0
%
% where
%
%     Z_t = [z_t,z_{t-1}, ....z_{t-N+2}]'    (nz*(N-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}
%              d*_t
%              y_t   
%              y*_t   ]      (2*N+7 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}
%              d*_{t-1}]    (6*N-3 x 1)
%
% The solution to the deterministic part of the problem is given by
%
%     Z_t = A Z_{t-1}
%
%

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

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

[V,L]    = eig(A2,-A1);
L        = diag(L);
[D,I]    = sort(abs(L));
if (~( (D(nX) <= 1.0000005) & (D(nX+1) > 1.0000005) ) )
  i      = sum(D<=1.0000005);
  disp(sprintf( ...
  'WARNING: program assumes %g eigenvalues inside unit circle, found %g',nX,i));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((N-1)*nz);
A(1:nz,1:nz-2)     = W(nA1-nz+1:nA1,[1,N,2*N-1,3*N-2,4*N-3:6*N-3]);
A(1:nz,nz+1:nz:nZ) = W(nA1-nz+1:nA1,2:N-1);
A(1:nz,nz+2:nz:nZ) = W(nA1-nz+1:nA1,N+1:2*N-2);
A(1:nz,nz+3:nz:nZ) = W(nA1-nz+1:nA1,2*N:3*N-3);
A(1:nz,nz+4:nz:nZ) = W(nA1-nz+1:nA1,3*N-1:4*N-4);
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+N-1} + a_1*Z_{t+N-2} .... + a_{N-1}*Z_t + a_N*Z_{t-1}
% + b_0*S_{t+N-1} + b_1*S_{t+N-2} .... + b_{N-1}*S_t               = 0
%
% where
%
%     S_t = [mu_t
%             :
%            mu_{t-N+1}
%            mu*_t
%             :
%            mu*_{t-N+1}]  (2N x 1)
%

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

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

ss             = [eye(nz);zeros(nZ-nz,nz)];
sa             = dR(:,1:nZ);
P              = kron(eye(2),[rhom,zeros(1,N-1);
                              eye(N-1),zeros(N-1,1)]);
sum1           = 0;
sum2           = 0;
aa             = zeros(nz,nZ);
bb             = zeros(nz,nS);
for i=1:N;
  sum1         = sum1 + kron(sa*ss,(P^(N-i))');
  aa(:)        = a0(:,i+1);
  sa           = sa*A + aa;
  bb(:)        = b0(:,i);
  sum2         = sum2 + (P^(N-i))'*bb';
end;
Psi            = [[rhom;zeros(N-2,1)],eye(N-1)];
dd             = [kron(eye(8),Psi), zeros(8*N-8,nz*nS-8*N)
                  zeros(nz*nS-8*N,8*N),eye(nz*nS-8*N)];
qq             = dd*sum1;
rr             = dd*vec(sum2);

nzero          = [2:N,N+2:2*N,2*N+2:3*N,3*N+2:4*N, ...
                  4*N+2:5*N,5*N+2:6*N,6*N+2:7*N,7*N+2:nS*nz];
B1             = zeros(nS,nz);
B1(nzero)      = -qq(:,nzero)\rr;
B              = ss*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,N-1);
Xbar      = log([tem(:);mu*ones(2*N,1)]);
G         = [A, B;
             zeros(nS,nZ), P];
G         = [G, (eye(nZ+nS)-G)*Xbar;
             zeros(1,nZ+nS),1];
H         = [zeros(nZ,2); 
             sigm,corm*sigm/sqrt(1-corm*corm);
             zeros(N-1,2); 
             0,sigm/sqrt(1-corm*corm); 
             zeros(N-1,2); 
             0,0];

if nti > 0;
  
  nimp     = N-1;
  nt       = nti+nimp;
  inputs   = [zeros(nimp-1,2);1,0;zeros(nti,2)];
  [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));
  dat      = Xt(:,2*N+5);
  datl     = Xtl(:,2*N+5);
  Yt       = exp(Xt(:,2*N+6));
  Yat      = exp(Xt(:,2*N+7));
  Mut      = exp(Xtl(:,nZ+1)); 
  Muat     = exp(Xtl(:,nZ+N+1));
  Mt       = cumprod(Mut);
  Mat      = cumprod(Muat);
  Pht      = pht.*Mt./Mut;
  Pft      = pft.*Mt./Mut;
  Pfat     = pfat.*Mat./Muat;
  Phat     = phat.*Mat./Muat;
  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']';
  Phit     = zeros(nt,N);
  Pfit     = zeros(nt,N);
  Pfait    = zeros(nt,N);
  Phait    = zeros(nt,N);
  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;
  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);    
  Pt       = (ome1^sup1*Phbt.^sup2+ome2^sup1*Pfbt.^sup2).^(1/sup2);
  Pat      = (ome1^sup1*Pfabt.^sup2+ome2^sup1*Phabt.^sup2).^(1/sup2);
  MPt      = Mt./Pt;
  MPat     = Mat./Pat;
  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);
  Lit      = ((Yhit+Yhait)./(Kitl(:,[N,1:N-1]).^alpha)).^(1/alph1);
  Lait     = ((Yfait+Yfit)./(Kaitl(:,[N,1:N-1]).^alpha)).^(1/alph1);
  Lt       = mean(Lit')';            
  Lat      = mean(Lait')';          
  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')';         
  Ct       = Yt-mean(Iit')';       
  Cat      = Yat-mean(Iait')';     
  RGDP     = Yt + mean(Yhait')' - mean(Yfit')';
  RGDPa    = Yat + mean(Yfit')' - mean(Yhait')';
  [U,UC,UM]= util(Ct,MPt,Lt,param);
  [Ua,UCa] = util(Cat,MPat,Lat,param);
  Dat      = dat.*Mt;
  Datl     = datl.*Mt./Mut;
  Vt       = 1-UM./UC;
  et       = (Vt.*Dat-Datl-mean(Pfit'.*Yfit')')./ ...
             (mean(Pfait'.*Yfait')'-Pat.*Yat);
  qt       = et.*Pat./Pt;
  pit      = Pt(2:nt) ./Pt(1:nt-1);  
  pit      = [pit(1);pit];
  rt       = UM ./(UC -UM)*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);
  Phit     = zeros(nts,N);
  Pfit     = zeros(nts,N);
  Pfait    = zeros(nts,N);
  Phait    = zeros(nts,N);

  for n = 1:nsim;

    inputs   = [mshocks(:,n),mashocks(:,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));
    dat      = Xt(:,2*N+5);
    datl     = Xtl(:,2*N+5);
    Yt       = exp(Xt(:,2*N+6));
    Yat      = exp(Xt(:,2*N+7));
    Mut      = exp(Xtl(:,nZ+1)); 
    Muat     = exp(Xtl(:,nZ+N+1));
    Mt       = cumprod(Mut);
    Mat      = cumprod(Muat);
    Pht      = pht.*Mt./Mut;
    Pft      = pft.*Mt./Mut;
    Pfat     = pfat.*Mat./Muat;
    Phat     = phat.*Mat./Muat;
    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']';
    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;
    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);    
    Pt       = (ome1^sup1*Phbt.^sup2+ome2^sup1*Pfbt.^sup2).^(1/sup2);
    Pat      = (ome1^sup1*Pfabt.^sup2+ome2^sup1*Phabt.^sup2).^(1/sup2);
    MPt      = Mt./Pt;
    MPat     = Mat./Pat;
    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);
    Lit      = ((Yhit+Yhait)./(Kitl(:,[N,1:N-1]).^alpha)).^(1/alph1);
    Lait     = ((Yfait+Yfit)./(Kaitl(:,[N,1:N-1]).^alpha)).^(1/alph1);
    Lt       = mean(Lit')';            
    Lat      = mean(Lait')';          
    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')';         
    Ct       = Yt-mean(Iit')';       
    Cat      = Yat-mean(Iait')';     
    RGDP     = Yt + mean(Yhait')' - mean(Yfit')';
    RGDPa    = Yat + mean(Yfit')' - mean(Yhait')';
    [U,UC,UM]= util(Ct,MPt,Lt,param);
    [Ua,UCa] = util(Cat,MPat,Lat,param);
    Dat      = dat.*Mt;
    Datl     = datl.*Mt./Mut;
    Vt       = 1-UM./UC;
    et       = (Vt.*Dat-Datl-mean(Pfit'.*Yfit')')./ ...
               (mean(Pfait'.*Yfait')'-Pat.*Yat);
    qt       = et.*Pat./Pt;
    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:nZ) = dR(:,N*nz+1:N*nz+N-1);
tem(:,2:nz:nZ) = dR(:,N*nz+N:N*nz+2*N-2);
tem(:,3:nz:nZ) = dR(:,N*nz+2*N-1:N*nz+3*N-3);
tem(:,4:nz:nZ) = dR(:,N*nz+3*N-2:N*nz+4*N-4);
tem(:,5:nz-2)  = dR(:,N*nz+4*N-3:N*nz+nX);
a0(:,N+1)      = vec(tem); 
