%This function is to replicate results for directional approach from Section 4.1 and Section 4.2
clear;clc;close all

addpath('data');
load('USdata');%X = [intercept, lagged h prices, lagged hh credit, detrended credit ratio,
               %    nfci, gdp gr, real mortgage rate, detrended hp/GDPper capita, detrended hhcredit/gdp]
               %Y = [h prices, hh credit]

%% Choose specification               



%Specifications from Section 4.1 and 4.2
%X = X(:,[1,2,3,4]);    %specifications with intercept, lagged prices and lagged credit, and credit ratio
%X = X(:,[1,2,3,5]);    %specifications with intercept, lagged prices and lagged credit, and nfci
%X = X(:,[1,2,3,10]); X = X(54:end,:); Y = Y(54:end,:); %specifications with intercept, lagged prices and lagged credit, and vix
%X = X(:,[1,2,3,8]);    %specifications with intercept, lagged prices and lagged credit, and price missalignment
%X = X(:,[1,2,3,9]);    %specifications with intercept, lagged prices and lagged credit, and hh leverage
%X = X(:,[1,2,3,6]);    %specifications with intercept, lagged prices and lagged credit, and gdp gr
%X = X(:,[1,2,3,7]);    %specifications with intercept, lagged prices and lagged credit, and real mortgage rate
%X = X(:,[1,2,3,11]); X = X(1:124,:); Y = Y(1:124,:); %specifications with intercept, lagged prices and lagged credit, and romer and romer mp shocks
%X = X(:,[1,2,3,12]); X = X(1:112,:); Y = Y(1:112,:); %specifications with intercept, lagged prices and lagged credit, and smets and wouters mp shocks

%Specification from Section 4.3
X = X(:,[1,2,3,4,6,8,9]); %credit ratio, gdp gr, hh leverage, hp misalignment

% Specification from Section 4.4
%X = X(:,[1,2,3,6,8]); %gdp growth, hp misalignment


%% %multicollinerity check
VIF = vif(X(:,2:end));   %vif for a predictor exceeding 10 is a problem
if any(VIF>10)
    disp('Multicollinearity issue');
end

%% Benchmark estimation
nboot=5000;
stdY = std(Y);
Yp = (1/sqrt(2))*Y(:,1) + (1/sqrt(2))*Y(:,2);

u_sqr = 0.05:0.05:0.95; 
for i=1:length(u_sqr)
    res_qr11 = rq(X,Yp,u_sqr(i)); 
    Beta(:,i) = res_qr11';
    Beta_mbb  = mbb(Yp,X,nboot,u_sqr(i)); 
    Beta_bsample(:,:,i) = Beta_mbb; 
    Beta90(:,:,i) = [quantile(Beta_mbb,0.05);quantile(Beta_mbb,0.95)];
    Beta80(:,:,i) = [quantile(Beta_mbb,0.10);quantile(Beta_mbb,0.90)];
    Beta68(:,:,i) = [quantile(Beta_mbb,0.16);quantile(Beta_mbb,0.84)];
    res_qr_cr = rq(X,Y(:,2),u_sqr(i)); 
    res_qr_hp = rq(X,Y(:,1),u_sqr(i)); 
    Beta_cr(:,i) = res_qr_cr';
    Beta_cr_mbb = mbb(Y(:,2),X,nboot,u_sqr(i)); 
    Beta_cr_bsample(:,:,i) = Beta_cr_mbb;
    Beta_hp(:,i) = res_qr_hp';
    Beta_hp_mbb = mbb(Y(:,1),X,nboot,u_sqr(i)); 
    Beta_hp_bsample(:,:,i) = Beta_hp_mbb;
end

%ols
res_ols = ols(Yp,X);

%% in-sample fit (quantile R2) - Koenker and Machado (1999), qBIC

q = 0.05;
ui=find(u_sqr(1,:)==q);
Yfit = X*Beta(:,ui);
SARq   = sum(check_fn(Yp - Yfit,q));
SARunc = sum(check_fn(Yp - quantile(Yp,q),q));
qR2 = 1 - (SARq/SARunc);
disp(' ')
fprintf([num2str(q) ': qR2: %.2f\n'], qR2)

n = size(Yp,1);k = size(X,2);
BICq = k*log(n) + n*log((1/n)*SARq);
fprintf([num2str(q) ': qBIC (no crossing imposed): %.3f\n'], BICq)


