// This file replicates the simulation results of the CIA model from 
// Thomas F. Cooley and Gary D. Hansen (1989) ``The Inflation Tax in a 
// Real Business Cycle Model'', The American Economic Review, 79(4), 
// pp.733--48.
// This file replicates Table 1 p.741: 
// -panel 2 when setting sig_e_g = 0
// -panel 3 when setting sig_e_g = 0.009 and gbar = 1.015
// -panel 4 when setting sig_e_g = 0.009 and gbar = 1.15
//
// Christophe Cahn, December 2009


var c x w r y k h phat mhat lam muu z g;
varexo e_z e_g;

parameters delt B bet thet gam alp sig_e_z sig_e_g gbar;
delt = 0.025 ;
bet = 0.99 ;
thet = 0.36 ;
B = 2.86 ;
gam = 0.95 ;
alp = 0.48 ;
gbar = 1.015 ;
sig_e_z = 0.00721 ;

// À choisir en fonction de la simulation
sig_e_g = 0.009;
// sig_e_g = 0;

model;
	// FOC on c, h, k, m
	1/c = lam + muu ;
	lam*w = B ;
	lam = bet*lam(1)*(r(1) + 1 - delt) ;
	lam = bet*(lam(1) + muu(1))/(phat(1)*g(1)/phat) ;
	// constraint for each multiplier muu, lam
	phat*c = (mhat(-1) + g - 1)/g ;
	c + x + mhat/phat = w*h + r*k(-1) + (mhat(-1)+g-1)/phat/g ;
	// investment x
	k = (1-delt)*k(-1) + x;
	// equilibrium condition for prices r w phat
	w = (1-thet)*y/h ;
	r = thet*y/k(-1) ;
	mhat = 1 ;
	// production y
	y = exp(z)*k(-1)^thet*h^(1-thet) ;
	// two shocks z, g
	z = gam*z(-1) + e_z ;
	log(g) = alp*log(g(-1)) + (1-alp)*log(gbar) +  e_g ;
end;

initval;
	g = gbar ;
	z = 0 ;
	mhat = 1 ;
	r = 1/bet - 1 + delt ;
	w = (1-thet)*(thet/r)^(thet/(1-thet)) ;
	lam = B/w ;
	c = bet/lam/gbar ;
	phat = 1/c ;
	h = c/(w + (thet/r)^(1/(1-thet))*(r-delt)) ;
	y = w*h/(1-thet) ;
	k = thet*y/r ;
	x = delt*k ;
	muu = 1/c - lam ;
end;

shocks;
	var e_z; stderr sig_e_z ;
	var e_g; stderr sig_e_g ;
end;

res = zeros(7,50) ;
correl = zeros(7,50) ;

// 5O iterations over 115-period simulations
for ii = 1:50
	stoch_simul(periods=115,order=1,nograph,noprint); 
	dlog_y = log(y)-hpfilter(log(y),1600);
	dlog_c = log(c)-hpfilter(log(c),1600);
	dlog_x = log(x)-hpfilter(log(x),1600);
	dlog_k = log(k)-hpfilter(log(k),1600);
	dlog_h = log(h)-hpfilter(log(h),1600);
	dlog_t = dlog_y - dlog_h ;
	dlog_p = log(phat.*cumprod(g))  - hpfilter(log(phat.*cumprod(g)),1600);
	res(1,ii) = std(dlog_y) ;
	res(2,ii) = std(dlog_c) ;
	res(3,ii) = std(dlog_x) ;
	res(4,ii) = std(dlog_k) ;
	res(5,ii) = std(dlog_h) ;
	res(6,ii) = std(dlog_t) ;
	res(7,ii) = std(dlog_p) ;
	correl(1,ii) = corr(dlog_y,dlog_y) ;
	correl(2,ii) = corr(dlog_c,dlog_y) ;
	correl(3,ii) = corr(dlog_x,dlog_y) ;
	correl(4,ii) = corr(dlog_k(1:end-1),dlog_y(2:end)) ;
	correl(5,ii) = corr(dlog_h,dlog_y) ;
	correl(6,ii) = corr(dlog_t,dlog_y) ;
	correl(7,ii) = corr(dlog_p,dlog_y) ;
end;

mat_output = [mean(res'*100)' std(res'*100)' mean(correl')' std(correl')'];
names = {'Output  ', 'Consumption', 'Investment', 'Capital Stock', 'Hours     ', 'Productivity', 'Price level'};
fprintf('\n')
fprintf('Results\n\n')
fprintf('Series\t\tStd Dev.\tCorr w. Y\n')
for	ii = 1:7
	fprintf('%s\t%1.2f(%1.2f)\t%1.2f(%1.2f)\n', char(names(ii)), mat_output(ii,:)) ;
end