%STOKEY2    This is the second part of a sequence of files used to
%           address Nancy Stokey's questions about the paper.
%           This file constructs two approximations to a sticky-price
%           economy (the benchmark economy in the paper with zero-income
%           effects preferences).  The first approximation is the 
%           linear approximation to the economy.  The second is
%           an approximation based on analytical calculations assuming
%           xi --> 1, eta --> 1, and beta=1 in the Euler equation for 
%           capital.  The difference between the economy here and that
%           in stokey1.m is the addition of adjustment costs which
%           are driven to infinity.
%
%           The economy has consumers, final goods producers, and 
%           intermediate goods producers.  First-order conditions
%           include:  
%
%           1.) the pricing equation for intermediate goods producers:
%               E_{t-1} (x_t-1/2 x_{t-1}-1/2 x_{t+1} - h_t - h_{t+1} ) = 0
%
%           2.) the money demand equation from consumer optimization:
%               U_mt - U_ct + beta E_t  U_{ct+1} Pbar_t/Pbar_{t+1} = 0
%
%           3.) the capital euler equation from consumer optimization:
%               U_ct = beta E_t  U_{ct+1}{ r_{t+1} + 1 - delta} = 0
%
%           4.) the resource constraint:
%               C_t + X_t = F(K_t,L_t) = Y_t
%
%           5.) the capital accumulation equation:
%               K_{t+1} = (1-delta) K_t + phi(X_t/K_t) K_t
%
% 
%           The approximate solution is:
%
%           log(y)
%            |
%            |            *                  
%            |             *          z[t]  = a * z[t-1] - mu[t]
%            |                        y[t]  = b * z[t-1] + c * mu[t]
%            |               *        (See notes for formulas for a,b, & c.) 
%            |                 *    
%            |                   *
%            |                     *     
%            |                         *  
%            |************-----------------******
%
%           with the capital stock not moving.
%

%                 

%           Ellen McGrattan, 12-15-97
%           Revised, 12-31-97, ERM


% Parameters:

% Utility  = (C^eta/eta - kappa*L^xi/xi + a*M^omega/omega)^(1-sigma)/(1-sigma)

a       = .04;
alpha   = .33;
b       =  50000;
beta    = .96;
delta   = .06;
eta     = .995;
kappa   = .9;
mu      = 1;
nu      = .995;
theta   = .9;
sigma   = 2.001;
xi      = 1.001;

% Steady state values for decision and state variables:

klratio = ( (1-beta*(1-delta))/(beta *alpha*theta) ) ^(1/(alpha-1));
% or, if beta=1 in Euler equation, then:
klratio = ( delta/(alpha*theta) ) ^(1/(alpha-1));
clratio = klratio^alpha-delta*klratio;
Lss     = (theta*(1-alpha)*klratio^alpha * clratio^(eta-1)/kappa)^(1/(xi-eta));
Css     = clratio*Lss;
Kss     = klratio*Lss;
Mss     = ( (1-beta)*Css^(eta-1)/ a) ^(1/(nu-1));
Yss     = Kss^alpha*Lss^(1-alpha);
v       = Css^eta/eta-kappa*Lss^xi/xi+a*Mss^nu/nu;

UC      = (Css^eta/eta-kappa*Lss^xi/xi+a*Mss^nu/nu)^(-sigma)*Css^(eta-1);
UL      =-(Css^eta/eta-kappa*Lss^xi/xi+a*Mss^nu/nu)^(-sigma)*Lss^(xi-1)*kappa;
UM      = (Css^eta/eta-kappa*Lss^xi/xi+a*Mss^nu/nu)^(-sigma)*Mss^(nu-1)*a;

X       = log([1/Mss;1/Mss;1/Mss;Yss;Yss;Kss;Kss;Kss;1;1;1]);
param   = [a;alpha;b;beta;delta;eta;kappa;nu;sigma;theta;xi];

del     = diag(max(abs(X)*1e-7,1e-8));
for i=1:length(X);
  dR(:,i)= (res_stokey2(X+del(:,i),param)- ...
            res_stokey2(X-del(:,i),param))/(2*del(i,i));
end;
dR1     = dR(1,:)';
dR2     = dR(2,:)';
dR3     = dR(3,:)';


% Coefficients for capital Euler equation if beta<1:

alp1    =-sigma/v*Css^eta-1+eta;       
alp2    = sigma/v*Css^eta+beta*(1-eta)*(1-delta);
alp3    = sigma/v*kappa*Lss^xi;
alp4    =-sigma/v*kappa*Lss^xi-xi*beta*(1/beta-1+delta);
alp5    =-sigma/v*a*Mss^nu+(nu-1)*(1-beta);
alp6    = sigma/v*a*Mss^nu;
alp7    =-b*delta;
alp8    = beta*(1/beta-1+delta)+beta*b*delta;
alp9    = b*delta;
alp10   =-beta*b*delta;


% Coefficients for capital Euler equation if beta=1:

alp1    =-sigma/v*Css^eta-1+eta;       
alp2    = sigma/v*Css^eta+(1-eta)*(1-delta);
alp3    = sigma/v*kappa*Lss^xi;
alp4    =-sigma/v*kappa*Lss^xi-xi*delta;
alp5    =-sigma/v*a*Mss^nu;
alp6    = sigma/v*a*Mss^nu;
alp7    =-b*delta;
alp8    = delta+b*delta;
alp9    = b*delta;
alp10   =-b*delta;

% Coefficients for money demand equation:

a1      =-sigma/v*Css^eta-(1-eta)/beta;
a2      = sigma/v*Css^eta+1-eta;
a3      = sigma/v*kappa*Lss^xi;
a4      =-sigma/v*kappa*Lss^xi;
a5      =-sigma/v*a*Mss^nu-(nu-1)*(1-beta)/beta;
a6      = sigma/v*a*Mss^nu;
a7      =-1;

% Approximation for a1,alp2,a5:

a1      = alp1;
alp2    = a2;
a5      = alp5;

% A check on derivatives:

dr1     = [-1/2;1;-1/2;
           -(alpha+xi-1)/(1-alpha)+(eta-1)*Yss/Css;
           -(alpha+xi-1)/(1-alpha)+(eta-1)*Yss/Css;
            alpha*xi/(1-alpha)+(eta-1)*(1-delta)*Kss/Css;
            alpha*xi/(1-alpha)+(eta-1)*(1-delta)*Kss/Css-(eta-1)*Kss/Css;
           -(eta-1)*Kss/Css;
            1/2;
            0;
            0];

dr2     = [1/2*(-a5+a7);
           1/2*(-a5-a6); 
           1/2*(-a6-a7); 
           a1*Yss/Css+a3/(1-alpha);
           a2*Yss/Css+a4/(1-alpha);
           a1*Kss/Css*(1-delta)-a3*alpha/(1-alpha);
          -a1*Kss/Css+a2*Kss/Css*(1-delta)-a4*alpha/(1-alpha);
          -a2*Kss/Css;
           1/2*(a5-a7);
           1/2*(2*a5+a6-a7);
           0]; 
dr3     = [-1/2*alp5;
           -1/2*(alp5+alp6); 
           -1/2*alp6;
           alp1*Yss/Css+alp3/(1-alpha);
           alp2*Yss/Css+alp4/(1-alpha);
           alp1*Kss/Css*(1-delta)-alp3*alpha/(1-alpha)+ ...
             alp7-(1-delta)*alp9/delta;
          -alp1*Kss/Css+alp2*Kss/Css*(1-delta)-alp4*alpha/(1-alpha)+ ...
             alp8+alp9/delta-(1-delta)*alp10/delta;
          -alp2*Kss/Css+alp10/delta;
           1/2*alp5;
           1/2*(2*alp5+alp6);
           0];


%dR1    = dr1;
%dR3    = dr2-dr3;
%dR2    = dr3;

phi11   =-dR1(6)/dR1(2);
phi12   =-dR1(8)/dR1(2);
phi13   =-dR1(4)/dR1(2);

phi21   =-dR3(6)/dR3(7);
phi22   =-dR3(8)/dR3(7);
phi23   =-dR3(1)/dR3(7);
phi24   =-dR3(2)/dR3(7);
phi25   =-dR3(3)/dR3(7);
phi26   =-dR3(4)/dR3(7);
phi27   =-dR3(5)/dR3(7);
phi28   =-dR3(9)/dR3(7);
phi29   =-dR3(10)/dR3(7);

phi31   =-dR2(5)/dR2(4);
phi32   =-dR2(1)/dR2(4);
phi33   =-dR2(2)/dR2(4);
phi34   =-dR2(3)/dR2(4);
phi35   =-dR2(6)/dR2(4);
phi36   =-dR2(7)/dR2(4);
phi37   =-dR2(8)/dR2(4);
phi38   =-dR2(9)/dR2(4);
phi39   =-dR2(10)/dR2(4);

A1      = [1,0,0,0,0;
           0,1,0,0,0;
           0,0,-1/2,-phi12,-phi13;
           0,0,-phi25,-phi22,-phi27;
           0,0,-phi34,-phi37,-phi31];
A2      = [0,0,-1,0,0;
           0,0,0,-1,0;
          -1/2,-phi11,1,-(phi11+phi12),-phi13;
          -phi23,-phi21,-phi24,1,-phi26;
          -phi32,-phi35,-phi33,-phi36,1];

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

a11     = A(1,1);
a12     = A(1,2);
a21     = A(2,1);
a22     = A(2,2);
a31     = A(3,1);
a32     = A(3,2);