q = 0.95;
ui=find(u_sqr(1,:)==q);
Yfit = X*Beta(:,ui);
SARq   = sum(check_fn(Yp - Yfit,q));
SARunc = sum(check_fn(Yp - quantile(Yp,q),q));
qR2 = 1 - (SARq/SARunc);
disp(' ')
fprintf([num2str(q) ': qR2: %.2f\n'], qR2)

n = size(Yp,1);k = size(X,2);
BICq = k*log(n) + n*log((1/n)*SARq);
fprintf([num2str(q) ': qBIC (no crossing imposed): %.3f\n'], BICq)



%% OUTPUT

%% Table 1a, 1b, 2 and 4 - one column for each specification
disp('QR coeffs:')
Beta(:,find(u_sqr(1,:)==0.05))
%% Table D1 and D2 - one column for each specification
disp('QR coeffs:')
Beta(:,find(u_sqr(1,:)==0.95))

%% Figure 3: index, estimated quantile (with confidence band) and bank sector crises 
% (for speficication with credit ratio, hh leverage, hp misalignment and real output growth)
q = 0.05;
ui=find(u_sqr(1,:)==q);
Yfit = X*Beta(:,ui);
for i=1:nboot
    Yfit_bsample(:,i) = X*Beta_bsample(i,:,ui)';
