// Title: The Financial Accelerator in an Estimated New Keynesian Model
// Authors: Christensen and Dib
// Publication: Review of Economic Dynamics 2008 (11), 155-178.
// This replication by Martin Harding.
//---------------------------------------------------------------------
// 1. Variable declaration 
//---------------------------------------------------------------------

var y //output (t)
    c //consumption (t)
    m //real balances (t)
    i //investment (t)
    h //hours (t)
    k //capital (t) (at the end of each period entrepreneurs purchase k(t+1)
    n //net worth (t)
    w //real wage (t)
    z //marginal productivity of capital (t)
    R //nominal interest rate (t)
    f //ex-post real return on capital held in (t), (t) ; [f_(t+1) in the paper refers to external funds rate]
    q //real price of capital to be used in period (t+1), (t)  
    lambdav //Lagrange multiplier of the budget constraint (t)
    xi //real marginal cost (t)
    mu //money growth rate (t)
    pi //inflation (t)
    e //preference shock (t)
    b //money-demand shock (t)
    A //technology shock (t)
    x //investment specific shock (t)
    g1 //aggregate price level
    g2 //price of optimizing firms
    pi_o //ratio between p(optimizing firms) and p(aggregate)
    S // external finance premium
    R_r //real interest rate
    rp
    log_y log_c log_m log_i log_h log_k log_n log_w log_z log_R log_f log_q log_lambdav log_xi log_mu log_pi log_pi_o log_S log_R_r log_rp
    y_hat c_hat m_hat i_hat h_hat k_hat n_hat w_hat z_hat R_hat f_hat q_hat lambdav_hat xi_hat mu_hat pi_hat pi_o_hat S_hat R_r_hat rp_hat
    ;


varexo eps_e eps_b eps_A eps_x eps_R;


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

parameters betap //discount factor
           gammap //elasticity of substitution between consumption and real balances
           etap //weight of leisure in the utility function
           nup //probability that an entrepreneur will survive 
           deltap //depreciation rate of capital
           alphap //share of labor in output
           psip //elasticity of the external finance premium w/r to a change in leverage of entrepreneurs
           chip //investment adjustment cost parameter
           phip //probability of not receiving calvo signal (indexing price with SS gross inflation)
           thetap // intermediate-goods elasticity of substitution
           kappap // parameter of the external finance premium function
           kappapp
           varrho_pi //inflation policy coefficient
           varrho_y //output policy coefficient
           varrho_mu //money-growth policy coefficient
           pi_ss //steady state gross inflation rate           
           e_ss  //constant associated to the preference shock
           x_ss  //constant associated to the investment specific shock
           A_ss // constant associated to the technology shock
           b_ss // constant associated to the money demand shock
           q_ss //steady state price of capital
           //R_ss
           rho_e
           rho_b
           rho_A
           rho_x
           pi_o_ss //steady state value of pi_o
           //iotap
           ;


betap=0.9854;
etap=1.315;
thetap=6;
deltap=0.025;
e_ss=1;
x_ss=1;
A_ss=1;
b_ss=0.0655;

//b_ss=1; 
pi_ss=1;
pi_o_ss=1;
q_ss=1;
nup=0.9781; 
kappap=0.3;
kappapp=-0.2;
psip=-0.06; // 0 in the NoFa

gammap=0.0598;
alphap=0.33;
varrho_pi=1.4059;
varrho_y=0.2947;
varrho_mu=0.6532;
chip=0.5882;
phip=0.7418;
rho_A=0.7625;
rho_e=0.6156;
rho_x=0.6562;
rho_b=0.7206;


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


model;

