% Runs the Metropolis Hastings RW chains for structural parameters and date breaks

% Note: The chain provided with the code in mhall_RBC.mat has 25000 draws for file size considerations. 
% To produce the full posterior distribution the chain length should be set to S_.Options.EST.N = 750000 
% in line 20 of run_mcmc.m

addpath functions

%% Load data and define posterior

warning('off','all') ; 

aus_data;

filename = 'hessmat.mat';
load(filename)

S_.Opt_.EST.kron = 0;   % for mcmc

post_ = @(x_,T_) loglikeF(x_,T_,S_) + prior(x_);

%% Chain options

S_.Options.EST.N = 750000;                             % number of draws for the chain. 
S_.Options.EST.burn_in = 0.25;                         % Cut x % of the chain to compute statistic
S_.Options.EST.numbersave = 4;                         % Number of times chain is saved 
S_.Options.EST.display_num = 10; 
S_.Options.EST.df = 12;                                % Degrees of freedom for multivariate t proposal
S_.Options.EST.kap_ = 0.075;                           % Tunning proposal density
S_.Options.EST.chains = 2;                             % number of chains to run
S_.Options.EST.maxT = find(t_ == 2008);                % Max date of break; 2008Q1 (pre-GFC)
S_.Options.EST.minT = find(t_ == 2002.25);             % Min date of break; 2002Q2
S_.Options.EST.upperT = 2;                             % maximum number of updated elements in T
S_.Options.EST.upper_xx = 10; % dim_x                  % maximum number of updated elements in xx

x0_ = x_mode;
dim_x = length(x_mode);
T0 = [find(t_ == 2003.25), find(t_ == 2003.25)];

c = S_.Options.EST.kap_*inv(hess_);
c = chol(c)';
S_.Options.EST.c = c;

S_.Options.EST.post_val = post_val;                    % loaded from hessmat.mat 

S_.param_names =  {'\sigma_{g}'
               '\sigma_{tauc}'
               '\sigma_{tauw}'
               '\sigma_{tauK}'
               '\sigma_{tau}'
               '\sigma_{z}'
               '\sigma_{rstar}'
               '\sigma_{zeta}'
               '\sigma_L'
               '\sigma_I'
               '\sigma_b'
               'Upsilon'
               'h'
               'psi_b'
               'gam_{gb}'
               'gam_{cb}'
               'gam_{wb}'
               'gam_{Kb}'
               'gam_{taub}'
               '\rho_g'
               '\rho2_g'
               '\rho_c'
               '\rho2_c'
               '\rho_w'
               '\rho2_w'
               '\rho_K'
               '\rho2_K'
               '\rho_{tau}'
               '\rho2_{tau}'
               '\rho_z'
               '\rho_r_{star}'
               '\rho_{zeta}'
               '\rho_{L}'
               '\rho_{I}'
               '\rho_b'
               'Delta_z'
               'mu_sig'
               '\Delta \tau_w'
               '\Delta \tau_K'};

%% MCMC PARFOR

maxproc = S_.Options.EST.chains ;
parpool(maxproc)


parfor runs = 1:maxproc
    
    warning('Off','all') ; 
    
    out = mcmc_func(x0_,T0,S_,runs)

    par_xx_out(:,:,runs) = out.xx_out ; 
    par_T_out(:,:,runs) = out.T_out ; 
    post_out(:,runs) = out.post_out; 
    acc_xx(:,runs) = out.acc_xx; 
    acc_T(:,runs) = out.acc_T;

    par_Prop_xx(:,:,runs) = out.Prop_xx ;
    par_T_prop(:,:,runs)  = out.T_prop ;

end

delete(gcp('nocreate'))

%% Stack the (trimmed) chain outputs

maxproc = S_.Options.EST.chains ;

params_xx_out   = [] ;
params_T_out   = [] ;

params_Prop_xx = [];

for runs=1:maxproc
    
    params_xx_out = [params_xx_out ; ...
        par_xx_out(round(length(par_xx_out)*S_.Options.EST.burn_in):end,:,runs)];
    
    params_T_out = [params_T_out ; ...
        par_T_out(round(length(par_T_out)*S_.Options.EST.burn_in):end,:,runs)] ;

    params_Prop_xx = [params_Prop_xx ; ...
        par_Prop_xx(round(length(par_xx_out)*S_.Options.EST.burn_in):end,:,runs)];
end


%% Compute Statistics from the Chain

compute_mode = 1;

params_T_out_mode = mode(params_T_out);
params_T_out_mean = mean(params_T_out);
params_T_out_std  = std(params_T_out);
params_T_out_sort = sort(params_T_out); 

params_xx_mode = zeros(1,dim_x);

if compute_mode
     pp_plots; % needed to compute mode of xx with kernel smoother.
end
  
params_xx_mean = mean(params_xx_out);

% Calculate probability intervals
params_xx_sort = sort(params_xx_out);
N_sort = size(params_xx_sort,1); 
params_ci_low = params_xx_sort(round(N_sort*0.05),:);
params_ci_high = params_xx_sort(round(N_sort*0.95),:);

disp('---------STDEV--OF--SHOCKS------------------')
disp('------PARAM------MODE---&------MEAN--&----FminMode--')
disp([(1:15)',(params_xx_mode(1:15))', (params_xx_mean(1:15))', (x_mode(1:15))])
disp('--------------------------------------------')
disp('---------AR(1)--OF--SHOCKS------------------')
disp('------PARAM------MODE---&------MEAN--&----FminMode--')
disp([(16:35)',(params_xx_mode(16:35))', (params_xx_mean(16:35))', (x_mode(16:35))])
disp('--------------------------------------------')
disp('---------PARAMETERS------------')
disp('------PARAM-----MODE--&-----MEAN--&----FminMode-&--')
disp([(36:dim_x)', params_xx_mode(36:end)',params_xx_mean(36:end)', x_mode(36:end)])

%% Save workspace
eval(['save mhall_RBC.mat;']);