//Endogenous variable
var l c t_w w r pi t_k r_k mc a k_g k pi_tilda t_c y i i_g c_g bd ge gr or d b m tl_g ksi y_t;

//Exogenous variable
varexo epsilon_a epsilon_m epsilon_i_g epsilon_c_g epsilon_t_c epsilon_t_w epsilon_t_k epsilon_or;

//Parameters
parameters sigma_l sigma_c beta delta_k alpha gamma rho_a omega delta_kg theta sigma_m rho_m eta_m_pi eta_m_y rho_ig rho_cg rho_tk rho_tw rho_tc rho_or k1 k3 k4 k5 k6 k7 k8 k9 k10 k11 k12 k13 k14 k15 k16 k17 k18 k19 k20 k21 k22 k23 k24;

sigma_l=2.22;
sigma_c=1.52;
beta=0.96;
delta_k=0.014;
alpha=0.4278;
gamma=0.19;
rho_a=0.84;
omega=0.58;
delta_kg=0.022;
theta=4.33;
sigma_m=2.24;
rho_m=0.56;
eta_m_pi=0.9898;
eta_m_y=2.967;
rho_ig=0.5;
rho_cg=0.5;
rho_tk=0.5;
rho_tw=0.5;
rho_tc=0.5;
rho_or=0.427;
//k1=t_w_ss/(1-t_w_ss)
k1=0.111;
//k3=t_c_ss/(1+t_c_ss)
k3=0.0244;
//k4=c_ss/y_ss
k4=0.515;
//k5=i_ss/y_ss
k5=0.244;
//k6=c_g_ss/y_ss
k6=0.146;
//k7=i_g_ss/y_ss
k7=0.094;
//k8=ge_ss/bd_ss
k8=2.107;
//k9=gr_ss/bd_ss
k9=1.107;
//k10=c_g_ss/ge_ss
k10=0.41;
//k11=i_g_ss/ge_ss
k11=0.39;
//k12=c_ss/gr_ss
k12=4.1;
//k13=c_ss/(gr_ss*(1+t_c_ss))
k13=4;
//k14=t_w_ss*w_ss*l_ss/gr_ss
k14=0.1;
//k15=t_k_ss*d_ss/gr_ss
k15=0.2;
//k16=or_ss/gr_ss
k16=0.6;
//k17=y_ss*mc_ss/d_ss
k17=7.2;
//k18=y_ss/y_t_ss
k18=0.866;
//k19=or_ss/y_t_ss
k19=0.134;
//k20=c_ss/(y_ss*(1+t_c))
k20=0.502;
//k21=b_ss/tl_g_ss
k21=0.969;
//k22=nd_g_ss/tl_g_ss
k22=0.031;
//k23=bd_ss/tl_g_ss
k23=0.868;
//k24=(beta*tl_g_ss+(1-beta)*b_ss)/(beta*tl_g_ss)
k24=1.041;


//The Model

model(linear);
sigma_l*l+sigma_c*c=-k1*t_w+w;
c=c(+1)-(1/sigma_c)*(r-pi(+1));
r_k(+1)=(1/(1-beta*(1-delta_k)))*(r-pi(+1));
mc=(1-alpha)*w+alpha*r_k-a-gamma*k_g(-1);
a=rho_a*a(-1)+epsilon_a;
l=r_k-w+k(-1);
pi_tilda=((1-omega)*(1-beta*omega)/omega)*(mc+k3*t_c)+beta*pi_tilda(+1);
pi=pi_tilda+k3*(t_c-t_c(-1));
y=a+alpha*k(-1)+(1-alpha)*l+gamma*k_g(-1);
k=(1-delta_k)*k(-1)+delta_k*i;
y_t=k4*c+k5*i+k6*c_g+k7*i_g-k3*k20*t_c;
k_g=(1-delta_kg)*k_g(-1)+delta_kg*i_g;
bd=k8*ge-k9*gr;
ge=k10*c_g+k11*i_g;
gr=k15*d+k3*k12*c+k3*k13*t_c+k14*(t_w+w+l)+k15*t_k+k16*or;
d=y-theta*k3*t_c-k17*mc;
m=(sigma_c/sigma_m)*c-(beta/((1-beta)*sigma_m))*r;
y_t=k18*y+k19*or;
tl_g=k21*b+k22*m;
tl_g=tl_g(-1)+(1/beta)*k21*r(-1)+((1-beta)/beta)*k21*b(-1)-k24*pi+k23*bd;
ksi=m-m(-1)+pi;
ksi=rho_m*ksi(-1)-eta_m_pi*pi(+1)-eta_m_y*y+epsilon_m;
i_g=rho_ig*i_g(-1)-(1-rho_ig)*tl_g(-1)+epsilon_i_g;
c_g=rho_cg*c_g(-1)-(1-rho_cg)*tl_g(-1)+epsilon_c_g;
t_c=rho_tc*t_c(-1)+(1-rho_tc)*tl_g(-1)+epsilon_t_c;
t_w=rho_tw*t_w(-1)+(1-rho_tw)*tl_g(-1)+epsilon_t_w;
t_k=rho_tk*t_k(-1)+(1-rho_tk)*tl_g(-1)+epsilon_t_k;
or=rho_or*or(-1)+epsilon_or;
end;

