var c pih pif q s y pi  r a rp es
//pi_1 pi_2 pi_3 pi_a
pistar ystar rstar;  

varexo eps_a eps_h eps_f eps_q eps_s eps_r eps_ystar eps_pistar eps_rstar;

parameters beta alpha h sigma phi eta deltah deltaf thetah thetaf 
rhor phi1 phi2 phi3 rhoa rhoq rhos rhoys rhopis rhors lambdah lambdaf; 

beta = 0.99;
alpha = 0.45;
h = 0.6;
sigma = 1;
phi = 1.5;
eta = 1;
deltah = 0.5;
deltaf = 0.5;
thetah = 0.5;
thetaf = 0.5;
rhor = 0.5;
phi1 = 1.5;
phi2 = 0.25;
phi3 = 0.25;
rhoa = 0.5;
rhoq = 0.5;
rhos = 0.5;
rhoys = 0.5;
rhopis = 0.5;
rhors = 0.5;
// replace paramters
lambdah = (1-beta*thetah)*(1-thetah)/thetah;
lambdaf=(1-beta*thetaf)*(1-thetaf)/thetaf;

model(linear);

        // consumption euler function
        c - h*c(-1) = c(+1) - h*c - (1-h)/sigma*(r-pi(+1));

        // Domestic inflation equation
        pih - deltah*pih(-1) = beta*( pih(+1) - deltah*pih ) + 
        lambdah*( phi*y - (1+phi)*a + alpha*s + sigma/(1-h)*(c-h*c(-1)) + eps_h);
        
        // Import inflation equation
        pif - deltaf*pif(-1) = beta*( pif(+1) - deltaf*pif ) +
        lambdaf*(q + (1-alpha)*s + eps_f);
        
        // Real UIP condition
        q(+1) - q = (r - pi(+1)) - (rstar - pistar(+1)) + rp;

        // terms of trade equation
        s - s(-1) = pif - pih + es;

        // market clearing condition
        y = (1-alpha)*c + alpha*eta*q + alpha*eta*s + alpha*ystar;

        // CPI inflation definition
        pi = (1-alpha)*pih + alpha*pif;
        
        // Taylor rule
        r = rhor*r(-1) + (1-rhor)*( phi1*pi + phi2*y + phi3*q) + eps_r;
    
        // annual cpi inflation
        //pi_1 = pi(-1);
        //pi_2 = pi(-2);
        //pi_3 = pi(-3);
        //pi_a = 400*(pi - pi_3(-1));
        
        // Exogenous shocks
        a = rhoa*a(-1) + eps_a;
        rp = rhoq*rp(-1) + eps_q;
        es = rhos*es(-1) + eps_s;
        ystar = rhoys*ystar(-1) + eps_ystar;
        pistar = rhopis*pistar(-1) + eps_pistar;
        rstar = rhors*rstar(-1) + eps_rstar;
end;

shocks;  // eps_a eps_h eps_f eps_q eps_s eps_r eps_ystar eps_pistar eps_rstar;
var eps_a; stderr 1.19;
var eps_h; stderr 2.66;
var eps_f; stderr 2.66;
var eps_q; stderr 0.53;
var eps_s; stderr 1.19;
var eps_r; stderr 0.40;
var eps_ystar; stderr 1.19;
var eps_pistar; stderr 1.19;
var eps_rstar; stderr 1.19;
end;

steady;
check;
//stoch_simul(irf = 30);

