

%Paper "Heterogeneity and the Macroeconomic Effects of Changes in
%Loan-to-Value Limits" by Emmanuel De Veirman and Jasper de Jong
%Journal of Money, Credit and Banking

%This program estimates the baseline VECM and computes impulse responses for 
%various scenarios including the Dutch scenario (i.e. the imposition and gradual
%reduction of the LTV limit to 100%)
 
%In its current settings, it produces the impulse responses in
%Figures 4 and 5 of the paper

%Key settings are:
%k : sets the VECM specification
%shockspec: sets the specification for the shock (series)
%bootstrap: determines whether bootstrapped confidence intervals are
%computed
%N : number of samples in bootstrap
%robustness: whether loop across specifications to examine robustness

addpath(genpath('JMCB_replic_deveirman-jong/Matlabtools'))

tic
clear; close all; clc;
rng('default');

cd('JMCB_replic_deveirman-jong\VECM estimation');

A = importdata('alldatarobu3.xlsx'); %Nine series June 2022.
%These are the processed data from Manu-April2022-VECM.prg.
A = A.data;
year = A(:,1); 
loghp = A(:,2); %log house prices 
loggdp = A(:,3); %log GDP
ltv_comb_nz = A(:,4); %our survey-based time series for the LTV
rh_ea1 = A(:,5); %ex ante real mortgage rate
logcons = A(:,6); %log consumption  
logires = A(:,7); %log residential investment
logibus = A(:,8); %log business investment

%Import the shocks for the Dutch scenario from Excel

B=importdata('matlab_dutchscenario2.xlsx'); %June 2022. %Update Oct 2022: add d1easp
colheaders=B.textdata;
B=B.data;

yearned=B(:,1); %years in the Dutch scenario
d1ep=B(:,2); %Baseline distributions: ex post with exceptions (as in original JMCB submission)
d1ea=B(:,3); %Bas. dist.: ex ante with exceptions 0.26. Baseline.
d1ea00=B(:,4); %Bas dist: ex ante without exceptions
d1ea20=B(:,5); %Bas dist: ex ante with exc 0.20
d1ea30=B(:,6); %Bas dist: ex ante with exc 0.30
d1ea50=B(:,7); %Bas dist: ex ante with exc 0.50
d2ea=B(:,8); %2010M4M12 as pre-cap dist: ex ante with exceptions 0.26
d3ea=B(:,9); %2010 as pre-cap dist: ex ante with exceptions 0.26
d1easp=B(:,10); %Bas. dist. ex ante with exceptions 0.26 and sprinkling.

%September 2022: compare 8-year (old way) to 7-year (new way) impulses
B2 = importdata('matlab_dutchscenario_long.xlsx');
colheaders = B2.textdata;
B2 = B2.data;
yearned8 = B2(:,1); %In all of these, the first row is 2011 instead of 2012
d1ea8 = B2(:,2); %This is different than d1ea: 2011 is the first row
d1ea8our = B2(:,3);
d1ea8ref = B2(:,4);

%October 2023: 9-year impulses.
B3 = importdata('matlab_dutchscenario_verylong.xlsx');
colheaders = B3.textdata;
B3 = B3.data;
yearned9 = B3(:,1); %Years 2011-2019
d1ea9our = B3(:,2); %Back to baseline where 5/12th to 2011. Exceptions 0.26  . "9" since it has 9 years, but 9th year is zero.
d1ea9del = B3(:,3); %Baseline + origination delays
d1easp9our = B3(:,4); %Baseline + sprinkling
d1easp9del = B3(:,5); %Baseline + sprinkling and origination delays
d1easp9bigdel = B3(:,6); %Baseline + sprinkling and long origination delays
d1ea9bigdel = B3(:,7); %Baseline + long origination delays
d2ea9our = B3(:,8); %Same as d1ea9our but with distribution(s) of last nine months
d3ea9our = B3(:,9); %Same as d1ea9our with with distribution for 2010

%Import the shocks for the Dutch scenario assuming various distributions
%differing (primarily) in terms of the standard deviation
%VERSION OCTOBER 2023
C = importdata('matlab_sigma.xlsx'); %up to version July 2022: .csv
colheaders = C.textdata; %C.colheaders; %.colheaders is for .csv file
C = C.data;
yearsigma=C(:,1); untr80=C(:,2); untr85=C(:,3); untr90=C(:,4); untr95=C(:,5);
untr100=C(:,6); tr80=C(:,7); tr85=C(:,8); tr90=C(:,9); tr95=C(:,10); tr100=C(:,11);
%yearsigma: years in Dutch scenario (sigma). untr: untruncated Pearson
%family distributions with stdevs ranging from 0.80 to 1 times the
%empirical stdev. tr: truncated Pearson family distributions.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Estimate a VECM(1) on annual data for 1979-2018 (+2 for diff&lags)
%with one cointegrating relation, using LeSage toolbox.
%EXCEPT for specifications with rh_ea1: 1981-2018

robustness = 0; %1 if do all robustness checks; 0 if only do baseline
%for all shockspec.  %When set this to 1, also set k to the 10-entry vector
%indicated in line 139.


%Create data matrices.
%Create array, which each cell containing a data matrix
yd = [loghp loggdp rh_ea1 ltv_comb_nz]; %Baseline
yl = [loghp logcons logires logibus rh_ea1 ltv_comb_nz];

%Alternative Cholesky orderings for baseline
%Version a)
% yd2 = [loggdp ltv_comb_nz loghp rh_ea1 ];
% yd3 = [ltv_comb_nz loggdp loghp rh_ea1];
% yd4 = [ltv_comb_nz loggdp rh_ea1 loghp];
%Version b) 
yd2 = [loggdp loghp rh_ea1 ltv_comb_nz];
yd3 = [loggdp rh_ea1 loghp ltv_comb_nz];
yd4 = [loggdp rh_ea1 ltv_comb_nz loghp];

yar = cell(17,1);
yar{4} = yd;
yar{8} = yd; yar{9} = yd; yar{10} = yd; %8-9-10 are 4 with 2-3-4 lags
yar{11} = yd; yar{12} = yd; %11-12 are 4 with 2-3 cointegration relationships
yar{15} = yd2; yar{16} = yd3; yar{17} = yd4; %Alternative Cholesky orderings for baseline
yar{23} = yl; %5, but without transactions