tem     = [1-a11/2-phi12*a21-phi13*a31,0, ...
             -a12/2-phi11-phi12-phi12*a22-phi13*a32, ...
              0,-phi13;
           -phi25,1-phi22*a22-phi25*a12-phi27*a32,-phi22,-phi26,-phi27;
           -phi22*a21-phi24-phi25*a11-phi27*a31,0, ...
              1-phi22*a22-phi25*a12-phi27*a32,0,-phi26;
           -phi34,-phi31*a32-phi34*a12-phi36-phi37*a22,-phi37,1,-phi31;
           -phi31*a31-phi33-phi34*a11-phi37*a21,0, ...
              -phi31*a32-phi34*a12-phi36-phi37*a22,0,1]\ ...
           [-1/2;phi29;phi28;phi39;phi38];
B       = [0,tem(1);tem(2:3)';tem(4:5)'];


psi     = -1/2*sigma/v*(1-beta)*Css^eta;
phi     = (sigma/v*Css^eta+1-eta)*Yss/Css-sigma/v*kappa*Lss^xi/(1-alpha);
h       = 2*alpha/(1-alpha)*(psi+1/2)/phi;
aa      = (1+h-2*sqrt(h))/(1-h);
dd      = -(1-alpha)/alpha * h * (1-sqrt(h))/(1-h);
ee      = (1-alpha)/alpha * h;

%[h, alpha/sigma*((2-beta-sigma*(1-beta))-kappa*Lss^xi/Css^eta)/ ...
%    (Yss/Css*(1-alpha)-kappa*Lss^xi/Css^eta), ...
%    alpha/sigma* ( (2-beta-sigma*(1-beta))*(1-alpha*theta)/(1-alpha)-theta)/ ...
%    (1-theta)]


A0      = [aa,0,0;
           0,1,0;
           dd,0,0];
B0      = [0,-aa;
           0,0;
           ee,-dd];

% Form of final solution:  Z[t] = A* Z[t-1] + B* S[t], where
%   Z[t] = [x[t],k[t+1],y[t]]
%   S[t] = [mu[t],mu[t-1]]

% Impulse responses for:
%
%       [ Z[t]   ] = [ A  B ]  [ Z[t-1] ]      [0 ] 
%       [ S[t+1] ]   [ 0  P ]  [ S[t]   ]   +  [S1] * eta[t+1]
%  or  
%         Xt[t+1]  = G * Xt[t] + H * eta[t+1]
%

% where eta[t] is N(0,1) and sig = [sig_e;0...;0].

P       = [0,0;1,0];
S1      = [.1;0];
Xbar    = log([1/Mss;Kss;Yss;mu;mu]);
G       = [A, B; zeros(2,3), P];
H       = [zeros(3,1); S1];
G       = [G (eye(5)-G)*Xbar;zeros(1,5),1];
H       = [H;0];
G0      = [A0, B0; zeros(2,3), P];
H0      = [zeros(3,1); S1];
G0      = [G0 (eye(5)-G0)*Xbar;zeros(1,5),1];
H0      = [H0;0];

nsim    = 1;
ny      = 25;
nimp    = 3;

imp     = [zeros(nimp-1,1);1;zeros(ny-nimp,1)];
%imp     = randn(ny,1);
[yt,xt] = dlsim(G,  H,  G,  H,  imp, [Xbar;1]);
[y0,x0] = dlsim(G0, H0, G0, H0, imp, [Xbar;1]);


Yt      = exp(yt(:,3));
Kt      = exp(xt(:,2));
Kft     = exp(yt(:,2));
Pt      = exp(yt(:,1));
Mut     = exp(xt(:,4));
mut     = xt(:,4);
mutl    = xt(:,5);
ytl     = xt(:,3);
ptl     = xt(:,1);
Xt      = (1+b*delta-sqrt(1+2*b*delta-2*b*(Kft./Kt-1+delta)))/b  .* Kt;
Ct      = Yt-Xt;
Lt      = (Yt./Kt.^alpha).^(1/(1-alpha));
kt      = xt(:,2);
kft     = yt(:,2);
ct      = log(Ct);
xt      = log(Xt);
yt      = log(Yt); 
pt      = log(Pt);

Y0      = exp(y0(:,3));
K0      = exp(x0(:,2));
Kf0     = exp(y0(:,2));
P0      = exp(y0(:,1));
X0      = Kf0-(1-delta)*K0;
C0      = Y0-X0;
L0      = (Y0./K0.^alpha).^(1/(1-alpha));
k0      = x0(:,2);
kf0     = y0(:,2);
c0      = log(C0);
x0      = log(X0);
y0      = log(Y0); 
p0      = log(P0);




return

bhld    = b;
a       = A(1,1);
b       = bhld*A(2,1);
c       = bhld*B(2,1);
d       = A(3,1);
e       = B(3,1);



[a-1/2-1/2*a^2-alpha*d*(1+a)/(1-alpha);
 1/2*a^2-1/2+delta*d*a/(1-alpha)-b+b*a;
 -1/2*a+1/2-delta*d/(1-alpha)-c-b;
 psi*(1-a^2)+phi*d*(1-a)+delta*d*a/(1-alpha)-b+b*a;
 psi*(a-1)+phi*(d+e)-delta*d/(1-alpha)-c-b]


