%----------------------------------------------------------
% Repeat estimation with simulated data using a New Keynesian DSGE model
% without zero lower bound on the nominal interest rate
%----------------------------------------------------------

clear all;
clc;

nrep = 500;           % number of replication

npara = 13;         % number of structural parameters (excluding s.d. of shocks)
nshock = 3;         % number of shocks

postmean1 = zeros(npara,nrep);
postmean2 = zeros(nshock,nrep);
posthpdinf1 = zeros(npara,nrep);
posthpdinf2 = zeros(nshock,nrep);
posthpdsup1 = zeros(npara,nrep);
posthpdsup2 = zeros(nshock,nrep);
irfmean = zeros(15,9,nrep);
irfhpdinf = zeros(15,9,nrep);
irfhpdsup = zeros(15,9,nrep);

load('simyy.mat');

for rep = 1:nrep
    
    y_obs = simyy(:,1,rep);
    pi_obs = simyy(:,2,rep);
    r_obs = simyy(:,3,rep);
    
    save('simdata','y_obs','pi_obs','r_obs');

    dynare NKmodel_est noclearall;
    
    postmean1(:,rep) = [
        oo_.posterior_mean.parameters.sig;
        oo_.posterior_mean.parameters.chi;
        oo_.posterior_mean.parameters.gam;
        oo_.posterior_mean.parameters.xi;
        oo_.posterior_mean.parameters.iota;
        oo_.posterior_mean.parameters.phi_pi;
        oo_.posterior_mean.parameters.phi_y;
        oo_.posterior_mean.parameters.phi_r;
        oo_.posterior_mean.parameters.piss;
        oo_.posterior_mean.parameters.rss;
        oo_.posterior_mean.parameters.rho_d;
        oo_.posterior_mean.parameters.rho_p;
        oo_.posterior_mean.parameters.rho_r];

    postmean2(:,rep) = [
        oo_.posterior_mean.shocks_std.e_d;
        oo_.posterior_mean.shocks_std.e_p;
        oo_.posterior_mean.shocks_std.e_r];
    
    posthpdinf1(:,rep) = [
        oo_.posterior_hpdinf.parameters.sig;
        oo_.posterior_hpdinf.parameters.chi;
        oo_.posterior_hpdinf.parameters.gam;
        oo_.posterior_hpdinf.parameters.xi;
        oo_.posterior_hpdinf.parameters.iota;
        oo_.posterior_hpdinf.parameters.phi_pi;
        oo_.posterior_hpdinf.parameters.phi_y;
        oo_.posterior_hpdinf.parameters.phi_r;
        oo_.posterior_hpdinf.parameters.piss;
        oo_.posterior_hpdinf.parameters.rss;
        oo_.posterior_hpdinf.parameters.rho_d;
        oo_.posterior_hpdinf.parameters.rho_p;
        oo_.posterior_hpdinf.parameters.rho_r];

    posthpdinf2(:,rep) = [
        oo_.posterior_hpdinf.shocks_std.e_d;
        oo_.posterior_hpdinf.shocks_std.e_p;
        oo_.posterior_hpdinf.shocks_std.e_r];

    posthpdsup1(:,rep) = [
        oo_.posterior_hpdsup.parameters.sig;
        oo_.posterior_hpdsup.parameters.chi;
        oo_.posterior_hpdsup.parameters.gam;
        oo_.posterior_hpdsup.parameters.xi;
        oo_.posterior_hpdsup.parameters.iota;
        oo_.posterior_hpdsup.parameters.phi_pi;
        oo_.posterior_hpdsup.parameters.phi_y;
        oo_.posterior_hpdsup.parameters.phi_r;
        oo_.posterior_hpdsup.parameters.piss;
        oo_.posterior_hpdsup.parameters.rss;
        oo_.posterior_hpdsup.parameters.rho_d;
        oo_.posterior_hpdsup.parameters.rho_p;
        oo_.posterior_hpdsup.parameters.rho_r];

    posthpdsup2(:,rep) = [
        oo_.posterior_hpdsup.shocks_std.e_d;
        oo_.posterior_hpdsup.shocks_std.e_p;
        oo_.posterior_hpdsup.shocks_std.e_r];
    
    irfmean(:,:,rep) = [
        oo_.PosteriorIRF.dsge.Mean.y_e_d ...
        oo_.PosteriorIRF.dsge.Mean.pi_e_d ...
        oo_.PosteriorIRF.dsge.Mean.r_e_d ...
        oo_.PosteriorIRF.dsge.Mean.y_e_p ...
        oo_.PosteriorIRF.dsge.Mean.pi_e_p ...
        oo_.PosteriorIRF.dsge.Mean.r_e_p ...
        oo_.PosteriorIRF.dsge.Mean.y_e_r ...
        oo_.PosteriorIRF.dsge.Mean.pi_e_r ...
        oo_.PosteriorIRF.dsge.Mean.r_e_r];

    irfhpdinf(:,:,rep) = [
        oo_.PosteriorIRF.dsge.HPDinf.y_e_d ...
        oo_.PosteriorIRF.dsge.HPDinf.pi_e_d ...
        oo_.PosteriorIRF.dsge.HPDinf.r_e_d ...
        oo_.PosteriorIRF.dsge.HPDinf.y_e_p ...
        oo_.PosteriorIRF.dsge.HPDinf.pi_e_p ...
        oo_.PosteriorIRF.dsge.HPDinf.r_e_p ...
        oo_.PosteriorIRF.dsge.HPDinf.y_e_r ...
        oo_.PosteriorIRF.dsge.HPDinf.pi_e_r ...
        oo_.PosteriorIRF.dsge.HPDinf.r_e_r];
    irfhpdsup(:,:,rep) = [
        oo_.PosteriorIRF.dsge.HPDsup.y_e_d ...
        oo_.PosteriorIRF.dsge.HPDsup.pi_e_d ...
        oo_.PosteriorIRF.dsge.HPDsup.r_e_d ...
        oo_.PosteriorIRF.dsge.HPDsup.y_e_p ...
        oo_.PosteriorIRF.dsge.HPDsup.pi_e_p ...
        oo_.PosteriorIRF.dsge.HPDsup.r_e_p ...
        oo_.PosteriorIRF.dsge.HPDsup.y_e_r ...
        oo_.PosteriorIRF.dsge.HPDsup.pi_e_r ...
        oo_.PosteriorIRF.dsge.HPDsup.r_e_r];
    
    save('results_temp');

    close all;
    clc;

end

save('results_all');