#R_ss = pi_ss/betap;
#S_ss = betap/nup;
#xi_ss = (thetap-1)/thetap;
#z_ss = S_ss*R_ss/pi_ss-1+deltap;
#k__y = alphap*(xi_ss/z_ss);
//#c__m = (((R_ss-1)/R_ss)^(gammap))/b_ss;
#i__y = deltap*k__y;
#c__y = 1-i__y;
#f_ss = 1/nup;
#h_ss = ((e_ss/(1+b_ss*((R_ss-1)/R_ss)^(1-gammap)))*(1-alphap)*xi_ss*(1-deltap*(alphap*(xi_ss/z_ss)))^(-1))/(etap + (e_ss/(1+b_ss*((R_ss-1)/R_ss)^(1-gammap)))*(1-alphap)*xi_ss*(1-deltap*(alphap*(xi_ss/z_ss)))^(-1));
#y_ss = ((alphap*(xi_ss/z_ss))^(alphap/(1-alphap)))*A_ss*h_ss;
#w_ss = (1-alphap)*xi_ss*A_ss*(alphap*(xi_ss/z_ss))^(alphap/(1-alphap));
#k_ss = y_ss*alphap*(xi_ss/z_ss);
#i_ss = deltap*k_ss;
#c_ss = y_ss-i_ss;
#n_ss = 0.5*k_ss;
#m_ss = b_ss*c_ss*(R_ss/(R_ss-1))^(gammap);
#lambdav_ss = (e_ss*c_ss^(-1/gammap))/(c_ss^((gammap-1)/gammap)+(b_ss^(1/gammap))*m_ss^((gammap-1)/gammap));
#mu_ss = pi_ss;
#g1_ss = (lambdav_ss*xi_ss*y_ss)/(1-phip*betap);
#g2_ss = (lambdav_ss*y_ss)/(1-phip*betap);
#R_r_ss=R_ss-pi_ss;
#rp_ss=f_ss-R_ss+pi_ss;


//eq.1
(e*c^(-1/gammap))/(c^((gammap-1)/gammap)+(b^(1/gammap))*m^((gammap-1)/gammap)) = lambdav;

//eq.2
((b*c)/m)^(1/gammap) = (R-1)/R;

//eq.3
etap/(1-h) = lambdav*w;

//eq.4
lambdav/R = betap*(lambdav(+1)/pi(+1));

//eq.5
z = alphap*xi*(y/k(-1));

//eq.6
w = (1-alphap)*xi*(y/h);

//eq.7
y = ((k(-1))^alphap)*(A*h)^(1-alphap);

//eq.8
y = c+i;

//eq.9
1 = phip*(pi_ss/pi)^(1-thetap)+(1-phip)*(pi_o)^(1-thetap);

//eq.10
f(+1) = S*(R/pi(+1));

//eq.11
S = S_ss + psip*(n/(q*k)-0.5) + (1/2)*kappap*(n/(q*k)-0.5)^2 + (1/3)*kappapp*(n/(q*k)-0.5)^3;
//S = n/(n+kappap*((n/(q(-1)*k))^(1/3)-1));
//S = S_ss;

//eq.12
f = (z+(1-deltap)*q)/q(-1);

//eq.13
n = nup*(f*q(-1)*k(-1)-f*(q(-1)*k(-1)-n(-1)));

//eq.14
k = x*i+(1-deltap)*k(-1);

//eq.15
q*x = 1+chip*(i/k(-1)-deltap);

//eq.16
R/(R_ss) = ((pi/pi_ss)^varrho_pi)*((y/y_ss)^varrho_y)*((mu/mu_ss)^varrho_mu)*exp(eps_R);

//eq.17
mu = (m*pi)/m(-1);

//eq.18
g1 = lambdav*xi*y+betap*phip*g1(+1);

//eq.19
g2 = lambdav*y*pi_o+betap*phip*(pi_ss/pi(+1))*(pi_o/pi_o(+1))*g2(+1);

//eq.20
thetap*g1 = (thetap-1)*g2;

//ex. extra.1
R_r=R-pi;

//ex. extra.2
rp=f(+1)-R+pi(+1);

//-------------------------------------------------
// 3.1 Exogenous processes for shocks 
//-------------------------------------------------

//nonlinear model (shocks in logs)

//eq.22
log(e)=rho_e*log(e(-1))+eps_e;

//eq.23
log(b)=(1-rho_b)*log(b_ss)+rho_b*log(b(-1))+eps_b;

//eq.24
log(A)=(1-rho_A)*log(A_ss)+rho_A*log(A(-1))+eps_A;

//eq.25
log(x)=rho_x*log(x(-1))+eps_x;



log_y=log(y);
log_c=log(c);
log_m=log(m);
log_i=log(i);
log_h=log(h); 
log_k=log(k);
log_n=log(n);
log_w=log(w); 
log_z=log(z); 
log_R=log(R); 
log_f=log(f); 
log_q=log(q); 
log_lambdav=log(lambdav); 
log_xi=log(xi); 
log_mu=log(mu); 
log_pi=log(pi); 
log_pi_o=log(pi_o); 
log_S=log(S); 
log_R_r=log(R_r); 
log_rp=log(rp);



