// This Dynare program generates the IRFs for the variables in Figure 1 & 2 
// in Den Haan and Kaltenbrunner (2009)
//
// It is not identical to the program used,
// but generates virtually identical results and is simpler than ours.
//
// In particular, we used the minimum subset of 1st-order perturbation solutions 
// and generated others from the model itself. For example, we obtain consumption 
// from the first-order solution and then solve investment from the (non-linear) 
// budget constraint. Here the IRFs of each variable is based on the first-order
// perturbation solution.
//
// In an earlier working paper version of the paper (CEPR discussion paper DP6063)
// we solved a slightly different model with a higer-order projection method and
// obtained similar results
//

var c, n, k, v, z, r, ygross, ynet, pbar, wage, firmv, inv, u, eta, inv2, 
z1, z2, z3, z4, z5, z6, z7, z8, z9, z10, z11, ns, profits, matches, rlevel, ulevel;
varexo e;

parameters beta, rho, alpha, gamma, delta, rox, nu, mu, psi, omega0, omega1, pbarss, kappa, phi, nstar;

rho   = 0.98;
beta  = 0.9966;
delta = 0.0084;
rox = 0.027289091487332;
nu = 0.5;
mu = 0.391729498506304;
nstar = 1.5938;

omega1 = 0.7547;

kappa = 2.652;
omega0 = 0.9772;

//phi   = 0.423140239077383 ;
phi =  0.433426865036823;
psi =  0.793664581158982;
pbarss =  1.151139658738700;
alpha = 0.317778687406195;

//gamma    = 0.45;
gamma    = 0.43;

// !!! The values of phi, psi, & pbarss are functions of other parameters and cannot be varied independently
// !!! To get the right values use the first part of "ourmodelviiisim.m" and use
// !!! psi=psi; phi=phiL; pbarss=log(pb);

model;
exp(n) = (1-rox)*exp(n(-1))+mu*exp((1-nu)*v+nu*u);
exp(-gamma*c)=beta*exp(-gamma*c(+1))*(alpha*exp(z)*exp((alpha-1)*(k-n))+1-delta);
exp(c)+exp(k)+psi*exp(v) = exp(z(-1))*exp(alpha*k(-1)+(1-alpha)*n(-1))+(1-delta)*exp(k(-1));
ygross  = z(-1) + alpha*k(-1) + (1-alpha)*n(-1);
exp(ynet) = exp(ygross) -  psi*exp(v);
pbar = log(1-alpha)+z(-1)+alpha*(k(-1)-n(-1));
r = log(alpha)+z(-1)+(alpha-1)*(k(-1)-n(-1));
exp(wage) = (1-omega1)*omega0*exp(pbarss)+omega1*omega0*exp(pbar);
exp(firmv) = beta*exp(-gamma*(c(+1)-c))*(exp(pbar(+1))-exp(wage(+1))+(1-rox)*exp(firmv(+1)));
exp(inv) = exp(k)-(1-delta)*exp(k(-1));
exp(inv2) = exp(inv)+psi*exp(v);
psi = exp(firmv)*mu*exp(nu*(u-v));
matches=log(mu)+nu*u+(1-nu)*v;
exp(eta)*mu*exp((1-nu)*(v-u))=phi*((nstar-exp(u)-exp(n(-1)))^(-kappa));
exp(eta)=beta*exp(-gamma*c(+1))*exp(wage(+1))+beta*exp(eta(+1))*(1-rox)-beta*phi*((nstar-exp(u(+1))-exp(n))^(-kappa));
exp(ns)=exp(u)+exp(n(-1));
exp(profits)=exp(pbar)-exp(wage);
rlevel=10000*exp(r);
//rlevel is in basis points (0.01 in model corresponds with 1% point & 100 basis points
ulevel=100*exp(u);
z = z1(-1);
z1 = z2(-1);
z2 = z3(-1);
z3 = z4(-1);
z4 = z5(-1);
z5 = z6(-1);
z6 = z7(-1);
z7 = z8(-1);
z8 = z9(-1);
z9 = z10(-1);
z10 = z11(-1);
z11 = rho*z11(-1)+e;
end;

initval;
c      =		 0.83611;
eta    =		 -0.174256;
firmv  =		 1.85494;
inv    =		 -0.518331;
inv2   =		 -0.274386;
k      =		 4.26119;
n      =		 -0.0583709;
ns     =		 1.66131e-013;
pbar   =		 0.860388;
u      =		 -2.86998;
v      =		 -2.57493;
wage   	=	 0.77374;
ygross 	=	 1.12083;
ynet   	=	 1.06571;
end;

shocks;
var e; stderr 0.0042;
end;


steady;
stoch_simul(order=1,nocorr,IRF=24) n ns inv2 inv c ygross;
stoch_simul(order=1,nocorr,IRF=24) profits firmv wage rlevel ulevel matches;