%RBC_ED   Use Vaughan's method to solve the neoclassical growth
%         model with elastic labor and distortions:
%
%            max sum_t beta^t u(C[t]/N[t],L[t]/N[t])*N[t]
%
%            subj. to C[t] + X[t] = r[t] K[t] + W[t] H[t] - Taxes + Transfers
%                     Taxes = tauc[t]*C[t] + taux[t]*X[t] + 
%                           + tauh[t]*w[t]*H[t]
%                           + tauk[t]*K[t]
%                           + taup[t](r[t]-delta-tauk[t])K[t]
%                           + taud[t]{r[t]K[t]-X[t]
%                                 - taup[t](r[t]-delta-tauk[t])K[t]
%                                 - taux[t] X[t]-tauk K[t]}
%                     K[t+1] = (1-delta) K[t] + X[t]
%                     1 = H[t]+L[t]
%
%                     Firms: Y[t] = F(K[t],Z[t]H[t])
%                     Resources: Y[t] = C[t]+X[t]+G[t]
%                     Utility: 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] (1+tauc[t])
%            u_l(C[t]/N[t],L[t]/N[t]) = mu[t] *W[t]*(1-tauh[t])
%
%            mu[t]*(1-taud[t])(1+taux[t]) = beta E_t (1-taud[t+1])*mu[t+1]*
%                 (1+(1-taup[t+1])(r[t+1]-delta-tauk[t+1])+(1-delta)taux[t+1])
%
%         Detrended first-order conditions:
%
%            u_c[t](1-taud[t])(1+taux[t]) = beta*(1+gz)^(-alpha) 
%               * E_t u_c[t+1]*(1+(1-taup[t+1])(r[t+1]-delta-tauk[t+1])
%                                +(1-delta)taux[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-tauh[t])/(1+tauc[t])
%            (1+gz)*(1+gn) k[t+1] = (1-delta)k[t]+x[t]
%            c[t]+x[t]+g[t] = y[t] = k[t]^theta *(z[t]h[t]) ^(1-theta) 
%            S[t] = S0+ P S[t-1] + Q eps[t], eps ~ N(0,Sigma), 
%
%            where c[t]=C[t]/(N[t]*(1+gz)^t), and similarly for y,k
%                  h[t]=H[t]/N[t]
%                  S = [ln z,tauc,taud,tauh,tauk,taup,taux]
%       

%       Ellen McGrattan, 11/24/13
%       Revised, ERM, 5/18/14

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

alpha      = 1.50;
beta       = 0.98;
delta      = 0.05;
gn         = 0.01;
gz         = 0.02;
psi        = 1.75;
theta      = 0.33;
zss        = 1.0;
taucss     = 0.06;
taudss     = 0.14;
tauhss     = 0.40;
taukss     = 0.00;
taupss     = 0.30;  % Compare hss with taup=0.2 case
tauxss     = 0.00;
gss        = 0.1;

P          = zeros(8);
Q          = P;
P(1,1)     = 0.99;
Q(1,1)     = 0.015;
T          = 250;
nsim       = 100;

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,tauc_t,taud_t,tauh_t,tauk_t,taup_t,taux_t,ln g_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(css);log(kss);log(css);Sbar;Sbar];
del        = max(abs(v)*1e-4,1e-8);
dR         = zeros(2,21);
for i=1:21;
  vp       = v;
  vm       = v;
  vp(i)    = vp(i)+del(i);
  vm(i)    = vm(i)-del(i);
  dR(:,i)  = (res_rbc_ed(vp,param)-res_rbc_ed(vm,param))/(2*del(i));
end;


%
% Step 2. Find 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},c_t]
%

a1         = dR(:,1);
a2         = dR(:,2);
a3         = dR(:,3);
a4         = dR(:,4);
a5         = dR(:,5);
a6         = dR(:,6:13);
a7         = dR(:,14:21);

nX         = 1;
nZ         = 2; 
nXZ        = 3;
A1         = [1,0,0;
              zeros(2,1),a4,a5];
A2         = [0,-1,0;
              a1,a2,a3];

[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:3);


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

resk       = a1+a2*A+a3*C(2)+a4*A*A+a5*C(2)*A;

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

I          = eye(8);
tem        =-[kron(I,a2+a4*A+a5*C(2))+kron(P',a4), ...
              kron(I,a3)+kron(P',a5)]\vec(a6+a7*P);
B          = tem(1:8)';
D          = [B;tem(9:16)'];

%
% Step 4. Simulate time series for x=[k,S,1],u=eps,y=[k',c]
%
%           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,[A0;log(css)-C(2)*log(kss)]] x[n] 
%
%                       (J)
%

A0         = (1-A)*log(kss)-B*Sbar;
P0         = (eye(8)-P)*Sbar;
G          = [A,B,A0;zeros(8,1),P,P0;zeros(1,9),1];
H          = [zeros(1,8);Q;zeros(1,8)];
J          = [C,D,[A0;log(css)-C(2)*log(kss)-D(2,:)*Sbar]];
X0         = [log(kss);Sbar;1];
%rng(5)
randn('seed',100)
vars       = zeros(nsim,6);
o          = ones(T,1);
for s=1:nsim
  eps      = randn(T,8);
  [Z,X]    = dlsim(G,H,J,zeros(2,8),eps,X0);
  ktp      = exp(Z(:,1));
  ct       = exp(Z(:,2));
  kt       = exp(X(:,1));
  zt       = exp(X(:,2));
  tauct    = X(:,3);
  taudt    = X(:,4);
  tauht    = X(:,5);
  taukt    = X(:,6);
  taupt    = X(:,7);
  tauxt    = X(:,8);
  gt       = exp(X(:,9));
  xt       = (1+gz)*(1+gn)*ktp-(1-delta)*kt;
  yt       = ct+xt+gt;
  ht0      = o./(1+psi*(1+tauct).*ct./(thet1*(1-tauht).*yt));
  ht       = (yt./(kt.^theta)).^(1/thet1) ./zt;
  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);

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


%
% Step 5. Print results
%
clc
disp('    RESULTS FOR THE RBC MODEL WITH ELASTIC LABOR & DISTORTIONS ')
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('    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(output)        %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,:)))


