function [imp,sim] = iscmtr1_ms(param,shocks);
%
%ISCMTR1_MS computes equilibria for the International-Sticky business
%           cycle model with Complete Markets and a Taylor Rule governing
%           interest rates in Chari, Kehoe, and McGrattan, ``Can Sticky 
%           Price Models Generate Volatile and Persistent Real Exchange 
%           Rates?'' This code assumes N==1.  It is an extension of 
%           iscmtr1.m to include government spending and technology shocks.
%
%           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})
%                   <=  P(s^t)*W(s^t)*L(s^t)) + M(s^{t-1}) + B(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.
%
%
%           INTEREST RATES:
%
%              r(s^t) = irule' [     r(s^{t-1})
%                                    r(s^{t-2})
%                                    r(s^{t-3})
%                                E_t log P(s^{t+1})-log P(s^t)
%                                    log P(s^t)    -log P(s^{t-1})
%                                    log P(s^{t-1})-log P(s^{t-2})
%                                    log P(s^{t-2})-log P(s^{t-3})
%                                    log gdp(s^t)
%                                    log gdp(s^{t-1})
%                                    log gdp(s^{t-2})  ]  + eps_r(s^t)
%

%           Ellen McGrattan, 9-15-99
%           Revised, ERM, 7-24-01

%______________________________________________________________________________
%
%PARAMETERS:
%
%  Staggered contracts:
%
%    T = fraction of year prices are fixed (T=1 ==> fixed for 1 year)
%    N = number of cohorts (= 1) 
%
%  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
%
%  Taylor rule for interest rates:
%
%    r[t]     = irule' * [rhs variables ] + constant + eps_r[t]
%    mu       = mean rate of inflation and money growth
%
%  Timing of shocks:
%
%    if timing = 1, monopolists observe only lagged shocks when choosing p[t-1]
%              = 2, monopolists observe lagged shocks and a[t]
%              = 3, monopolists observe lagged shocks and g[t],a[t]
%______________________________________________________________________________
%

T        = param(1);
N        = 1;
Utype    = param(4);
beta     = param(11);
ome1     = param(12);
rho      = param(13);
theta    = param(14);
alpha    = param(15);
delta    = param(16);
b        = param(17);
mu       = param(19);
rhoe     = param(33);
sige     = param(34);
core     = param(35);
gs       = param(36);
rhog     = param(37);
sigg     = param(38);
corg     = param(39);
rhoa     = param(40);
siga     = param(41);
cora     = param(42);
timing   = param(43);
nti      = param(47);
nts      = param(48);
nsim     = param(49);
rimp     = 1;
gimp     = 0;
aimp     = 0;

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

%______________________________________________________________________________



%______________________________________________________________________________
%
% STEADY STATE FOR THE ECONOMY
%

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

%______________________________________________________________________________



%______________________________________________________________________________
%
% LINEARIZED FIRST ORDER CONDITIONS
%

% The ordering of variables to the residuals is
%
%     ZZ  = [z_{t+1},z_t,ph_{t-2},ph_{t-3},ph_{t-4},
%            pf_{t-2},pf_{t-3},pf_{t-4},pf*_{t-2},pf*_{t-3},pf*_{t-3},
%            ph*_{t-2},ph*_{t-3},ph*_{t-4},k_{t-1},k*_{t-1},
%            y_{t-1},y_{t-2},y*_{t-1},y*_{t-2},
%            r_{t-1},r_{t-2},r_{t-3},r*_{t-1},r*_{t-2},r*_{t-3},
%            eps_{r,t},eps*_{r,t},g_{t+1},g_t,g*_{t+1},g*_t, 
%            a_{t+1},a_t,a*_{t+1},a*_t]'
%
%     z_t = [ph_{t-1}          NOTES: (1) variables are logged and demeaned
%            pf_{t-1}                 (2) prices are normalized by mu^t
%            pf*_{t-1}                     
%            ph*_{t-1}
%            k_t
%            k*_t
%            y_t   
%            y*_t   
%            r_t   
%            r*_t   ]      (10 x 1)
%

