//%--------------------------------------------
//% Declare endogenous and exogenous variables  
//%--------------------------------------------

var c l w pai d A q x y_L k_L l_L n_L z_L ksi_L y_S k_S l_S n_S z_S ksi_S y k i ksi rl bl rs bs rf f rd rr;

varexo eps_r eps_A eps_x;


//%-----------------------------------------------
//% Declare model parameters
//%-----------------------------------------------
parameters beta sigma phi chi delta phi_up z_L_ss rl_ss nju k_L_ss n_L_ss psi_L alpha z_S_ss rs_ss
gamma k_S_ss n_S_ss psi_S rho y_S_ss y_L_ss l_S_ss l_L_ss ksi_S_ss ksi_L_ss c_ss i_ss
sigma1 bs_ss rf_ss eta d_ss rhor rhoy phirl rho_A bl_ss f_ss rho_x;

//%----------------------------------------------
//% Setting the values of calibrated parameters
//%----------------------------------------------
beta=0.99;
sigma=1;
phi=2;
chi=0.3;
delta=0.035;
phi_up=0.5;
z_L_ss=1;
rl_ss=0.00135;
nju=0.98;
k_L_ss=0.65;
n_L_ss=0.25;
psi_L=0.9;
alpha=0.4;
z_S_ss=1;
rs_ss=0.1;
gamma=0.6;
k_S_ss=0.35;
n_S_ss=0.15;
psi_S=0.8;
rho=6;
y_S_ss=0.3859;
y_L_ss=0.6141;
l_S_ss=0.7;
l_L_ss=0.3;
ksi_S_ss=0.5;
ksi_L_ss=0.5;
c_ss=0.4;
i_ss=0.6;
sigma1=0.3;
bs_ss=0.4;
rf_ss=0.009;
eta=0.15;
d_ss=2;
rhor=1.5;
rhoy=0.125;
phirl=0.7;
rho_A=0.9;
bl_ss=0.5;
f_ss=0.1;
rho_x=0.5;




//%---------------------------------------------
//% model equations
//%---------------------------------------------
model(linear);

//% 1
c=c(+1)-(rd-pai)/sigma;

//% 2
w=sigma*c+phi*l;

//% 3
q=chi*delta*(i-k)-x;


//% 4
k(+1)=delta*i+delta*x+(1-delta)*k;

//% 5
pai=beta*pai(+1)+(1-beta*phi_up)*(1-phi_up)/phi_up*ksi;


//% 6
rl=z_L_ss/rl_ss*z_L(+1)+(1-delta)/rl_ss*q(+1)-q;


//% 7
n_L(+1)/nju*rl_ss=k_L_ss/n_L_ss*rl-(k_L_ss/n_L_ss-1)*(rl(-1)-pai)-psi_L*(k_L_ss/n_L_ss-1)*(k_L+q(-1))+(psi_L*(k_L_ss/n_L_ss-1)+1)*n_L;

//% 8
z_L=y_L+ksi-k_L;

//% 9
w=y_L+ksi-l_L;

//% 10
y_L=alpha*k_L+(1-alpha)*l_L+(1-alpha)*A;

//% 11
rs(+1)=z_S_ss/rs_ss*z_S(+1)+(1-delta)/rs_ss*q(+1)-q;

//% 12
n_S(+1)/(nju*(1-gamma)*rs_ss)=k_S_ss/n_S_ss*rs-(k_S_ss/n_S_ss-1)*(rs(-1)-pai)-psi_S*(k_S_ss/n_S_ss-1)*(k_S+q(-1))+(psi_S*(k_S_ss/n_S_ss-1)+1)*n_S+gamma/(1-gamma)*(rs+n_S);

//% 13
z_S=y_S+ksi-k_S;

//% 14
w=y_S+ksi-l_S;

//% 15
y_S=alpha*k_S+(1-alpha)*l_S+(1-alpha)*A;

//% 16
y=(y_S_ss)*y_S+(y_L_ss)*y_L;

//% 17
k=k_S_ss*k_S+k_L_ss*k_L;

//% 18
l=l_S_ss*l_S+l_L_ss*l_L;

//% 19
y=(ksi_S_ss)*ksi_S+(ksi_S_ss)*ksi_L;

//% 20
y=c_ss*c+i_ss*i;

//% 21
rs_ss*rs=rl_ss*rl+sigma1*bs_ss*bs;

//% 22
rf_ss*rf=rl_ss*rl;

//% 23
bl_ss*bl+bs_ss*bs+f_ss*f=(1-eta)*d;

//% 24
rd=rhor*pai+rhoy*y+eps_r;

//% 25
rl=(1-phirl)*rd+phirl*rl(-1);

//% 26
rr=rd;

//% 27
bl=q+k_L(+1)-n_L(+1);

//% 28
bs=q+k_S(+1)-n_S(+1);

//% 29
A=rho_A*A(-1)+eps_A;

//% 30
f=alpha*d;

//%31
x=rho_x*x(-1)+eps_x;

//%32
c+d=rd*d(-1)+w*l;

end;


//%-------------------------------------------------
//% one standerr shock
//%-------------------------------------------------
shocks;
var eps_A=0.01;
var eps_r=0.01;
var eps_x=0.01;
end;

%----------------------------------------------------------------
% generate IRFs
%----------------------------------------------------------------
stoch_simul(order = 1,irf=15) y c i ;