var a rn pi yg;   %unobservables
var y int p;% pie;    %observables
var lambda kappa; %Auxilary variables
var r_f r_b; %Auxilary variables
% var psi_yr psi_yp psi_pr psi_pp; %Auxilary variables

varexo e_a e_y e_i e_p;% e_e;

parameters sigma eta theta bet zet;
parameters phi_pi phi_y;
parameters rho_a psi;
parameters mupi;%mupie;

sigma = 1; % risk aversion, from Rudebusch 2002
eta = 1; % disutility from labor supply, elasticity of labor supply?
theta = 0.75; % price rigidity
bet = 0.99; % subject discount rate
zet = 0.9; % indexation ratio
phi_pi = 1.5;             % Sensitivity of the central bank with respect to inflation.
phi_y = 0.5/4;            % Sensitivity of the central bank with respect to the output gap.
rho_a = 0.9;            % Persistence of the technology shock.
psi = (1 + eta) / (sigma + eta); % coefficient of a (in yn equation), set alpha = 0 (CRS production technology)
mupi = 0.0076; % drift of inflation
% mupie = 0.0018; % drift of inflation expectation, measurement error


model;
% Auxilary variables

lambda = (1 - theta) * (1 - theta * bet) / theta;
kappa = lambda * (sigma + eta);
r_f = bet/(1 + zet * bet);
r_b = zet/(1 + zet * bet);

%(1 - r_f * psi_pp - rho_a * r_f) * psi_pr = kappa * psi_yr;
%(1 - r_f * psi_pp) * psi_pp = r_b + kappa * psi_yp;
%(sigma - sigma * rho_a + phi_y) * psi_yr - 1 = (sigma * psi_yp - phi_pi + rho_a + psi_pp) * psi_pp;
%(sigma + phi_y) * psi_yp = (sigma * psi_yp - phi_pi - psi_pi) * psi_pp;
%Lambda_a = 1/((1 - bet * rho_a) * (sigma * (1 - rho_a) + phi_y) + kappa * (phi_pi + phi_y));
%ksi = - (1 - sigma * (1 - rho_a) * (1 - bet * rho_a) * Lambda_a)  / (sigma * (1 - rho_a) * (1 - bet * rho_a) * Lambda_a);
    
% State transition comes from the following

a = rho_a * a(-1) + e_a;                        % Eq. 7: Technology shocks follow an AR(1) process with persistence rho_a.
rn =  sigma * psi * (rho_a - 1)* a;              % Eq. 3: The evolution of the natural rate of interest.
yg = yg(+1) - 1 / sigma * ( int - pi(+1) - rn );    % Eq. 1: The Dynamic IS equation.
pi = r_f * pi(+1) + r_b * pi(-1) + kappa * yg;      % Eq. 2: The New Keynesian Philips Curve.

//yg = (1 - bet * rho_a) * Lambda_a * rn;
//pi = kappa * Lambda_a * rn;

% Measurement equation
 
y = 4 * (psi * (a - a(-1)) + yg - yg(-1) + e_y); 
int = phi_y * yg + phi_pi * pi + e_i;
p = 4 * (pi + mupi + e_p);
% pie = 4 * (-kappa/r_f * yg + 1/r_f * pi - r_b/r_f * pi(-1) + mupi + mupie + e_e);

end;


%--------------------------------------------------------------------------
% 3. Initial Values
%--------------------------------------------------------------------------
 
resid(1) // display the residuals of the static equations of the model
         // using the values given for the endogenous in the last initval or endval block 
         // (or the steady state file if you provided one, see section Steady state)
steady (solve_algo=2);
check;

%--------------------------------------------------------------------------
% 4. Shocks and solution
%--------------------------------------------------------------------------
 
 
shocks;
var e_a; stderr 0.01;
var e_y; stderr 0.0012;
var e_i; stderr 0.00;
var e_p; stderr 0.0047;
%var e_e; stderr 0.00;

end;

%--------------------------------------------------------------------------
% 5. Estimation
%--------------------------------------------------------------------------

estimated_params;

mupi, normal_pdf, 0.0076, 0.001;
zet, gamma_pdf, 0.5, 0.2;

stderr e_a, inv_gamma_pdf, 0.0085, 0.001;
stderr e_y, inv_gamma_pdf, 0.005, 0.001;
stderr e_i, inv_gamma_pdf, 0.005, 0.001;
stderr e_p, inv_gamma_pdf, 0.0076, 0.002;
%stderr e_e, inv_gamma_pdf, 0.005, 0.001;

corr e_a,e_y, normal_pdf, 0, 0.0025;
corr e_a,e_i, normal_pdf, 0, 0.0025;
corr e_a,e_p, normal_pdf, 0, 0.0025;
corr e_y,e_i, normal_pdf, 0, 0.0025;
corr e_y,e_p, normal_pdf, 0, 0.0025;
corr e_i,e_p, normal_pdf, 0, 0.0025;

end;


varobs y int p; %pie;

estimation(mode_check, datafile=data_DNKWO, mh_nblocks = 2, mode_compute=4, mh_replic=20000, mh_jscale=0.5);