chind = nan(4,4); %Cholesky index. Leave first row empty, so row 2-3-4 are Cholesky specification 2-3-4.
%Version a
%chind(2,:) = [2 4 1 3];
%chind(3,:) = [4 2 1 3];
%chind(4,:) = [4 2 3 1];
%Version b
chind(2,:) = [2 1 3 4];
chind(3,:) = [2 3 1 4];
chind(4,:) = [2 3 4 1];

 
%for shockspec = 5 (sigma): use k = 4.
for k = 4:4  % k = [4 8 9 10 11 12 15 16 17 23] 
%To check for robustness: set robustness =1 and set k to the 10-entry vector that is currently commented out
    
    %Define data matrix
    y = yar{k};
      
    y = y(3:end,:); %Sample 1979-2018 (because of availability rh_ea1)
      
 
    if k == 4 | k == 15 | k == 16 | k == 17 | k == 23 %Specifications with 1 lag and 1 cointegration relation
        result = ecm_manu1(y,1,1); %Estimate VECM with 1 lag and 1 cointegration relation
    elseif k == 8 %Baseline but with 2 lags
        result = ecm_manu1(y,2,1);
    elseif k == 9 %Baseline but with 3 lags
        result = ecm_manu1(y,3,1);
    elseif k == 10 %Baseline butwith 4 lags
        result = ecm_manu1(y,4,1);
    elseif k == 11 %Baseline but with 2 cointegration relations
        result = ecm_manu1(y,1,2);
    elseif k == 12 %Baseline but with 3 cointegration relations
        result = ecm_manu1(y,1,3);  
    end
    
    if k == 4 | k == 8 | k == 9 | k == 10 | k == 11 | k == 12
        vnames = strvcat('loghp','loggdp','r','ltv')
    elseif k == 15
        vnames = strvcat('loggdp','loghp','r','ltv') 
    elseif k == 16
        vnames = strvcat('loghp','loggdp','ltv','r')  
    elseif k == 17
        vnames = strvcat('loggdp','loghp','ltv','r')  
    elseif k == 23
        vnames = strvcat('loghp','logcons','logires','logibus','r','ltv')
    end
    
    prt(result,vnames) %print results in command window
     
    result(1).cointvec %Cointegration vector(s)
    result(1).data = y; %Put data matrix in result structure
    result(1).err = nan(size(result(1).resid,1),size(y,2));  
    %The above line pre-defines error matrix:
    %as many rows as one residual series and as many columns as data matrix
    for i = 1:size(y,2)
        result(1).err(:,i) = result(i).resid; %filling in the error matrix
    end
       
    if k >= 15 & k <= 17 %Alternative Cholesky orderings
        result(1).chind = chind(k - 13,:); %Take the relevant row of Cholesky indices into impulsovecm_manu.m
    end
    result(1).k = k; %Robustness index, to take into impulsovecm_manu.m
          
    bootstrap = 1; %1 if generate N sets of artificial data from
    %re-sampled residuals so as to construct confidence intervals; 0 otherwise
     
    if k == 15 | k == 16 | k == 17
        bootstrap = 0; %No bootstrap for Cholesky orderings
    end
    
    if bootstrap == 1
        %Do bootstrapping for all VECM specifications
        
        %Generate N artificial samples for bootstrapping
        N = 100000; %Number of artificial samples
        strap = bootstrap2(result,N); %Generate N samples from VECM estimates and errors
        neqs = strap.neqs; nobs = strap.nobs; nlag = strap.nlag;
        yhp_boot = nan(size(strap.Xmat,1),1); ygdp_boot = nan(size(strap.Xmat,1),1);
        yltv_boot = nan(size(strap.Xmat,1),1); yr_boot = nan(size(strap.Xmat,1),1);
        ycons_boot = nan(size(strap.Xmat,1),1); yires_boot = nan(size(strap.Xmat,1),1);
        yibus_boot = nan(size(strap.Xmat,1),1); yrent_boot = nan(size(strap.Xmat,1),1);
        ypbouw_boot = nan(size(strap.Xmat,1),1); ytransac_boot = nan(size(strap.Xmat,1),1);
        if k == 4 | k == 8 | k == 9 | k == 10 | k == 11 | k == 12
            for i = 1:N
                yhp_boot(:,i) = strap.Xmat(:,i*neqs-3);
                ygdp_boot(:,i) = strap.Xmat(:,i*neqs-2);
                yr_boot(:,i) = strap.Xmat(:,i*neqs-1);
                yltv_boot(:,i) = strap.Xmat(:,i*neqs);
            end
        elseif k == 23
            for i = 1:N
                yhp_boot(:,i) = strap.Xmat(:,i*neqs-5);
                ycons_boot(:,i) = strap.Xmat(:,i*neqs-4);
                yires_boot(:,i) = strap.Xmat(:,i*neqs-3);
                yibus_boot(:,i) = strap.Xmat(:,i*neqs-2);
                yr_boot(:,i) = strap.Xmat(:,i*neqs-1);
                yltv_boot(:,i) = strap.Xmat(:,i*neqs);
            end
        end
        
        %Estimate VECM on each of N samples and store IRFs for each of N
        %samples
        allstrapresult = cell(N,1);
        for j = 1:N
            if k == 4 | k == 8 | k == 9 | k == 10 | k == 11 | k == 12
                yboot = [yhp_boot(:,j) ygdp_boot(:,j) yr_boot(:,j) yltv_boot(:,j)];
            elseif k == 23
                yboot = [yhp_boot(:,j) ycons_boot(:,j) yires_boot(:,j) yibus_boot(:,j) yr_boot(:,j) yltv_boot(:,j)];
            end
            
            if k == 4 | k == 23
                strapresult = ecm_manu1(yboot,1,1);
            elseif k == 8
                strapresult = ecm_manu1(yboot,2,1);
            elseif k == 9
                strapresult = ecm_manu1(yboot,3,1);
            elseif k == 10
                strapresult = ecm_manu1(yboot,4,1);
            elseif k == 11
                strapresult = ecm_manu1(yboot,1,2);
            elseif k == 12
                strapresult = ecm_manu1(yboot,1,3);           
            end
            
           strapresult(1).cointvec; %cointvec if there is at least one
            strapresult(1).data = yboot; %Put data matrix in result structure
            strapresult(1).err = nan(size(strapresult(1).resid,1),size(y,2));
            %The above line pre-defines error matrix:
            %as many rows as one residual series and as many columns as data matrix
            for i = 1:size(y,2)
                strapresult(1).err(:,i) = strapresult(i).resid; %filling in the error matrix
            end
            strapresult(1).k = k; %Robustness index, to take into impulsovecm_manu.m
            %strapresult(1).ltvmat = ltvmat;
            %strapresult(1).ltimat = ltimat; %Store these two matrices as I will use their size in impulsovecm_manu.m
            allstrapresult{j}=strapresult; %jth sheet of array contains bootstrap results for jth sample
        end
    end
 
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%
    %Produce impulse responses to the LTV shock

    shockspec = 3; 
    %3 if Dutch scenario;
    %5 for ten heterogeneity (sigma) impulse paths;
   
    T = 10; %impulse response horizon (for long-run shocks), defined as T years after the shock
    Tplot = 10; %impulse response horizon to plot (for heterogeneity)
    impuls = zeros(T+1,1); %T+1 IRFs: shock and T years.
    
    if shockspec == 3 %Dutch scenario
        
        %Case d1ep: baseline distribution, ex post with exceptions
        impuls = zeros(T+1,1); %T+1 IRFs: shock and T years.
        impuls(1:size(d1ep,1)) = d1ep/100; %Dutch scenario, case d1ep
        irfres_d1ep=impulsovecm_manu(result,T,impuls,shockspec);
         
        if bootstrap == 1 
            
            %Prepare to fetch responses of all variables that could be in the VECM
            %(some will stay empty for now)
            strapirf_ltv_d1ep = nan(T+1,N); strapirf_hp_d1ep = nan(T+1,N);
            strapirf_gdp_d1ep = nan(T+1,N); strapirf_r_d1ep = nan(T+1,N);
            strapirf_rent_d1ep = nan(T+1,N); strapirf_pbouw_d1ep = nan(T+1,N);
            strapirf_cons_d1ep = nan(T+1,N); strapirf_ires_d1ep = nan(T+1,N);
            strapirf_ibus_d1ep = nan(T+1,N); strapirf_transac_d1ep = nan(T+1,N);
            strapirf_dltv_d1ep = nan(T+1,N); strapirf_dhp_d1ep = nan(T+1,N);
            strapirf_dgdp_d1ep = nan(T+1,N); strapirf_dr_d1ep = nan(T+1,N);

            for j = 1:N
                strapresult = allstrapresult{j};
                strapirfres_d1ep=impulsovecm_manu(strapresult,T,impuls,shockspec);
                
                if k == 4 | k == 8 | k == 9 | k == 10 | k == 11 | k == 12
                    strapirf_hp_d1ep(:,j) = strapirfres_d1ep.resplev(:,1);
                    strapirf_gdp_d1ep(:,j) = strapirfres_d1ep.resplev(:,2);
                    strapirf_r_d1ep(:,j) = strapirfres_d1ep.resplev(:,3);
                    strapirf_ltv_d1ep(:,j) = strapirfres_d1ep.resplev(:,4);  
                elseif k == 23
                    strapirf_hp_d1ep(:,j) = strapirfres_d1ep.resplev(:,1);
                    strapirf_cons_d1ep(:,j) = strapirfres_d1ep.resplev(:,2);
                    strapirf_ires_d1ep(:,j) = strapirfres_d1ep.resplev(:,3);
                    strapirf_ibus_d1ep(:,j) = strapirfres_d1ep.resplev(:,4);
                    strapirf_r_d1ep(:,j) = strapirfres_d1ep.resplev(:,5);
                    strapirf_ltv_d1ep(:,j) = strapirfres_d1ep.resplev(:,6);
                end
                clear strapirfres_d1ep %This is only for one of N specifications
            end

            pctstrap_ltv_d1ep = prctile(strapirf_ltv_d1ep,[2.5 5 50 95 97.5],2);
            pctstrap_hp_d1ep = prctile(strapirf_hp_d1ep,[2.5 5 50 94 97.5],2);
            if k == 4  | k == 8 | k == 9 | k == 10 | k == 11 | k == 12 
                pctstrap_gdp_d1ep = prctile(strapirf_gdp_d1ep,[2.5 5 50 95 97.5],2);
            end
            if k == 4 | k == 8 | k == 9 | k == 10 | k == 11 | k == 12 | k == 23
                pctstrap_r_d1ep = prctile(strapirf_r_d1ep,[2.5 5 50 95 97.5],2);
            end
            if k == 23
                pctstrap_cons_d1ep = prctile(strapirf_cons_d1ep,[2.5 5 50 95 97.5],2);
                pctstrap_ires_d1ep = prctile(strapirf_ires_d1ep,[2.5 5 50 95 97.5],2);
                pctstrap_ibus_d1ep = prctile(strapirf_ibus_d1ep,[2.5 5 50 95 97.5],2); 
            end
        end        
        
        %Case d1ea: baseline distribution, ex ante with exceptions 0.26
        %(with imposition fully attributed to 2012.)
        impuls = zeros(T+1,1); %T+1 IRFs: shock and T years.
        impuls(1:size(d1ea,1)) = d1ea/100; %Dutch scenario, case d1ea
        irfres_d1ea=impulsovecm_manu(result,T,impuls,shockspec);
        
        if bootstrap == 1 
            
            %Prepare to fetch responses of all variables that could be in the VECM
            %(some will stay empty for now)
            strapirf_ltv_d1ea = nan(T+1,N); strapirf_hp_d1ea = nan(T+1,N);
            strapirf_gdp_d1ea = nan(T+1,N); strapirf_r_d1ea = nan(T+1,N);
            strapirf_rent_d1ea = nan(T+1,N); strapirf_pbouw_d1ea = nan(T+1,N);
            strapirf_cons_d1ea = nan(T+1,N); strapirf_ires_d1ea = nan(T+1,N);
            strapirf_ibus_d1ea = nan(T+1,N); strapirf_transac_d1ea = nan(T+1,N);
            strapirf_dltv_d1ea = nan(T+1,N); strapirf_dhp_d1ea = nan(T+1,N);
            strapirf_dgdp_d1ea = nan(T+1,N); strapirf_dr_d1ea = nan(T+1,N);
            for j = 1:N
                strapresult = allstrapresult{j};
                strapirfres_d1ea=impulsovecm_manu(strapresult,T,impuls,shockspec);
                
                if k == 4 | k == 8 | k == 9 | k == 10 | k == 11 | k == 12
                    strapirf_hp_d1ea(:,j) = strapirfres_d1ea.resplev(:,1);
                    strapirf_gdp_d1ea(:,j) = strapirfres_d1ea.resplev(:,2);
                    strapirf_r_d1ea(:,j) = strapirfres_d1ea.resplev(:,3);
                    strapirf_ltv_d1ea(:,j) = strapirfres_d1ea.resplev(:,4);
                elseif k == 23
                    strapirf_hp_d1ea(:,j) = strapirfres_d1ea.resplev(:,1);
                    strapirf_cons_d1ea(:,j) = strapirfres_d1ea.resplev(:,2);
                    strapirf_ires_d1ea(:,j) = strapirfres_d1ea.resplev(:,3);
                    strapirf_ibus_d1ea(:,j) = strapirfres_d1ea.resplev(:,4);
                    strapirf_r_d1ea(:,j) = strapirfres_d1ea.resplev(:,5);
                    strapirf_ltv_d1ea(:,j) = strapirfres_d1ea.resplev(:,6);
                end
                clear strapirfres_d1ea %This is only for one of N specifications
            end
            pctstrap_ltv_d1ea = prctile(strapirf_ltv_d1ea,[2.5 5 50 95 97.5],2);
            pctstrap_hp_d1ea = prctile(strapirf_hp_d1ea,[2.5 5 50 94 97.5],2);
            if k == 4 | k == 8 | k == 9 | k == 10 | k == 11 | k == 12
                pctstrap_gdp_d1ea = prctile(strapirf_gdp_d1ea,[2.5 5 50 95 97.5],2);
            end
            if k == 4 | k == 8 | k == 9 | k == 10 | k == 11 | k == 12 | k == 23
                pctstrap_r_d1ea = prctile(strapirf_r_d1ea,[2.5 5 50 95 97.5],2);
            end   
            if k == 23
                pctstrap_cons_d1ea = prctile(strapirf_cons_d1ea,[2.5 5 50 95 97.5],2);
                pctstrap_ires_d1ea = prctile(strapirf_ires_d1ea,[2.5 5 50 95 97.5],2);
                pctstrap_ibus_d1ea = prctile(strapirf_ibus_d1ea,[2.5 5 50 95 97.5],2);
            end
        end
        
        %Case d1ea00: baseline distribution, ex ante without exceptions
        impuls = zeros(T+1,1); %T+1 IRFs: shock and T years.
        impuls(1:size(d1ea00,1)) = d1ea00/100;
        irfres_d1ea00=impulsovecm_manu(result,T,impuls,shockspec);
        
        %Case d1ea20: baseline distribution, ex ante with exceptions 0.20
        impuls = zeros(T+1,1); %T+1 IRFs: shock and T years.
        impuls(1:size(d1ea20,1)) = d1ea20/100;
        irfres_d1ea20=impulsovecm_manu(result,T,impuls,shockspec);
        
        %Case d1ea30: baseline distribution, ex ante with exceptions 0.30
        impuls = zeros(T+1,1); %T+1 IRFs: shock and T years.
        impuls(1:size(d1ea30,1)) = d1ea30/100;
        irfres_d1ea30=impulsovecm_manu(result,T,impuls,shockspec);
        
        %Case d1ea50: baseline distribution, ex ante with exceptions 0.50
        impuls = zeros(T+1,1); %T+1 IRFs: shock and T years.
        impuls(1:size(d1ea50,1)) = d1ea50/100;
        irfres_d1ea50=impulsovecm_manu(result,T,impuls,shockspec);
        
        %Case d2ea: 2010M4-M12 as pre-cap, ex ante with exceptions 0.26
        impuls = zeros(T+1,1); %T+1 IRFs: shock and T years.
        impuls(1:size(d2ea,1)) = d2ea/100;
        irfres_d2ea=impulsovecm_manu(result,T,impuls,shockspec);
        
        %Case d3ea: 2010 as pre-cap, ex ante with exceptions 0.26
        impuls = zeros(T+1,1); %T+1 IRFs: shock and T years.
        impuls(1:size(d3ea,1)) = d3ea/100;
        irfres_d3ea=impulsovecm_manu(result,T,impuls,shockspec);
        
        %Case d1easp: baseline distribution, ex ante with exceptions 0.26
        %and sprinkling
        impuls = zeros(T+1,1); %T+1 IRFs: shock and T years.
        impuls(1:size(d1easp,1)) = d1easp/100;
        irfres_d1easp=impulsovecm_manu(result,T,impuls,shockspec);

        %Case d1ea8: baseline distribution, ex ante with exceptions 0.26
        %assign entire imposition  to 2012, i.e. baseline as of Sep 2022
        %This is still the overall baseline, but written with 8 annual shocks so 2011 is first shock
        impuls = zeros(T+1,1); %T+1 IRFs: shock and T years.
        impuls(1:size(d1ea8,1)) = d1ea8/100; %Dutch scenario, case d1ea8
        irfres_d1ea8=impulsovecm_manu(result,T,impuls,shockspec);
        
        %Case d1ea8our: baseline distribution, ex ante with exceptions 0.26
        %our way: assign 5/12th to 2011 and 7/12th to 2012
        %written with 8 annual shocks so 2011 is first shock 
        %BASELINE OCT 2023
        impuls = zeros(T+1,1); %T+1 IRFs: shock and T years.
        impuls(1:size(d1ea8our,1)) = d1ea8our/100; 
        irfres_d1ea8our=impulsovecm_manu(result,T,impuls,shockspec);
                
        %Case d1ea9our: baseline distribution, ex ante with exceptions 0.26
        %assign 5/12th to 2011 and 7/12th to 2012
        %written with 8 annual shocks so 2011 is first shock, but writing
        %9th year as zero (this should be the same as the above)
        impuls = zeros(T+1,1); %T+1 IRFs: shock and T years.
        impuls(1:size(d1ea9our,1)) = d1ea9our/100; 
        irfres_d1ea9our=impulsovecm_manu(result,T,impuls,shockspec);

        %Case d2ea9our: distribution(s) last nine months, ex ante with exceptions 0.26
        %our way: assign 5/12th to 2011 and 7/12th to 2012
        %written with 8 annual shocks so 2011 is first shock, but writing
        %9th year as zero
        impuls = zeros(T+1,1); %T+1 IRFs: shock and T years.
        impuls(1:size(d2ea9our,1)) = d2ea9our/100; 
        irfres_d2ea9our=impulsovecm_manu(result,T,impuls,shockspec);

        %Case d3ea9our: 2010 distribution, ex ante with exceptions 0.26
        %our way: assign 5/12th to 2011 and 7/12th to 2012
        %written with 8 annual shocks so 2011 is first shock, but writing
        %9th year as zero
        impuls = zeros(T+1,1); %T+1 IRFs: shock and T years.
        impuls(1:size(d3ea9our,1)) = d3ea9our/100;
        irfres_d3ea9our=impulsovecm_manu(result,T,impuls,shockspec);

        if bootstrap == 1
            
            %Prepare to fetch responses of all variables that could be in the VECM
            %(some will stay empty for now)
            strapirf_ltv_d1ea9our = nan(T+1,N); strapirf_hp_d1ea9our = nan(T+1,N);
            strapirf_gdp_d1ea9our = nan(T+1,N); strapirf_r_d1ea9our = nan(T+1,N);
            strapirf_rent_d1ea9our = nan(T+1,N); strapirf_pbouw_d1ea9our = nan(T+1,N);
            strapirf_cons_d1ea9our = nan(T+1,N); strapirf_ires_d1ea9our = nan(T+1,N);
            strapirf_ibus_d1ea9our = nan(T+1,N); strapirf_transac_d1ea9our = nan(T+1,N);
            strapirf_dltv_d1ea9our = nan(T+1,N); strapirf_dhp_d1ea9our = nan(T+1,N);
            strapirf_dgdp_d1ea9our = nan(T+1,N); strapirf_dr_d1ea9our = nan(T+1,N);
            for j = 1:N
                strapresult = allstrapresult{j};
                strapirfres_d1ea9our=impulsovecm_manu(strapresult,T,impuls,shockspec);
                
                if k == 4 | k == 8 | k == 9 | k == 10 | k == 11 | k == 12
                    strapirf_hp_d1ea9our(:,j) = strapirfres_d1ea9our.resplev(:,1);
                    strapirf_gdp_d1ea9our(:,j) = strapirfres_d1ea9our.resplev(:,2);
                    strapirf_r_d1ea9our(:,j) = strapirfres_d1ea9our.resplev(:,3);
                    strapirf_ltv_d1ea9our(:,j) = strapirfres_d1ea9our.resplev(:,4);
                elseif k == 23
                    strapirf_hp_d1ea9our(:,j) = strapirfres_d1ea9our.resplev(:,1);
                    strapirf_cons_d1ea9our(:,j) = strapirfres_d1ea9our.resplev(:,2);
                    strapirf_ires_d1ea9our(:,j) = strapirfres_d1ea9our.resplev(:,3);
                    strapirf_ibus_d1ea9our(:,j) = strapirfres_d1ea9our.resplev(:,4);
                    strapirf_r_d1ea9our(:,j) = strapirfres_d1ea9our.resplev(:,5);
                    strapirf_ltv_d1ea9our(:,j) = strapirfres_d1ea9our.resplev(:,6);
                end
                clear strapirfres_d1ea9our %This is only for one of N specifications
            end
            pctstrap_ltv_d1ea9our = prctile(strapirf_ltv_d1ea9our,[2.5 5 50 95 97.5],2);
            pctstrap_hp_d1ea9our = prctile(strapirf_hp_d1ea9our,[2.5 5 50 94 97.5],2);
            if k == 4 | k == 8 | k == 9 | k == 10 | k == 11 | k == 12
                pctstrap_gdp_d1ea9our = prctile(strapirf_gdp_d1ea9our,[2.5 5 50 95 97.5],2);
            end 
            if k == 4 | k == 8 | k == 9 | k == 10 | k == 11 | k == 12 | k == 23
                pctstrap_r_d1ea9our = prctile(strapirf_r_d1ea9our,[2.5 5 50 95 97.5],2);
            end
            if k == 23
                pctstrap_cons_d1ea9our = prctile(strapirf_cons_d1ea9our,[2.5 5 50 95 97.5],2);
                pctstrap_ires_d1ea9our = prctile(strapirf_ires_d1ea9our,[2.5 5 50 95 97.5],2);
                pctstrap_ibus_d1ea9our = prctile(strapirf_ibus_d1ea9our,[2.5 5 50 95 97.5],2);
            end
        end
    
        %Case d1ea9del: baseline distribution, ex ante with exceptions 0.26
        %account for delays in mortgage origination
        %9 years' worth of impulses
        impuls = zeros(T+1,1); %T+1 IRFs: shock and T years.
        impuls(1:size(d1ea9del,1)) = d1ea9del/100;
        irfres_d1ea9del=impulsovecm_manu(result,T,impuls,shockspec);

        %Case d1easp9our: baseline distribution, ex ante with exceptions 0.26
        %sprinkling
        impuls = zeros(T+1,1); %T+1 IRFs: shock and T years.
        impuls(1:size(d1easp9our,1)) = d1easp9our/100; 
        irfres_d1easp9our=impulsovecm_manu(result,T,impuls,shockspec);

        %Case d1easp9del: baseline distribution, ex ante with exceptions 0.26
        %sprinkling and origination delays
        impuls = zeros(T+1,1); %T+1 IRFs: shock and T years.
        impuls(1:size(d1easp9del,1)) = d1easp9del/100; 
        irfres_d1easp9del=impulsovecm_manu(result,T,impuls,shockspec);

        %Case d1easp9bigdel: baseline distribution, ex ante with exceptions 0.26
        %sprinkling and long origination delays
        impuls = zeros(T+1,1); %T+1 IRFs: shock and T years.
        impuls(1:size(d1easp9bigdel,1)) = d1easp9bigdel/100; 
        irfres_d1easp9bigdel=impulsovecm_manu(result,T,impuls,shockspec);

        %Case d1ea9bigdel: baseline distribution, ex ante with exceptions 0.26
        %and long origination delays
        impuls = zeros(T+1,1); %T+1 IRFs: shock and T years.
        impuls(1:size(d1ea9bigdel,1)) = d1ea9bigdel/100; 
        irfres_d1ea9bigdel=impulsovecm_manu(result,T,impuls,shockspec);

        %Case d1ea8ref: baseline distribution, ex ante with exceptions 0.26
        impuls = zeros(T+1,1); %T+1 IRFs: shock and T years.
        impuls(1:size(d1ea8ref,1)) = d1ea8ref/100; 
        irfres_d1ea8ref=impulsovecm_manu(result,T,impuls,shockspec);
       
        if k == 4 | k == 8 | k == 9 | k == 10 | k == 11 | k == 12 %k== 4: Baseline  [loghp loggdp rh_ea1 ltv_comb_nz].
               
            %GRAPH FOR PAPER
            %Plot IRF with ex ante and exceptions (5/12 to 2011): baseline
            %vs. baseline + sprinkling and/OR origination delays
            time=[0:T]';
            zeroline=zeros(size(time));
            figure; % We represent  100*ln(x) or 100*x, so it is in pct or pps
            subplot(2,2,4);
            h2=plot(time,100*irfres_d1ea9our(1).resplev(:,1),time,100*irfres_d1ea9del(1).resplev(:,1),time,100*irfres_d1easp9our(1).resplev(:,1),time,100*irfres_d1easp9del(1).resplev(:,1),time,zeroline,'k')%, title('Real House Prices'), ylabel('%'), xlabel('years since first shock')
            set(h2(1),'Color','b','Linewidth',2,'LineStyle','--'); set(h2(2),'Color','[0.0 0.7 0.0]','LineWidth',2); set(h2(3),'Color','[0.95,0.025,0.098]','LineWidth',2); set(h2(4),'Color','k','LineWidth',2); set(h2(5),'LineWidth',1);
            grid on; axis([0 T -4 0]); yticks([-4:1:0])
            set(gca,'FontSize',17)
            title('Real House Prices','FontSize',22),ylabel('%','FontSize',18), xlabel('years since first shock','FontSize',18)
            
            subplot(2,2,3);
            h4=plot(time,100*irfres_d1ea9our(1).resplev(:,2),time,100*irfres_d1ea9del(1).resplev(:,2),time,100*irfres_d1easp9our(1).resplev(:,2),time,100*irfres_d1easp9del(1).resplev(:,2),time,zeroline,'k')%, title({'Real Gross Domestic Product'}), ylabel('%'), xlabel('years since first shock')
            set(h4(1),'Color','b','Linewidth',2,'LineStyle','--'); set(h4(2),'Color','[0.0 0.7 0.0]','LineWidth',2); set(h4(3),'Color','[0.95,0.025,0.098]','LineWidth',2); set(h4(4),'Color','k','LineWidth',2); set(h4(5),'LineWidth',1);
            grid on; axis([0 T -2 0]); yticks([-2:1:0])
            set(gca,'FontSize',17)
            title('Real Gross Domestic Product','FontSize',22),ylabel('%','FontSize',18), xlabel('years since first shock','FontSize',18)
            
            subplot(2,2,2);
            h3=plot(time,100*irfres_d1ea9our(1).resplev(:,3),time,100*irfres_d1ea9del(1).resplev(:,3),time,100*irfres_d1easp9our(1).resplev(:,3),time,100*irfres_d1easp9del(1).resplev(:,3),time,zeroline,'k')%, title({'','Real Mortgage Rate'}), ylabel('pps'), xlabel('years since first shock'), legend('ex ante','ex post','Location','NorthEast')
            set(h3(1),'Color','b','Linewidth',2,'LineStyle','--'); set(h3(2),'Color','[0.0 0.7 0.0]','LineWidth',2); set(h3(3),'Color','[0.95,0.025,0.098]','LineWidth',2); set(h3(4),'Color','k','LineWidth',2); set(h3(5),'LineWidth',1);
            grid on; axis([0 T -1 1]); yticks([-1:0.5:1])
            set(gca,'FontSize',17)
            title('Real Mortgage Rate','FontSize',22),ylabel('pps','FontSize',18), xlabel('years since first shock','FontSize',18)   
            
            subplot(2,2,1);
            h1=plot(time,100*irfres_d1ea9our(1).resplev(:,4),time,100*irfres_d1ea9del(1).resplev(:,4),time,100*irfres_d1easp9our(1).resplev(:,4),time,100*irfres_d1easp9del(1).resplev(:,4),time,zeroline,'k')%,title('Loan-to-Value Ratio'),ylabel('pps'), xlabel('years since first shock')
            set(h1(1),'Color','b','Linewidth',2,'LineStyle','--'); set(h1(2),'Color','[0.0 0.7 0.0]','LineWidth',2); set(h1(3),'Color','[0.95,0.025,0.098]','LineWidth',2); set(h1(4),'Color','k','LineWidth',2); set(h1(5),'LineWidth',1);
            grid on; axis([0 T -2 0]); yticks([-2:1:0])
            set(gca,'FontSize',17)
            title('Loan-to-Value Ratio','FontSize',22),ylabel('pps','FontSize',18), xlabel('years since first shock','FontSize',18)
            legend('baseline','with origination delays','with rounding','with delays & rounding','Location','NorthEast','FontSize',16)
                             
                               
            %Plot IRFs for 3 different distributions
            cd('JMCB_replic_deveirman-jong\VECM estimation');
            time=[0:T]';
            zeroline=zeros(size(time));
            figure; % We represent  100*ln(x) or 100*x, so it is in pct or pps
            subplot(2,2,4);
            h2=plot(time,100*irfres_d1ea(1).resplev(:,1),time,100*irfres_d2ea(1).resplev(:,1),time,100*irfres_d3ea(1).resplev(:,1),time,zeroline,'k')%, title('Real House Prices'), ylabel('%'), xlabel('years since first shock')
            set(h2(1),'LineWidth',3); set(h2(2),'LineWidth',3); set(h2(3),'LineWidth',3); grid on; axis([0 T -5 0]); yticks([-5:1:0])
            set(gca,'FontSize',17)
            title('Real House Prices','FontSize',22),ylabel('%','FontSize',18), xlabel('years since first shock','FontSize',18)
                           
            subplot(2,2,3);
            h4=plot(time,100*irfres_d1ea(1).resplev(:,2),time,100*irfres_d2ea(1).resplev(:,2),time,100*irfres_d3ea(1).resplev(:,2),time,zeroline,'k')%, title({'','Real Gross Domestic Product'}), ylabel('%'), xlabel('years since first shock')
            set(h4(1),'LineWidth',3); set(h4(2),'LineWidth',3); set(h4(3),'LineWidth',3); grid on; axis([0 T -2 0]); yticks([-2:1:0])
            set(gca,'FontSize',17)
            title('Real Gross Domestic Product','FontSize',22),ylabel('%','FontSize',18), xlabel('years since first shock','FontSize',18)
            
            subplot(2,2,2);
            h3=plot(time,100*irfres_d1ea(1).resplev(:,3),time,100*irfres_d2ea(1).resplev(:,3),time,100*irfres_d3ea(1).resplev(:,3),time,zeroline,'k')%, title({'Real Mortgage Rate'}), ylabel('pps'), xlabel('years since first shock'), legend('2010M8-2011M7','2010M4-2010M12','2010','Location','NorthEast')
            set(h3(1),'LineWidth',3); set(h3(2),'LineWidth',3); set(h3(3),'LineWidth',3); grid on; axis([0 T -1 1]); yticks([-1:0.5:1])
            set(gca,'FontSize',17)
            title('Real Mortgage Rate','FontSize',22),ylabel('pps','FontSize',18), xlabel('years since first shock','FontSize',18), legend('Baseline distributions','Nine-month distributions','Pre-cap distribution 2010','Location','NorthEast', 'FontSize',16,'NumColumns',1)  
            
            subplot(2,2,1);
            h1=plot(time,100*irfres_d1ea(1).resplev(:,4),time,100*irfres_d2ea(1).resplev(:,4),time,100*irfres_d3ea(1).resplev(:,4),time,zeroline,'k')%,title('Loan-to-Value Ratio'),ylabel('pps'), xlabel('years since first shock')
            set(h1(1),'LineWidth',3); set(h1(2),'LineWidth',3); set(h1(3),'LineWidth',3); grid on; axis([0 T -3 0]); yticks([-3:1:0])
            set(gca,'FontSize',17)
            title('Loan-to-Value Ratio','FontSize',22),ylabel('pps','FontSize',18), xlabel('years since first shock','FontSize',18)
            
            %Plot IRFs with different fractions of exceptions
            time=[0:T]';
            zeroline=zeros(size(time));
            figure; % We represent  100*ln(x) or 100*x, so it is in pct or pps
            subplot(2,2,4);
            h2=plot(time,100*irfres_d1ea(1).resplev(:,1),time,100*irfres_d1ea00(1).resplev(:,1),time,100*irfres_d1ea20(1).resplev(:,1),time,100*irfres_d1ea30(1).resplev(:,1),time,100*irfres_d1ea50(1).resplev(:,1),time,zeroline,'k')%, title('Real House Prices'), ylabel('%'), xlabel('years since first shock')
            set(h2(1),'LineWidth',3); set(h2(2),'LineWidth',3); set(h2(3),'LineWidth',3); set(h2(4),'LineWidth',3); set(h2(5),'LineWidth',3); grid on; axis([0 T -5 0]); yticks([-5:1:0])
            set(gca,'FontSize',17)
            title('Real House Prices','FontSize',22),ylabel('%','FontSize',18), xlabel('years since first shock','FontSize',18)
                       
            subplot(2,2,3);
            h4=plot(time,100*irfres_d1ea(1).resplev(:,2),time,100*irfres_d1ea00(1).resplev(:,2),time,100*irfres_d1ea20(1).resplev(:,2),time,100*irfres_d1ea30(1).resplev(:,2),time,100*irfres_d1ea50(1).resplev(:,2),time,zeroline,'k')%, title({'','Real Gross Domestic Product'}), ylabel('%'), xlabel('years since first shock')
            set(h4(1),'LineWidth',3); set(h4(2),'LineWidth',3); set(h4(3),'LineWidth',3); set(h4(4),'LineWidth',3); set(h4(5),'LineWidth',3); grid on; axis([0 T -2 0]); yticks([-2:1:0])
            set(gca,'FontSize',17)
            title('Real Gross Domestic Product','FontSize',22),ylabel('%','FontSize',18), xlabel('years since first shock','FontSize',18)
            
            subplot(2,2,2);
            h3=plot(time,100*irfres_d1ea(1).resplev(:,3),time,100*irfres_d1ea00(1).resplev(:,3),time,100*irfres_d1ea20(1).resplev(:,3),time,100*irfres_d1ea30(1).resplev(:,3),time,100*irfres_d1ea50(1).resplev(:,3),time,zeroline,'k')%, title({'','Real Mortgage Rate'}), ylabel('pps'), xlabel('years since first shock'), legend('0.26','0.00','0.20','0.30','0.50','Location','NorthEast')
            set(h3(1),'LineWidth',3); set(h3(2),'LineWidth',3); set(h3(3),'LineWidth',3); set(h3(4),'LineWidth',3); set(h3(5),'LineWidth',3); grid on; axis([0 T -1 1]); yticks([-1:0.5:1])
            set(gca,'FontSize',17)
            title('Real Mortgage Rate','FontSize',22),ylabel('pps','FontSize',18), xlabel('years since first shock','FontSize',18), legend('26% exceptions','no exceptions','20%','30%','50%','Location','NorthEast','FontSize',16,'NumColumns',3)  
            
            subplot(2,2,1);
            h1=plot(time,100*irfres_d1ea(1).resplev(:,4),time,100*irfres_d1ea00(1).resplev(:,4),time,100*irfres_d1ea20(1).resplev(:,4),time,100*irfres_d1ea30(1).resplev(:,4),time,100*irfres_d1ea50(1).resplev(:,4),time,zeroline,'k')%,title('Loan-to-Value Ratio'),ylabel('pps'), xlabel('years since first shock')
            set(h1(1),'LineWidth',3); set(h1(2),'LineWidth',3); set(h1(3),'LineWidth',3); set(h1(4),'LineWidth',3); set(h1(5),'LineWidth',3); grid on; axis([0 T -3 0]); yticks([-3:1:0])
            set(gca,'FontSize',17)
            title('Loan-to-Value Ratio','FontSize',22),ylabel('pps','FontSize',18), xlabel('years since first shock','FontSize',18)                             
            
        end
        
        
        if bootstrap ==1
                        
            if k == 4  %k=4 is baseline [loghp loggdp rh_ea1 ltv_comb_nz].
                
                %Plot IRFs for baseline: baseline distribution, ex ante with 0.26
                %exceptions. 
                %bootstrapped 90% and 95% confidence intervals
                cd('JMCB_replic_deveirman-jong\VECM estimation');
                time=[0:T]'; zeroline=zeros(size(time));
                figure; % We represent  100*ln(x) or 100*x, so it is in pct or pps
                
                subplot(2,2,1)
                lo95ltv = 100*pctstrap_ltv_d1ea9our(:,1); lo90ltv = 100*pctstrap_ltv_d1ea9our(:,2);
                up90ltv = 100*pctstrap_ltv_d1ea9our(:,4); up95ltv = 100*pctstrap_ltv_d1ea9our(:,5);
                x = [lo95ltv lo90ltv-lo95ltv up90ltv-lo90ltv up95ltv-up90ltv];
                hold on
                harea = area(time,x,-3.5,'EdgeColor','none');
                set(harea(1),'Facecolor',[1 1 1],'FaceAlpha',0)
                set(harea(2),'Facecolor',[0.3 0.3 0.6],'FaceAlpha',0.5)
                set(harea(3),'Facecolor',[0 0 1],'FaceAlpha',0.5)
                set(harea(4),'Facecolor',[0.3 0.3 0.6],'FaceAlpha',0.5)
                plot(time,100*irfres_d1ea9our(1).resplev(:,4),'Color',[0 0 1],'LineWidth',3)
                grid on; ax=gca; ax.GridAlpha = 0.45;
                axis([0 T -3.5 1]); yticks([-3:1:1])
                plot(time,zeroline,'k','LineWidth',1.5)
                hold off
                set(gca,'FontSize',17)
                title('Loan-to-Value Ratio','FontSize',22),ylabel('pps','FontSize',18), xlabel('years since first shock','FontSize',18)
                subplot(2,2,4);
                lo95hp = 100*pctstrap_hp_d1ea9our(:,1); lo90hp = 100*pctstrap_hp_d1ea9our(:,2);
                up90hp = 100*pctstrap_hp_d1ea9our(:,4); up95hp = 100*pctstrap_hp_d1ea9our(:,5);
                x = [lo95hp lo90hp-lo95hp up90hp-lo90hp up95hp-up90hp];
                hold on
                harea = area(time,x,-10,'EdgeColor','none');
                set(harea(1),'Facecolor',[1 1 1],'FaceAlpha',0)
                set(harea(2),'Facecolor',[0.3 0.3 0.6],'FaceAlpha',0.5)
                set(harea(3),'Facecolor',[0 0 1],'FaceAlpha',0.5)
                set(harea(4),'Facecolor',[0.3 0.3 0.6],'FaceAlpha',0.5)
                plot(time,100*irfres_d1ea9our(1).resplev(:,1),'Color',[0 0 1],'Linewidth',3)
                grid on; ax=gca; ax.GridAlpha = 0.45;
                axis([0 T -8 2]); yticks([-8:2:2])
                plot(time,zeroline,'k','LineWidth',1.5)
                hold off
                set(gca,'FontSize',17)
                title('Real House Prices','FontSize',22),ylabel('%','FontSize',18), xlabel('years since first shock','FontSize',18)
                subplot(2,2,3);
                %axes(ha(3));
                lo95gdp = 100*pctstrap_gdp_d1ea9our(:,1); lo90gdp = 100*pctstrap_gdp_d1ea9our(:,2);
                up90gdp = 100*pctstrap_gdp_d1ea9our(:,4); up95gdp = 100*pctstrap_gdp_d1ea9our(:,5);
                x = [lo95gdp lo90gdp-lo95gdp up90gdp-lo90gdp up95gdp-up90gdp];
                hold on
                harea = area(time,x,-3,'EdgeColor','none');
                set(harea(1),'Facecolor',[1 1 1],'FaceAlpha',0)
                set(harea(2),'Facecolor',[0.3 0.3 0.6],'FaceAlpha',0.5)
                set(harea(3),'Facecolor',[0 0 1],'FaceAlpha',0.5)
                set(harea(4),'Facecolor',[0.3 0.3 0.6],'FaceAlpha',0.5)
                plot(time,100*irfres_d1ea9our(1).resplev(:,2),'Color',[0 0 1],'Linewidth',3)
                grid on; ax=gca; ax.GridAlpha = 0.45;
                axis([0 T -3 1]); yticks([-3:1:1])
                plot(time,zeroline,'k','LineWidth',1.5)
                hold off
                set(gca,'FontSize',17)
                title('Real Gross Domestic Product','FontSize',22),ylabel('%','FontSize',18), xlabel('years since first shock','FontSize',18)
                subplot(2,2,2);
                lo95r = 100*pctstrap_r_d1ea9our(:,1); lo90r = 100*pctstrap_r_d1ea9our(:,2);
                up90r = 100*pctstrap_r_d1ea9our(:,4); up95r = 100*pctstrap_r_d1ea9our(:,5);
                x = [lo95r lo90r-lo95r up90r-lo90r up95r-up90r];
                hold on
                harea = area(time,x,-1,'EdgeColor','none');
                set(harea(1),'Facecolor',[1 1 1],'FaceAlpha',0)
                set(harea(2),'Facecolor',[0.3 0.3 0.6],'FaceAlpha',0.5)
                set(harea(3),'Facecolor',[0 0 1],'FaceAlpha',0.5)
                set(harea(4),'Facecolor',[0.3 0.3 0.6],'FaceAlpha',0.5)
                plot(time,100*irfres_d1ea9our(1).resplev(:,3),'Color',[0 0 1],'Linewidth',3)
                grid on; ax=gca; ax.GridAlpha = 0.45;
                axis([0 T -1 1]); yticks([-1:0.5:1])
                plot(time,zeroline,'k','LineWidth',1.5)
                set(gca,'FontSize',17)
                legend([harea(2) harea(3)],'95% confidence interval','90% confidence interval','Location','SouthEast','FontSize',17)
                title('Real Mortgage Rate','FontSize',22),ylabel('pps','FontSize',18), xlabel('years since first shock','FontSize',18)
                hold off                           
                 
            elseif k == 8 | k == 9 | k == 10 | k == 11 | k == 12
                %Plot IRFs for baseline: baseline distribution, ex ante with 0.26
                %exceptions
                %bootstrapped 90% and 95% confidence intervals
                cd('JMCB_replic_deveirman-jong\VECM estimation');
                time=[0:T]'; zeroline=zeros(size(time));
                figure; % We represent  100*ln(x) or 100*x, so it is in pct or pps
                
                subplot(2,2,1)
                lo95ltv = 100*pctstrap_ltv_d1ea(:,1); lo90ltv = 100*pctstrap_ltv_d1ea(:,2);
                up90ltv = 100*pctstrap_ltv_d1ea(:,4); up95ltv = 100*pctstrap_ltv_d1ea(:,5);
                x = [lo95ltv lo90ltv-lo95ltv up90ltv-lo90ltv up95ltv-up90ltv];
                hold on
                harea = area(time,x,-9,'EdgeColor','none');
                set(harea(1),'Facecolor',[1 1 1],'FaceAlpha',0)
                set(harea(2),'Facecolor',[0.3 0.3 0.6],'FaceAlpha',0.5)
                set(harea(3),'Facecolor',[0 0 1],'FaceAlpha',0.5)
                set(harea(4),'Facecolor',[0.3 0.3 0.6],'FaceAlpha',0.5)
                plot(time,100*irfres_d1ea(1).resplev(:,4),'Color',[0 0 1],'LineWidth',3)
                grid on; ax=gca; ax.GridAlpha = 0.45;
                axis([0 T -9 0.5]); yticks([-9:3:0])
                plot(time,zeroline,'k','LineWidth',1.5)
                hold off
                set(gca,'FontSize',17)
                title('Loan-to-Value Ratio','FontSize',22),ylabel('pps','FontSize',18), xlabel('years since first shock','FontSize',18)
                subplot(2,2,4);
                %axes(ha(4));
                lo95hp = 100*pctstrap_hp_d1ea(:,1); lo90hp = 100*pctstrap_hp_d1ea(:,2);
                up90hp = 100*pctstrap_hp_d1ea(:,4); up95hp = 100*pctstrap_hp_d1ea(:,5);
                x = [lo95hp lo90hp-lo95hp up90hp-lo90hp up95hp-up90hp];
                hold on
                harea = area(time,x,-30,'EdgeColor','none');
                set(harea(1),'Facecolor',[1 1 1],'FaceAlpha',0)
                set(harea(2),'Facecolor',[0.3 0.3 0.6],'FaceAlpha',0.5)
                set(harea(3),'Facecolor',[0 0 1],'FaceAlpha',0.5)
                set(harea(4),'Facecolor',[0.3 0.3 0.6],'FaceAlpha',0.5)
                plot(time,100*irfres_d1ea(1).resplev(:,1),'Color',[0 0 1],'Linewidth',3)
                grid on; ax=gca; ax.GridAlpha = 0.45;
                axis([0 T -30 10]); yticks([-30:10:10])
                plot(time,zeroline,'k','LineWidth',1.5)
                hold off
                set(gca,'FontSize',17)
                title('Real house prices','FontSize',22),ylabel('%','FontSize',18), xlabel('years since first shock','FontSize',18)
                subplot(2,2,3);
                lo95gdp = 100*pctstrap_gdp_d1ea(:,1); lo90gdp = 100*pctstrap_gdp_d1ea(:,2);
                up90gdp = 100*pctstrap_gdp_d1ea(:,4); up95gdp = 100*pctstrap_gdp_d1ea(:,5);
                x = [lo95gdp lo90gdp-lo95gdp up90gdp-lo90gdp up95gdp-up90gdp];
                hold on
                harea = area(time,x,-8,'EdgeColor','none');
                set(harea(1),'Facecolor',[1 1 1],'FaceAlpha',0)
                set(harea(2),'Facecolor',[0.3 0.3 0.6],'FaceAlpha',0.5)
                set(harea(3),'Facecolor',[0 0 1],'FaceAlpha',0.5)
                set(harea(4),'Facecolor',[0.3 0.3 0.6],'FaceAlpha',0.5)
                plot(time,100*irfres_d1ea(1).resplev(:,2),'Color',[0 0 1],'Linewidth',3)
                grid on; ax=gca; ax.GridAlpha = 0.45;
                axis([0 T -8 3]); yticks([-8:2:3])
                plot(time,zeroline,'k','LineWidth',1.5)
                hold off
                set(gca,'FontSize',17)
                title('Real Gross Domestic Product','FontSize',22),ylabel('%','FontSize',18), xlabel('years since first shock','FontSize',18)
                subplot(2,2,2);
                %axes(ha(2));
                lo95r = 100*pctstrap_r_d1ea(:,1); lo90r = 100*pctstrap_r_d1ea(:,2);
                up90r = 100*pctstrap_r_d1ea(:,4); up95r = 100*pctstrap_r_d1ea(:,5);
                x = [lo95r lo90r-lo95r up90r-lo90r up95r-up90r];
                hold on
                harea = area(time,x,-3,'EdgeColor','none');
                set(harea(1),'Facecolor',[1 1 1],'FaceAlpha',0)
                set(harea(2),'Facecolor',[0.3 0.3 0.6],'FaceAlpha',0.5)
                set(harea(3),'Facecolor',[0 0 1],'FaceAlpha',0.5)
                set(harea(4),'Facecolor',[0.3 0.3 0.6],'FaceAlpha',0.5)
                plot(time,100*irfres_d1ea(1).resplev(:,3),'Color',[0 0 1],'Linewidth',3)
                grid on; ax=gca; ax.GridAlpha = 0.45;
                axis([0 T -3 3]); yticks([-3:1:3])
                plot(time,zeroline,'k','LineWidth',1.5)
                set(gca,'FontSize',17)
                legend([harea(2) harea(3)],'95% confidence interval','90% confidence interval','Location','SouthEast','FontSize',17)
                title('Real mortgage rate','FontSize',22),ylabel('pps','FontSize',18), xlabel('years since first shock','FontSize',18)
                hold off         

            elseif k == 23   %[loghp logcons logires logibus rh_ea1 ltv_comb_nz]
                
                %ASSIGN 5/12TH OF IMPOSITION TO 2011.  
                %Plot IRFs for baseline: baseline distribution, ex ante with 0.26
                %exceptions  
                %bootstrapped 90% and 95% confidence intervals
                cd('JMCB_replic_deveirman-jong\VECM estimation');
                time=[0:T]'; zeroline=zeros(size(time));
                figure; % We represent  100*ln(x) or 100*x, so it is in pct or pps
                
                subplot(2,3,1)
                lo95ltv = 100*pctstrap_ltv_d1ea9our(:,1); lo90ltv = 100*pctstrap_ltv_d1ea9our(:,2);
                up90ltv = 100*pctstrap_ltv_d1ea9our(:,4); up95ltv = 100*pctstrap_ltv_d1ea9our(:,5);
                x = [lo95ltv lo90ltv-lo95ltv up90ltv-lo90ltv up95ltv-up90ltv];
                hold on
                harea = area(time,x,5,'EdgeColor','none');
                set(harea(1),'Facecolor',[1 1 1],'FaceAlpha',0)
                set(harea(2),'Facecolor',[0.3 0.3 0.6],'FaceAlpha',0.5)
                set(harea(3),'Facecolor',[0 0 1],'FaceAlpha',0.5)
                set(harea(4),'Facecolor',[0.3 0.3 0.6],'FaceAlpha',0.5)
                plot(time,100*irfres_d1ea9our(1).resplev(:,6),'Color',[0 0 1],'LineWidth',3)
                grid on; ax=gca; ax.GridAlpha = 0.45;
                axis([0 T -3.5 1]); yticks([-3:1:1])
                plot(time,zeroline,'k','LineWidth',1.5)
                hold off
                set(gca,'FontSize',17)
                title('Loan-to-Value Ratio','FontSize',22),ylabel('pps','FontSize',18), xlabel('years since first shock','FontSize',18)
                subplot(2,3,2);
                lo95hp = 100*pctstrap_hp_d1ea9our(:,1); lo90hp = 100*pctstrap_hp_d1ea9our(:,2);
                up90hp = 100*pctstrap_hp_d1ea9our(:,4); up95hp = 100*pctstrap_hp_d1ea9our(:,5);
                x = [lo95hp lo90hp-lo95hp up90hp-lo90hp up95hp-up90hp];
                hold on
                harea = area(time,x,-12,'EdgeColor','none');
                set(harea(1),'Facecolor',[1 1 1],'FaceAlpha',0)
                set(harea(2),'Facecolor',[0.3 0.3 0.6],'FaceAlpha',0.5)
                set(harea(3),'Facecolor',[0 0 1],'FaceAlpha',0.5)
                set(harea(4),'Facecolor',[0.3 0.3 0.6],'FaceAlpha',0.5)
                plot(time,100*irfres_d1ea9our(1).resplev(:,1),'Color',[0 0 1],'Linewidth',3)
                grid on; ax=gca; ax.GridAlpha = 0.45;
                axis([0 T -10 4]); yticks([-9:3:3])         
                plot(time,zeroline,'k','LineWidth',1.5)
                hold off
                set(gca,'FontSize',17)
                title('Real house prices','FontSize',22),ylabel('%','FontSize',18), xlabel('years since first shock','FontSize',18)
                subplot(2,3,4)
                lo95cons = 100*pctstrap_cons_d1ea9our(:,1); lo90cons = 100*pctstrap_cons_d1ea9our(:,2);
                up90cons = 100*pctstrap_cons_d1ea9our(:,4); up95cons = 100*pctstrap_cons_d1ea9our(:,5);
                x = [lo95cons lo90cons-lo95cons up90cons-lo90cons up95cons-up90cons];
                hold on
                harea = area(time,x,-6,'EdgeColor','none');
                set(harea(1),'Facecolor',[1 1 1],'FaceAlpha',0)
                set(harea(2),'Facecolor',[0.3 0.3 0.6],'FaceAlpha',0.5)
                set(harea(3),'Facecolor',[0 0 1],'FaceAlpha',0.5)
                set(harea(4),'Facecolor',[0.3 0.3 0.6],'FaceAlpha',0.5)
                plot(time,100*irfres_d1ea9our(1).resplev(:,2),'Color',[0 0 1],'LineWidth',3)
                grid on; ax=gca; ax.GridAlpha = 0.45;
                axis([0 T -3.6 0.5]); yticks([-3:1:0.5])
                plot(time,zeroline,'k','LineWidth',1.5)
                hold off
                set(gca,'FontSize',17)
                title('Consumption','FontSize',22),ylabel('%','FontSize',18), xlabel('years since first shock','FontSize',18)
                subplot(2,3,5);
                lo95ires = 100*pctstrap_ires_d1ea9our(:,1); lo90ires = 100*pctstrap_ires_d1ea9our(:,2);
                up90ires = 100*pctstrap_ires_d1ea9our(:,4); up95ires = 100*pctstrap_ires_d1ea9our(:,5);
                x = [lo95ires lo90ires-lo95ires up90ires-lo90ires up95ires-up90ires];
                hold on
                harea = area(time,x,-12,'EdgeColor','none');
                set(harea(1),'Facecolor',[1 1 1],'FaceAlpha',0)
                set(harea(2),'Facecolor',[0.3 0.3 0.6],'FaceAlpha',0.5)
                set(harea(3),'Facecolor',[0 0 1],'FaceAlpha',0.5)
                set(harea(4),'Facecolor',[0.3 0.3 0.6],'FaceAlpha',0.5)
                plot(time,100*irfres_d1ea9our(1).resplev(:,3),'Color',[0 0 1],'Linewidth',3)
                grid on; ax=gca; ax.GridAlpha = 0.45;
                axis([0 T -12 10]); yticks([-10:5:10])
                plot(time,zeroline,'k','LineWidth',1.5)
                hold off
                set(gca,'FontSize',17)
                title('Residential investment','FontSize',22),ylabel('%','FontSize',18), xlabel('years since first shock','FontSize',18)
                subplot(2,3,6);
                lo95ibus = 100*pctstrap_ibus_d1ea9our(:,1); lo90ibus = 100*pctstrap_ibus_d1ea9our(:,2);
                up90ibus = 100*pctstrap_ibus_d1ea9our(:,4); up95ibus = 100*pctstrap_ibus_d1ea9our(:,5);
                x = [lo95ibus lo90ibus-lo95ibus up90ibus-lo90ibus up95ibus-up90ibus];
                hold on
                harea = area(time,x,-12,'EdgeColor','none');
                set(harea(1),'Facecolor',[1 1 1],'FaceAlpha',0)
                set(harea(2),'Facecolor',[0.3 0.3 0.6],'FaceAlpha',0.5)
                set(harea(3),'Facecolor',[0 0 1],'FaceAlpha',0.5)
                set(harea(4),'Facecolor',[0.3 0.3 0.6],'FaceAlpha',0.5)
                plot(time,100*irfres_d1ea9our(1).resplev(:,4),'Color',[0 0 1],'Linewidth',3)
                grid on; ax=gca; ax.GridAlpha = 0.45;
                axis([0 T -10 3]); yticks([-10:5:5])
                plot(time,zeroline,'k','LineWidth',1.5)
                hold off
                set(gca,'FontSize',17)
                title('Business investment','FontSize',22),ylabel('%','FontSize',18), xlabel('years since first shock','FontSize',18)
                subplot(2,3,3);
                lo95r = 100*pctstrap_r_d1ea9our(:,1); lo90r = 100*pctstrap_r_d1ea9our(:,2);
                up90r = 100*pctstrap_r_d1ea9our(:,4); up95r = 100*pctstrap_r_d1ea9our(:,5);
                x = [lo95r lo90r-lo95r up90r-lo90r up95r-up90r];
                hold on
                harea = area(time,x,-12,'EdgeColor','none');
                set(harea(1),'Facecolor',[1 1 1],'FaceAlpha',0)
                set(harea(2),'Facecolor',[0.3 0.3 0.6],'FaceAlpha',0.5)
                set(harea(3),'Facecolor',[0 0 1],'FaceAlpha',0.5)
                set(harea(4),'Facecolor',[0.3 0.3 0.6],'FaceAlpha',0.5)
                plot(time,100*irfres_d1ea9our(1).resplev(:,5),'Color',[0 0 1],'Linewidth',3)
                grid on; ax=gca; ax.GridAlpha = 0.45;
                axis([0 T -1 1]); yticks([-1:0.5:1])
                plot(time,zeroline,'k','LineWidth',1.5)
                hold off
                set(gca,'FontSize',17)
                title('Real Mortgage Rate','FontSize',22),ylabel('%','FontSize',18), xlabel('years since first shock','FontSize',18)
                legend([harea(2) harea(3)],'95% confidence interval','90% confidence interval','Location','SouthEast','FontSize',15)
                hold off
                     
            end
        end

        %Store impulses for d1ea so I can compare them across VECM specifications
        if robustness == 1 %If loop for all robustness checks
            if k == 4
                yres{4} = irfres_d1ea8our.resplev;  %5/12th of imposition in 2011, 7/12th in 2012.
            elseif k == 8
                yres{8} = irfres_d1ea.resplev;
            elseif k == 9
                yres{9} = irfres_d1ea.resplev;
            elseif k == 10
                yres{10} = irfres_d1ea.resplev;
            elseif k == 11
                yres{11} = irfres_d1ea.resplev;
            elseif k == 12
                yres{12} = irfres_d1ea.resplev;
            elseif k == 15
                yres{15} = irfres_d1ea8our.resplev;  
            elseif k == 16
                yres{16} = irfres_d1ea8our.resplev;   
            elseif k == 17
                yres{17} = irfres_d1ea8our.resplev;  
            end
        end
        
    end
        
    if shockspec == 5  & k == 4 %Influence of heterogeneity. Ten IRFs. %Only for baseline VECM specification.
        respLTVhet = nan(T+1,size(C,2)); respHPhet = nan(T+1,size(C,2));
        respGDPhet = nan(T+1,size(C,2)); respRhet = nan(T+1,size(C,2));
        time = [0:T]';
        respLTVhet(:,1) = time; respHPhet(:,1) = time;
        respGDPhet(:,1) = time; respRhet(:,1) = time;
        
        %Little bootstrap precursor: embryo of three-dimensional arrays
        %(this part will actually stay empty, only other sheets will be
        %filled)
        if bootstrap == 1
            strapirf_ltv(:,:) = nan(T+1,N); strapirf_hp(:,:) = nan(T+1,N);
            strapirf_gdp(:,:) = nan(T+1,N); strapirf_r(:,:) = nan(T+1,N);
            pctstrap_ltv(:,:) = nan(T+1,5); pctstrap_hp(:,:) = nan(T+1,5);
            pctstrap_gdp(:,:) = nan(T+1,5); pctstrap_r(:,:) = nan(T+1,5);
        end
        
        for i = 2:size(C,2)
            %Back to main (non-bootstrap) analysis
            impuls = zeros(T+1,1);
            impuls(1:size(C,1)) = C(:,i)./100; %fetch LTV impulse
            irfres = impulsovecm_manu(result,T,impuls,shockspec);
            respHPhet(:,i) = irfres(1).resplev(:,1);%store IRF path
            respGDPhet(:,i) = irfres(1).resplev(:,2);
            respRhet(:,i) = irfres(1).resplev(:,3);
            respLTVhet(:,i) = irfres(1).resplev(:,4);
            
            %Do bootstrapped confidence intervals for some of these 10
            %impulse paths
            if (i==2|i==6|i==7|i==11) & bootstrap == 1
                %untruncated 0.75
                %and 1.25 stdev and truncated 0.75 and
                %1.25 times stdev 
                %Define the rest of the arrays, to be filled
                strapirf_ltv(:,:,i) = nan(T+1,N); strapirf_hp(:,:,i) = nan(T+1,N);
                strapirf_gdp(:,:,i) = nan(T+1,N); strapirf_r(:,:,i) = nan(T+1,N);
                pctstrap_ltv(:,:,i) = nan(T+1,5); pctstrap_hp(:,:,i) = nan(T+1,5);
                pctstrap_gdp(:,:,i) = nan(T+1,5); pctstrap_r(:,:,i) = nan(T+1,5);
                for j = 1:N
                    strapresult = allstrapresult {j};
                    strapirfres=impulsovecm_manu(strapresult,T,impuls,shockspec);
                    
                    %Fill in IRFs. For each variable, jth column is
                    %jth simulation and ith sheet is ith impulse
                    strapirf_hp(:,j,i) = strapirfres.resplev(:,1);
                    strapirf_gdp(:,j,i) = strapirfres.resplev(:,2);
                    strapirf_r(:,j,i) = strapirfres.resplev(:,3);
                    strapirf_ltv(:,j,i) = strapirfres.resplev(:,4);
                end
                clear strapirfres %This is only for one of N specifications
                
                %IRF percentiles for ith shock specification
                pctstrap_ltv(:,:,i) = prctile(strapirf_ltv(:,:,i),[2.5 5 50 95 97.5],2);
                pctstrap_hp(:,:,i) = prctile(strapirf_hp(:,:,i),[2.5 5 50 95 97.5],2);
                pctstrap_gdp(:,:,i) = prctile(strapirf_gdp(:,:,i),[2.5 5 50 95 97.5],2);
                pctstrap_r(:,:,i) = prctile(strapirf_r(:,:,i),[2.5 5 50 95 97.5],2);
            end
            
        end         
        
       cd('JMCB_replic_deveirman-jong\VECM estimation');
   
    end
    %Back to initial directory
    cd('JMCB_replic_deveirman-jong\VECM estimation');
    if robustness == 1
        %close all %close figures when checking robustness
    end
    %end %end of shockspec loop (ONLY WHEN COMPARE IRFs PEARSON SYSTEM VS. EMPIRICAL DISTRIBUTIONS)