end
Yfit90 = [quantile(Yfit_bsample',0.05);quantile(Yfit_bsample',0.95)];
ind_ = find(Yp(:,1)<Yfit90(1,:)');
length(ind_)/size(Yp,1)

h=1:size(Y,1);
figure
set(gcf, 'Color', 'w');
hold on
plot(Yfit90(1,:),'--','Color',[0.88 0.88 0.88]);
plot(Yfit90(2,:),'--','Color',[0.88 0.88 0.88]);
x2 = [h, fliplr(h)];
inBetween = [Yfit90(1,:),fliplr(Yfit90(2,:))];
a1 = fill(x2,inBetween,[0.88 0.88 0.88],'LineStyle','none');
a2 = plot(Yp(:,1),'k-','Linewidth',1.5);
a3 = plot(Yfit(:,1),'k--','Linewidth',1);
ylimits = ylim;
plot([45,45],[ylimits(1),ylimits(2)],'r-','LineWidth',1.5);
plot([121,121],[ylimits(1),ylimits(2)],'r-','LineWidth',1.5);
plot(Yp(:,1),'k-','Linewidth',1.5);
plot(Yfit(:,1),'k--','Linewidth',1);
ylim([ylimits(1) ylimits(2)]);
legend([a2,a3],'Index of house price growth and credit growth',[num2str(q) '-th quantile of conditional distribution of index']);
xlim([1 size(Y,1)]);
set(gca,'XTick',1:28:size(Y,1))
set(gca,'XTickLabel',dates(1:28:size(Y,1))) 

%% Figure 4: contributions to one year ahead index change (for speficication with credit ratio, hh leverage, hp misalignment and real output growth)
q = 0.05;
ui=find(u_sqr(1,:)==q);
Ybar = X.*Beta(:,ui)';
Yfit = X*Beta(:,ui);

figure
set(gcf, 'Color', 'w');
hold on
bar(Ybar(:,2:end),'stacked')
plot(Yfit(:,1)-Beta(1,ui),'k-','Linewidth',1.5);
xlim([1 size(Y,1)]);
set(gca,'XTick',1:20:size(Y,1))
set(gca,'XTickLabel',dates(1:20:size(Y,1))) 
legend('House price growth','Household credit growth','Credit ratio','GDP growth','House price misalignment','Household leverage')



%% Figure C2: quantile function for credit ratio (for specification with credit ratio)
figure
set(gcf, 'Color', 'w');
hold on
a1 = plot(Beta(4,:),'Linewidth',1.5);
a2 = plot(squeeze(Beta90(1,4,:)),'k--');
plot(squeeze(Beta90(2,4,:)),'k--');
a3 = plot(squeeze(Beta68(1,4,:)),'k-.');
plot(squeeze(Beta68(2,4,:)),'k-.');
a4 = plot(ones(1,19)*res_ols.beta(4,1),'r--','LineWidth',1.5);
xlim([1,length(u_sqr)]);
set(gca,'XTick',1:3:length(u_sqr))
set(gca,'XTickLabel',u_sqr(1:3:length(u_sqr))) 
xlabel('Quantile');
legend([a1,a2,a3,a4],'Quantile function','90% confidence band','68% confidence band','OLS');

%% Figure D1: two fitted geometric quantiles (robustness check)
q1 = 0.05;
ui=find(u_sqr(1,:)==q1);
Yfit1 = X*Beta(:,ui);
q2 = 0.10;
ui=find(u_sqr(1,:)==q2);
Yfit2 = X*Beta(:,ui);
 
h=1:size(Y,1);
figure
set(gcf, 'Color', 'w');
hold on
plot(Yfit1(:,1),'r-');
plot(Yfit2(:,1),'k--');
ylim([-17,12]);
xlim([1 size(Y,1)]);
set(gca,'XTick',1:28:size(Y,1))
set(gca,'XTickLabel',dates(1:28:size(Y,1))) 
yyaxis right
plot(Yfit2-Yfit1,'-.');
plot(zeros(1,size(Yfit1,1)));
legend('0.05-th conditional quantile','0.10-th conditional quantile','difference (right axis)');

%% Figure C4: 5th percentile of index, 5th perc of cred gr and 5th perc of h prices
%(for speficication with credit ratio, hh leverage, hp misalignment and real output growth)
q = 0.05;
ui=find(u_sqr(1,:)==q);
Yfit = X*Beta(:,ui);
ind = find(Yp(:,1)<Yfit(:,1));
Yfit_cr = X*Beta_cr(:,ui);
ind_cr = find(Y(:,2)<Yfit_cr(:,1));
Yfit_hp = X*Beta_hp(:,ui);
ind_hp = find(Y(:,1)<Yfit_hp(:,1));

h=1:size(Y,1);
figure
set(gcf, 'Color', 'w');
 
subplot(1,3,1)
 hold on
 plot(Yp(:,1),'k-','Linewidth',1.5);
 plot(Yfit(:,1),'k--','Linewidth',1);
 ylimits = ylim;
 for i=1:size(Yp,1)
    if ismember(i,ind)
        h = area([i-0.5,i+0.5],[ylimits(2)+0.01,ylimits(2)+0.01],'basevalue',ylimits(1),'FaceColor',[0.85 0.85 0.85],'EdgeColor',[0.85 0.85 0.85]);        
        dates{i,1}
    end
 end
 plot(Yp(:,1),'k-','Linewidth',1.5);
 plot(Yfit(:,1),'k--','Linewidth',1);
 ylim([ylimits(1) ylimits(2)]);
 xlim([1 size(Y,1)]);
 set(gca,'XTick',1:28:size(Y,1))
 set(gca,'XTickLabel',dates(1:28:size(Y,1))) 
 title('Index');
 legend('Observed data',[num2str(q) '-th conditional quantile'],'Data below quantile');
subplot(1,3,2)
 hold on
 plot(Y(:,1),'k-','Linewidth',1.5);
 plot(Yfit_hp(:,1),'k--','Linewidth',1);
 ylimits = ylim;
 for i=1:size(Yp,1)
    if ismember(i,ind_hp)
        h = area([i-0.5,i+0.5],[ylimits(2)+0.01,ylimits(2)+0.01],'basevalue',ylimits(1),'FaceColor',[0.85 0.85 0.85],'EdgeColor',[0.85 0.85 0.85]);        
        dates{i,1}
    end
 end
 plot(Y(:,1),'k-','Linewidth',1.5);
 plot(Yfit_hp(:,1),'k--','Linewidth',1);
 ylim([ylimits(1) ylimits(2)]);
 xlim([1 size(Y,1)]);
 set(gca,'XTick',1:28:size(Y,1))
 set(gca,'XTickLabel',dates(1:28:size(Y,1))) 
 title('House prices');
subplot(1,3,3)
 hold on
 plot(Y(:,2),'k-','Linewidth',1.5);
 plot(Yfit_cr(:,1),'k--','Linewidth',1);
 ylimits = ylim;
 for i=1:size(Yp,1)
    if ismember(i,ind_cr)
        h = area([i-0.5,i+0.5],[ylimits(2)+0.01,ylimits(2)+0.01],'basevalue',ylimits(1),'FaceColor',[0.85 0.85 0.85],'EdgeColor',[0.85 0.85 0.85]);        
        dates{i,1}
    end
 end
 plot(Y(:,2),'k-','Linewidth',1.5);
 plot(Yfit_cr(:,1),'k--','Linewidth',1);
 ylim([ylimits(1) ylimits(2)]);
 xlim([1 size(Y,1)]);
 set(gca,'XTick',1:28:size(Y,1))
 set(gca,'XTickLabel',dates(1:28:size(Y,1))) 
 title('Household credit');

 disp('both credit and prices')
 for i=1:size(Yp,1)
    if ismember(i,ind_hp) && ismember(i,ind_cr)
        dates{i,1}
    end
 end
 


%% (Pseudo) out-of-sample validation
% {
T = size(Yp,1);
qx = 0.05;
T_start = 100;%100;
Yfitx = NaN(T-T_start,1);
R2_nom = 0;
R2_denom = 0;
for iT = T_start:(T-1)
    iT
    Yx = Y(1:iT,:);
    Xx = X(1:iT,:);
    %Ypx = (1/stdY(1,1))*Yx(:,1) + (1/stdY(1,2))*Yx(:,2);
    Ypx = (1/(1/sqrt(2)))*Yx(:,1) + (1/(1/sqrt(2)))*Yx(:,2);

    res_qr11x = rq(Xx,Ypx,qx); 
    Betax(:,iT-T_start+1) = res_qr11x;
    Yfitx(iT+1-T_start,1) = X(iT+1,:)*Betax(:,iT-T_start+1);
    
    Betax_boot = mbb(Ypx,Xx,1000,qx);
    Yfitx_boot = X(iT+1,:)*Betax_boot';
    Yfit90(iT+1-T_start,:) = [quantile(Yfitx_boot,0.05);quantile(Yfitx_boot,0.95)];
    res_qr11x_ = quantile(Ypx,qx);
    resids(iT+1-T_start,1) = Yp(iT+1,1)-Yfitx(iT+1-T_start,1);
    R2_nom = R2_nom + check_fn(Yp(iT+1,1)-Yfitx(iT+1-T_start,1),qx);
    R2_denom = R2_denom + check_fn(Yp(iT+1,1)-res_qr11x_,qx);
end
R2 = 1 - (R2_nom)/(R2_denom);
%}

%% Figure 5: pseudo out of sample forecasting performance
load('Yfit_insample')   %fitted quantile from estimation based on whole sample
ui=find(u_sqr(1,:)==qx);
h=1:size(Y,1);
hx=(T_start+1):size(Y,1);

figure
set(gcf, 'Color', 'w');
hold on
plot(h,Yfit_insample(:,1),'b-.','Linewidth',1.5);
plot(hx,Yfitx(:,1),'r--','Linewidth',1.5);
plot(h,Yp,'k-','Linewidth',1.5);

plot(hx,Yfit90(:,1),'--','Color',[255/255 183/255 183/255]);
plot(hx,Yfit90(:,2),'--','Color',[255/255 183/255 183/255]);
x2 = [hx, fliplr(hx)];
inBetween = [Yfit90(:,1)',fliplr(Yfit90(:,2)')];
a4 = fill(x2,inBetween,[255/255 183/255 183/255],'LineStyle','none');

a2 = plot(h,Yfit_insample(:,1),'b-.','Linewidth',1.5);
a3 = plot(hx,Yfitx(:,1),'r--','Linewidth',1.5);
a1 = plot(h,Yp,'k-','Linewidth',1.5);

ylim([-60,35]);

legend([a1,a2,a3,a4],'Index of house price growth and credit growth','0.05-th quantile of conditional distribution of index: in-sample','0.05-th quantile of conditional distribution of index: out-of-sample','90% confidence band');
xlim([1 size(Y,1)]);
set(gca,'XTick',1:28:size(Y,1))
set(gca,'XTickLabel',dates(1:28:size(Y,1))) 


%% Figure B1: Data
load('USdata');
tit = {'House price growth','Household credit growth','Credit ratio','NFCI','GDP growth','Real mortgage rate','House price misalignment','Household leverage'};
 
figure
set(gcf, 'Color', 'w');
for i=2:9
   subplot(2,4,i-1)
   hold on
   plot(X(:,i));
   title(tit{1,i-1})
   xlim([1 size(X,1)]);
   set(gca,'XTick',1:38:size(X,1))
   set(gca,'XTickLabel',dates(1:38:size(X,1))) 
end