var c w n m r pi h epsilon_h q_h y_h k_h q_l l n_h k_g n_c k_c y_c z i_h i_c t b g i_g y i g_b epsilon_r;
var g_obs m_obs y_obs y_h_obs i_obs i_g_obs; 
varexo e_h e_l e_z e_r e_g e_b ;
parameters sigma nu psi beta phi delta_h phi_h mu tau gamma delta epsilon phi_p alpha eta phi_b phi_g gamma_r gamma_pi gamma_y
rho_h rho_l  rho_z rho_r rho_g rho_b ;

%fixed parameters
beta = 0.98;
delta_h = 0.01;
delta = 0.025;
epsilon = 4.61;
gamma_r = 0.76;
gamma_pi = 1.37;
gamma_y = 0.125;

%parameters to be estimated
sigma = 0.8;
nu = 1.5;
psi = 0.8;
phi = 0.8;
phi_h = 0.1;
mu = 0.3;
tau = 0.3;
gamma = 0.1;
phi_p = 10;
alpha = 0.6;
eta = 0.03;
phi_b = 0.1;
phi_g = 0.1;
rho_h = 0.85;
rho_l = 0.85;
rho_z = 0.85;
rho_r = 0.85;
rho_g = 0.85;
rho_b = 0.85;

model(linear);

%steady-state conditions
#t_star = 0.8061;
#m_star = 2.7115;
#b_star= 1.3729;
#i_star = 0.5953;
#g_star = 0.2223;
#y_hstar = 0.0535;
#i_hstar = 0.1962;
#i_gstar = 0.1404;
#n_hstar = 0.0138;
#ql_star=1-t_star;
#c_star=1-i_star-g_star;
#y_cstar=1-y_hstar;
#i_cstar=1-i_hstar;
#g_bstar=1-i_gstar;
#n_cstar=1-n_hstar;

%log-linearized equilibrium conditions
nu*n=-sigma*c+w;

m=sigma/psi*c-beta/(psi*(1-beta))*r;

r=-sigma*(c-c(+1))+pi(+1);

(phi*(1-beta*(1-delta_h))+phi_h*(1+beta))*h-(1-beta*(1-delta_h))*epsilon_h=sigma*c-q_h+phi_h*(h(-1)+beta*h(+1))+beta*(1-delta_h)*(q_h(+1)-sigma*c(+1));

y_h(+1)=k_h-q_h(+1)+sigma/(beta*(delta-1)+1)*(c(+1)-c);

q_l=q_h+y_h-l;

q_h+y_h-n_h=w;

y_h=mu*k_h(-1)+tau*l+(1-mu-tau)*n_h+gamma*k_g(-1);

i_h=1/delta*k_h-(1-delta)/delta*k_h(-1);

sigma*c=sigma*beta*(1-delta)*c(+1)-(1-beta*(1-delta))*((1+nu)*n_c(+1)-k_c);

pi=beta/(1+beta)*pi(+1)+pi(-1)/(1+beta)-(1-epsilon)/(phi_p*(1+beta))*((1+nu)*n_c-y_c+sigma*c);

y_c=z+alpha*k_c(-1)+(1-alpha)*n_c+eta*k_g(-1);

i_c=1/delta*k_c-(1-delta)/delta*k_c(-1);

t_star*t+b_star*(b-b(-1))+m_star*(m-m(-1))+ql_star*(q_l+l)+(b_star+m_star)*pi=g;

i_g=1/delta*k_g-(1-delta)/delta*k_g(-1);

t=phi_b*b(-1)+phi_g*g;

y_c=c_star*c+i_star*i+g_star*g;

y_h=1/delta_h*h-(1-delta_h)/delta_h*h(-1);

y=y_cstar*y_c+y_hstar*y_h;

i=i_cstar*i_c+i_hstar*i_h;

g=g_bstar*g_b+i_gstar*i_g;

n=n_hstar*n_h+n_cstar*n_c;

r=gamma_r*r(-1)+(1-gamma_r)*(gamma_pi*pi(-1)+gamma_y*y(-1))+epsilon_r;

%shocks
epsilon_h=rho_h*epsilon_h(-1)+e_h;
l=rho_l*l(-1)+e_l;
z=rho_z*z(-1)+e_z;
epsilon_r=rho_r*epsilon_r(-1)+e_r;
i_g=rho_g*i_g(-1)+e_g;
g_b=rho_b*g_b(-1)+e_b;

%observation equations
g_obs = g;
m_obs = m;
y_obs = y;
y_h_obs = y_h; 
i_obs = i;
i_g_obs = i_g;

end;

check;
steady;
model_diagnostics ;
shocks;
var e_h ; stderr 0.33;
var e_l  ; stderr 0.55;
var e_z  ;   stderr 0.15;
var e_r  ;  stderr 0.45;
var e_g ; stderr 0.95;
var e_b  ;  stderr 0.11;

end;

estimated_params;
sigma, normal_pdf, 1, 0.5;
nu, normal_pdf, 1, 0.5;
psi, normal_pdf, 1, 0.5;
phi, normal_pdf, 1, 0.5;
phi_h, inv_gamma_pdf, 0.01, 2;
gamma, beta_pdf, 0.1, 0.1;
eta, beta_pdf, 0.03, 0.1;
phi_p, inv_gamma_pdf, 10, 2;
mu, beta_pdf, 0.3, 0.01;
tau, beta_pdf, 0.3, 0.01;
alpha, beta_pdf, 0.627, 0.01;
phi_b, inv_gamma_pdf, 0.3, 1;
phi_g, inv_gamma_pdf, 0.2, 1;

rho_h,  beta_pdf, 0.7, 0.1;
rho_l,  beta_pdf, 0.7, 0.1;
rho_z,  beta_pdf, 0.7, 0.1;
rho_r,  beta_pdf, 0.7, 0.1;
rho_g,  beta_pdf, 0.7, 0.1;
rho_b,  beta_pdf, 0.7, 0.1;

stderr e_h,  inv_gamma_pdf, 0.2, 1;
stderr e_l,  inv_gamma_pdf, 0.2, 1;
stderr e_z,  inv_gamma_pdf, 0.2, 1;
stderr e_r,  inv_gamma_pdf, 0.2, 1;
stderr e_g,  inv_gamma_pdf, 0.2, 1;
stderr e_b,  inv_gamma_pdf, 0.2, 1;
end;

varobs g_obs m_obs y_obs y_h_obs i_obs i_g_obs;

estimated_params_init(use_calibration);
end;

estimation(datafile = da,
%mode_file=newmodel_mode, 
%moments_varendo, 
mode_check,
mode_compute =4,
mcmc_jumping_covariance=identity_matrix,
%presample=4,
%prefilter=0,
mh_replic = 10000,
mh_nblocks=1,
mh_jscale=0.40,
mh_drop=0.1);
identification;