//estimated_params;
    //          STARTING VALUES  &  PRIORS
    //h,          0.6, 1E-5, 0.999,   BETA_PDF,   0.60, 0.20;
    //sigma,      1, 1E-5, 10,        GAMMA_PDF,  1.00, 0.50;
    //phi,        1.5, -10, 10,       NORMAL_PDF, 1.50, 0.25;
    //eta,        1, 1E-5, 10,        GAMMA_PDF,  1.00, 0.50; 
    //deltah,     0.7, 1E-5, 0.999,   BETA_PDF,   0.70, 0.20;
    //deltaf,     0.7, 1E-5, 0.999,   BETA_PDF,   0.70, 0.20;
    //thetah,     0.5, 1E-5, 0.999,   BETA_PDF,   0.50, 0.20;
    //thetaf,     0.5, 1E-5, 0.999,   BETA_PDF,   0.50, 0.20;
    //rhoys,      0.5, 1E-5, 0.999,   BETA_PDF,   0.50, 0.20;
    //rhopis,     0.5, 1E-5, 0.999,   BETA_PDF,   0.50, 0.20;
    //rhors,      0.5, 1E-5, 0.999,   BETA_PDF,   0.50, 0.20;
    //rhoa,       0.5, 1E-5, 0.999,   BETA_PDF,   0.50, 0.20;
    //rhoq,       0.9, 1E-5, 0.999,   BETA_PDF,   0.90, 0.20;
    //rhos,       0.25, 1E-5, 0.999,  BETA_PDF,   0.25, 0.20;
    //rhor,       0.5, 1E-5, 0.999,   BETA_PDF,   0.50, 0.20;
    //phi1,       1.5, 1.0001, 10,    GAMMA_PDF,  1.50, 0.30;
    //phi2,       0.25, 1E-5, 10,     GAMMA_PDF,  0.25, 0.10;
    //phi3,       0.25, 1E-5, 10,     GAMMA_PDF,  0.25, 0.10;
    
    //stderr eps_a, 1,1E-8,20,        INV_GAMMA_PDF, 1.00, 0.40;
    //stderr eps_h, 1,1E-8,20,        INV_GAMMA_PDF, 1.00, 0.40;
    //stderr eps_f, 1,1E-8,20,        INV_GAMMA_PDF, 1.00, 0.40;
    //stderr eps_q, 1,1E-8,20,        INV_GAMMA_PDF, 2.00, 0.50;
    //stderr eps_s, 1,1E-8,20,        INV_GAMMA_PDF, 1.00, 0.40;
    //stderr eps_r, 1,1E-8,20,        INV_GAMMA_PDF, 1.00, 0.40;
    //stderr eps_ystar, 1,1E-8,20,    INV_GAMMA_PDF, 1.00, 0.40;
    //stderr eps_pistar, 1,1E-8,20,  INV_GAMMA_PDF, 1.00, 0.40;
    //stderr eps_rstar, 1,1E-8,20,    INV_GAMMA_PDF, 1.00, 0.40;

estimated_params;
    //          STARTING VALUES  &  PRIORS
    h,             BETA_PDF,   0.60, 0.20;
    sigma,             GAMMA_PDF,  1.00, 0.50;
    phi,              NORMAL_PDF, 1.50, 0.25;
    eta,               GAMMA_PDF,  1.00, 0.50; 
    deltah,       BETA_PDF,   0.70, 0.20;
    deltaf,       BETA_PDF,   0.70, 0.20;
    thetah,       BETA_PDF,   0.50, 0.20;
    thetaf,      BETA_PDF,   0.50, 0.20;
    rhoys,        BETA_PDF,   0.50, 0.20;
    rhopis,       BETA_PDF,   0.50, 0.20;
    rhors,        BETA_PDF,   0.50, 0.20;
    rhoa,        BETA_PDF,   0.50, 0.20;
    rhoq,        BETA_PDF,   0.90, 0.20;
    rhos,      BETA_PDF,   0.25, 0.20;
    rhor,         BETA_PDF,   0.50, 0.20;
    phi1,          GAMMA_PDF,  1.50, 0.30;
    phi2,            GAMMA_PDF,  0.25, 0.10;
    phi3,           GAMMA_PDF,  0.25, 0.10;
    
    stderr eps_a,       INV_GAMMA_PDF, 1.00, 0.40;
    stderr eps_h,        INV_GAMMA_PDF, 1.00, 0.40;
    stderr eps_f,        INV_GAMMA_PDF, 1.00, 0.40;
    stderr eps_q,        INV_GAMMA_PDF, 2.00, 0.50;
    stderr eps_s,       INV_GAMMA_PDF, 1.00, 0.40;
    stderr eps_r,       INV_GAMMA_PDF, 1.00, 0.40;
    stderr eps_ystar,     INV_GAMMA_PDF, 1.00, 0.40;
    stderr eps_pistar,  INV_GAMMA_PDF, 1.00, 0.40;
    stderr eps_rstar,     INV_GAMMA_PDF, 1.00, 0.40;
end;
  
varobs pif q s y pi r pistar ystar rstar;

estimation(datafile=nz_data,mode_compute=4,mh_nblocks = 2,mh_replic=1000,
mh_jscale=0.25);

// mode_compute=3
// load_mh_file
// nodiagnostic
// mode_file=MODEL1_mode