% Code from Mertens and Ravn (2018) modified by Dias and Duarte
% All mistakes are our own. Please email joao.duarte@novasbe.pt for any
% mistakes or bugs found.
% last update: 8/10/2018


clear all;

% =========================================================================
% DESCRIPTION 
% This script:
%   1- loads in quarterly data, 
%   2- Identifies monetary shocks using High-frequency data.
%   3- Plots the IRFs that are computed one at a time

addpath('auxiliary files')


%% Do IRS one by one

varlist = {'HP', 'HR', 'RVR', 'HOMEOWN'};


for i=1:4
% Create DATASET
%%%%%%%%%%%%%%%%
mdata = xlsread('datasets/data/quarterly_hfi_var_data.xlsx');
%mdata(:,15) = mdata(:,15)./mdata(:,3);
mdata(:,3:4) = log(mdata(:,3:4))*100;
mdata(:,9:10) = log(mdata(:,9:10))*100;
mdata(:,9:10) = log(mdata(:,9:10))*100;
mdata(:,14:end) = log(mdata(:,14:end))*100;
mdata(isnan(mdata))=0;
mdata = mdata(7:end,:); % Activate to start the sample in 1981 
% 1 year	2 quarter	3 cpi 4 ip 5 ff	6 gs1	7 gs2	8 ebp  

pdata = xlsread('datasets/data/quarterly_factors.xlsx');
ppdata = zeros(length(mdata),size(pdata,2));
ppdata(end-length(pdata)+1:end,:) = pdata;
ppdata(isnan(ppdata))=0;
% 1 year	2 quarter	3 mp1_tc 4 ff4_tc	5 ed2_tc	6 ed3_tc	7 ed4_tc

DATASET.TSERIES = [mdata ppdata(:,3:end)]; 
DATASET.LABEL   = {'YEAR','MONTH','CPI','IP','FF','GS1','GS2','EBP','HP','HR','HOMEOWN','HVR','RVR','REALHP','REALHR','OWNEQ','MP1','FF4','ED2','ED3','ED4'};
DATASET.UNIT    = [0         0      1     1    2    2     2          2       2    2      2     2 2 2 2 2 2 2 2 2  ]; 
DATASET.FIGLABELS= [{'YEAR','MONTH','Price Level','Industrial Production','Federal Funds Rate','One-Year Rate','Two-Year Rate','Excess Bond Premium','HP','HR','HOMEOWN','HVR','RVR','REALHP','REALHR','OWNEQ','MP1','FF4','ED2','ED3','ED4'};];  
DATASET.MAP = containers.Map(DATASET.LABEL,1:size(DATASET.TSERIES,2));

% VAR specification
%%%%%%%%%%%%%%%%%%%%
VAR.p      = 4;                                % Number of Lags
VAR.irhor  = 20;                                % Impulse Response Horizon
VAR.select_vars      = {'GS1','IP','EBP','CPI',varlist{i}};
VAR.vars             = DATASET.TSERIES(:,cell2mat(values(DATASET.MAP,VAR.select_vars)));
VAR.MAP              = containers.Map([VAR.select_vars],[1:size(VAR.vars,2)]);
VAR.proxies          = DATASET.TSERIES(:,cell2mat(values(DATASET.MAP,{'FF4'})));
VAR.DET              = ones(length(VAR.vars),1); % Deterministic Terms

VAR        = doProxySVAR(VAR);

% Inference:
%%%%%%%%%%%%
% method  1: Mertens and Ravn (2013) wild bootstrap 
%         2: Montiel-Olea Stock Watson (2016) parametric bootstrap
%         3: Delta Method 
%         4: Montiel-Olea Stock Watson (2016) asy weak IV
%         5: Jentsch and Lunsford Moving Block Bootstrap
%         6: Jentsch and Lunsford Moving Block Bootstrap (adjusted to allow non zero-mean proxies)

nboot     = 5000;         % Number of Bootstrap Samples (Paper does 5000)
clevel    = 68;          % Bootstrap Percentile Shown
BlockSize = floor(5.03*length(DATASET.TSERIES).^0.25); % size of blocks in the MBB bootstrap
seed      = 2;           % seed for random number generator
rng(seed);               % iniate the random number generator

%VARci_wildbs   = doProxySVARci(VAR,clevel,1,nboot); 
%VARci_mswbs    = doProxySVARci(VAR,clevel,2,nboot);
%VARci_delta    = doProxySVARci(VAR,clevel,3);
% VARci_msw_wiv  = doProxySVARci(VAR,clevel,4);
VARci_mbb      = doProxySVARci(VAR,clevel,6,nboot,BlockSize);

shocksize = -0.25;
VAR.irs   = shocksize*VAR.irs;

legendflag = 1;
shock = 1;
%VARci.irsL  = shocksize*VARci_delta.irsL(:,:,shock);
%VARci.irsH  = shocksize*VARci_delta.irsH(:,:,shock);
%VARci.irsL2 = shocksize*VARci_wildbs.irsL(:,:,shock);
%VARci.irsH2 = shocksize*VARci_wildbs.irsH(:,:,shock);
%VARci.irsL3 = shocksize*VARci_mswbs.irsL(:,:,shock);
%VARci.irsH3 = shocksize*VARci_mswbs.irsH(:,:,shock);
VARci.irs(:,i) = shocksize*VARci_mbb.irs(:,5,shock);
VARci.irsL4(:,i) = shocksize*VARci_mbb.irsL(:,5,shock);
VARci.irsH4(:,i) = shocksize*VARci_mbb.irsH(:,5,shock);
%legendlabels = {'Delta Method','Wild Bootstrap','Par. Bootstrap','Block Bootstrap'};
name = 'Figure2';
LEGENDPANEL = 1;
xlab_text='horizon (quarters)';

end

plotdisplay = {'HP', 'HR', 'RVR', 'HOMEOWN'};
DATASET.FIGLABELS= {'Housing Prices', 'Housing Rents', 'Housing Stock for Renting Vacancy Rate', 'Homeownership Rate'};
do_IRS_Figure2DD


