var
W Q R N I Y C K B m M P Gamma G_Gamma G_Z//--- 13 Endogenous state variables
Psi Phi bb A Lamdaa Lamdaz Z Ksi Theta //--- 9 variables appearing in shock pricesses
DLogGDP  DLogC  DLogI  DLogP  DLogM //--- 5 Observed data
;
varexo Eps_psi Eps_phi Eps_b Eps_a Eps_lamdaa Eps_lamdaz Eps_z Eps_ksi Eps_theta;
parameters
h omega  //---  estimated parameters
alpha beta  g_gamma g_z deta r theta_bar iy  //---  calibrated parameters
rho_psi rho_phi rho_b  rho_a rho_lamdaa rho_lamdaz rho_z rho_ksi rho_theta;//---  parameters related to shock processes
rho_psi = 0;
rho_phi = 0;
rho_b = 0;
rho_a = 0;
rho_lamdaa = 0;
rho_lamdaz = 0;
rho_z = 0;
rho_ksi = 0;
rho_theta = 0;
h = 0;
omega = 2;
alpha = 0.41;
beta = 0.9702;
deta = 0.025;
g_gamma = 1.02477375;
g_z = 1.73261;
theta_bar = 1;
r = 1.01626;
iy = 0.352347;
model(linear);
# phiMgamma = 1-beta/g_gamma;
# Q_bar = 1/beta;
# ksi_bar = ((g_gamma*g_z*Q_bar)/beta-r-(1-deta)*Q_bar-(Q_bar-1)*r)/Q_bar;
# yk = r/(alpha*g_gamma*g_z);
# ik = iy*yk;
# cy = 1- iy;
# bi = 1-r/ik-ksi_bar*Q_bar/ik;

Gamma =(g_gamma/(g_gamma-beta*h))*((h*(C(-1) - G_Gamma)/(g_gamma - h)) -
g_gamma*C/g_gamma - h) - beta*h/(g_gamma - beta*h)*(h*C/(g_gamma - h) - g_gamma*(C(+1) + G_Gamma(+1))/g_gamma - h);
Psi = Gamma + Y - N;
Gamma + G_Gamma(+1)/g_gamma = phiMgamma*(Phi - M) + beta*(Gamma(+1)) ;
Q_bar*g_gamma*g_z*(Q + Gamma + G_Gamma(+1) + G_Z(+1)) = beta*(r*R(+1) + (1 - deta)*Q_bar*Q(+1) +
(Q_bar - 1)*(r*R(+1) + ksi_bar*Q_bar*(Ksi(+1) + Q(+1))) + Gamma(+1)*(r + (1 - deta)*Q_bar + (Q_bar - 1)*(r + ksi_bar*Q_bar)));
B = beta*Q_bar*theta_bar*(Gamma(+1) + Q(+1) - Gamma - P(+1) + B(+1) + m - m(-1) + Theta);
P = (1 + beta)*omega*g_gamma^2*g_z^2*I + omega*g_gamma^2*g_z^2*(G_Gamma + G_Z) -
omega*g_gamma^2*g_z^2*I(-1) - beta*omega*g_gamma^2*g_z^2*(I(+1) + G_Gamma(+1) + G_Z(+1));
W = Y - N;
Y = alpha*K + (1 - alpha)*N;
K(+1) = (1 - deta)*(K - G_Gamma(+1) - G_Z(+1))/g_gamma*g_z + ik*(I - G_Gamma(+1) - G_Z(+1))/g_gamma*g_z;
Y = cy*C + iy*I;
P + I = r*(R + K)/yk + ksi_bar*Q_bar*(Ksi + Q + K)/yk + bi*B;
m = theta_bar*(Theta(-1) + m(-1));
bb = M - M(-1);
G_Gamma = (alpha/(1 - alpha))*(Lamdaz + Z - Z(-1))+(1/(1-alpha))*(A + Lamdaa - Lamdaa(-1));
G_Z = Lamdaz + Z - Z(-1);

Psi = rho_psi*Psi(-1) + Eps_psi;
Phi = rho_phi*Phi(-1) + Eps_phi;
bb = rho_b*bb(-1) + Eps_b;
A = rho_a*A(-1) + Eps_a;
Lamdaa = rho_lamdaa*Lamdaa(-1) + Eps_lamdaa;
Lamdaz = rho_lamdaz*Lamdaz(-1) + Eps_lamdaz;
Z = rho_z*Z(-1) + Eps_z;
Ksi = rho_ksi*Ksi(-1) + Eps_ksi;
Theta = rho_theta*Theta(-1) + Eps_theta;
DLogP =  (P -P(-1)+ log(g_z)+G_Z);
DLogI = (I - I(-1) + log(g_gamma) + G_Gamma);
DLogC = (C-C(-1) + log(g_gamma) + G_Gamma);
DLogM = (M - M(-1) + log(g_gamma) + G_Gamma);
DLogGDP = iy*DLogI + (1-iy)*DLogC;
end;
shocks;
var Eps_psi; stderr 0;
var Eps_phi; stderr 0;
var Eps_b; stderr 0;
var Eps_a; stderr 0;
var Eps_lamdaa; stderr 0;
var Eps_lamdaz; stderr 0;
var Eps_z; stderr 0;
var Eps_ksi; stderr 0;
var Eps_theta; stderr 0;
end;
check;
estimated_params;
stderr Eps_psi,INV_GAMMA2_PDF,0.01,inf;
stderr Eps_phi,INV_GAMMA2_PDF,0.01,inf;
stderr Eps_b,INV_GAMMA2_PDF,0.01,inf;
stderr Eps_a,INV_GAMMA2_PDF,0.01,inf;
stderr Eps_lamdaa,INV_GAMMA2_PDF,0.01,inf;
stderr Eps_lamdaz,INV_GAMMA2_PDF,0.01,inf;
stderr Eps_z,INV_GAMMA2_PDF,0.01,inf;
stderr Eps_ksi,INV_GAMMA2_PDF,0.01,inf;
stderr Eps_theta,INV_GAMMA2_PDF,0.01,inf;
h,BETA_PDF,0.333333333333333,0.235702260395516;
omega,GAMMA_PDF,2,2;
rho_psi,BETA_PDF,0.85,0.1;
rho_phi,BETA_PDF,0.85,0.1;
rho_b,BETA_PDF,0.85,0.1;
rho_a,BETA_PDF,0.85,0.1;
rho_lamdaa,BETA_PDF,0.85,0.1;
rho_lamdaz,BETA_PDF,0.85,0.1;
rho_z,BETA_PDF,0.85,0.1;
rho_ksi,BETA_PDF,0.85,0.1;
rho_theta,BETA_PDF,0.85,0.1;
end;
varobs DLogP DLogI DLogC DLogM DLogGDP ;
estimation(bayesian_irf,irf=40,datafile=data_hyy,first_obs=1,nobs=83,mh_replic=2000,mh_nblocks=4,mh_jscale=1);

  error('The steady state contains NaN or Inf')