(* ::Package:: *)

%%
% This code takes the estiamted decopose HF MP shocks, and regresses them
% on a number of varialbes to orthoganlise them
% For Hambur and Haque (2023)

%% Set inital setting
clear
close all

addpath('Code\Decomposing shocks\additional');
% 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

q_data _start = "2001-06-30"; % including missing year
% q_data _start = "2002-12-31"; % exclude missing yearr
% q_data _end = "2019-12-31"; % End at end of shockc
q_data _end = "2021-03-31"; % End at end of data

gfc_start _prev = "2008-08-31";
gfc_end _prev ="2009-04-30";
% full list vars
output_vars _m = {'twi', 'money_market _spread', 'us_cp _spread', 'pc22', 'commodity_prices', 'asx200', 'baa', 'vix','employment', 'unemployment'};
output_vars _m _l = {'asx200' }; % those in levels that need taking logs. employment, commodity, twi and vix already logged.

reg_vars= {'employment', 'unemployment', 'twi', 'money_market _spread', 'us_cp _spread', 'pc22', 'commodity_prices', 'asx200', 'baa', 'vix'};
reg_vars _ 12m= {'employment', 'unemployment'}; % annual vars
reg_vars _ 3m= {'twi', 'money_market _spread', 'us_cp _spread', 'pc22', 'commodity_prices', 'asx200', 'baa', 'vix'}; % 3m change vars


shocks = {'exp_pc1', 'exp_pc2', 'tp_pc1', 'exp_dir _pc1', 'exp_dir _pc2', 'tp_dir _pc1'};


opts = detectImportOptions('\Input data\VAR data\orthog_series.csv');
opts = setvartype(opts,output_vars _q,'double'); 
% VAR data
var_data = readtable('G\Input data\VAR data\orthog_series.csv',opts);
var_data= var_data(:,2:size(var_data,2));
var_table = table2timetable(var_data);

% Table F17 for yields for PC
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 VAR and PC data together,ake 3m and 12m changes
var_table = synchronize(var_table, pc_table, 'intersection');
var_table{:,output_vars _m _l} =log(var_table{:,output_vars _m _l});
var_table{:,reg_vars _ 12m} = lag(var_table{:,reg_vars _ 12m},2)-lag(var_table{:,reg_vars _ 12m},14); % EXTRA LAG TO REFLECT TIMING OF lfs
var_table{:,reg_vars _ 3m} = lag(var_table{:,reg_vars _ 3m},1)-lag(var_table{:,reg_vars _ 3m},4);

% Shock data
opts = detectImportOptions('Input data\VAR data\pc_yield_old _baseline.csv'); % change this if want to bring in bootstrapped or all events
opts = setvaropts(opts,{'dates'},'InputFormat','dd/MM/uuuu'); 
shock_data = readtable('Input data\VAR data\pc_yield_old _baseline.csv',opts); % change this if want to bring in bootstrapped or all events
shock_table = table2timetable(shock_data);

shock_table.dates= eomdate(shock_table.dates);
% bring together
data_table = synchronize(var_table, shock_table, 'intersection');
data_table _m = data_table(:, ismember(data_table.Properties.VariableNames, [output_vars _m shocks]));
data_table _m = data_table _m(timerange(month_data _start,month_data _end),:);
data_table _m{:,:}(isnan(data_table _m{:,:}))=0;
data_table _m{:,:}(isnan(data_table _m{:,:}))=0;

% table as for regs
reg_table = data_table _m(:, ismember(data_table _m.Properties.VariableNames, reg_vars));
% Do regs
reg.y = table2array(data_table _m(timerange(month_data _start _instr,month_data _end),ismember(data_table _m.Properties.VariableNames, shocks)));
reg.x = table2array(reg_table(timerange(month_data _start _instr,month_data _end),:));
reg.x=[ones(size(reg.x,1),1), reg.x].';
reg.a = reg.y'*reg.x.'*inv(reg.x*reg.x.');
reg.fitted = (reg.a*reg.x).';

reg.residuals= reg.y-(reg.a*reg.x).'; 
reg.residuals2 = reg.residuals.^2;
reg.SSE=sum(reg.residuals2,1);
reg.y2 = (reg.y-mean(reg.y,1)).^2;
reg.SSY=sum(reg.y2,1);
reg.r2=1-reg.SSE./reg.SSY;

% now do a version removing GFC periosd
reg_table _pre _gfc = reg_table(timerange(month_data _start _instr,gfc_start _prev),:);
reg_table _post _gfc =reg_table(timerange(gfc_end _prev,month_data _end),:);
reg_table _no _gfc= [reg_table _pre _gfc; reg_table _post _gfc];

shock_pre _gfc = data_table _m(timerange(month_data _start _instr,gfc_start _prev),ismember(data_table _m.Properties.VariableNames, shocks));
shock_post _gfc =data_table _m(timerange(gfc_end _prev,month_data _end),ismember(data_table _m.Properties.VariableNames, shocks));
shock_no _gfc= [shock_pre _gfc; shock_post _gfc];
% stats for share of variation explained
regnogfc.y = table2array(shock_no _gfc);
regnogfc.x = table2array(reg_table _no _gfc);
regnogfc.x=[ones(size(regnogfc.x,1),1), regnogfc.x].';
regnogfc.a = regnogfc.y'*regnogfc.x.'*inv(regnogfc.x*regnogfc.x.');
regnogfc.fitted = (regnogfc.a*regnogfc.x).';

regnogfc.residuals= regnogfc.y-(regnogfc.a*regnogfc.x).'; 
regnogfc.residuals2 = regnogfc.residuals.^2;
regnogfc.SSE=sum(regnogfc.residuals2,1);
regnogfc.y2 = (regnogfc.y-mean(regnogfc.y,1)).^2;
regnogfc.SSY=sum(regnogfc.y2,1);
regnogfc.r2=1-regnogfc.SSE./regnogfc.SSY;


% export
export_table = array2table(reg.residuals);
 export_table.Properties.VariableNames=shocks;
 export_table = [data_table _m(timerange(month_data _start _instr,month_data _end),ismember(data_table _m.Properties.VariableNames, date)), export_table];
writetimetable(export_table, '\Input data\VAR data\pc_yield_orthog.csv')


export_table _no _gfc = array2table(regnogfc.residuals);
 export_table _no _gfc.Properties.VariableNames=shocks;
 dates_pre _gfc=reg_table(timerange(month_data _start _instr,gfc_start _prev),ismember(data_table _m.Properties.VariableNames, date));
 dates_post _gfc=reg_table(timerange(gfc_end _prev,month_data _end),ismember(data_table _m.Properties.VariableNames, date));
dates_no _gfc = [dates_pre _gfc; dates_post _gfc];
 export_table _no _gfc = [dates_no _gfc, export_table _no _gfc];
writetimetable(export_table _no _gfc, '\Input data\VAR data\pc_yield_orthog _no _gfc.csv')
