(* ::Package:: *)

%% Code to run IV SVAR for quarterly VAR
% Hambur and Haque (2023)
clear
close all
addpath(genpath('Code/IV_SVAR'))                       % set path (needs to be modified)
outputs = "U/Outputs/VAR outputs/VAR_Quarterly";                       % set path and file name for output (needs to be modified)
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

q_data _start = "1994-03-31";                      % including missing year
q_data _start _instr = "2002-09-30";                % exclude missing yearr
q_data _end = "2019-12-31";                        % End at end of data

modelSpec.nLags                 =4;               % VAR lags
modelSpec.nHorizons             =16;              % max horizon for IRFs
modelSpec.cLevel                =[68 90];         % size of error bands
modelSpec.bootSize              =1000;            % size of bootstrap sample for error bands


output_vars _q = {'unemployment', 'aus_gdp _deflated _cpi', 'aus_cpi', 'rtwi', 'pc21', 'pc22', 'm_ 120', 'aus_hp','retail_sales','dwelling_approvals','cash_rate', 'mortgage_spread'};
output_vars _q _fix = {'unemployment', 'aus_gdp _deflated _cpi', 'aus_cpi', 'rtwi'};

shocks = {'exp_pc1'};

% VAR data
opts = detectImportOptions('/Input data/VAR data/var_data_m.csv');                      % path for the VAR data (needs to be modified)
opts = setvartype(opts,output_vars _q _fix,'double');

var_data = readtable('/Input data/VAR data/var_data_m.csv',opts);                       % path for the VAR data (needs to be modified)
var_data= var_data(:,2:size(var_data,2));
var_table = table2timetable(var_data);

opts = detectImportOptions('/Input data/VAR data/Full_hist_month.csv');                 % path for yield curve data (needs to be modified)
opts = setvaropts(opts,{'Date'},'InputFormat','dd/MM/uuuu');
pc_data = readtable('/Input data/VAR data/Full_hist_month.csv',opts);                   % path for yield curve data (needs to be modified)
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);


% Shock data
opts = detectImportOptions('/Input data/VAR data/pc_yield_baseline.csv');              % path for the shocks (needs to be modified)
opts = setvaropts(opts,{'dates'},'InputFormat','dd/MM/uuuu');
shock_data = readtable('/Input data/VAR data/pc_yield_baseline.csv',opts);             % path for the shocks (needs to be modified)
shock_table = table2timetable(shock_data);

shock_table.dates= eomdate(shock_table.dates);

data_table = synchronize(var_table, pc_table, shock_table);

% Bring together
data_table _q = data_table(:, ismember(data_table.Properties.VariableNames, [output_vars _q shocks]));
data_table _q = convert2quarterly(data_table _q, 'Aggregation', ["mean" "mean" "mean" "mean" "mean" "mean" "lastvalue" "lastvalue" "lastvalue" "mean" "mean" "mean" "sum" ]);

data_table _q = data_table _q(timerange(q_data _start,q_data _end),:);
data_table _q{:,:}(isnan(data_table _q{:,:}))=0;

%-2-PSVAR GK--------------------------------------------------------------%
var_table _q=data_table _q(:,ismember(data_table _q.Properties.VariableNames, output_vars _q ));
shock_table _q=data_table _q(timerange(month_data _start _instr,month_data _end),ismember(data_table _q.Properties.VariableNames, shocks ));
dataStructure.data        = data_table _q.Variables;
dataStructure.dates       = data_table _q.date;
dataStructure.varname     = data_table _q.Properties.VariableNames;
dataStructure.varLongName = data_table _q.Properties.VariableNames;

IV.data        = shock_table _q.Variables;
IV.dates       = shock_table _q.date;
IV.labels     = shock_table _q.Properties.VariableNames;
IV.description = shock_table _q.Properties.VariableNames;

modelSpec.identification        ='PSVAR';

psvarSet                        ={'unemployment'; 'rtwi'; 'aus_cpi'; 'aus_hp'; 'aus_gdp _deflated _cpi'; 'dwelling_approvals'; 'mortgage_spread'; 'pc21';'pc22'};
[~,dataSelection]               =ismember(psvarSet,dataStructure.varname);


