var y c m i h k n w z R f q lamda xi pi mu; 

%Financial Accelerators variables: 

varexo eps_R A B e x;
   
parameters beta eta theta delta b piss v S kn psi chi alpha gamma phi varrho_pi 
varrho_y varrho_mu Ass xiss Rss fss zss lamdass css mss kss yss wss hss iss;

%Parameter calibrations, see Table 1 and Table 2 in paper
beta       =   0.9928;
eta        =   1.315;
theta      =   6;
delta      =   0.025;
b          =   0.062;
piss       =   1.0079;
v          =   0.9728;
S          =   1.0075;
kn         =   2;

%Parameter estimates using maximum likelihood
psi        =   0.0420;
chi        =   0.5882;
alpha      =   0.3384;
gamma      =   0.0598;
phi        =   0.7418;
varrho_pi  =   1.4059;
varrho_y   =   0.2947;
varrho_mu  =   0.6532;
Ass        =   800; %Steady state technology level

%steady-state equilibrium using for the log-linearized equations

xiss=(theta-1)/theta;
Rss=piss/beta;
fss=S*Rss/piss;
zss=fss-1+delta;
lamdass_css=(1+b*(piss/(piss-beta))^(gamma-1))^(-1);
lamdass_mss=lamdass_css*b*(piss/(piss-beta))^gamma;
kss_yss=alpha *(xiss/zss); %kss/yss
css_yss=1-delta*(kss_yss); %css/yss
wss_hss_lamdass=(1-alpha)*(lamdass_css)*xiss/(css_yss);
hss=wss_hss_lamdass / (eta+wss_hss_lamdass);
yss=Ass*hss*(kss_yss)^(alpha/(1-alpha));
kss=kss_yss*yss;
css=css_yss*yss;
iss=delta*kss;
lamdass=lamdass_css/css;
mss=lamdass_mss/lamdass;
wss=wss_hss_lamdass/(hss*lamdass);

%The log-linearized equilibrium system
model(linear); 
    
    ((1-gamma)*(lamdass*css)-1)*c = gamma*lamda + (lamdass*mss)*(Rss-1)*(B+(gamma-1)*m)/Rss - gamma*e;
    (gamma*R)/(Rss-1) = B + c - m;
    hss*h = (1-hss)*(w+lamda);
    y = alpha*k + (1-alpha)*h + (1-alpha)*A;
    yss*y = css*c + iss*i;
    w = y + xi - h;
    z = y + xi - k;
    mu = m - m(-1) + pi;
    R = varrho_pi*pi + varrho_mu*mu + varrho_y*y+eps_R;
    f = zss*z/fss + (1-delta)*q/fss - q(-1);
    q = chi*(i-k) - x;
    
    %New Keynesian Phillips Curve:
    pi = beta*pi(+1) + (1-beta*phi)*(1-phi)*xi/phi;
    
    lamda(+1) = lamda - R + pi(+1);
    k(+1) = delta*i + delta*x + (1-delta)*k;
    f(+1) = R - pi(+1) + psi*(q+k(+1)-n(+1));
    n(+1)/(v*fss) = (kn*f) - (kn-1)*(R(-1) - pi) - psi*(kn-1)*(k+q(-1)) + (psi*(kn-1)+1)*n;
    
end;

initval;
y = 0; 
c = 0;
R = 0;
h = 0;
k = 0;
i = 0; 
pi = 0;
f = 0; 
n = 0;
end;

steady;
resid;
check;

shocks;
var eps_R;      
stderr 0.0058;
var A;       
stderr 0.0096;
var B;       
stderr 0.0103;
var e;       
stderr 0.0073;
var x;       
stderr 0.0331;
end;

estimated_params; 
psi,         beta_pdf, 0.0420, 0.0137; 
chi,         beta_pdf, 0.5882, 0.1742; 
alpha,       beta_pdf, 0.3384, 0.0259;
gamma,       beta_pdf, 0.0598, 0.0039;
phi,         beta_pdf, 0.7418, 0.0118;
varrho_pi,   beta_pdf, 1.4059, 0.0788;
varrho_y,    beta_pdf, 0.2947, 0.0690;
varrho_mu,   beta_pdf, 0.6532, 0.0783;
end;

stoch_simul(order=1,irf=100) y i c h R pi k n f;