phs      = ss(1);
ys       = ss(2);
ks       = ss(14);
rs       = mu/beta-1;
z        = [phs*ones(4,1);ks;ks;ys;ys;exp(rs);exp(rs)];
nz       = 10;
ZZ       = log([z;z;phs*ones(12,1);ks;ks;ys;ys;ys;ys;exp(rs*ones(6,1));
                1;1;gs;gs;gs;gs;1;1;1;1]);
nZZ      = 2*nz+34;
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_iscmtr1_ms(ZZp,param)-res_iscmtr1_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+1}]' + A2*[X_t,z_t]' + shock terms = 0.
%
% where 
%
%     z_t = [ph_{t-1}
%            pf_{t-1}
%            pf*_{t-1}
%            ph*_{t-1}
%            k_t
%            k*_t
%            y_t
%            y*_t
%            r_t
%            r*_t   ]      (10 x 1)
%
%     X_t = [ph_{t-2},
%              :
%            ph_{t-4},
%            pf_{t-2},
%              :
%            pf_{t-4}
%            pf*_{t-2}
%              :
%            pf*_{t-4}
%            ph*_{t-2}
%              :
%            ph*_{t-4}
%            k_{t-1}
%            k*_{t-1}
%            y_{t-1}
%            y_{t-2}
%            y*_{t-1}
%            y*_{t-2}
%            r_{t-1}
%            r_{t-2}
%            r_{t-3}
%            r*_{t-1}
%            r*_{t-2}
%            r*_{t-3}]   (24 x 1)
%
% The solution to the deterministic part of the problem is given by
%
%     Z_t = A Z_{t-1}
%
% where Z_t = [z_t,z_{t-1},z_{t-2}]'    (3*nz x 1)
%

nX       = 24;
nZ       = 3*nz;
nA1      = nX+nz;

A1       = [eye(nX,nX),zeros(nX,nz);
            zeros(nz,nX),dR(:,1:nz)];

I1       = [kron(eye(4),[zeros(1,3);
                         eye(2),zeros(2,1)]),zeros(12,12);
            zeros(2,nX);
            zeros(4,14),kron(eye(2),[0,0;1,0]),zeros(4,6);
            zeros(6,18),kron(eye(2),[0,0,0;1,0,0;0,1,0])];
I2       = [kron(eye(4),[1;zeros(2,1)]),zeros(12,6);
            zeros(2,4),eye(2),zeros(2,4);
            zeros(4,6),kron(eye(2),[1;0]),zeros(4,2);
            zeros(6,8),kron(eye(2),[1;0;0])];