y_hat=log(y)-log(y_ss);
c_hat=log(c)-log(c_ss);
m_hat=log(m)-log(m_ss);
i_hat=log(i)-log(i_ss);
h_hat=log(h)-log(h_ss); 
k_hat=log(k)-log(k_ss);
n_hat=log(n)-log(n_ss);
w_hat=log(w)-log(w_ss); 
z_hat=log(z)-log(z_ss); 
R_hat=log(R)-log(R_ss); 
f_hat=log(f)-log(f_ss); 
q_hat=log(q)-log(q_ss); 
lambdav_hat=log(lambdav)-log(lambdav_ss); 
xi_hat=log(xi)-log(xi_ss); 
mu_hat=log(mu)-log(mu_ss); 
pi_hat=log(pi)-log(pi_ss); 
pi_o_hat=log(pi_o)-log(pi_o_ss); 
S_hat=log(S)-log(S_ss); 
R_r_hat=log(R_r)-log(R_r_ss); 
rp_hat=log(rp)-log(rp_ss);

end;

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

steady_state_model;

e = e_ss;
x = x_ss;
A = A_ss;
b = b_ss;
pi = pi_ss;
pi_o = pi_o_ss;
q = q_ss;
R = pi/betap;
xi = (thetap-1)/thetap;
f = 1/nup;
z = f-1+deltap;
h = ((e/(1+b*((R-1)/R)^(1-gammap)))*(1-alphap)*xi*(1-deltap*(alphap*(xi/z)))^(-1))/(etap + (e/(1+b*((R-1)/R)^(1-gammap)))*(1-alphap)*xi*(1-deltap*(alphap*(xi/z)))^(-1));
y = ((alphap*(xi/z))^(alphap/(1-alphap)))*A*h;
w = (1-alphap)*xi*A*(alphap*(xi/z))^(alphap/(1-alphap));
k = y*alphap*(xi/z);
i = deltap*k;
c = y-i;
n = 0.5*k;
m = b*c*(R/(R-1))^(gammap);
lambdav = (e*c^(-1/gammap))/(c^((gammap-1)/gammap)+(b^(1/gammap))*m^((gammap-1)/gammap));
mu = pi;
g1 = (lambdav*xi*y)/(1-phip*betap);
g2 = (lambdav*y)/(1-phip*betap);
S = betap/nup;
R_r=R-pi;
rp=f-R+pi;

log_y=log(y);
log_c=log(c);
log_m=log(m);
log_i=log(i);
log_h=log(h); 
log_k=log(k);
log_n=log(n);
log_w=log(w); 
log_z=log(z); 
log_R=log(R); 
log_f=log(f); 
log_q=log(q); 
log_lambdav=log(lambdav); 
log_xi=log(xi); 
log_mu=log(mu); 
log_pi=log(pi); 
log_pi_o=log(pi_o); 
log_S=log(S); 
log_R_r=log(R_r); 
log_rp=log(rp);


end;

shocks;

var eps_e ; stderr 0.0073; // 0; // 
var eps_b ; stderr 0.0103; // 0; //
var eps_A ; stderr 0.0093; // 0; // 
var eps_x ; stderr 0.0331; // 0; //
var eps_R ; stderr 0.0058; // 0; // 
 


/*
var eps_R ; stderr 0.029;
var eps_e ; stderr 0.0365; 
var eps_b ; stderr 0.0515; 
var eps_A ; stderr 0.048; 
var eps_x ; stderr 0.1655;
*/

end;

steady ;
check ;

stoch_simul(order=1,pruning,irf=20) y,i,c,h,pi,R,S,q,n;
//stoch_simul(order=1,pruning,irf=20) log_y,log_i,log_c, log_h, log_R, log_pi, log_q, log_n log_S; // log_y,log_i,log_S,log_n ; // y_hat,i_hat,c_hat,h_hat,pi_hat,R_hat,S_hat,n_hat,A,x,e,b; //
//stoch_simul(order=2,pruning,replic=200,irf=20) log_y,log_i,log_S,log_n; // y_hat,i_hat,c_hat,h_hat,pi_hat,R_hat,S_hat,n_hat,A,x,e,b; //
//stoch_simul(order=3,pruning,replic=200,periods=500,irf=20) log_y,log_i,log_S,log_n; // y_hat,i_hat,c_hat,h_hat,pi_hat,R_hat,S_hat,n_hat,A,x,e,b; //