%---------------------------------------------------------------------
% 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];
P0      = (eye(4)-P)*Sbar;

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

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

%
% 7a. Steady state:
%
x     = (eye(4)-P)\P0;
z     = exp(x(1));
taul  = x(2);
taux  = x(3);
g     = exp(x(4));
betah = beta*(1+gz)^(-sigma);
kl    = ((1+taux)*(1-betah*(1-delta))/(betah*theta))^(1/(theta-1))*z;
yk    = (kl/z)^(theta-1);
xi1   = yk-(1+gz)*(1+gn)+1-delta;
xi2   = (1-taul)*(1-theta)*(kl)^theta*z^(1-theta)/psi;
xi3   = xi2/kl;
k     = (xi2+g)/(xi1+xi3);
c     = xi1*k-g;
l     = k/kl;
y     = yk*k;
x     = y-c-g;
r     = theta*y/k;
X0    = [log(k);log(z);taul;taux;log(g);1];
Y0    = [log(y);log(x);log(l);log(g)];
xbar  = X0;

%
% 7b. Approximations for l,y,c,x
%
philh  =-(psi*y*(1-theta)+(1-theta)*(1-taul)*y*(1-l)/l*theta+ ...
         (1-theta)*(1-taul)*y);
philk  = (psi*y*theta+psi*(1-delta)*k-(1-theta)*(1-taul)*y*(1-l)/l *theta)/philh;
philz  = (psi*y*(1-theta)-(1-theta)^2*(1-taul)*y*(1-l)/l)/philh;
phill  = ((1-theta)*(1-taul)*y*(1-l)/l *(1/(1-taul)))/philh;
philx  = 0;
philg  = (-psi*g)/philh;
philkp = (-psi*(1+gz)*(1+gn)*k)/philh;
phiyk  = theta+(1-theta)*philk;
phiyz  = (1-theta)*(1+philz);
phiyl  = (1-theta)*phill;
phiyx  = 0;
phiyg  = (1-theta)*philg;
phiykp = (1-theta)*philkp;
phixk  = -k/x*(1-delta);
phixkp = k/x*(1+gz)*(1+gn);
phixz  = 0;
phixl  = 0;
phixx  = 0;
phixg  = 0;
phick  = y/c*phiyk-x/c*phixk;
phicz  = y/c*phiyz;
phicl  = y/c*phiyl;
phicx  = 0;
phicg  = y/c*phiyg-g/c;
phickp = y/c*phiykp-x/c*phixkp;

%
% 7c. Compute gamma_k
%
coef1 = betah*r*(1-theta);
coef2 = -(1+taux)*psi*(1-sigma)*l/(1-l);
coef3 = -(1+taux)*sigma;

q      = roots([coef1*philkp+coef2*philkp+coef3*phickp, ...
                coef1*(philk-1)+coef2*(philk-philkp)+coef3*(phick-phickp), ...
               -coef2*philk-coef3*phick]);
i      = find(abs(q)<1);
gammak = q(i);

%
% 7d. Compute other gamma's
%
kap0  = coef2*[philz;phill;philx;philg]+ ...
        coef3*[phicz;phicl;phicx;phicg]+[0;0;1;0];
kap1  = coef2*philkp+coef3*phickp-(coef1+coef2)*(philk+philkp*gammak)- ...
        coef3*(phick+phickp*gammak)+coef1;
coef  = -coef1-coef2;
zet0  = -(coef1+coef2)*[philz;phill;philx;philg]- ...
        coef3*[phicz;phicl;phicx;phicg]-[coef1;0;betah*(1-delta);0];
zet1  =  coef*philkp-coef3*phickp;

gamma = -(kap1*eye(4)+zet1*P')\(kap0+P'*zet0);
gammaz= gamma(1);
gammal= gamma(2);
gammax= gamma(3);
gammag= gamma(4);
gamma0= (1-gammak)*log(k)-gammaz*log(z)-gammal*taul-gammax*taux-gammag*log(g);
gamma = [gammak;gamma;gamma0];

%
% 7e. 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)
%
A     = [ gammak, gammaz, gammal, gammax, gammag, 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,phiyx,phiyg]+phiykp*gamma(1:5)';
          [phixk,phixz,phixl,phixx,phixg]+phixkp*gamma(1:5)';
          [philk,philz,phill,philx,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,xbar);
