close all;

var eg_hat y_hat c_hat u_hat Ncap_hat k_hat kp_hat mu_hat z_hat i_hat
 lam_hat eb_hat R_hat pi_hat upsilon_hat rk_hat N_hat xi_hat eta_star_hat
theta_star_hat V_hat q_hat U_hat thetacap_hat zeta_hat
s_hat m_hat x_hat w_hat M_hat emp_hat gy_hat gc_hat gi_hat gw_hat;


varexo emu eb ez ezeta eg emp ethetastar eetastar;

parameters delta h phiu2 HC phiI alpha chi phiP theta sigma rho phiw e
rhor rhopi rhoy rhomu rhoeb rhozeta rhoz
 rhog rhomp tau phiv rhoetastar rhothetastar 
Uss qss piss Rss zss gy ;

// fixed parameters
delta=0.025 ; % capital depreciation rate
alpha=0.33;   % capital share
theta=6;    % elasticity of substituation btw goods
chi=0.01;   % backward-looking price setting
tau=0.4;   % replacement ratio
rho=0.085;   %job destruction rate
sigma=0.65;   % elasticity of match to unemp
qss=0.7;     % prob. to fill a vacancy within a quater
Uss=0.0578;   % unemp. rate
zss=1.0039;  % quarterly groos growth rate
piss=1.0088; % quaterly gross inflation rate
Rss=1.0139;  % quarterly gross nominal interest rate
gy=0.2;       %gy=g/y
// estimated parameters initialisation
phiv=0.5;   
HC=0.005;
h=0.7;
phiI=5;
phiu2=0.5;
phiP=60;
phiw=150;
e=0.5;
rhor=0.7;
rhopi=1.5;
rhoy=0.5;
rhoz=0.3;
rhomp=0.5;
rhomu=0.5;
rhoeb=0.5;
rhozeta=0.5;
rhothetastar=0.5;
rhoetastar=0.5;
rhog=0.5;


model(linear); 
#Nss=1-Uss;
#Vss=rho*Nss/qss;
#mss=qss*Vss;
#beta=zss*piss/Rss;
#khi=1-rho;
#egss=1/(1-gy);
#xiss=(theta-1)/theta;
#rkss=(zss/beta)-1+delta;
#phiu1=rkss;
#kys=alpha*xiss/rkss;        %%kys=kss/yss
#iks=zss-1+delta;           %% iks=iss/kss
#iys=iks*kys;                %iys=iss/yss
#cys=(1/egss)-(HC)-iys;     %cys=css/yss
#Sss=1-khi*Nss;
#zeta=qss*(Vss/Sss)^(sigma);
#s=zeta*(Vss/Sss)^(1-sigma);  % wys=wss/yss
#wys=(1/Nss)*(xiss*(1-alpha)-(1-beta)*khi*2*(HC)/rho);
#upsilonss=(xiss*(1-alpha)+(1+beta*khi*s/rho)*2*(HC))/(wys*Nss);
#etass=(1-tau)/(upsilonss-tau);
#Mss=etass/(1-etass);



