close all;
addpath mfiles_jh;
var c_t i_t  y_t u_t k_t n_t lam_t x_t mu_t eta_t rk_t a_t z_t a1_t z1_t na1_t  nz1_t ;
varexo e_a1  e_z1 ne_a1  ne_z1 ;%e_a2 ne_a2  e_z2 ne_z2;

parameters beta theta alpha gamma phi uParam omega delta sigma rho_a1 rho_a2 rho_z1 rho_z2 psi;

%
% Parameters:
%

beta = 0.985; % subjective discount factor
theta = 1.4; % labor elasticity
alpha = 1-0.64; % labor share in production
gamma = 0.001; % The KEY Parameter
phi = 1.3; % Adjustment cost param
uParam = 1/0.15; %Util Param.
psi = 2.3; %CHECK CHECK
omega = 1 + uParam;
delta=0.025;
sigma=1;

rho_a1=0;
rho_a2=0;
rho_z1=0;
rho_z2=0;

nkratio=((1/beta-1+delta)/alpha)^(1/(1-alpha));
ckratio=(nkratio^(1-alpha)-delta);
cnratio=ckratio/(nkratio);
Para10=((1-alpha)*nkratio^(-alpha));
Para20=(theta*cnratio-Para10*gamma/((beta*(1-gamma)-1)))*psi;
N_ss=(Para10/Para20)^(1/theta);
%%#N_ss=(Para3/psi)^(1/theta);
%psi=alpha*nkratio^(alpha-1)/((theta*cnratio+((gamma/(beta*(1-gamma)-1))*alpha*nkratio^(alpha-1)))*N_ss^theta);
C_ss=N_ss*cnratio;
K_ss=N_ss/nkratio;
I_ss=K_ss*delta;
Y_ss=K_ss^alpha*N_ss^(1-alpha);
X_ss=C_ss;
Para1=(C_ss-psi*N_ss^theta*X_ss)^(-sigma);
Para2=(C_ss-psi*N_ss^theta*X_ss);
Mu_ss=Para1*psi*N_ss^theta/(beta*(1-gamma)-1);
Lam_ss=Para1+gamma*Mu_ss;

%Check the Computation of SteadyState
C=C_ss;N=N_ss;X=X_ss;I=I_ss; K=K_ss; Mu=Mu_ss;Lam=Lam_ss;U=1;
funval(1)=(C-psi*N^theta*X)^(-sigma)+gamma*Mu-Lam;
funval(2)=(C-psi*N^theta*X)^(-sigma)*theta*psi*N^(theta-1)*X-(1-alpha)*Lam*N^(-alpha)*(U*K)^(alpha);
funval(3)=(C-psi*N^theta*X)^(-sigma)*psi*N^(theta)+Mu-beta*(1-gamma)*Mu;
funval(4)=1-beta*((1-delta)+alpha*U^alpha*K^(alpha-1)*N^(1-alpha));
funval(5)=Y_ss-C-I;
if sum(abs(funval)>1e-14)>0;disp(funval);error('Check the SteadyState');end;

model(linear);
#scalar=10000;
#nkratio=(((1/beta-1+delta)/alpha)^(1/(1-alpha)));
#ckratio=((nkratio)^(1-alpha)-delta);
#cnratio=(ckratio)/(nkratio);
#Para10=((1-alpha)*nkratio^(-alpha));
#Para20=(theta*cnratio-Para10*gamma/((beta*(1-gamma)-1)))*psi;
#N_ss=(Para10/Para20)^(1/theta);
//%%#N_ss=(Para3/psi)^(1/theta);
//%psi=alpha*nkratio^(alpha-1)/((theta*cnratio+((gamma/(beta*(1-gamma)-1))*alpha*nkratio^(alpha-1)))*N_ss^theta);
#C_ss=N_ss*cnratio;
#K_ss=N_ss/nkratio;
#I_ss=K_ss*delta;
#Y_ss=K_ss^alpha*N_ss^(1-alpha);
#X_ss=C_ss;
#Para1=(C_ss-psi*N_ss^theta*X_ss)^(-sigma);
#Para2=(C_ss-psi*N_ss^theta*X_ss);
#Mu_ss=Para1*psi*N_ss^theta/(beta*(1-gamma)-1);
#Lam_ss=Para1+gamma*Mu_ss;
#rkss=alpha*nkratio^(1-alpha);
//1. National Income
//y_t=(alpha)*(u_t+k_t(-1)-1/(1-alpha)*a_t-1/(1-alpha)*z_t)+(1-alpha)*n_t;
y_t=alpha*u_t+alpha*k_t(-1)+(1-alpha)*n_t-alpha/(1-alpha)*a_t-alpha/(1-alpha)*z_t;