% Something to store F Stata
fstats = zeros(1,2);
robfstats = zeros(1,2);
Ls =zeros(1,2);


for i = 13
    modelSpec.selectedInstrument=dataStructure.varname(i);
    for q = 11:12

        modelSpec.dataSelection         =dataSelection;

        shockVar                        =dataStructure.varname(q);

        % declare shock variable and shock size
        modelSpec.shockVar              =ismember(psvarSet,shockVar);
        modelSpec.shockSize             =1*double(ismember(psvarSet,shockVar));

        % compute responses
        irfSVARIV_ 1                     = computeIRFs_boot _cr _ 10(dataStructure,modelSpec, IV,W2);

        % SAVE OUTPUT
        runid =outputs +'IRFs_q _' +string(modelSpec.selectedInstrument) +'_' + q ;
        save(runid,...
            'modelSpec', 'irfSVARIV_ 1', 'dataStructure', 'W2')



        % SAVE charting OUTPUT
        runid =outputs + 'IRFsg_' +string(modelSpec.selectedInstrument) +'_' + q + '.png' ;

        irfSVARIV_ 1.irfs = [irfSVARIV_ 1.irfs; ([irfSVARIV_ 1.irfs(end-1,:); irfSVARIV_ 1.irfs(end,:)]'*W2(:,end))'];
        irfSVARIV_ 1.irfs = [irfSVARIV_ 1.irfs; ([irfSVARIV_ 1.irfs(end-2,:); irfSVARIV_ 1.irfs(end-1,:)]'*W2(:,2))'];

        irfs_temp = zeros(size(irfSVARIV_ 1.irfs_u,1)+2,size(irfSVARIV_ 1.irfs_u,2), size(irfSVARIV_ 1.irfs_u,3));
        irfs_temp(1:size(irfSVARIV_ 1.irfs_u,1), :, :) = irfSVARIV_ 1.irfs_u;
        irfs_temp(size(irfSVARIV_ 1.irfs_u,1)+1,: , 1) = ([irfSVARIV_ 1.irfs_u(end-1,:,1); irfSVARIV_ 1.irfs_u(end,:,1)]'*W2(:,end))';
        irfs_temp(size(irfSVARIV_ 1.irfs_u,1)+2, :, 1) = ([irfSVARIV_ 1.irfs_u(end-1,:,1); irfSVARIV_ 1.irfs_u(end,:,1)]'*W2(:,2))';
        irfs_temp(size(irfSVARIV_ 1.irfs_u,1)+1,:, 2) = ([irfSVARIV_ 1.irfs_u(end-1,:,2); irfSVARIV_ 1.irfs_u(end,:,2)]'*W2(:,end))';
        irfs_temp(size(irfSVARIV_ 1.irfs_u,1)+2, :, 2) = ([irfSVARIV_ 1.irfs_u(end-1,:,2); irfSVARIV_ 1.irfs_u(end,:,2)]'*W2(:,2))';
        irfSVARIV_ 1.irfs_u=irfs_temp;


        irfs_temp = zeros(size(irfSVARIV_ 1.irfs_l,1)+2,size(irfSVARIV_ 1.irfs_l,2), size(irfSVARIV_ 1.irfs_l,3));
        irfs_temp(1:size(irfSVARIV_ 1.irfs_l,1), :, :) = irfSVARIV_ 1.irfs_l;
        irfs_temp(size(irfSVARIV_ 1.irfs_l,1)+1,: , 1) = ([irfSVARIV_ 1.irfs_l(end-1,:,1); irfSVARIV_ 1.irfs_l(end,:,1)]'*W2(:,end))';
        irfs_temp(size(irfSVARIV_ 1.irfs_l,1)+2, :, 1) = ([irfSVARIV_ 1.irfs_l(end-1,:,1); irfSVARIV_ 1.irfs_l(end,:,1)]'*W2(:,2))';
        irfs_temp(size(irfSVARIV_ 1.irfs_l,1)+1,:, 2) = ([irfSVARIV_ 1.irfs_l(end-1,:,2); irfSVARIV_ 1.irfs_l(end,:,2)]'*W2(:,end))';
        irfs_temp(size(irfSVARIV_ 1.irfs_l,1)+2, :, 2) = ([irfSVARIV_ 1.irfs_l(end-1,:,2); irfSVARIV_ 1.irfs_l(end,:,2)]'*W2(:,2))';
        irfSVARIV_ 1.irfs_l=irfs_temp;



        %-plot irfs---------------------------------------------------------------%
        Ccolor  =[ 1 .6 .2];
        GKcolor =[.0 .0 .4];
        Scolor  =[ 1 .0 .0];    bandFillColor1 =[.85 .85 .85];  bandFillColor2 =[.7 .7 .7];

        % plot labels
        varname                 =dataStructure.varLongName;


        figure;
        plotColumns =3; pln=1;
        nHorizon    =modelSpec.nHorizons;
        n           =length(psvarSet)+2;                                     % added for 10 year year rate and cash rate
        modelSpec.dataSelection=[modelSpec.dataSelection; 10; 3];            % added for 10 year year rate and cash rate

        for j=1:n

            subplot(ceil(n/plotColumns),plotColumns,pln)

            hold on

            % bands 90
            fill([0:1:nHorizon, fliplr(0:1:nHorizon)],...
                [squeeze(irfSVARIV_ 1.irfs_u(j,:,2)) fliplr(squeeze(irfSVARIV_ 1.irfs_l(j,:,2)))],...
                bandFillColor1,'EdgeColor','none');

            % bands 68
            fill([0:1:nHorizon, fliplr(0:1:nHorizon)],...
                [squeeze(irfSVARIV_ 1.irfs_u(j,:,1)) fliplr(squeeze(irfSVARIV_ 1.irfs_l(j,:,1)))],...
                bandFillColor2,'EdgeColor','none');


            % zero line
            plot(0:nHorizon,zeros(size(0:nHorizon)),'k')
            hold on

            % irfs
            plot(0:nHorizon,irfSVARIV_ 1.irfs(j,:),'-.','LineWidth',1.5,'color',GKcolor)

            hold off; axis tight

            grid on; xlim([0 nHorizon]);
            set(gca,'XTick',0:6:nHorizon,'XTickLabel',cellstr(num2str((0:6:nHorizon)')),'layer','top')
            title(varname{modelSpec.dataSelection(j)},'FontSize',9,'FontWeight','normal')

            if j==1

                ylabel('% points','FontSize',9)
            end


            pln=pln+1;
            dataid = outputs +'IRFs_q _' +string(modelSpec.selectedInstrument) +'_' + q + '.xls' ;

            output_irf = [irfSVARIV_ 1.irfs(j,:); squeeze(irfSVARIV_ 1.irfs_l(j,:,2)); squeeze(irfSVARIV_ 1.irfs_u(j,:,2)); squeeze(irfSVARIV_ 1.irfs_l(j,:,1)); squeeze(irfSVARIV_ 1.irfs_u(j,:,1))];
            output_irf = array2table(output_irf');
            row_labels = {'irfs', 'lower_ 90', 'upper_ 90', 'lower_ 68','upper_ 68'};
            output_irf.Properties.VariableNames=row_labels;
            if j<n-1
                writetable(output_irf,dataid,'Sheet',char(varname{modelSpec.dataSelection(j)}));
            end
            if j==n
                writetable(output_irf,dataid,'Sheet','Short rate');
            end
            if j==n-1
                writetable(output_irf,dataid,'Sheet','10 year yield');
            end

        end

        saveas(gcf, runid)



        fstats(i-12, q-10) = irfSVARIV_ 1.fstat;
        robfstats(i-12, q-10) = irfSVARIV_ 1.robfstat;
        Ls(i-12, q-10) = irfSVARIV_ 1.L;
    end
end
fstats_tab= array2table(fstats);
robfstats_tab= array2table(robfstats);
runid =outputs + 'IRFs_fstats.xls' ;

writetable(fstats_tab,runid,'Sheet','F_stat');
writetable(robfstats_tab,runid,'Sheet','robust_F _stat')