%(1) yt
(cys+iys)*y_hat=cys*c_hat+iys*i_hat+(phiu1*kys)*u_hat+(1/egss)*eg_hat+2*HC*Ncap_hat;
%(2) kt
k_hat=u_hat+kp_hat(-1)-z_hat;
%(3) kp
zss*kp_hat=(1-delta)*kp_hat(-1)-(1-delta)*z_hat+(zss-1+delta)*mu_hat+(zss-1+delta)*i_hat;
%(4) lam
lam_hat=eb_hat+R_hat+lam_hat(+1)-pi_hat(+1)-z_hat(+1);
%(5) ct
lam_hat=((beta*h*zss)/((zss-beta*h)*(zss-h)))*c_hat(+1)-(((zss^2)+beta*h^2)/((zss-beta*h)*(zss-h)))*c_hat+((h*zss)/((zss-beta*h)*(zss-h)))*c_hat(-1)
+((beta*h*zss)/((zss-beta*h)*(zss-h)))*z_hat(+1)-((h*zss)/((zss-beta*h)*(zss-h)))*z_hat;
%(6) rk
rk_hat=(phiu2/phiu1)*u_hat;
%(7) i
upsilon_hat=((1+beta)*(phiI*zss^2))*i_hat+(phiI*zss^2)*z_hat-
(phiI*zss^2)*i_hat(-1)-mu_hat-(beta*phiI*zss^2)*i_hat(+1)-(beta*phiI*zss^2)*z_hat(+1);
%(8) upsilon
upsilon_hat=lam_hat(+1)-lam_hat-z_hat(+1)+((1-delta)*beta*zss^(-1))*upsilon_hat(+1)+(beta*(zss^(-1))*rkss)*rk_hat(+1);
%(9) u
y_hat=alpha*k_hat+(1-alpha)*N_hat;
%(10) xi
xi_hat=rk_hat-y_hat+k_hat;
%(11) pi
pi_hat=(chi/(1+beta*chi))*pi_hat(-1)+(beta/(1+beta*chi))*pi_hat(+1)+(1/(1+beta*chi))*((theta-1)/phiP)*xi_hat-theta_star_hat;
%(12) N
N_hat=khi*N_hat(-1)+(1-khi)*q_hat+(1-khi)*V_hat;
%(13) U
U_hat=-(Nss/(1-Nss))*N_hat;
%(14) Thetacap
thetacap_hat=V_hat+(khi*Nss/Sss)*N_hat(-1);
%(15) q
q_hat=zeta_hat-sigma*thetacap_hat;
%(16) s
s_hat=zeta_hat+(1-sigma)*thetacap_hat;
%(17) Ncap_hat
Ncap_hat=((phiv*Vss)/(phiv*Vss+(1-phiv)*mss))*V_hat+(((1-phiv)*mss)/(phiv*Vss+(1-phiv)*mss))*m_hat-N_hat;
%(18) m_hat
m_hat=q_hat+V_hat;
%(19) x_hat
x_hat=m_hat-N_hat;
%(20) V
(2*2*HC*khi*Ncap_hat/rho)=(2*HC/rho)*x_hat+(1-alpha)*xiss*xi_hat+(wys*Nss-(beta*khi/rho)*2*HC)*(y_hat-N_hat)
-(wys*Nss)*w_hat+(beta*khi*2*HC/rho)*(lam_hat(+1)-lam_hat+y_hat(+1)-N_hat(+1)+2*Ncap_hat(+1)-x_hat(+1));
%(21) w_hat
(wys*Nss/etass)*w_hat=(1-alpha)*xiss*xi_hat+((1-alpha)*xiss+2*HC)*(y_hat-N_hat)+2*2*HC*Ncap_hat
+((beta*khi*s*2*HC)/rho)*(s_hat(+1)+lam_hat(+1)-lam_hat+2*Ncap_hat(+1)+y_hat(+1)-m_hat(+1))
-(wys*Nss-(1-alpha)*xiss-(1+(beta*khi/rho))*2*HC)*M_hat-beta*khi*(1-s)*2*HC*M_hat(+1)/rho;
%(22) M_hat
M_hat=eta_star_hat+(beta*khi*phiw*(1/wys))*z_hat(+1)+(beta*khi*phiw*(1/wys))*pi_hat(+1)+(beta*khi*phiw*(1/wys))*w_hat(+1)
-((phiw*(1/wys))*(1+beta*khi))*w_hat
-((phiw*(1/wys))*(1+beta*khi*e))*pi_hat-(phiw*(1/wys))*z_hat+(phiw*(1/wys))*w_hat(-1)+(phiw*(1/wys))*e*pi_hat(-1);
%(23) R
R_hat=rhor*R_hat(-1)+((1-rhor)*rhopi/4)*(pi_hat+pi_hat(-1)+pi_hat(-2)+pi_hat(-3))
+((1-rhor)*rhoy/4)*(gy_hat+gy_hat(-1)+gy_hat(-2)+gy_hat(-3))+emp_hat;
%(24)gy_hat
gy_hat=y_hat-y_hat(-1)+z_hat;
%(25) gc_hat
gc_hat=c_hat-c_hat(-1)+z_hat;
%(26) gi_hat
gi_hat=i_hat-i_hat(-1)+z_hat;
%(27) gw_hat
gw_hat=w_hat-w_hat(-1)+z_hat;
%%shocks
%(28) mu_hat
mu_hat=rhomu*mu_hat(-1)+emu;
%(29) eb_hat
eb_hat=rhoeb*eb_hat(-1)+eb;
%(30) z_hat
z_hat=rhoz*z_hat(-1)+ez;
%(31) zeta
zeta_hat=rhozeta*zeta_hat(-1)+ezeta;