//2. Capital Accumulation
//k_t=(1-delta)*(k_t(-1)-1/alpha*a_t-1/alpha*z_t)+delta*i_t-rkss*u_t;
k_t=(1-delta)*k_t(-1)-(1-delta)/(1-alpha)*a_t-(1-delta)/(1-alpha)*z_t+delta*i_t-rkss*u_t;

//3. F.O.C on conumption;
Lam_ss*lam_t=+gamma*Mu_ss*mu_t+gamma*Mu_ss*(gamma-1)*(c_t-x_t(-1)+1/(1-alpha)*a_t+alpha/(1-alpha)*z_t)-sigma*Para1/Para2*(C_ss*c_t-psi*N_ss^(theta)*X_ss*(theta*n_t+x_t));

//4. F.O.C on habit;
(beta*(1-gamma)-1)*theta*n_t-sigma*(beta*(1-gamma)-1)/Para2*(C_ss*c_t-psi*N_ss^theta*X_ss*(theta*n_t+x_t))+1*mu_t=beta*(1-gamma)*mu_t(+1)+beta*(1-gamma)*gamma*(c_t(+1)-x_t)+beta*(1-gamma)*(gamma-sigma)/(1-alpha)*(a_t(+1)+alpha*z_t(+1));

//5. F.O.C on labor;
(theta-1)*n_t+x_t-sigma/Para2*C_ss*c_t+sigma/Para2*psi*N_ss^(theta)*X_ss*(theta*n_t+x_t)=lam_t+(alpha)*u_t+alpha*k_t(-1)-alpha/(1-alpha)*a_t-alpha/(1-alpha)*z_t-alpha*n_t;

//6. F.O.C on utilization;
  uParam*u_t+eta_t=rk_t;

//7. Asset Pricing;
//eta_t=beta*(alpha)*nkratio^(1-alpha)*(lam_t(+1)-sigma/alpha*(a_t+(1-alpha)*z_t)+a_t(+1)+alpha*(n_t(+1)-k_t-u_t(+1)))+beta*(1-delta)*(eta_t(+1)-sigma/alpha*(a_t+(1-alpha)*z_t)-z_t(+1));
rk_t=a_t+z_t+(1-alpha)*(n_t)-(1-alpha)*u_t-(1-alpha)*k_t(-1);

lam_t+eta_t=lam_t(+1)+beta*(1-delta)*(eta_t(+1))-sigma/(1-alpha)*(a_t(+1)+alpha*z_t(+1))-z_t(+1)+(1-beta*(1-delta))*rk_t(+1);

//8. F.O.C on I_t
eta_t=phi*(i_t-i_t(-1)+1/(1-alpha)*a_t+1/(1-alpha)*z_t)-beta*phi*(i_t*(+1)-i_t+1/(1-alpha)*a_t(+1)+1/(1-alpha)*z_t(+1));

//9. accumulation of x_t;
x_t=gamma*c_t+(1-gamma)*x_t-(1-gamma)/(1-alpha)*a_t-alpha*(1-gamma)/(1-alpha)*z_t;

Y_ss*y_t=C_ss*c_t+I_ss*i_t;
//exp(y_t)=exp(c_t)+exp(i_t);
//*********************Shock Series***********************************//
a1_t=rho_a1*a1_t(-1)+e_a1+na1_t(-4);
//a2_t=rho_a2*a2_t(-1)+e_a2;
z1_t=rho_z1*z1_t(-1)+e_z1+nz1_t(-4);
//z2_t=rho_z2*z2_t(-1)+e_z2;
na1_t=ne_a1;
//na2_t=ne_a2;
nz1_t=ne_z1;
//nz2_t=ne_z2;

//a_t=a1_t+a2_t-a2_t(-1)+na1_t(-4)+na2_t(-4)-na2_t(-5);
a_t=a1_t;
//z_t=z1_t+z2_t-z2_t(-1)+nz1_t(-4)+nz2_t(-4)-nz2_t(-5);
z_t=z1_t;
end;

check;

shocks;
var e_a1;stderr 1;
//var e_a2;stderr 1;
var e_z1;stderr 1;
//var e_z2;stderr 1;
var ne_a1;stderr 1;
//var ne_a2;stderr 1;
var ne_z1;stderr 1;
//var ne_z2;stderr 1;
end;


nsteps=10;
nperiod=4;
stoch_simul(periods=1000,irf=10)c_t i_t n_t y_t k_t z_t a_t;

