// This is the Ramsey model with adjustment costs.  Jermann(1998),JME 41, pages 257-275
// Olaf Weeken
// Bank of England, 13 June, 2005
// modified January 20, 2006 by Michel Juillard

//---------------------------------------------------------------------
// 1. Variable declaration
//---------------------------------------------------------------------

var c, d, erp1, i, k, m1, r1, rf1, w, y, z, mu; 
varexo ez; 

//---------------------------------------------------------------------
// 2. Parameter declaration and calibration
//---------------------------------------------------------------------

parameters alf, chihab, xi, delt, tau, g, rho, a1, a2, betstar, bet;

alf        = 0.36;    // capital share in production function
chihab     = 0.819;   // habit formation parameter
xi         = 1/4.3;   // capital adjustment cost parameter
delt       = 0.025;   // quarterly deprecition rate
g          = 1.005;    //quarterly growth rate (note zero growth =>g=1)
tau        = 5;       // curvature parameter with respect to c
rho        = 0.95;    // AR(1) parameter for technology shock

a1         = (g-1+delt)^(1/xi);             
a2         = (g-1+delt)-(((g-1+delt)^(1/xi))/(1-(1/xi)))*((g-1+delt)^(1-(1/xi))); 
betstar    = g/1.011138;
bet        = betstar/(g^(1-tau));             

//---------------------------------------------------------------------
// 3. Model declaration
//---------------------------------------------------------------------

model;  
g*k  = (1-delt)*k(-1) + ((a1/(1-1/xi))*(g*i/k(-1))^(1-1/xi)+a2)*k(-1);
d    = y - w - i; 
w    = (1-alf)*y;
y    = z*g^(-alf)*k(-1)^alf;
c    = w + d; 
mu   = (c-chihab*c(-1)/g)^(-tau)-chihab*bet*(c(+1)*g-chihab*c)^(-tau);
mu   = (betstar/g)*mu(+1)*(a1*(g*i/k(-1))^(-1/xi))*(alf*z(+1)*g^(1-alf)*(k^(alf-1))+
       (1-delt+(a1/(1-1/xi))*(g*i(+1)/k)^(1-1/xi)+a2)/
       (a1*(g*i(+1)/k)^(-1/xi))-g*i(+1)/k);
log(z) = rho*log(z(-1)) + ez;

m1   = (betstar/g)*mu(+1)/mu;
rf1  = 1/m1;
r1   = (a1*(g*i/k(-1))^(-1/xi))*(alf*z(+1)*g^(1-alf)*(k^(alf-1))+
       (1-delt+(a1/(1-1/xi))*(g*i(+1)/k)^(1-1/xi)+a2)/
       (a1*(g*i(+1)/k)^(-1/xi))-g*i(+1)/k);
erp1 = r1 - rf1;
end;

//---------------------------------------------------------------------
// 4. Initial values and steady state
//---------------------------------------------------------------------

initval;
m1     = betstar/g;
rf1    = (1/m1);
r1     = (1/m1);
erp1   = r1-rf1;

z      = 1;
k      = (((g/betstar)-(1-delt))/(alf*g^(1-alf)))^(1/(alf-1));
y      = (g^(1-alf))*k^alf;
w      = (1-alf)*y;
i      = (1-(1/g)*(1-delt))*k;
d      = y - w - i;
c      = w + d;

mu     = ((c-(chihab*c/g))^(-tau))-chihab*bet*((c*g-chihab*c)^(-tau));

ez     = 0;
end;

steady;        
check;

//---------------------------------------------------------------------
// 5. Shock declaration  
//                       
//---------------------------------------------------------------------

shocks;
var ez; stderr 0.01;  
end;

// No pruning case

stoch_simul(order=3, irf=40, nograph);

nShocks = M_.exo_nbr;                                             % number of shocks 
nIrf = 40;                                                        % length of impulse responses
burnin = 10000;
ex_=zeros(burnin, nShocks);                                       % shock innovation matrix
sim_path_0 = simult_(oo_.dr.ys,oo_.dr, ex_, 3); 
ergodicmean = sim_path_0(:,end);                                  % actually ergodic mean in the absence of shocks     
initial_condition_states = repmat(ergodicmean,1,M_.maximum_lag);

%% Calculating impulse response of a zero shock (starting from the egodic mean in the absence of shocks)

ex_=zeros(nIrf, nShocks);                                         % generate the shock innovation matrix
sim_path = simult_(initial_condition_states,oo_.dr, ex_, 3);      

toez_irf          = sim_path(:,M_.maximum_lag+1:end) - repmat(ergodicmean,1,nIrf);
ctoez_irf_noprune = toez_irf(strmatch('c',M_.endo_names,'exact'),:);


// Pruning case

stoch_simul(order=3, irf=40, pruning, nograph);

nShocks = M_.exo_nbr;                                             % number of shocks 
nIrf = 40;                                                        % length of impulse responses
burnin = 10000;
ex_=zeros(burnin, nShocks);                                       % shock innovation matrix
sim_path_0   = simult_(oo_.dr.ys,oo_.dr, ex_, 3); 
ergodicmean  = sim_path_0(:,end);                                 % actually ergodic mean in the absence of shocks     
initial_condition_states = repmat(ergodicmean,1,M_.maximum_lag);

%% Calculating impulse response of a zero shock (starting from the egodic mean in the absence of shocks)

ex_=zeros(nIrf, nShocks);                                         % generate the shock innovation matrix
sim_path = simult_(initial_condition_states,oo_.dr, ex_, 3); 

toez_irf        = sim_path(:,M_.maximum_lag+1:end) - repmat(ergodicmean,1,nIrf);
ctoez_irf_prune = toez_irf(strmatch('c',M_.endo_names,'exact'),:);

%% Plot IRF to zero shock without pruning and with pruning

figure;
tt = 1 : nIrf;
plot(tt,ctoez_irf_noprune,tt,ctoez_irf_prune,'LineWidth',2);
title('Response of consumption to zero shock');
legend('No pruning','Pruning');