var y y_c y_h c_s c_b h_s h_b k i k_c k_h r n_s n_b pi w r_k q eta eta_l l v z_c z_h mu_h mu_c;
varexo e_v e_zc e_zh e_mu_h e_mu_c e_r;

parameters beta gamma sigma_s sigma_b phi_s phi_b phi_h psi_s psi_b psi delta alpha epsilon_c epsilon_h tao tao_c tao_h gamma_r gamma_pi gamma_y rho_zc rho_zh rho_mu_c rho_mu_h rho_v 
g_hs g_hb g_c g_h g_cs g_cb g_ci g_kc g_kh sigma_v sigma_zc sigma_zh sigma_mu_c sigma_mu_h sigma_r; 


//calibrations
beta=0.99;
gamma=0.97;
sigma_s=0.5;
sigma_b=0.6;
phi_s=0.5;
phi_b=0.8; 
phi_h=0;
psi_s=1.01; 
psi_b=1.01; 
psi=0.15;  
delta=0.025;
alpha=0.5;
epsilon_c=0.75; 
epsilon_h=0.5; 
tao=0.8; 
tao_c=0.5; 
tao_h=0.7;
gamma_r=0.7;
gamma_pi=1.5; 
gamma_y=0.2;
g_hs=0.3; 
g_hb=0.7;//g_hb=1-g_hs;  
g_c=0.6; 
g_h=0.4;//g_h=1-g_c; 
g_cs=0.4; 
g_cb=0.3; 
g_ci=0.3;//g_ci=1-g_cc-g_cb;
g_kc=0.7; 
g_kh=0.3;//g_kh=1-g_kc; 
rho_zc=0.8;
rho_zh=0.8; 
rho_mu_c=0.8;
rho_mu_h=0.8;
rho_v=0.75;
sigma_zc=0.01;
sigma_zh=0.01; 
sigma_mu_c=0.01; 
sigma_mu_h=0.01; 
sigma_v=0.01;
sigma_r=0.01;

model(linear);
//consumption equation
c_s=c_s(+1)-(r-pi(+1))/sigma_s;
c_b=c_b(+1)-(r-pi(+1))/sigma_b;

//capital vlaue 
eta=beta*(1-delta)*eta(+1)-r+pi(+1)+(1-beta*(1-delta))*r_k(+1);

//investment
i=i(-1)/(1+beta)+beta*i(+1)/(1+beta)+psi*eta/(1+beta);

//labor supply
psi_s*n_s=-sigma_s*c_s+w;
psi_b*n_b=-sigma_b*c_b+w;
//n_b=n_s;

k=(1-delta)*k(-1)+delta*i;

//house demand
(phi_h*(1+beta)+phi_s*(1-beta))*h_s=beta*(q(+1)-sigma_s*c_s(+1))-q+sigma_s*c_s+phi_h*(h_s(-1)+beta*h_s(+1))+(1-beta)*v;
(phi_h*(1+gamma)+phi_b*(1-gamma-tao*(beta-gamma)))*h_b=-gamma*sigma_b*c_b(+1)+(gamma+tao*(beta-gamma))*q(+1)-q+sigma_b*c_b+phi_h*(h_b(-1)+gamma*h_b(+1))+tao*(beta-gamma)*(eta_l+pi(+1)-r)+(1-gamma-tao*(beta-gamma))*v;

//product functon
y_c=z_c+tao_c*k_c(-1)+(1-tao_c)*(alpha*n_s+(1-alpha)*n_b);
y_h=z_h+tao_h*k_h(-1)+(1-tao_h)*(alpha*n_s+(1-alpha)*n_b);

l=q(+1)+pi(+1)+h_b-r;

n_b=k_c(-1)+r_k-w;
n_b=k_h(-1)+r_k-w;

//inflation
pi=beta*pi(+1)+(1-beta*epsilon_c)*(1-epsilon_c)*((1-tao_c)*w+tao_c*r_k-z_c+mu_c)/epsilon_c;
q=(q(-1)-pi)/(1+beta)+beta*(q(+1)+pi(+1))/(1+beta)+(1-beta*epsilon_h)*(1-epsilon_h)*(mu_h+(1-tao_h)*w+tao_h*r_k-q)/((1+beta)*epsilon_h); 

y_h=g_hs*(h_s-h_s(-1))+g_hb*(h_b-h_b(-1));
y_c=g_cs*c_s+g_cb*c_b+g_ci*i;
y=g_c*y_c+g_h*y_h;
k=g_kc*k_c+g_kh*k_h;

//monetary policy
r=gamma_r*r(-1)+(1-gamma_r)*(gamma_pi*pi(-1)+gamma_y*y(-1))+sigma_r*e_r;

//shocks
z_c=rho_zc*z_c(-1)+sigma_zc*e_zc;
z_h=rho_zh*z_h(-1)+sigma_zh*e_zh;
v=rho_v*v(-1)+sigma_v*e_v;
mu_c=rho_mu_c*mu_c(-1)+sigma_mu_c*e_mu_c;
mu_h=rho_mu_h*mu_h(-1)+sigma_mu_h*e_mu_h;

end;

shocks;              
var e_zc;
stderr 1;

var e_zh;
stderr 1;

var e_v;
stderr 1;

var e_mu_c;
stderr 1;

var e_mu_h;
stderr 1;

end;

check;

stoch_simul(irf=30,periods=1000,replic=1000) y y_c y_h c_s c_b h_s h_b k i k_c k_h r n_s n_b pi w q eta eta_l l;