Z_t_e_a1=cumsum(z_t_e_a1);
//Z_t_e_a2=cumsum(z_t_e_a2);
Z_t_e_z1=cumsum(z_t_e_z1);
//Z_t_e_z2=cumsum(z_t_e_z2);
Z_t_ne_a1=cumsum(z_t_ne_a1);
//Z_t_ne_a2=cumsum(z_t_ne_a2);
Z_t_ne_z1=cumsum(z_t_ne_z1);
//Z_t_ne_z2=cumsum(z_t_ne_z2);

A_t_e_a1=cumsum(a_t_e_a1);
//A_t_e_a2=cumsum(a_t_e_a2);
A_t_e_z1=cumsum(a_t_e_z1);
//A_t_e_z2=cumsum(a_t_e_z2);
A_t_ne_a1=cumsum(a_t_ne_a1);
//A_t_ne_a2=cumsum(a_t_ne_a2);
A_t_ne_z1=cumsum(a_t_ne_z1);
//A_t_ne_z2=cumsum(a_t_ne_z2);

 figure;
 subplot(2,1,1);
 hold on;
 plot(n_t_e_a1(1:nsteps))
 plot(n_t_ne_a1(1:nsteps),'r--');
 plot(nperiod,n_t_ne_a1(nperiod),'g*');
 title('Permanent N-shock');
 //subplot(2,2,2);
 //hold on;
 //plot(n_t_e_a2(1:nsteps));
 //plot(n_t_ne_a2(1:nsteps),'r--');
 //plot(nperiod,n_t_ne_a2(nperiod),'g*');
//title('Temporary N-shock');
 subplot(2,1,2);
 hold on;
 plot(n_t_e_z1(1:nsteps));
 plot(n_t_ne_z1(1:nsteps),'r--');
plot(nperiod,n_t_ne_z1(nperiod),'g*');
 title('Permanent I-shock');
//subplot(2,2,4);
// hold on;
// plot(n_t_e_z2(1:nsteps));
// plot(n_t_ne_z2(1:nsteps),'r--');
//plot(nperiod,n_t_ne_z2(nperiod),'g*');
// title('Temporary I-shock');
legend('Tech Shock','News Shock (Two quarters advance)');
suptitle('Response of Labor');

figure;
 subplot(2,1,1);
 hold on;
 tmp0=c_t_e_a1+1/(1-alpha)*A_t_e_a1+(alpha)/(1-alpha)*Z_t_e_a1;
 plot(tmp0(1:nsteps));

 tmp=c_t_ne_a1+1/(1-alpha)*A_t_ne_a1+alpha/(1-alpha)*Z_t_ne_a1;
 plot(tmp(1:nsteps),'r--');

 plot(nperiod,tmp(nperiod),'g*');
 title('Permanent N-shock');
 C_t_e_a1=tmp0;
 C_t_ne_a1=tmp; 

 //subplot(2,2,2);
 //hold on;
 //tmp0=c_t_e_a2+1/alpha*A_t_e_a2+(1-alpha)/alpha*Z_t_e_a2;
 //plot(tmp0(1:nsteps));

 //tmp=c_t_ne_a2+1/alpha*A_t_ne_a2+(1-alpha)/alpha*Z_t_ne_a2;
 //plot(tmp(1:nsteps),'r--');
 //plot(nperiod,tmp(nperiod),'g*');
 //title('Temporary N-shock');
 //C_t_e_a2=tmp0;
 //C_t_ne_a2=tmp; 


 subplot(2,1,2);
 hold on;
 tmp0=(c_t_e_z1+1/(1-alpha)*A_t_e_z1+(alpha)/(1-alpha)*Z_t_e_z1);
 plot(tmp0(1:nsteps));

 tmp=c_t_ne_z1+1/(1-alpha)*A_t_ne_z1+(alpha)/(1-alpha)*Z_t_ne_z1;
 plot(tmp(1:nsteps),'r--');

 plot(nperiod,tmp(nperiod),'g*');
 title('Permanent I-shock');
 C_t_e_z1=tmp0;
 C_t_ne_z1=tmp; 

 //subplot(2,2,4);
 //hold on;
 //tmp0=(c_t_e_z2+1/alpha*A_t_e_z2+(1-alpha)/alpha*Z_t_e_z2);
 //plot(tmp0(1:nsteps));
//% plot(c_t_ne_z2+1/alpha*A_t_ne_z2+(1-alpha)/alpha*Z_t_ne_z2,'r--');
//tmp=c_t_ne_z2+1/alpha*A_t_ne_z2+(1-alpha)/alpha*Z_t_ne_z2;
//plot(tmp(1:nsteps),'r--');
//plot(nperiod,tmp(nperiod),'g*');
//legend('Tech Shock','News Shock (Four quarters advance)');
//C_t_e_z2=tmp0;
//C_t_ne_z2=tmp; 

