close all;

var 		e pi piw q r y yp yw z
            obs_e obs_pi obs_q obs_r obs_y;
varexo 		eps_piw eps_q eps_r eps_yw eps_z;
parameters 	alpha kapa psi_1 psi_2 psi_3 r_ss tau  
			rho_r rho_q rho_piw rho_yw rho_z
			sigma_piw sigma_r sigma_q sigma_yw sigma_z;

//LS2007 - Canada			
alpha 			= 0.11;
kapa 			= 0.32;
psi_1			= 1.30;
psi_2			= 0.23;
psi_3			= 0.14;
tau 			= 0.31;

r_ss            = 2.52;

rho_r			= 0.69;  
rho_q 			= 0.31;
rho_piw         = 0.46;
rho_yw          = 0.97;
rho_z           = 0.42;

sigma_piw       = 2.00;
sigma_q         = 1.25;
sigma_r 		= 0.36;
sigma_yw        = 1.29;
sigma_z			= 0.84;

/* NZ paper
alpha 			= 0.211;
kapa 			= 1.578;
psi_1			= 3.719;
psi_2			= 0.189;
psi_3			= 0.075;
tau 			= 0.531;

r_ss            = 3.794;

rho_r			= 0.383;  
rho_q 			= 0.203;
rho_piw         = 0.619;
rho_yw          = 0.907;
rho_z           = 0.312;

sigma_piw       = 2.376;
sigma_q         = 0.857;
sigma_r 		= 0.838;
sigma_yw        = 2.339;
sigma_z			= 1.647;
*/
model(linear);
    #beta = exp(-r_ss/400);

	//eqn 1 - IS curve
	y = y(+1) - (tau + alpha*(2 - alpha)*(1 - tau))*(r - pi(+1)) - rho_z*z - alpha*(tau + alpha*(2 - alpha)*(1 - tau))*(q(+1) - q) + alpha*(2 - alpha)*((1 - tau)/tau)*(yw(+1) - yw);

	//eqn 2 - phillips curve
	pi = beta*pi(+1) + alpha*beta*(q(+1) - q) - alpha*(q - q(-1)) + (kapa/(tau + alpha*(2 - alpha)*(1 - tau)))*(y - yp);

	//eqn 3 - ppp
	pi = (e - e(-1)) + (1 - alpha)*(q - q(-1)) + piw;

	//eqn 4 - monetary policy
	r = rho_r*r(-1) + (1 - rho_r)*(psi_1*pi + psi_2*y + psi_3*(e - e(-1))) + eps_r;

	//eqn 5 - terms of trade law of motion
	q = q(-1) + rho_q*(q(-1) - q(-2)) + eps_q;

	//eqn 6 - world output
	yw = rho_yw*yw(-1) + eps_yw;

	//eqn 7 - world inflation
	piw = rho_piw*piw(-1) + eps_piw;

	//eqn 8 - growth rate of world technology (following LS2005)
	z = rho_z*z(-1) + eps_z;
    
	//eqn 9 - potential output
	yp = -alpha*(2 - alpha)*((1 - tau)/tau)*yw; //CHECK YWORLD!!!

    //observed variables
    obs_r   = 4*r;
    obs_pi  = 4*pi;
    obs_y   = y - y(-1) + z; //CHECK!!!
    obs_e   = e - e(-1);
    obs_q   = q - q(-1);
end;

initval;
e 			= 0;
pi			= 0;
piw         = 0;
q			= 0;
r			= 0;
y			= 0;
yp          = 0;
yw          = 0;
z			= 0;
end;
steady;

check;

shocks; 


var eps_piw; stderr sigma_piw;
var eps_q; stderr sigma_q;
var eps_r; stderr sigma_r;
var eps_yw; stderr sigma_yw;
var eps_z; stderr sigma_z;


/*
var eps_piw; stderr 2.00;
var eps_q; stderr 1.25;
var eps_r; stderr 0.36;
var eps_yw; stderr 1.29;
var eps_z; stderr 0.84;
*/

/*
var eps_piw = sigma_piw^2;
var eps_q = sigma_q^2;
var eps_r = sigma_r^2;
var eps_yw = sigma_yw^2;
var eps_z = sigma_z^2;
*/

/*
var eps_piw = sigma_piw;
var eps_q = sigma_q;
var eps_r = sigma_r;
var eps_yw = sigma_yw;
var eps_z = sigma_z;
*/
end;

stoch_simul(periods = 5000, irf=12);

%===priors===%
estimated_params;
%   parameter name      initival    lb          ub          prior shape     param1      param2      param3      param4      jscale
    alpha,                                                  beta_pdf,       0.20,       0.05;
    kapa,                                                   gamma_pdf,      0.50,       0.25;
    psi_1,                                                  gamma_pdf,      1.50,       0.50; 
    psi_2,                                                  gamma_pdf,      0.25,       0.13;
    psi_3,                                                  gamma_pdf,      0.25,       0.13;
    r_ss,                                                   gamma_pdf,      2.50,       1.00; 
    tau,                                                    beta_pdf,       0.50,       0.20;
    rho_r,                                                  beta_pdf,       0.50,       0.20; //CHECK ALTERNATIVE
    rho_q,                                                  beta_pdf,       0.40,       0.20;
    rho_piw,                                                beta_pdf,       0.80,       0.10;
    rho_yw,                                                 beta_pdf,       0.90,       0.05;
    rho_z,                                                  beta_pdf,       0.20,       0.05;
	sigma_piw,                                              inv_gamma_pdf,  0.55,       4.00;
    sigma_r,                                                inv_gamma_pdf,  0.50,       4.00;
    sigma_q,                                                inv_gamma_pdf,  1.50,       4.00;
    sigma_yw,                                               inv_gamma_pdf,  1.50,       4.00;
    sigma_z,                                                inv_gamma_pdf,  1.00,       4.00;
end;

%====observed variables====%
varobs obs_e obs_pi obs_q obs_r obs_y;

%===estimation===%
//estimation(datafile=data_canada, mh_replic=50000, mh_jscale=0.40, mh_nblocks=2, bayesian_irf);