var c h k b Y A w g tau I rk rb pistar pi D lamda y omega epsilon f1 f2 ;
varexo e_A e_g e_tau e_omega e_epsilon;
parameters habit alpha beta delta phi theta Bphi site rho_A rho_g rho_tau rho_omega rho_epsilon c_ss h_ss lamda_ss w_ss I_ss k_ss y_ss Y_ss pistar_ss pi_ss b_ss rb_ss g_ss tau_ss rk_ss D_ss f1_ss f2_ss ;

habit=0.9;
alpha=0.45;
beta=0.99;
delta=0.025;
phi=0.276;
theta=4;
Bphi=0.75;
rho_A=0.9;
rho_g=0.8;
rho_tau=0.8;
rho_omega=0.33;
rho_epsilon=0.74;

h_ss=1/3;
tau_ss=0.2;
omega_ss=0.1;
epsilon_ss=0.1;
rb_ss=(1/beta)-1;
rk_ss=((1/beta)-1+delta)/(1-tau_ss);
pi_ss=1;
pistar_ss=1;
w_ss=0.032329;
k_ss=h_ss*(alpha/(1-alpha))*w_ss/rk_ss;
I_ss=k_ss*delta;
y_ss=(k_ss^alpha)*(h_ss^(1-alpha));
Y_ss=y_ss;
g_ss=0.2*Y_ss;
c_ss=Y_ss-g_ss-I_ss;
lamda_ss=((1-habit*beta)*omega_ss/c_ss)/(1-habit);
site=lamda_ss*w_ss*(1-tau_ss)*epsilon_ss/(h_ss^phi);
D_ss=1;
b_ss=0.4*Y_ss;
f1_ss=(lamda_ss*((theta-1)/theta)*Y_ss)/(1-beta*Bphi);
f2_ss=(lamda_ss*Y_ss)/(1-beta*Bphi);

model;

#comsumer' problem
lamda=(c(+1)+habit*habit*beta*c(-1)-habit*(1+beta)*c+(1-habit)*omega+beta*habit*(habit-1)*omega(+1))/(1+habit*habit*beta-habit-habit*beta)-((1-2*habit+habit*habit)*c+(habit*habit-habit)*c(-1)+(1-habit)*c(+1))/((1-habit)*(1-habit));
phi*h=lamda+w+epsilon-tau*tau_ss/(1-tau_ss);
0=lamda(+1)-lamda+((1-tau_ss)*rk_ss*rk(+1)-rk_ss*tau_ss*tau(+1))/(rk_ss-rk_ss*tau_ss+1-delta);
0=lamda(+1)-lamda+rb_ss*rb(+1)/(1+rb_ss);
I_ss*I=k_ss*k-(1-delta)*k_ss*k(-1);

#firm's problem
w=A+alpha*k(-1)-alpha*h;
rk=A+(alpha-1)*k(-1)+(1-alpha)*h;
Y=A+alpha*k(-1)+(1-alpha)*h;

#NKPC
pistar=f1-f2;
f1=((lamda_ss*((theta-1)/theta)*Y_ss)*(lamda+((1-alpha)^(alpha-1)*(alpha)^(-alpha)*(w_ss)^(1-alpha)*(rk_ss)^alpha*((1-alpha)*w+alpha*rk-A))+Y)+beta*Bphi*f1_ss*(theta*pi(+1)+f1(+1))/(lamda_ss*((theta-1)/theta)*Y_ss+beta*Bphi*f1_ss));
f2=((lamda_ss*Y_ss)*(lamda+Y)+beta*Bphi*f2_ss*((theta-1)*pi(+1)+f2(+1)))/(lamda_ss*Y_ss+beta*Bphi*f2_ss);
(1-Bphi)*(1-theta)*pistar+Bphi*(theta-1)*pi=0;

#inflation process
D=(1-Bphi)*(-theta)*pistar+Bphi*(theta*pi+D(-1));
Y+D=y;

#goverment
g_ss*g+b_ss*rb_ss*rb+(1+rb)*b_ss*b(-1)=b_ss*b+h_ss*tau_ss*w_ss*(w+h+tau)+k_ss*rk_ss*tau_ss*(rk+k(-1)+tau);

#external shock
tau=rho_tau*tau(-1)+e_tau;
A=rho_A*A(-1)+e_A;
g=rho_g*g(-1)+e_g;
omega=rho_omega*omega(-1)+e_omega;
epsilon=rho_epsilon*epsilon(-1)+e_epsilon;

#market clearing
Y_ss*Y=c_ss*c+I_ss*I+g_ss*g;

end;

initval;
h=0;
tau=0;
omega=0;
epsilon=0;
rb=0;
rk=0;
pi=0;
pistar=0;
w=0;
k=0;
I=0;
y=0;
Y=0;
D=0;
g=0;
c=0;
lamda=0;
b=0;
A=0;
f1=0;
f2=0;
e_tau=0;  e_A=0;  e_g=0; e_omega=0; e_epsilon=0;
end;
steady;
check;
shocks;
var e_tau;stderr 0.01;
var e_A;stderr 0.01;
var e_g;stderr 0.01;
var e_omega;stderr 0.01;
var e_epsilon;stderr 0.01;
end;
stoch_simul(order=1,irf=40,hp_filter=1600,periods=0) Y I c h w;