// title('Temporary I-shock');
suptitle('Response of Consumption');

 figure;
 subplot(2,2,1);
 hold on;
 tmp0=(i_t_e_a1+1/(1-alpha)*A_t_e_a1+(1)/(1-alpha)*Z_t_e_a1);
 plot(tmp0(1:nsteps));
 %plot(i_t_ne_a1+1/(1-alpha)*A_t_ne_a1+(1-alpha)/(1-alpha)*Z_t_ne_a1,'r--');
 tmp=i_t_ne_a1+1/(1-alpha)*A_t_ne_a1+1/(1-alpha)*Z_t_ne_a1;
 plot(tmp(1:nsteps),'r--');
 plot(nperiod,tmp(nperiod),'g*');
 title('Permanent N-shock');
I_t_e_a1=tmp0;
I_t_ne_a1=tmp; 

 //subplot(2,2,2);
 //hold on;
 //tmp0=(i_t_e_a2+1/(1-alpha)*A_t_e_a2+(1)/(1-alpha)*Z_t_e_a2);
 //plot(tmp0(1:nsteps));
 //%plot(i_t_ne_a2+1/(1-alpha)*A_t_ne_a2+(1)/(1-alpha)*Z_t_ne_a2,'r--');
 //tmp=i_t_ne_a2+1/(1-alpha)*A_t_ne_a2+(1)/(1-alpha)*Z_t_ne_a2;
 //plot(tmp(1:nsteps),'r--');
 //plot(nperiod,tmp(nperiod),'g*');
 //title('Temporary N-shock');

 //I_t_e_a2=tmp0;
 //I_t_ne_a2=tmp; 
 subplot(2,2,3);
 hold on;
 tmp0=(i_t_e_z1+1/(1-alpha)*A_t_e_z1+(1)/(1-alpha)*Z_t_e_z1);
 plot(tmp0(1:nsteps));
% plot(i_t_ne_z1+1/(1-alpha)*A_t_ne_z1+(1)/(1-alpha)*Z_t_ne_z1,'r--');
 tmp=i_t_ne_z1+1/(1-alpha)*A_t_ne_z1+(1)/(1-alpha)*Z_t_ne_z1;
 plot(tmp(1:nsteps),'r--');

 plot(nperiod,tmp(nperiod),'g*');
 title('Permanent I-shock');

 I_t_e_z1=tmp0;
 I_t_ne_z1=tmp; 

 //subplot(2,2,4);
 //hold on;
 //tmp0=(i_t_e_z2+1/(1-alpha)*A_t_e_z2+1/(1-alpha)*Z_t_e_z2);
 //plot(tmp0(1:nsteps));
 //%plot(i_t_ne_z2+1/(1-alpha)*A_t_ne_z2+(1)/(1-alpha)*Z_t_ne_z2,'r--');
 //tmp=(i_t_ne_z2+1/(1-alpha)*A_t_ne_z2+(1)/(1-alpha)*Z_t_ne_z2);
//plot(tmp(1:nsteps),'r--');
// plot(nperiod,tmp(nperiod),'g*');
// title('Temporary I-shock');
legend('Tech Shock','News Shock (Four quarters advance)');
suptitle('Response of Investment');

//I_t_e_z2=tmp0;
//I_t_ne_z2=tmp; 

%% Output
figure;
 subplot(2,2,1);
 hold on;
 tmp0=y_t_e_a1+1/(1-alpha)*A_t_e_a1+(1-alpha)/(1-alpha)*Z_t_e_a1;
 
 plot(tmp0(1:nsteps));
 tmp=y_t_ne_a1+1/(1-alpha)*A_t_ne_a1+(1-alpha)/(1-alpha)*Z_t_ne_a1;
 plot(tmp(1:nsteps),'r--');
 tmp=y_t_ne_a1+1/(1-alpha)*A_t_ne_a1+(1-alpha)/(1-alpha)*Z_t_ne_a1;
 plot(nperiod,tmp(nperiod),'g*');
 title('Permanent N-shock');
 Y_t_e_a1=tmp0;
 Y_t_ne_a1=tmp; 

