var w_sc w_sh w_bc w_bh l_sc l_sh l_bc l_bh q pi r_s c_s h_s u_s i_s r_b c_b h_b u_b i_b 
    b_b m_b y_h m_h b_h a_h y_c a_c y x;

varexo e, u;

parameters eta tau beta1 beta2 beta3 delta phi_i alpha j theta rho x4 phi_r phi_pi phi_y phi_h phi_b
           x1 q_star rb_star mb_star mh_star  x_star ib_cstar bb_cstar wl_cstar hb_cstar yc_ystar yh_ystar
           cs_ycstar cb_ycstar is_yhstar ib_yhstar ls_lbstar gamma gamma1;             
           
eta  = 2;
tau = 1;
beta1 = 0.99;
beta2 = 0.98;
beta3 = 0.97;
delta = 0.025;
phi_i = 2;
alpha = 0.55;
j = 0.1;
theta = 0.64;
rho = 0.75;
x4 = 0.5;
phi_r = 0.7;
phi_pi = 0.5;
phi_y = 0.5;
phi_h = 0.5;
phi_b = 0.5;
gamma = 0.95;
gamma1 = 0.95;
q_star=2;
rb_star=1.025;
mb_star=0.56;
mh_star=0.68;
x_star=1.2;


model(linear);
(x1+1)*(w_sc-c_s)=(x1*(eta-1)+tau)*l_sc+(eta-1-tau)*l_sh;
(x1+1)*(w_sh-c_s)=x1*(eta-1-tau)*l_sc+(x1*tau+eta-1)*l_sh;
r_s+c_s=pi(+1)+c_s(+1);
h_s=beta1*(1-delta)/(1-beta1*(1-delta))*u_s(+1)-u_s/(1-beta1*(1-delta));
q-c_s=u_s-phi_i*(i_s-i_s(-1))+beta1*phi_i*(i_s(+1)-i_s);
i_s=1/delta*h_s-(1-delta)/delta*h_s(-1);

(x1+1)*(w_bc-c_b)=(x1*(eta-1)+tau)*l_bc+(eta-1-tau)*l_bh;
(x1+1)*(w_bh-c_b)=x1*(eta-1-tau)*l_bc+(x1*tau+eta-1)*l_bh;

c_b+q_star*ib_cstar*(q+i_b)+rb_star*bb_cstar*(r_b(-1)+b_b(-1)-pi)=bb_cstar*b_b+alpha*wl_cstar*(w_bc+l_bc)+(1-alpha)*wl_cstar*(w_bh+l_bh);
-j/(hb_cstar*q_star)*h_b-u_b+(1-delta)*beta2*u_b(+1)=(1/rb_star-beta2)*mb_star*(m_b+q(+1))+mb_star/rb_star*(pi(+1)-c_b-r_b)+beta2*mb_star*c_b(+1);
 q-c_b=u_b-phi_i*(i_b-i_b(-1))+beta2*phi_i*(i_b(+1)-i_b);
i_b=1/delta*h_b-(1-delta)/delta*h_b(-1);
b_b=m_b+q(+1)+h_b+pi(+1)-r_b;

y_h=a_h+theta*l_sh+(1-theta)*l_bh;
(1-beta3*mh_star)*(w_sh+l_sh-q-y_h)=mh_star/rb_star*(b_h-w_sh-l_sh)-beta3*mh_star*m_h;
(1-beta3*mh_star)*(w_bh+l_bh-q-y_h)=mh_star/rb_star*(b_h-w_bh-l_bh)-beta3*mh_star*m_h;
b_h=m_b+q+y_h+pi(+1)-r_b;

y_c=a_c+theta*l_sc+(1-theta)*l_bc;
w_sc=y_c-x-l_sc;
w_bc=y_c-x-l_bc;
pi=beta3*pi(+1)+(1-rho)*(1-beta3*rho)/rho*((1-theta)*w_bc+theta*w_sc-a_c);

r_s=phi_r*r_s(-1)+(1-phi_r)*((1+phi_pi)*pi+phi_y*(y-y(-1)));
y=yc_ystar*y_c+yh_ystar*y_h;
y_c=x4*cs_ycstar*c_s+(1-x4)*cb_ycstar*c_b;
y_h=x4*is_yhstar*i_s+(1-x4)*ib_yhstar*i_b;
m_h=-phi_h*(b_h-b_h(-1));
m_b =-phi_b*(b_b-b_b(-1));
x4*ls_lbstar*(alpha*l_sc+(1-alpha)*l_sh)=(x4-1)*(alpha*l_bc+(1-alpha)*l_bh);
a_h = gamma*a_h(-1) + e;
a_c = gamma1*a_c(-1) + u;
end;

steady_state_model;
x1=alpha/(1-alpha);
hb_cstar=j/(q_star*(mb_star/rb_star*(1-beta2*rb_star)-((1-delta)*beta2-1)));
ib_cstar=delta*hb_cstar;
bb_cstar=mb_star*q_star/rb_star*hb_cstar;
wl_cstar=1+(rb_star-1)*bb_cstar;
yh_ystar=1/(1+(q_star*alpha*x_star*(1-beta3*mh_star+mh_star/rb_star))/(1-alpha));
yc_ystar=1-yh_ystar;
cb_ycstar=((1-theta)/(alpha*x_star))/(1+(rb_star-1)*bb_cstar);
cs_ycstar=(1-(1-x4)*cb_ycstar)/x4;
ib_yhstar=(1/yh_ystar-1)*cb_ycstar*ib_cstar;
is_yhstar=(1-(1-x4)*ib_yhstar)/x4;
ls_lbstar=(x4/(1/cb_ycstar+x4-1)*(theta/(1-theta)))^(1/eta);

end;

shocks;
var e; stderr 0.009;
var u; stderr 0.009;

end;

stoch_simul;