%RBC    Use Vaughan's method to solve the neoclassical growth
%       model with elastic labor, but without distortions:
%
%          max sum_t beta^t u(C[t]/N[t],L[t]/N[t])*N[t]
%
%          subj. to C[t] + X[t] = Y[t] = K[t]^theta (Z[t]H[t])^(1-theta)
%                   K[t+1] = (1-delta) K[t] + X[t]
%                   1 = H[t]+L[t]
%                   u(c,l) = {[c l^psi]^(1-alpha)-1}/(1-alpha)
%
%       with Z[t] = z[t](1+gz)^t and N[t]=(1+gn)^t.
%
%       First-order conditions:
%
%          u_c(C[t]/N[t],L[t]/N[t]) = mu[t]
%          u_l(C[t]/N[t],L[t]/N[t]) = mu[t] *(1-theta) Y[t]/H[t]
%          mu[t] = beta E_t mu[t+1](1-delta+theta Y[t+1]/K[t+1])
%
%       Detrended first-order conditions:
%
%          u_c[t] = beta*(1+gz)^(-alpha) E_t u_c[t+1]*(1-delta+y[t+1]/k[t+1])
%          u_c[t] = c[t]^(-alpha)*l[t]^(psi*(1-alpha)) 
%          psi c[t] = (1-theta) y[t]/h[t] *(1-h[t])
%          (1+gz)*(1+gn) k[t+1] = (1-delta)k[t]+x[t]
%          c[t]+x[t] = y[t] = k[t]^theta *(z[t]h[t]) ^(1-theta) 
%          ln z[t] = rho ln z[t-1] + eps[t], eps ~ N(0,sigma^2)
%
%          where c[t]=C[t]/(N[t]*(1+gz)^t), and similarly for y,k
%                h[t]=H[t]/N[t]
%       

%       Ellen McGrattan, 11/24/13
%       Revised, ERM, 12/5/13

%---------------------------------------------------------------------
% PARAMETERS
%    

alpha      = 1.0001;
beta       = 0.995;
delta      = 0.008;
gn         = 0.0025;
gz         = 0.005;
psi        = 1.365;
theta      = 0.33;
zss        = 1.00;
taucss     = 0.065;
taudss     = 0.144;
tauhss     = 0.382;
taukss     = 0.003;
taupss     = 0.330;
tauxss     = 0.00;
gss        = 0.001;

P          = zeros(8);
Q          = P;
P(1,1)     = 0.9334;
Q(1,1)     = 1;

T          = 120;
%
% If iplot==1, load shocks in eps.dat and plot hours
%
iplot      = 1;
if iplot==1;
  load eps.dat
  nsim     = 1;
else  
  Q(1,1)   = .01;
  nsim     = 100;
end;
param      = [alpha;beta;delta;gn;gz;psi;theta];

%---------------------------------------------------------------------
% COMPUTE EQUILIBRIUM 
%    
%     X_{t+1} = A X_t + B S_t
%     Z_t     = C X_t + D S_t
%
%     where X_t = ln k_t           (all deviations from steady state)
%           Z_t = [ln k_{t+1},c_t]
%           S_t = ln z_t
%---------------------------------------------------------------------
%
% Step 1. Steady state 
%

thet1      = 1-theta;
betg       = beta*(1+gz)^(-alpha);
Sbar       = [log(zss);taucss;taudss;tauhss;taukss;taupss;tauxss;log(gss)];
rss        = ( (1+tauxss)-betg*(1+(1-taupss)*(-delta-taukss)+(1-delta)* ...
             tauxss) )/( betg*(1-taupss) );
kdivy      = theta/rss;
xdivk      = (1+gz)*(1+gn)-1+delta;
kdivh      = kdivy^(1/thet1)*zss;
cdivl      = (1-theta)*(1-tauhss)*kdivh^theta*zss^thet1/(psi*(1+taucss));
kss        = (cdivl+gss)/(1/kdivy-xdivk+cdivl/kdivh);
yss        = kss/kdivy;
hss        = kss/kdivh;
xss        = xdivk*kss;
css        = yss-xss-gss;

