(* ::Package:: *)

%% Code to run Bayesian Local Projection on identified shocks
%% for Hambur and Haque (2023) 
clear
close all
addpath(genpath('\Code\LP\'))
addpath('Code\Nom ATSM') % Contains some of the sub-functions
% month_data _start = "2001-04-30"; % including missing year
month_data _start _instr = "2002-09-30"; % exclude missing yearr
month_data _start = "1994-01-31"; % exclude missing yearr

month_data _end = "2019-12-31"; % End at end of shockc
% month_data _end = "2021-6-30"; % End at end of data

init_reg = {'pc21', 'pc22', 'retail_sales', 'unemployment', 'twi', 'mortgage_spread' };
adit_reg = {'dwelling_approvals', 'aus_hp' };
%% Get data loaded
output_vars _m = {'employment', 'unemployment', 'unemployment_gap', 'mortgage_spread', 'mortgage_rate', 'cash_rate', 'aus_hp', 'dwelling_approvals', 'retail_sales', 'bbsw_ 3m',  'ebp', 'gz_spread',  'vix', 'twi' };
% output_vars _m _l = { 'asx200' };

shocks = {'exp_pc1', 'exp_pc2', 'tp_pc1'};

opts = detectImportOptions('\Input Data\Var data\var_data_m.csv');
opts = setvartype(opts,output_vars _m,'double'); 
var_data = readtable('\Input Data\Var data\var_data_m.csv'',opts);
var_data= var_data(:,2:size(var_data,2));
var_table = table2timetable(var_data);
% Bring in yield data for PCs
opts = detectImportOptions('\Input Data\Var data\Full_hist_month.csv');
opts = setvaropts(opts,{'Date'},'InputFormat','dd/MM/uuuu'); 
pc_data = readtable('\Input Data\Var data\Full_hist_month.csv',opts);
pc_input = standardise(table2array(pc_data(:,2:end)));
[pc2,W2]=extract(pc_input,2);
pc2(:,2)=-pc2(:,2);
W2(2,:)=-W2(2,:);
pc_table = [pc_data(:,1), pc_data(:,42),array2table(pc2) ];
pc_table= table2timetable(pc_table);
% Bring in shocks
opts = detectImportOptions('\Input Data\Var data\pc_yield_old _adj.csv');
opts = setvaropts(opts,{'dates'},'InputFormat','dd/MM/uuuu'); 
shock_data = readtable('\Input Data\Var data\pc_yield_old _adj.csv',opts);
shock_table = table2timetable(shock_data);

shock_table.dates= eomdate(shock_table.dates);

data_table = synchronize(var_table, pc_table, shock_table);

data_table _m = data_table(:, ismember(data_table.Properties.VariableNames, [output_vars _m init_reg shocks]));
data_table _m = data_table _m(timerange(month_data _start,month_data _end),:);
% data_table _m{:,output_vars _m _l} =log(data_table _m{:,output_vars _m _l});
data_table _m{:,:}(isnan(data_table _m{:,:}))=0;
data_table _m{:,:}(isnan(data_table _m{:,:}))=0;

names=data_table _m.Properties.VariableNames;
data0=data_table _m.Variables;

[~,id0]               =ismember(init_reg,names);
[~,id1]               =ismember(adit_reg,names);
[~,instr_num]               =ismember(shocks,names);
NN=length(id1);
T0=size(data_table(timerange(month_data _start,month_data _start _instr),:),1); % length of pre-sample
SS=length(instr_num);
% Run Bayesian LP over each instrument looking at
for i = 1:SS
    id_m=instr_num(i); 
    shock_name=names(id_m);
    for q = 1:2
        shockVAR=q;

        [irf,irf1,irf2,irfyy,irf_all,irfBLP]=run_LP(data0,id0,T0,id_m,shockVAR,names,W2,12,1,2);
        relstat=irfBLP.reliability;
        FF=irfBLP.Fstatistic;



        parfor j=1:NN
            id=[id0 id1(j)];


            [irf,irf1,irf2,irfyy,irf_all]=run_LP(data0,id,T0,id_m,shockVAR,names,W2,12,1,2);

            save_irfyy2(:,:,j)=irfyy (:,:,end)';
            save_irfyy1(:,:,j)=irfyy (:,:,1)';
            save_irf(:,:,j)=irf_all;

        end
    names1=names(id1);
% save results
results = 'Outputs\LP outputs \r esults_LP _' + string(shock_name) +'_' + q;
save(results,...
    'relstat','FF','save_irf','save_irfyy1','save_irfyy2','names1','irf_all','irfyy','irf','irf1','irf2','names', 'W2')
% save('results_ 2_sl',...
%    'relstat','FF','irf_all','irfyy','irf','irf1','irf2','names')
init_reg = {'pc21', 'pc22', 'retail_sales', 'unemployment', 'mortgage_spread', 'twi' };
adit_reg = {'aus_hp', 'dwelling_approvals'};
% Graphing
namesx{3}='Factor 1';
namesx{4}='Factor 2';
namesx{5}='Retail_sales';
namesx{6}='Unemployment';
namesx{7}='Mortgage_spread';
namesx{8}='TWI';
namesx{8}='House prices';
namesx{9}='Dwelling approvals';


namesx{1}='10 year yield';
namesx{2}='Short rate';
Horz=36;
HH=0:Horz;

irf_ 10 = [irf(:,1) irf(:,2)]*W2(:,end);
irf_ 10_u = [irf1(:,1) irf1(:,2)]*W2(:,end);
irf_ 10_l = [irf2(:,1) irf2(:,2)]*W2(:,end);

irf_sr = [irf(:,1) irf(:,2)]*W2(:,2);
irf_sr _u = [irf1(:,1) irf1(:,2)]*W2(:,2);
irf_sr _l = [irf2(:,1) irf2(:,2)]*W2(:,2);

figure(2)
subplot(4,3,1)
plotx3(HH', [irf_ 10 irf_ 10_u irf_ 10_l]);
title(namesx{1},'Interpreter','Latex','FontSize',12);



subplot(4,3,2)
plotx3(HH',[irf_sr irf_sr _u irf_sr _l]);
title(namesx{2},'Interpreter','Latex','FontSize',12);

for j=1:length(id0)
    subplot(4,3,j+2)
plotx3(HH',[irf(:,j) irf1(:,j) irf2(:,j)]);

title(namesx{j+2},'Interpreter','Latex','FontSize',12);
end
jj=1;
for j=length(id0)+1:length(id0)+length(id1)
    subplot(4,3,j+2)
plotx3(HH',save_irf (:,:,jj)');

title(names1{jj},'Interpreter','Latex','FontSize',12);
jj=jj+1;
end

chart_name = 'Outputs\LP outputs\chart_orthog_'+results;
saveas(figure(2), chart_name);
close Figure 2

% output
mean = [irf_ 10, irf_sr, irf, save_irf (1,:,1)', save_irf (1,:,2)' ];
upper = [irf_ 10_u, irf_sr _u, irf1, save_irf (3,:,1)', save_irf (3,:,2)'];
lower = [irf_ 10_l, irf_sr _l, irf2, save_irf (2,:,1)', save_irf (2,:,2)' ];
    end 
end