end

%close all  

if robustness == 1 & shockspec == 3
 
    %IRFs from k = 4 [loghp loggdp rh_ea1 ltv_comb_nz]
    yresd = yres{4};
    irfd_hp = yresd(:,1); irfd_gdp = yresd(:,2); irfd_r = yresd(:,3); irfd_ltv = yresd(:,4);

    %IRFs from k = 8 [loghp loggdp rh_ea1 ltv_comb_nz] with 2 lags
    yresh = yres{8};
    irfh_hp = yresh(:,1); irfh_gdp = yresh(:,2); irfh_r = yresh(:,3); irfh_ltv = yresh(:,4);

    %IRFs from k = 9 [loghp loggdp rh_ea1 ltv_comb_nz] with 3 lags
    yresi = yres{9};
    irfi_hp = yresi(:,1); irfi_gdp = yresi(:,2); irfi_r = yresi(:,3); irfi_ltv = yresi(:,4);

    %IRFs from k = 10 [loghp loggdp rh_ea1 ltv_comb_nz] with 4 lags
    yresj = yres{10};
    irfj_hp = yresj(:,1); irfj_gdp = yresj(:,2); irfj_r = yresj(:,3); irfj_ltv = yresj(:,4);

    %IRFs from k = 11 [loghp loggdp rh_ea1 ltv_comb_nz] with 2
    %cointegration relations
    yresk = yres{11};
    irfk_hp = yresk(:,1); irfk_gdp = yresk(:,2); irfk_r = yresk(:,3); irfk_ltv = yresk(:,4);

    %IRFs from k = 12 [loghp loggdp rh_ea1 ltv_comb_nz] with 3
    %cointegration relations
    yresl = yres{12};
    irfl_hp = yresl(:,1); irfl_gdp = yresl(:,2); irfl_r = yresl(:,3); irfl_ltv = yresl(:,4);

    %IRFs from k = 15. Alternative Cholesky ordering:
    %Version a: [loggdp ltv loghp rh_ea1]
    %Version b: [loggdp loghp rh_ea1 ltv_comb_nz]
    yreso = yres{15};
    %Version a
    %irfo_hp = yreso(:,3); irfo_gdp = yreso(:,1); irfo_r = yreso(:,4); irfo_ltv = yreso(:,2);
    %Version b
    irfo_hp = yreso(:,2); irfo_gdp = yreso(:,1); irfo_r = yreso(:,3); irfo_ltv = yreso(:,4);

    %IRFs from k = 16. Alternative Cholesky ordering:
    %Version a: [ltv loggdp loghp rh_ea1]
    %Version b: [loggdp rh_ea1 loghp ltv]
    yresp = yres{16};
    %Version a
    %irfp_hp = yresp(:,3); irfp_gdp = yresp(:,2); irfp_r = yresp(:,4); irfp_ltv = yresp(:,1);
    %Version b
    irfp_hp = yresp(:,3); irfp_gdp = yresp(:,1); irfp_r = yresp(:,2); irfp_ltv = yresp(:,4);

    %IRFs from k = 17. Alternative Cholesky ordering:
    %Version a: [ltv loggdp rh_ea1 loghp]
    %Version b: [loggdp rh_ea1 ltv loghp]
    yresq = yres{17};
    %Version a
    %irfq_hp = yresq(:,4); irfq_gdp = yresq(:,2); irfq_r = yresq(:,3); irfq_ltv = yresq(:,1);
    %Version b
    irfq_hp = yresq(:,4); irfq_gdp = yresq(:,1); irfq_r = yresq(:,2); irfq_ltv = yresq(:,3);


    %Combine k=4, k=15, k=16, k=17. This is the baseline with different Cholesky orders
    time=[0:T]';
    zeroline=zeros(size(time));
    figure; % We represent  100*ln(x) or 100*x, so it is in pct or pps

    subplot(2,2,4);
    h2=plot(time,100*irfd_hp,time,100*irfo_hp,time,100*irfp_hp,time,100*irfq_hp,time,zeroline,'k')%, title('Real house price'), ylabel('%'), xlabel('years since first shock')
    set(h2(1),'LineWidth',2); set(h2(2),'LineWidth',2); set(h2(3),'LineWidth',2); set(h2(4),'LineWidth',2); grid on; axis([0 T -5 0]); yticks([-5:1:0])
    set(gca,'FontSize',17)
    title('Real House Prices','FontSize',22),ylabel('%','FontSize',18), xlabel('years since first shock','FontSize',18)

    subplot(2,2,3);
    h3=plot(time,100*irfd_gdp,time,100*irfo_gdp,time,100*irfp_gdp,time,100*irfq_gdp,time,zeroline,'k')%, title({'','GDP'}), ylabel('%'), xlabel('years since first shock')
    set(h3(1),'LineWidth',2); set(h3(2),'LineWidth',2); set(h3(3),'LineWidth',2); set(h3(4),'LineWidth',2); grid on; axis([0 T -2 0]); yticks([-2:1:0])
    set(gca,'FontSize',17)
    title('Real Gross Domestic Product','FontSize',22),ylabel('%','FontSize',18), xlabel('years since first shock','FontSize',18)

    subplot(2,2,1);
    h1=plot(time,100*irfd_ltv,time,100*irfo_ltv,time,100*irfp_ltv,time,100*irfq_ltv,time,zeroline,'k')%,title('LTV'),ylabel('%ge points'), xlabel('years since first shock')
    set(h1(1),'LineWidth',2); set(h1(2),'LineWidth',2); set(h1(3),'LineWidth',2); set(h1(4),'LineWidth',2); grid on; axis([0 T -3 0]); yticks([-3:1:0])
    %legend('HP GDP r LTV','GDP HP r LTV','HP GDP LTV r', 'GDP HP LTV r')
    set(gca,'FontSize',17)
    title('Loan-to-Value Ratio','FontSize',22),ylabel('pps','FontSize',18), xlabel('years since first shock','FontSize',18)

    subplot(2,2,2);
    h4=plot(time,100*irfd_r,time,100*irfo_r,time,100*irfp_r,time,100*irfq_r,time,zeroline,'k')%, title('Real mortgage rate'), ylabel('%ge points'), xlabel('years since first shock')
    set(h4(1),'LineWidth',2); set(h4(2),'LineWidth',2);set(h4(3),'LineWidth',2); set(h4(4),'LineWidth',2); grid on; axis([0 T -1 1]); yticks([-1:0.5:1])
    set(gca,'FontSize',17)
    title('Real Mortgage Rate','FontSize',22),ylabel('pps','FontSize',18), xlabel('years since first shock','FontSize',18), legend('baseline ordering','alternative 1', 'alternative 2', 'alternative 3','Location','NorthEast', 'FontSize',16,'NumColumns',2)
    % 'HP GDP r LTV','GDP HP r LTV','HP GDP LTV r', 'GDP HP LTV r'
    %'baseline ordering','alternative 1', 'alternative 2', 'alternative 3'


    %Combine k=4, k=8, k=9. This is the baseline with 1-3 LAGS
    time=[0:T]';
    zeroline=zeros(size(time));
    figure; % We represent  100*ln(x) or 100*x, so it is in pct or pps

    subplot(2,2,4);
    h2=plot(time,100*irfd_hp,time,100*irfh_hp,time,100*irfi_hp,time,zeroline,'k')%, title('Real house price'), ylabel('%'), xlabel('years since first shock')
    set(h2(1),'LineWidth',2); set(h2(2),'LineWidth',2); set(h2(3),'LineWidth',2); grid on; axis([0 T -5 0]); yticks([-5:1:0])
    set(gca,'FontSize',17)
    title('Real House Prices','FontSize',22),ylabel('%','FontSize',18), xlabel('years since first shock','FontSize',18)

    subplot(2,2,3);
    h3=plot(time,100*irfd_gdp,time,100*irfh_gdp,time,100*irfi_gdp,time,zeroline,'k')%, title({'','GDP'}), ylabel('%'), xlabel('years since first shock')
    set(h3(1),'LineWidth',2); set(h3(2),'LineWidth',2); set(h3(3),'LineWidth',2); grid on; axis([0 T -2 0]); yticks([-2:1:0])
    set(gca,'FontSize',17)
    title('Real Gross Domestic Product','FontSize',22),ylabel('%','FontSize',18), xlabel('years since first shock','FontSize',18)

    subplot(2,2,1);
    h1=plot(time,100*irfd_ltv,time,100*irfh_ltv,time,100*irfi_ltv,time,zeroline,'k')%,title('LTV'),ylabel('%ge points'), xlabel('years since first shock')
    set(h1(1),'LineWidth',2); set(h1(2),'LineWidth',2); set(h1(3),'LineWidth',2); grid on; axis([0 T -4 0]); yticks([-4:1:0])
    set(gca,'FontSize',17)
    title('Loan-to-Value Ratio','FontSize',22),ylabel('pps','FontSize',18), xlabel('years since first shock','FontSize',18)

    subplot(2,2,2);
    h4=plot(time,100*irfd_r,time,100*irfh_r,time,100*irfi_r,time,zeroline,'k')%, title('Real mortgage rate'), ylabel('%ge points'), xlabel('years since first shock')
    set(h4(1),'LineWidth',2); set(h4(2),'LineWidth',2);set(h4(3),'LineWidth',2); grid on; axis([0 T -1 1]); yticks([-1:0.5:1])
    set(gca,'FontSize',17)
    title('Real Mortgage Rate','FontSize',22),ylabel('pps','FontSize',18), xlabel('years since first shock','FontSize',18), legend('1 lag','2 lags','3 lags','Location','NorthEast', 'FontSize',16,'NumColumns',2)

end


toc