%
% Step 2. Linearize first-order conditions
%

v          = [log(kss);log(kss);log(kss);Sbar;Sbar;
              log(css);log(css);log(hss);log(hss)];
del        = max(abs(v)*1e-4,1e-8);
dR         = zeros(5,23);
for i=1:23;
  vp       = v;
  vm       = v;
  vp(i)    = vp(i)+del(i);
  vm(i)    = vm(i)-del(i);
  dR(:,i)  = (res_rbc(vp,param)-res_rbc(vm,param))/(2*del(i));
end;
i1         = 1:4;
i2         = 5;
j1         = 1:19;
j2         = 20:23;
Coef       = -inv(dR(i1,j2))*dR(i1,j1);
dRnew      = dR(i2,j1)+dR(i2,j2)* Coef;

%
% Step 3. Find the deterministic part of solution
%
%         Residuals are A1*[X_{t+1}] + A2*[X_t]  + shock terms = 0
%                          [Z_{t+1}]      [Z_t]
%
%                           X_t  = k_t, Z_t  = [k_{t+1}]
%

nX         = 1;
nZ         = 1;
nXZ        = 2;

a1         = dRnew(1);
a2         = dRnew(2);
a3         = dRnew(3);
a4         = dRnew(4:11);
a5         = dRnew(12:19);


A1         = [1,0;
              0,a3];
A2         = [0,-1;
              a1,a2];

[V,L]      = eig(A2,-A1);
L          = diag(L);
[D,I]      = sort(abs(L));
if (~( (D(nX) <= 1) & (D(nX+1) > 1) ) )
  i      = sum(D<=1);
  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:nXZ,I(1:nX));
W          = real([V11*diag(L(I(1:nX)))/V11;
                   V21/V11]);
A          = W(1);
C          = W(2);

