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

%Financial Accelerators variables: 

varexo eps_R e_eps_A e_eps_b e_eps_e e_eps_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 rhoA
rhob rhoe rhox;

%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

%shock parameters
rhoA = 0.7625;
rhob = 0.7206;
rhoe = 0.6156;
rhox = 0.6562;

%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/f + (1-delta)*q/fss - q(-1);
    q = chi*(i-k) - x;
    pi = beta*pi(+1) + (1-beta*phi)*(1-phi)*xi/phi;
    lamda = lamda(+1) + R - pi(+1);
    (1-delta)*k = k(+1) - delta*i - delta*x;
    f(+1) = R - pi(+1) + psi*(q+k(+1)-n(+1));
    (psi*(kn-1) +1)*n = n(+1)/(v*f) - (kn*f) + (kn-1)*(R(-1) + pi) + psi*(kn-1)*(k+q(-1));
    %shock process
    log(A) = (1-rhoA)*log(Ass) + rhoA*log(A(-1)) + e_eps_A;
    log(B) = (1-rhob)*log(b) + rhob*log(B(-1)) + e_eps_b;
    log(e) = rhoe*log(e(-1)) + e_eps_e;
    log(x) = rhox*log(x(-1)) + e_eps_x;
end;

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

check;
steady;

shocks;
var eps_R;      
stderr 0.0058;
var e_eps_A;       
stderr 0.0096;
var e_eps_b;       
stderr 0.0103;
var e_eps_e;       
stderr 0.0073;
var e_eps_x;       
stderr 0.0331;
end;


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