%---------------------------------------------------------------------
% 1. Default parameters for 
%    gn, gz, beta, delta, psi, sigma, theta, ghat, a_0, a, b, D, R
%
gn      = (1.015)^(1/4)-1;
gz      = (1.016)^(1/4)-1;
beta    = .9722^(1/4);
delta   = 1-(1-.0464)^(1/4);
psi     = 2.24;
sigma   = 1.000001;
theta   = .35;

Sbar = [ ...
  -1.2218e-01
   2.3456e-01
   2.1814e-01
  -1.7104e+00];
P = [ ...
   7.6890e-01   4.7057e-02   4.3179e-01  -4.1890e-02
  -2.3290e-02   9.9352e-01   4.1715e-02  -4.5134e-03
  -1.0129e-01   2.5306e-02   1.1832e+00  -1.7991e-02
  -3.0636e-02   4.2339e-02  -1.6011e-03   9.9730e-01];
Q = [ ...
   1.0814e-02            0            0            0
  -2.5674e-03  -6.2297e-03            0            0
   3.7148e-03  -8.8764e-04   1.9564e-03            0
  -5.0087e-03  -5.0538e-03   1.4790e-02  -2.0158e-06];
Q       = chol(Q*Q')';
P0      = (eye(4)-P)*Sbar;

param   = [gn;gz;beta;delta;psi;sigma;theta];

%---------------------------------------------------------------------
% 7.  Compute equilibrium
%---------------------------------------------------------------------

%
% 7a. Steady state:
%
tem        = (eye(4)-P)\P0;
zs         = exp(tem(1));
tauls      = tem(2);
tauxs      = tem(3);
gs         = exp(tem(4));
beth       = beta*(1+gz)^(-sigma);
kls        = ((1+tauxs)*(1-beth*(1-delta))/(beth*theta))^(1/(theta-1))*zs;
A          = (zs/kls)^(1-theta)-(1+gz)*(1+gn)+1-delta;
B          = (1-tauls)*(1-theta)*kls^theta*zs^(1-theta)/psi;
ks         = (B+gs)/(A+B/kls);
cs         = A*ks-gs;
ls         = ks/kls;
ys         = ks^theta*(zs*ls)^(1-theta);
xs         = ys-cs-gs;
X0         = [log(ks);log(zs);tauls;tauxs;log(gs);1];
Y0         = [log(ys);log(xs);log(ls);log(gs)];
                                                                                
%
% 7b. Call subroutine with residuals:
%
                                                                                
Z          = [log(ks);log(ks);log(ks);log(zs);log(zs);tauls;tauls;
              tauxs;tauxs;log(gs);log(gs)];
del        = max(abs(Z)*1e-5,1e-8);
for i=1:11;
  Zp       = Z;
  Zm       = Z;
  Zp(i)    = Z(i)+del(i);
  Zm(i)    = Z(i)-del(i);
  dR(i,1)  = (res_wedge(Zp,param)-res_wedge(Zm,param))/(2*del(i));
end;
                                                                                
%
% 7c. Solution:  log k[t+1] = gamma0 + gammak* log k[t] + gamma* S[t]
%
                                                                                
a0         = dR(1);
a1         = dR(2);
a2         = dR(3);
b0         = dR(4:2:11)';
b1         = dR(5:2:11)';
tem        = roots([a0,a1,a2]);
gammak     = tem(find(abs(tem)<1));
gamma      = -((a0*gammak+a1)*eye(4)+a0*P')\(b0*P+b1)';
gamma0     = (1-gammak)*log(ks)-gamma'*[log(zs);tauls;tauxs;log(gs)];
Gamma      = [gammak;gamma;gamma0];
                                                                                
%
% 7c. State-space system:   X[t+1] = A X[t] + B eps[t+1]
%                           Y[t]   = C X[t] + ome[t]
%                           ome[t] = D ome[t-1] + eta[t],  E eta eta' ~ N(0,R)
%
philh      =-(psi*ys*(1-theta)+(1-theta)*(1-tauls)*ys*(1-ls)/ls*theta+ ...
             (1-theta)*(1-tauls)*ys);
philk      = (psi*ys*theta+psi*(1-delta)*ks- ...
             (1-theta)*(1-tauls)*ys*(1-ls)/ls *theta)/philh;
philz      = (psi*ys*(1-theta)-(1-theta)^2*(1-tauls)*ys*(1-ls)/ls)/philh;
phill      = ((1-theta)*(1-tauls)*ys*(1-ls)/ls *(1/(1-tauls)))/philh;
philg      = (-psi*gs)/philh;
philkp     = (-psi*(1+gz)*(1+gn)*ks)/philh;
phiyk      = theta+(1-theta)*philk;
phiyz      = (1-theta)*(1+philz);
phiyl      = (1-theta)*phill;
phiyg      = (1-theta)*philg;
phiykp     = (1-theta)*philkp;
phixk      = -ks/xs*(1-delta);
phixkp     = ks/xs*(1+gz)*(1+gn);

                                                                                
A     = [ gammak,                         gamma', gamma0;
         [0;0;0;0],            P,                     P0;
               0,      0,      0,      0,      0,      1];
B     = [ 0,0,0,0;
            Q;
          0,0,0,0];
C     = [ [phiyk,phiyz,phiyl,    0,phiyg]+phiykp*Gamma(1:5)';
          [phixk,    0,    0,    0,    0]+phixkp*Gamma(1:5)';
          [philk,philz,phill,    0,philg]+philkp*Gamma(1:5)';
          0,0,0,0,1];
phi0  = Y0-C*X0(1:5);
C     = [C,phi0];
shks  = randn(180,4);
[Yt,Xt] =  dlsim(A,B,C,zeros(4),shks,X0);