// subplot(2,2,2);
// hold on;
// tmp0=y_t_e_a2+1/(1-alpha)*A_t_e_a2+(1-alpha)/(1-alpha)*Z_t_e_a2;
// plot(tmp0(1:nsteps));
// tmp=y_t_ne_a2+1/(1-alpha)*A_t_ne_a2+(1-alpha)/(1-alpha)*Z_t_ne_a2;
// plot(tmp(1:nsteps),'r--');
// plot(nperiod,tmp(nperiod),'g*');
// title('Temporary N-shock');
// Y_t_e_a2=tmp0;
// Y_t_ne_a2=tmp; 


 subplot(2,2,3);
 hold on;
 tmp0=(y_t_e_z1+1/(1-alpha)*A_t_e_z1+(1-alpha)/(1-alpha)*Z_t_e_z1);
 plot(tmp0(1:nsteps));
 tmp=y_t_ne_z1+1/(1-alpha)*A_t_ne_z1+(1-alpha)/(1-alpha)*Z_t_ne_z1;
 plot(tmp(1:nsteps),'r--');
 plot(nperiod,tmp(nperiod),'g*');
 title('Permanent I-shock');
 Y_t_e_z1=tmp0;
 Y_t_ne_z1=tmp; 

 //subplot(2,2,4);
 //hold on;
 //tmp0=(y_t_e_z2+1/(1-alpha)*A_t_e_z2+(1-alpha)/(1-alpha)*Z_t_e_z2);
 //plot(tmp0(1:nsteps));
//% plot(y_t_ne_z2+1/(1-alpha)*A_t_ne_z2+(1-alpha)/(1-alpha)*Z_t_ne_z2,'r--');
//tmp=y_t_ne_z2+1/(1-alpha)*A_t_ne_z2+(1-alpha)/(1-alpha)*Z_t_ne_z2;
//plot(tmp(1:nsteps),'r--');
//plot(nperiod,tmp(nperiod),'g*');
//legend('Tech Shock','News Shock (Four quarters advance)');
//Y_t_e_z2=tmp0;
//Y_t_ne_z2=tmp; 
// title('Temporary I-shock');
suptitle('Response of Output');

%% Replicate JR2008 Figure 2
figure;
subplot(2,3,1);
hold on;
plot([0:nsteps],[0; Y_t_ne_a1(1:nsteps)]);
plot([0:nsteps],[0;Y_t_ne_z1(1:nsteps)],'r.-');
plot(nperiod,Y_t_ne_z1(nperiod),'g*');
title('Output');

subplot(2,3,2);
hold on;
plot([0:nsteps],[0; C_t_ne_a1(1:nsteps)]);
plot([0:nsteps],[0; C_t_ne_z1(1:nsteps)],'r.-');
plot(nperiod,C_t_ne_z1(nperiod),'g*');
title('Consumption');

subplot(2,3,3);
hold on;
plot([0:nsteps],[0; I_t_ne_a1(1:nsteps)]);
plot([0:nsteps],[0; I_t_ne_z1(1:nsteps)],'r.-');
plot(nperiod,I_t_ne_z1(nperiod),'g*');
title('Investment');

subplot(2,3,4);
hold on;
plot([0:nsteps],[0; n_t_ne_a1(1:nsteps)]);
plot([0:nsteps],[0; n_t_ne_z1(1:nsteps)],'r.-');
plot(nperiod,n_t_ne_z1(nperiod),'g*');
legend1=legend('N-shock','I-shock');
set(legend1,'Position',[0.6274 0.1742 0.1893 0.1373]);

title('Labor Hours');


suptitle('Response to News on Permanent Shocks');

figure;
subplot(2,3,1);
hold on;
plot([0:nsteps],[0; Y_t_e_a1(1:nsteps)]);
plot([0:nsteps],[0;Y_t_e_z1(1:nsteps)],'r.-');
plot(nperiod,Y_t_e_z1(nperiod),'g*');
title('Output');

subplot(2,3,2);
hold on;
plot([0:nsteps],[0; C_t_e_a1(1:nsteps)]);
plot([0:nsteps],[0; C_t_e_z1(1:nsteps)],'r.-');
plot(nperiod,C_t_e_z1(nperiod),'g*');
title('Consumption');

subplot(2,3,3);
hold on;
plot([0:nsteps],[0; I_t_e_a1(1:nsteps)]);
plot([0:nsteps],[0; I_t_e_z1(1:nsteps)],'r.-');
plot(nperiod,I_t_e_z1(nperiod),'g*');
title('Investment');

subplot(2,3,4);
hold on;
plot([0:nsteps],[0; n_t_e_a1(1:nsteps)]);
plot([0:nsteps],[0; n_t_e_z1(1:nsteps)],'r.-');
plot(nperiod,n_t_e_z1(nperiod),'g*');
legend1=legend('N-shock','I-shock');
set(legend1,'Position',[0.6274 0.1742 0.1893 0.1373]);

title('Labor Hours');

suptitle('Response to Permanent Shocks');