%(34) eg_hat
eg_hat=rhog*eg_hat(-1)+eg;
%(35) emp
emp_hat=rhomp*emp_hat(-1)+emp;

%(36) theta_star_hat
%theta_star_hat=(1/((1+beta*chi)*phiP))*th;
%(37) eta_star_hat
%eta_star_hat=(1/(1-etass))*eta_hat;

theta_star_hat=rhothetastar*theta_star_hat(-1)-ethetastar;
eta_star_hat=rhoetastar*eta_star_hat(-1)+eetastar;
end;

steady_state_model;
eg_hat=0;
 y_hat=0;
 c_hat=0;
 u_hat=0;
 Ncap_hat=0;
 k_hat=0;
 kp_hat=0;
 mu_hat=0;
 z_hat=0;
 i_hat=0;
 lam_hat=0;
 eb_hat=0;
 R_hat=0;
 pi_hat=0;
upsilon_hat=0;
 rk_hat=0;
 N_hat=0;
 xi_hat=0;
 eta_star_hat=0;
theta_star_hat =0;
V_hat =0;
q_hat =0;
U_hat=0;
 thetacap_hat=0;
 zeta_hat=0;
s_hat=0;
 m_hat=0;
 x_hat=0;
 w_hat=0;
 M_hat=0;
 emp_hat=0;
 gy_hat=0;
 gc_hat=0;
 gi_hat=0;
 gw_hat=0;

end;

shocks;

var ez;
stderr 0.0127;
var emp;
stderr 0.0021;
var emu;
stderr 0.0617;
var eb;
stderr 0.0082;
var ezeta;
stderr 0.0231;
var ethetastar;
stderr 0.0010;
var eetastar;
stderr 0.0010;
var eg;
stderr 0.0047;
end;
steady ;
check ;
stoch_simul(order=2,nograph,pruning,irf=20,hp_filter=1600);


varobs gy_hat gc_hat gi_hat gw_hat V_hat U_hat pi_hat R_hat;

estimated_params;
phiv, beta_pdf, 0.5, 0.2;
h, beta_pdf, 0.7, 0.1;
phiI,gamma_pdf,5,1;
phiu2,gamma_pdf,0.5,0.1;
phiP,gamma_pdf,60,10;
phiw,gamma_pdf,150,25;
e,beta_pdf,0.5,0.2;
rhor,beta_pdf,0.7,0.1;
rhopi,gamma_pdf,1.5,0.1;
rhoy,gamma_pdf,0.5,0.1;
HC,gamma_pdf,5,1;

rhoz, beta_pdf, 0.3, 0.1;
rhomp, beta_pdf, 0.5, 0.2;
rhomu,beta_pdf, 0.5, 0.2;
rhoeb, beta_pdf, 0.5, 0.2;
rhozeta, beta_pdf, 0.5, 0.2;
rhothetastar, beta_pdf, 0.5, 0.2;
rhoetastar, beta_pdf, 0.5, 0.2;
rhog, beta_pdf, 0.7, 0.1;


stderr ez, inv_gamma_pdf, 0.1,3;
stderr emp, inv_gamma_pdf, 0.1,3;
stderr emu, inv_gamma_pdf, 0.1,3;
stderr eb, inv_gamma_pdf, 0.1,3;
stderr ezeta, inv_gamma_pdf, 0.1,3;
stderr ethetastar, inv_gamma_pdf, 0.1,3;
stderr eetastar, inv_gamma_pdf, 0.1,3;
stderr eg, inv_gamma_pdf, 0.1,3;

end;



estimation(datafile=data_furlanetto, mh_jscale=0.2,mh_drop=0.5, mode_compute=4,mh_replic=500000,conf_sig=.95);