//Steady state
initval;
l=0;
c=0;
t_w=0;
w=0;
r=0;
pi=0;
t_k=0;
r_k=0;
mc=0;
a=0;
k_g=0;
k=0;
pi_tilda=0;
t_c=0;
y=0;
i=0;
i_g=0;
c_g=0;
bd=0;
ge=0;
gr=0;
or=0;
d=0;
b=0;
m=0;
tl_g=0;
ksi=0;
y_t=0;
end;

//Blanchard kan condition check
check;
steady;

//The Exogenous shocks
shocks;
var epsilon_a;
stderr 0.01;
var epsilon_m;
stderr 0.01;
var epsilon_i_g;
stderr 0.01; 
var epsilon_c_g;
stderr 0.01;
var epsilon_t_c;
stderr 0.01;
var epsilon_t_w;
stderr 0.01;
var epsilon_t_k;
stderr 0.01;
var epsilon_or;
stderr 0.01;
end;

//Observation variables
varobs c_g i_g m or t_c t_k t_w y_t;

// DECLARE ESTIMATED PARAMETERS
estimated_params;
sigma_l, gamma_pdf, 2.22, 0.05;
sigma_c, gamma_pdf, 1.52, 0.05;
beta, beta_pdf, 0.96, 0.025, 0, 1;
alpha, beta_pdf, 0.4278, 0.025, 0,1; 
gamma, beta_pdf, 0.19, 0.025, 0,1;
rho_a, beta_pdf, 0.84, 0.025, -1,1;
omega, beta_pdf, 0.58, 0.025, 0,1;
theta, gamma_pdf, 4.33, 0.05, 1; 
sigma_m, gamma_pdf, 2.24, 0.05;
rho_m, beta_pdf, 0.56, 0.025, -1,1;
eta_m_pi, normal_pdf, 0.9898, 0.15;
eta_m_y, normal_pdf, 2.967, 0.15;
rho_ig, beta_pdf, 0.5, 0.025, -1,1;
rho_cg, beta_pdf, 0.5, 0.025, -1,1;
rho_tw, beta_pdf, 0.5, 0.025, -1,1;
rho_tk, beta_pdf, 0.5, 0.025, -1,1;
rho_tc, beta_pdf, 0.5, 0.025, -1,1;
rho_or, beta_pdf, 0.427, 0.025, -1,1;


stderr epsilon_a, INV_GAMMA_PDF,0.01,inf;
stderr epsilon_m, INV_GAMMA_PDF,0.01,inf;
stderr epsilon_i_g, INV_GAMMA_PDF,0.01,inf;
stderr epsilon_c_g, INV_GAMMA_PDF,0.01,inf;
stderr epsilon_t_c, INV_GAMMA_PDF,0.01,inf;
stderr epsilon_t_w, INV_GAMMA_PDF,0.01,inf;
stderr epsilon_t_k, INV_GAMMA_PDF,0.01,inf;
stderr epsilon_or, INV_GAMMA_PDF,0.01,inf;
end;

/*
dynare_sensitivity(lik_init=3);//, pvalue_ks=0.05, pvalue_corr=0.05);
identification(advanced=1,max_dim_cova_group=3,prior_mc=250, lik_init=4);
*/


irf_calibration;
y(1:10), epsilon_t_w, -;
y(1:10), epsilon_t_k, -;
end;

dynare_sensitivity;

/*
estimation (datafile = siasat_data, mh_jscale=0.6, mh_replic=100 ,mh_nblocks=5, mode_compute=6, lik_init=2, plot_priors=0, graph_format=pdf, bayesian_irf)l c t_w w r pi t_k r_k mc a k_g k pi_tilda t_c y i i_g c_g bd ge gr or d b m tl_g ksi y_t;
*/