A2       = [-I1,-I2;
            dR(:,2*nz+1:2*nz+nX),dR(:,nz+1:2*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(nZ);
A(1:nz,1:nz)       = W(nA1-nz+1:nA1,[1,4,7,10,13:15,17,19,22]);
A(1:nz,nz+1:nz:nZ) = W(nA1-nz+1:nA1,2:3);
A(1:nz,nz+2:nz:nZ) = W(nA1-nz+1:nA1,5:6);
A(1:nz,nz+3:nz:nZ) = W(nA1-nz+1:nA1,8:9);
A(1:nz,nz+4:nz:nZ) = W(nA1-nz+1:nA1,11:12);
A(1:nz,2*nz-3)     = W(nA1-nz+1:nA1,16);
A(1:nz,2*nz-2)     = W(nA1-nz+1:nA1,18);
A(1:nz,2*nz-1)     = W(nA1-nz+1:nA1,20);
A(1:nz,2*nz)       = W(nA1-nz+1:nA1,23);
A(1:nz,3*nz-1)     = W(nA1-nz+1:nA1,21);
A(1:nz,3*nz)       = W(nA1-nz+1:nA1,24);
A(nz+1:nZ,1:nZ-nz) = eye(nZ-nz);

%______________________________________________________________________________


%______________________________________________________________________________
%
% SOLUTION OF COEFFICIENTS ON POLICY SHOCKS (B)
%

% The residuals can be written as:
%
%   a_0*Z_{t+1} + a_1*Z_t + a_2*Z_{t-1}
% + b_0*S_{t+1} + b_1*S_t = 0
%
% where
%
%     S_t = [eps_{r,t}
%            eps_{r,t-1}
%            eps*_{r,t}
%            eps*_{r,t-1}
%            g_t
%            g_{t-1}
%            g*_t
%            g*_{t-1}
%            a_t
%            a_{t-1}
%            a*_t
%            a*_{t-1}]     (12 x 1)
%

nS             = 12;
a0             = [dR(:,1:2*nz),zeros(nz,nz)];
a1             = zeros(nz,nZ);
               % Note: a1=0 and a2 not needed to compute B
b0             = [zeros(nz,1),dR(:,2*nz+nX+1),zeros(nz,1),dR(:,2*nz+nX+2), ...
                  dR(:,2*nz+nX+3:2*nz+nX+10)];
b1             = zeros(nz,nS);

P              = [rhoe,0,0,0,0,0,0,0,0,0,0,0;
                  1,0,0,0,0,0,0,0,0,0,0,0;
                  0,0,rhoe,0,0,0,0,0,0,0,0,0;
                  0,0,1,0,0,0,0,0,0,0,0,0;
                  0,0,0,0,rhog,0,0,0,0,0,0,0;
                  0,0,0,0,1,0,0,0,0,0,0,0;
                  0,0,0,0,0,0,rhog,0,0,0,0,0;
                  0,0,0,0,0,0,1,0,0,0,0,0;
                  0,0,0,0,0,0,0,0,rhoa,0,0,0;
                  0,0,0,0,0,0,0,0,1,0,0,0;
                  0,0,0,0,0,0,0,0,0,0,rhoa,0;
                  0,0,0,0,0,0,0,0,0,0,1,0];

if timing == 2;
  Psi          = [rhoe,1,0,0,0,0,0,0,0,0,0,0;
                  0,0,rhoe,1,0,0,0,0,0,0,0,0;
                  0,0,0,0,rhog,1,0,0,0,0,0,0;
                  0,0,0,0,0,0,rhog,1,0,0,0,0;
                  0,0,0,0,0,0,0,0,1,0,0,0;
                  0,0,0,0,0,0,0,0,0,1,0,0;
                  0,0,0,0,0,0,0,0,0,0,1,0;
                  0,0,0,0,0,0,0,0,0,0,0,1];
  dd           = [kron(eye(4),Psi), zeros(4*(nS-4),nz*nS-4*nS)
                  zeros(nz*nS-4*nS,4*nS),eye(nz*nS-4*nS)];
  nzero        = [2,4,6,8:12,14,16,18,20:24,26,28,30,32:36,38,40,42,44:nz*nS];
elseif timing == 3;
  Psi          = [[rhoe,1,0,0;
                   0,0,rhoe,1],zeros(2,8);
                   zeros(8,4),eye(8)];
  dd           = [kron(eye(4),Psi), zeros(4*(nS-2),nz*nS-4*nS)
                  zeros(nz*nS-4*nS,4*nS),eye(nz*nS-4*nS)];
  nzero        = [2,4:12,14,16:24,26,28:36,38,40:nz*nS];
else;
  Psi          = [rhoe,1,0,0,0,0,0,0,0,0,0,0;
                  0,0,rhoe,1,0,0,0,0,0,0,0,0;
                  0,0,0,0,rhog,1,0,0,0,0,0,0;
                  0,0,0,0,0,0,rhog,1,0,0,0,0;
                  0,0,0,0,0,0,0,0,rhoa,1,0,0;
                  0,0,0,0,0,0,0,0,0,0,rhoa,1];
  dd           = [kron(eye(4),Psi), zeros(4*(nS-6),nz*nS-4*nS)
                  zeros(nz*nS-4*nS,4*nS),eye(nz*nS-4*nS)];
  nzero        = [2:2:48,49:nz*nS];
end;
qq             = dd*(kron(a0*A+a1,eye(nS))+kron(a0,P'));
rr             = dd*vec(b1'+P'*b0');
B1             = zeros(nS,nz);
B1(nzero)      = -qq(:,nzero)\rr;
B              = [eye(nz);zeros(nZ-nz,nz)]*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)

Xbar      = log([z;z;z;1;1;1;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);
             sige,core*sige/sqrt(1-core*core),0,0,0,0;
             0,0,0,0,0,0;
             0,sige/sqrt(1-core*core),0,0,0,0;
             0,0,0,0,0,0;
             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];


cs        = ss(3);
MPs       = ss(4);
ls        = ss(18);
[Us,UCs,UMs,ULs,UCCs,UMMs,ULLs,UCMs,UCLs,ULMs] = util(cs,MPs,ls,param);

if nti > 0;
  
  nt       = nti+1;
  inputs   = [rimp,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));
  Kt       = exp(Xt(:,5));
  Kat      = exp(Xt(:,6));
  Ktl      = exp(Xtl(:,5));
  Katl     = exp(Xtl(:,6));
  Yt       = exp(Xt(:,7));
  Yat      = exp(Xt(:,8));
  rt       = Xt(:,9);
  rat      = Xt(:,10);
  Gt       = exp(Xtl(:,nZ+5));
  Gat      = exp(Xtl(:,nZ+7));
  At       = exp(Xtl(:,nZ+9));
  Aat      = exp(Xtl(:,nZ+11));
  Pht      = pht .* mu.^[0:nt-1]';
  Pft      = pft .* mu.^[0:nt-1]';
  Pfat     = pfat.* mu.^[0:nt-1]';
  Phat     = phat.* mu.^[0:nt-1]';
  Pt       = (ome1^sup1*Pht.^sup2+ome2^sup1*Pft.^sup2).^(1/sup2);
  Pat      = (ome1^sup1*Pfat.^sup2+ome2^sup1*Phat.^sup2).^(1/sup2);
  Yht      = (ome1*Pt./Pht).^sup1.*Yt;
  Yft      = (ome2*Pt./Pft).^sup1.*Yt;
  Yfat     = (ome1*Pat./Pfat).^sup1.*Yat;
  Yhat     = (ome2*Pat./Phat).^sup1.*Yat;
  Lt       = (((Yht+Yhat)./(Ktl.^alpha)).^(1/alph1))./At;
  Lat      = (((Yfat+Yft)./(Katl.^alpha)).^(1/alph1))./Aat;
  KK       = Kt./Ktl;
  KKa      = Kat./Katl;
  It       = (1+b*delta-sqrt(1+2*b-2*b*KK))/b .* Ktl;
  Iat      = (1+b*delta-sqrt(1+2*b-2*b*KKa))/b.* Katl;
  Ct       = Yt-It-Gt;
  Cat      = Yat-Iat-Gat;
  MPt      = exp((-(UCMs/UMs-UCCs/UCs)*cs*log(Ct/cs)- ...
                   (ULMs/UMs-UCLs/UCs)*ls*log(Lt/ls)+ ...
                    beta/mu/rs*(rt-rs))/((UMMs/UMs-UCMs/UCs)*MPs))*MPs;
  MPat     = exp((-(UCMs/UMs-UCCs/UCs)*cs*log(Cat/cs)- ...
                   (ULMs/UMs-UCLs/UCs)*ls*log(Lat/ls)+ ...
                    beta/mu/rs*(rat-rs))/((UMMs/UMs-UCMs/UCs)*MPs))*MPs;
  Mt       = MPt.*Pt;
  RGDP     = Yt + Yhat - Yft;
  RGDPa    = Yat + Yft - Yhat;
  [U,UC,UM]= util(Ct,MPt,Lt,param);
  [Ua,UCa] = util(Cat,MPat,Lat,param);
  qt       = UCa./UC;
  et       = qt.*Pt./Pat;
  pit      = Pt(2:nt) ./Pt(1:nt-1);  
  pit      = [pit(1);pit];
  rrt      = ((1+rt(1:nt-1)) ./pit(2:nt)-1)*100;
  rrt(nt)  = rrt(nt-1);
  rt       = rt*100;    

  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(2:1+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);

  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));
    Kt       = exp(Xt(:,5));
    Kat      = exp(Xt(:,6));
    Ktl      = exp(Xtl(:,5));
    Katl     = exp(Xtl(:,6));
    Yt       = exp(Xt(:,7));
    Yat      = exp(Xt(:,8));
    rt       = Xt(:,9);
    rat      = Xt(:,10);
    Gt       = exp(Xtl(:,nZ+5));
    Gat      = exp(Xtl(:,nZ+7));
    At       = exp(Xtl(:,nZ+9));
    Aat      = exp(Xtl(:,nZ+11));
    Pht      = pht .* mu.^[0:nts-1]';
    Pft      = pft .* mu.^[0:nts-1]';
    Pfat     = pfat.* mu.^[0:nts-1]';
    Phat     = phat.* mu.^[0:nts-1]';
    Pt       = (ome1^sup1*Pht.^sup2+ome2^sup1*Pft.^sup2).^(1/sup2);
    Pat      = (ome1^sup1*Pfat.^sup2+ome2^sup1*Phat.^sup2).^(1/sup2);
    Yht      = (ome1*Pt./Pht).^sup1.*Yt;
    Yft      = (ome2*Pt./Pft).^sup1.*Yt;
    Yfat     = (ome1*Pat./Pfat).^sup1.*Yat;
    Yhat     = (ome2*Pat./Phat).^sup1.*Yat;
    Lt       = (((Yht+Yhat)./(Ktl.^alpha)).^(1/alph1))./At;
    Lat      = (((Yfat+Yft)./(Katl.^alpha)).^(1/alph1))./Aat;
    KK       = Kt./Ktl;
    KKa      = Kat./Katl;
    It       = (1+b*delta-sqrt(1+2*b-2*b*KK))/b .* Ktl;
    Iat      = (1+b*delta-sqrt(1+2*b-2*b*KKa))/b.* Katl;
    Ct       = Yt-It-Gt;
    Cat      = Yat-Iat-Gat;
    MPt      = exp((-(UCMs/UMs-UCCs/UCs)*cs*log(Ct/cs)- ...
                     (ULMs/UMs-UCLs/UCs)*ls*log(Lt/ls)+ ...
                      beta/mu/rs*(rt-rs))/((UMMs/UMs-UCMs/UCs)*MPs))*MPs;
    MPat     = exp((-(UCMs/UMs-UCCs/UCs)*cs*log(Cat/cs)- ...
                     (ULMs/UMs-UCLs/UCs)*ls*log(Lat/ls)+ ...
                      beta/mu/rs*(rat-rs))/((UMMs/UMs-UCMs/UCs)*MPs))*MPs;
    Mt       = MPt.*Pt;
    RGDP     = Yt + Yhat - Yft;
    RGDPa    = Yat + Yft - Yhat;
    [U,UC,UM]= util(Ct,MPt,Lt,param);
    [Ua,UCa] = util(Cat,MPat,Lat,param);
    qt       = UCa./UC;
    et       = qt.*Pt./Pat;
    nex1     = et.*Phat.*Yhat;
    nex2     = Pht.*Yhat;
    nim1     = Pft.*Yft;
    nim2     = et.*Pfat.*Yft;
    rex1     = Yhat;
    rim1     = Yft;
    tot1     = Pft./(et.*Phat);
    tot2     = et.*Pfat./Pht;
    tot3     = Pft./Pht;
    nnxt     = nex1-nim1;
    rnxt     = rex1-rim1;
    NGDP     = Pt .*Yt +nnxt;
    nnxy     = nnxt ./NGDP;
    rnxy     = rnxt ./RGDP;
    phi      = b/2*(It./Ktl-delta).^2 .*Ktl;
    adj      = phi./It * 100;

    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 a1 = ...
%   (2) check that a0*A*A+a1*A+a2 is roughly 0.

a2            = zeros(nz,nZ);
a2(:,1:nz:nZ) = dR(:,2*nz+1:2*nz+3);
a2(:,2:nz:nZ) = dR(:,2*nz+N+3:2*nz+6);
a2(:,3:nz:nZ) = dR(:,2*nz+2*N+5:2*nz+9);
a2(:,4:nz:nZ) = dR(:,2*nz+3*N+7:2*nz+12);
a2(:,5:nz-3)  = dR(:,2*nz+4*N+9:2*nz+15);
a2(:,nz-2)    = dR(:,2*nz+17);
a2(:,nz-1)    = dR(:,2*nz+19);
a2(:,nz)      = dR(:,2*nz+22);
a2(:,2*nz-3)  = dR(:,2*nz+16);
a2(:,2*nz-2)  = dR(:,2*nz+18);
a2(:,2*nz-1)  = dR(:,2*nz+20);
a2(:,2*nz)    = dR(:,2*nz+23);
a2(:,3*nz-1)  = dR(:,2*nz+21);
a2(:,3*nz)    = dR(:,2*nz+24);
a0*A*A+a1*A+a2