%
% Step 4. Find the stochastic part of the solution (ie., matrices B,D)
%
%            E_t[a1* ln k[t] + a2* ln k[t+1] + a3* ln k[t+2] 
%                            + a4 * S[t]  + a5* S[t+1] ] = 0
%
%           where S =[ln(z),tauc,taud,tauh,tauk,taup,taux,ln(g)]
%
%           ==>  (a1+a2*A+a3*A*A)*k
%              + (a2*B+a3*A*B+a3*B*P+a4+a5*P)*S = 0
%
%           Note: need vec(A+B)=vec(A)+vec(B), vec(ABC)=(C' kron A)*vec(B)
%

resk       = a1+a2*A+a3*A*A;

if (abs(resk)>1e-8);
  disp('Check deterministic part of the solution')
end;

I          = eye(8);
tem        =-([kron(I,a2+a3*A)+kron(P',a3)])\vec(a4+a5*P);
B          = zeros(1,8);
B(:)       = tem;
D          = B;

ii         = [1,3];   % indices for c[t],h[t]
C2         = Coef(ii,1)+Coef(ii,2)*A;
D2         = Coef(ii,4:11)+Coef(ii,2)*B;

%
% Step 5. Simulate time series for x=[k,z,1],u=eps,y=[k']
%
%           x[n+1] = [A,   B, A0] x[n] + [ 0 ] u[n]
%                    [0,   P, P0]        [ Q ]
%                    [0,   0,  1]        [ 0 ]
%
%                       (G)               (H)
%
%           y[n]   = [C, D, C0] x[n] 
%
%                       (J)
%

A0         = (1-A)*log(kss)-B*Sbar;
P0         = (I-P)*Sbar;
C0         = log([css;hss])-C2*log(kss)-D2*Sbar;
G          = [A,B,A0;zeros(8,1),P,P0;zeros(1,9),1];
H          = [zeros(1,8);Q;zeros(1,8)];
J          = [C2,D2,C0;A,B,A0];
X0         = [log(kss);Sbar;1];
%rng(5)
randn('seed',100)
vars       = zeros(nsim,7);
for s=1:nsim
  epst     = randn(T,8);
  if iplot==1;
    espt     = zeros(T,8);
    epst(:,1)= eps;
  end;
  [Z,X]    = dlsim(G,H,J,zeros(3,8),epst,X0);
  ct       = exp(Z(1:T,1));
  ht       = exp(Z(1:T,2));
  ktp      = exp(Z(1:T,3));
  kt       = exp(X(1:T,1));
  zt       = exp(X(1:T,2));
  tauct    = X(1:T,3);
  taudt    = X(1:T,4);
  tauht    = X(1:T,5);
  taukt    = X(1:T,6);
  taupt    = X(1:T,7);
  tauxt    = X(1:T,8);
  gt       = exp(X(1:T,9));
  xt       = (1+gz)*(1+gn)*ktp-(1-delta)*kt;
  yt       = ct+xt+gt;
  lwt      = (psi*(1+tauct).*ct./(1-ht))./(thet1*(1-tauht).*yt./ht);
  lyt      = log(yt)-hptrend(log(yt),1600);
  lct      = log(ct)-hptrend(log(ct),1600);
  lxt      = log(xt)-hptrend(log(xt),1600);
  lkt      = log(kt)-hptrend(log(kt),1600);
  lht      = log(ht)-hptrend(log(ht),1600);
  llwt     = log(lwt)-hptrend(log(lwt),1600);
  clwy     = corrcoef([llwt,lyt]);

  vars(s,:)= [[std(lyt), std(lct), std(lxt), std(lkt), std(lht),  ...
               std(llwt)]*100, clwy(1,2)];
end;
if nsim>1;
  tab      = [mean(vars)',std(vars)'];
else
  tab      = [vars',0*vars'];
end;

%
% Step 6. Print results
%
clc
disp('    RESULTS FOR THE RBC MODEL WITH ELASTIC LABOR ')
disp('    -------------------------------------------- ')
disp(' Inputs: ')
disp(' ')
disp(sprintf('    Utility curvature, alpha        %7.3f',alpha))
disp(sprintf('    Discount factor, beta           %7.3f',beta))
disp(sprintf('    Depreciation rate, delta        %7.3f',delta))
disp(sprintf('    Population growth, gn           %7.3f',gn))
disp(sprintf('    Technological growth, gz        %7.3f',gz))
disp(sprintf('    Shock persistence, rho          %7.3f',P(1,1)))
disp(sprintf('    Shock serial correlation, sigma %7.3f',Q(1,1)))
disp(sprintf('    Capital share, theta            %7.3f',theta))
disp(' ')
disp(' Outputs: ')
disp(' ')
disp('    Optimal decision rule: ')
disp(' ')
disp(sprintf('    ln k[t+1] = %6.2f +%6.3f ln k[t] +%6.3f ln z[t]',[A0,A,B(1)]))
disp(sprintf('                       %7.3f tauc[t] +%6.3f taud[t]',[B(2:3)]))
disp(sprintf('                       %7.3f tauh[t]  %6.3f tauk[t]',[B(4:5)]))
disp(sprintf('                       %7.3f taup[t]  %6.3f taux[t]  %6.3f ln g[t]',[B(6:8)]))
disp(' ')
disp('    Time series statistics (logged and detrended variables): ')
disp(' ')
disp(sprintf('    std(gdp)           %7.3f (%7.3f)',tab(1,:)))
disp(sprintf('    std(consumption)   %7.3f (%7.3f)',tab(2,:)))
disp(sprintf('    std(investment)    %7.3f (%7.3f)',tab(3,:)))
disp(sprintf('    std(capital)       %7.3f (%7.3f)',tab(4,:)))
disp(sprintf('    std(hours)         %7.3f (%7.3f)',tab(5,:)))
disp(sprintf('    std(labor wedge)   %7.3f (%7.3f)',tab(6,:)))
disp(sprintf('    cor(wedge, gdp)    %7.3f (%7.3f)',tab(7,:)))

%
% Step 7. Plot hours
%
if (iplot==1);
  load ../MLEdata/ldhpctot.dat
  t        = ptime(1985.2,119,4);
  hact     = ldhpctot(149:267);
  plot(t,exp(hact-mean(hact))*hss,t,ht(2:120))
  title('Per Capita Hours of Work, 1985:1-2014:4')
  legend('Actual','Predicted')
end;

