%This script is to replicate results for direct approach from Section 4.1 and Section 4.2
%It estimates multivariate quantile linear regression following Chakraborty (2003). See Appendix A for the estimation procedure.
%Note that notation in the code below closely follows Chakraborty (2003)
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
%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

%% Set parameters
TR = 1;           %indicator of transformation retransformation procedure (1.. TR switched on)

%% Preliminaries
n = size(X,1);
k = size(X,2);
ij_last=[1,2,3,4,5,6];t_old = 1000;

%% Step 1: Finding Sigma_hat (unbiased estimate)
res_ols11=ols(Y(:,1),X);
res_ols11.bint90 = [res_ols11.beta-1.645*res_ols11.bstd,res_ols11.beta+1.645*res_ols11.bstd];
res_ols12=ols(Y(:,2),X);
res_ols12.bint90 = [res_ols12.beta-1.645*res_ols12.bstd,res_ols12.beta+1.645*res_ols12.bstd];
Sigma_hat = cov(res_ols11.resid,res_ols12.resid);

if TR
    
%% Steps 2-4: Computing transformation matrix E(alpha), V(alpha) and t(alpha)
%             and minimizing t(alpha) wrt alpha
Tol_t = 1+1e-9;             %tolerance when looking for transformation matrix (Step 2-4)

disp('Finding transformation matrix E(alpha)...');
ij = 0;flag = 0;
hh = waitbar(0,'Minimization t(alpha)...');
tic;
for i1=1:n-5
    if ij_last(1,1)>i1
        continue
    end
    for i2=(i1+1):n-4
        if ij_last(1,2)>i2
            continue
        end
        for i3=(i2+1):n-3
            if ij_last(1,3)>i3
                continue
            end
            for i4=(i3+1):n-2         
                for j1=(i4+1):n-1
                    for j2=(j1+1):n
                        alpha = [i1,i2,i3,i4,j1,j2];
                                            
                        W = [X(i1,:)',X(i2,:)',X(i3,:)',X(i4,:)'];          %k x k
                        Z = [Y(i1,:)',Y(i2,:)',Y(i3,:)',Y(i4,:)'];          %2 x k
                        E = [Y(j1,:)'-Z*(W\X(j1,:)'),...
                                Y(j2,:)'-Z*(W\X(j2,:)')];
                        V = E'*(Sigma_hat\E);
                        if det(V)<0
                            t = t_old;  %just to filter out cases when det(V)<0
                            disp('negative det');
                        else
                            t = (trace(V)/2)/(det(V)^(1/2)); %note that t(alpha)>=1
                        end
                        if t<t_old
                            alpha_hat = alpha
                            E_hat = E;
                            t_old = t;
                            if t<Tol_t
                                flag = 1;
                            break
                            end
                        end
                        ij=ij+1;
                    
                        if flag
                            break
                        end
                    end
                    if flag
                        break
                    end
                end
                if flag
                   break
                end           
            end
            if flag
                break
            end
        end
        if flag
            break
        end
    end
    if flag
        break
    end
waitbar(i1/(n-(k+1)),hh);    
end
estimation_time_alpha_hat=toc/(60*60);
fprintf('t(alpha) minimization (hours): %.2f\n',estimation_time_alpha_hat);
close(hh);
ij_last = [i1,i2,i3,i4,j1,j2];
if i1 == (n-5)
   disp('all combinations checked for minimum'); 
   all_combinations = 1;
   t_last = t_old;
else
   t_last = t; 
end

else %TR=0
   E_hat = eye(2); 
   alpha_hat = [];
end %end of TR=0,1

%% Step 5: Compute z^(alpha_hat)_i...z,v(alpha_hat)...v  
for i=1:n 
    if ~ismember(i,alpha_hat)
        z(i,:)=(E_hat\Y(i,:)')';
    else
        z(i,:) = [NaN,NaN];
    end
end

%prepare set of all 2-quantiles
k0 = 1;clear U;
for k1 = -5:5
    for k2 = -5:5
        U(:,k0) = [k1/10;k2/10];
        k0=k0+1;
    end
end

Newton1 = NaN(1,size(U,2));
clear Betas Beta1 BetaSQR U_sqr
for iu = 1:size(U,2)
    u = U(:,iu);
    if sum(u==0) == 2
        v = [0;0];
    else
        v = ((E_hat\u)*norm(u))/norm(E_hat\u);
    end

    %% Step 6: Degeneracy condition
    flag = 0;
    %{
    disp('Checking non-degeneracy condition...');
    hh = waitbar(0,'Checking non-degeneracy condition...');
    for i1=1:(n-3)
        for i2=(i1+1):(n-2)
            for i3=(i2+1):(n-1)
                for i4=(i3+1):n
                    alpha_dc = [i1,i2,i3,i4];
                    if ~any(ismember(alpha_dc,alpha_hat))
                        Wdc = [X(i1,:)',X(i2,:)',X(i3,:)',X(i4,:)'];   
                        Zdc = [z(i1,:)',z(i2,:)',z(i3,:)',z(i4,:)']; 
                
                        DCleft = 0;
                        for j=1:n
                            if ~ismember(j,[alpha_hat,i1,i2,i3,i4])   
                                DCleft = DCleft + ...
                                         kron(((z(j,:)' - Zdc*(Wdc\X(j,:)'))/norm(z(j,:)' - Zdc*(Wdc\X(j,:)')) + v),X(j,:)');
                            end      
                        end
                        DCleft = norm(DCleft);
                        DCright = 0;
                        for j=alpha_dc
                            DCright = DCright + norm(X(j,:)');
                        end
                        DCright = (1+norm(v))*DCright;
                
                        %Degeneracy condition
                        if  DCleft<=DCright
                            disp('Degeneracy Condition satisfied.');
                            flag = 1;
                            break
                        end
                    end
                end
                if flag
                    break
                end
            end
            if flag
                break
            end
        end
        if flag
            break
        end
        waitbar(i1/(n-(k+1)),hh);
    end
    close(hh);
    %}

    %% Step 7: Iteration to get Gamma(u) - Newton's method
    % the method converges quandratically - the number of correct digits
    % doubles each iteration => not necessarry to employ too many tries

    %various init values of Gamma
    indices = find(~isnan(z(:,1)));
    res_ols1 = ols(z(indices,1),X(indices,:));
    res_ols2 = ols(z(indices,2),X(indices,:));
    res_qr1 = rq(X(indices,:),z(indices,1),(v(1,1)+1)/2); %single output QR
    res_qr2 = rq(X(indices,:),z(indices,2),(v(2,1)+1)/2); %single output QR
    
    %usual Newton method can be used:
    %[Gamma1,Newton1(1,iu)] = newton([res_qr1';res_qr2'],1e-13,alpha_hat,z,X,v);

    %or fminsearch
    % G is 4x2 
    fun = @(G)sum(    sqrt(sum((z-X*G).^2,2)) +   (z-X*G)*v, 'omitnan');
    options = optimset('MaxIter',50000,'MaxFunEvals',50000);
    [Gamma4,~,Der_free(1,iu)] = fminsearch(fun,[res_qr1';res_qr2']',options);


    %% Results
    if flag == 1
        Betas(:,:,iu) = E_hat*(Zdc/Wdc);
    else
        %Beta1(:,:,iu) = E_hat*Gamma1;
        Beta4(:,:,iu) = E_hat*Gamma4';
        if Der_free(1,iu) == 1
            Betas(:,:,iu) = Beta4(:,:,iu);
        else
            Betas(:,:,iu) = NaN(2,k);
        end
    end   
end

%% confidence bands - moving block bootstrap
nboot = 5000;                                       %number of bootstrap samples, 0 means no bootstrap
qboot_set = [-0.5,-0.5;-0.4,-0.4;0,0;0.5,0.5]';     %geometric quantiles for which confidence bands are computed
%qboot_set = [-0.5,-0.5]';
%qboot_set = [0.5,0.5]';
if nboot>0
    clear Betas_boot Beta90_left Beta90_right Beta80_left Beta80_right ...
          Beta68_left Beta68_right Gamma1_boot Gamma2_boot Gamma3_boot
    Newton1boot = NaN(size(qboot_set,2),nboot);
    Der_free_boot = NaN(size(qboot_set,2),nboot);

    b = floor(n^(1/3));
    m = floor(n/b);
    q = n-b+1;
    
    for iq = 1:size(qboot_set,2)
        Betas_boot = NaN(nboot,2,k);
        qboot = qboot_set(:,iq);
        
        if sum(qboot==0) == 2
            vboot = [0;0];
        else
            vboot = ((E_hat\qboot)*norm(qboot))/norm(E_hat\qboot);
        end
    
        zi = NaN(q,b,2);
        Xi = NaN(q,b,k);
        for i=1:q   
            zi(i,:,:) = z(i:i+b-1,:); 
            Xi(i,:,:) = X(i:i+b-1,:);
        end
        
        hh = waitbar(0,['Bootstrapping:',num2str(iq),'/',num2str(size(qboot_set,2)),'...']);
        i = 1;
        while i<=nboot
            z_boot=[];X_boot=[];i_boot=[];
            for im=1:m       
                i_boot = randi([1 q]);
                z_boot = [z_boot;squeeze(zi(i_boot,:,:))];
                X_boot = [X_boot;squeeze(Xi(i_boot,:,:))];
            end
                    
            indices = find(~isnan(z_boot(:,1)));
            indices_nan = find(isnan(z_boot(:,1)));
            res_qr1 = rq(X_boot(indices,:),z_boot(indices,1),(vboot(1,1)+1)/2); %single output QR
            res_qr2 = rq(X_boot(indices,:),z_boot(indices,2),(vboot(2,1)+1)/2); %single output QR            

            %[Gamma1_boot,Newton1boot(iq,i)] = newton([res_qr1';res_qr2'],Tol_Gamma,indices_nan',z_boot,X_boot,vboot);
            
            fun = @(G)sum(    sqrt(sum((z_boot-X_boot*G).^2,2)) +   (z_boot-X_boot*G)*vboot, 'omitnan');
            options = optimset('MaxIter',50000,'MaxFunEvals',50000);
            [Gamma4_boot,~,Der_free_boot(iq,i)] = fminsearch(fun,[res_qr1';res_qr2']',options);

            if Der_free_boot(iq,i) == 1
                Betas_boot(i,:,:) = E_hat*Gamma4_boot';
                i = i + 1;
            end
        waitbar(i/nboot,hh);
        end    
    
        Betas_boot(isnan(Betas_boot(:,1,1)),:,:)=[];
        Beta90_left(:,:,iq) = [quantile(squeeze(Betas_boot(:,1,:)),0.05);quantile(squeeze(Betas_boot(:,2,:)),0.05)];
        Beta90_right(:,:,iq) = [quantile(squeeze(Betas_boot(:,1,:)),0.95);quantile(squeeze(Betas_boot(:,2,:)),0.95)];
        Beta80_left(:,:,iq) = [quantile(squeeze(Betas_boot(:,1,:)),0.10);quantile(squeeze(Betas_boot(:,2,:)),0.10)];
        Beta80_right(:,:,iq) = [quantile(squeeze(Betas_boot(:,1,:)),0.90);quantile(squeeze(Betas_boot(:,2,:)),0.90)];       
        Betas_boot_all(:,:,:,iq) = Betas_boot;
        close(hh);
    end
end

%% single QR and OLS
u_sqr = 0.05:0.05:0.95; 
for i=1:length(u_sqr)
    res_qr11 = rq(X,Y(:,1),u_sqr(i)); %single output QR
    res_qr12 = rq(X,Y(:,2),u_sqr(i)); %single output QR
    BetaSQR(:,:,i) = [res_qr11';res_qr12'];

    %confidence bands for marginal quantile regressions
    Beta_mbb  = mbb(Y(:,1),X,nboot,u_sqr(i)); 
    BetaSQR1_90(:,:,i) = [quantile(Beta_mbb,0.05);quantile(Beta_mbb,0.95)];
    BetaSQR1_68(:,:,i) = [quantile(Beta_mbb,0.16);quantile(Beta_mbb,0.84)];

    Beta_mbb  = mbb(Y(:,2),X,nboot,u_sqr(i)); 
    BetaSQR2_90(:,:,i) = [quantile(Beta_mbb,0.05);quantile(Beta_mbb,0.95)];
    BetaSQR2_68(:,:,i) = [quantile(Beta_mbb,0.16);quantile(Beta_mbb,0.84)];
end

[BetaOLS1,BetaOLS1_90] = regress(Y(:,1),X,0.9);  
[BetaOLS1,BetaOLS1_80] = regress(Y(:,1),X,0.8); 
[BetaOLS2,BetaOLS2_90] = regress(Y(:,2),X,0.9);  
[BetaOLS2,BetaOLS2_80] = regress(Y(:,2),X,0.8); 

%% OUTPUT

%analysis of convergence of minimization procedure
fprintf('ratio of converging fminsearch %.2f\n',sum(Der_free>0)/size(Der_free,2));
fprintf('ratio of converging fminsearch in bootstrap %.2f\n',sum(Der_free>0)/size(Der_free,2));


%% Table 1a, 1b, 2 and 4 - one column for each specification
disp('QR coeffs:')
Betas(:,:,intersect(find(U(1,:)==-0.5),find(U(2,:)==-0.5)))'
%% Table D1 and D2 - one column for each specification
disp('QR coeffs:')
Betas(:,:,intersect(find(U(1,:)==0.5),find(U(2,:)==0.5)))'
%% Table D3 and D4 - one column for each specification
disp('QR coeffs:')
Betas(:,:,intersect(find(U(1,:)==-0.4),find(U(2,:)==-0.4)))'

%% Figure B2: scatter plot of data
figure
set(gcf, 'Color', 'w');
hold on
scatter(Y(:,1),Y(:,2));
ax = gca;
ax.XAxisLocation = 'origin';
ax.YAxisLocation = 'origin';
xlabel('House price growth (%)');
ylabel('Household credit growth (%)');

%% Figure C1: 3D plot of coeffs at NFCI for various index vectors (for specification with NFCI)
axe_length = 11;
qq  = [-0.5:0.1:0.5];
figure
set(gcf, 'Color', 'w');
subplot(1,2,1)
    mesh(qq,qq,reshape(squeeze(Betas(1,4,:)),axe_length,axe_length));
    xlabel('u_{1}');
    xlim([qq(1,1),qq(1,end)]);
    xticks([qq(1,1):0.2:qq(1,end)]);
    ylabel('u_{2}');
    ylim([qq(1,1),qq(1,end)]);
    yticks([qq(1,1):0.2:qq(1,end)]);
    zlabel('Coefficient at NFCI');
    title('Equation for house prices');
subplot(1,2,2)
    mesh(qq,qq,reshape(squeeze(Betas(2,4,:)),axe_length,axe_length));
    xlabel('u_{1}');    
    xlim([qq(1,1),qq(1,end)]);
    xticks([qq(1,1):0.2:qq(1,end)]);
    ylabel('u_{2}');
    ylim([qq(1,1),qq(1,end)]);
    yticks([qq(1,1):0.2:qq(1,end)]);
    zlabel('Coefficient at NFCI');
    title('Equation for household credit');


%% Figure C3: Fitted quantile and univariate marginal quantiles (for specification with credit ratio)
Yfit = [mean(X)*Betas(1,:,1)',mean(X)*Betas(2,:,1)'];
for i=1:nboot
    Yfit_boot(i,:) = [mean(X)*squeeze(Betas_boot_all(i,1,:,1)),mean(X)*squeeze(Betas_boot_all(i,2,:,1))]'; 
end

indSQR = 5;
YfitSQR1 = mean(X)*BetaSQR(1,:,indSQR)';
YfitSQR2 = mean(X)*BetaSQR(2,:,indSQR)';

YfitSQR1_90 = [mean(X)*BetaSQR1_90(1,:,indSQR)',mean(X)*BetaSQR1_90(2,:,indSQR)'];
YfitSQR2_90 = [mean(X)*BetaSQR2_90(1,:,indSQR)',mean(X)*BetaSQR2_90(2,:,indSQR)'];

tmp1_center = YfitSQR1_90(1,1) + (YfitSQR1_90(1,2)-YfitSQR1_90(1,1))/2;
tmp2_center = YfitSQR2_90(1,1) + (YfitSQR2_90(1,2)-YfitSQR2_90(1,1))/2;

figure
set(gcf, 'Color', 'w');
hold on
a1 = scatter(Yfit(1,1),Yfit(1,2),80,'ko','MarkerFaceColor','k');
a2 = scatter(Yfit_boot(:,1),Yfit_boot(:,2),'bo');
scatter(Yfit(1,1),Yfit(1,2),80,'ko','MarkerFaceColor','k');
a4 = plot([tmp1_center,tmp1_center],[YfitSQR2_90(1,1),YfitSQR2_90(1,2)],'r','LineWidth',1.5);
plot([YfitSQR1_90(1,1),YfitSQR1_90(1,2)],[tmp2_center,tmp2_center],'r','LineWidth',1.5)
a3 = scatter(tmp1_center,YfitSQR2,80,'ro','MarkerFaceColor','r');
scatter(YfitSQR1,tmp2_center,80,'ro','MarkerFaceColor','r');
ax = gca;
ax.XAxisLocation = 'origin';
ax.YAxisLocation = 'origin';

xlabel('House price growth (%)');
ylabel('Household credit growth (%)');

title('[-0.5,-0.5]-th geometric quantile');
legend([a1,a2,a3,a4],'Geometric quantile','Bootstrapped values',...
    '0.25-th quantile of household credit growth/house price growth','90% confidence interval')



Yfit = [mean(X)*Betas(1,:,121)',mean(X)*Betas(2,:,121)'];
for i=1:nboot
    Yfit_boot(i,:) = [mean(X)*squeeze(Betas_boot_all(i,1,:,4)),mean(X)*squeeze(Betas_boot_all(i,2,:,4))]'; 
end

indSQR = 15;
YfitSQR1 = mean(X)*BetaSQR(1,:,indSQR)';
YfitSQR2 = mean(X)*BetaSQR(2,:,indSQR)';

YfitSQR1_90 = [mean(X)*BetaSQR1_90(1,:,indSQR)',mean(X)*BetaSQR1_90(2,:,indSQR)'];
YfitSQR2_90 = [mean(X)*BetaSQR2_90(1,:,indSQR)',mean(X)*BetaSQR2_90(2,:,indSQR)'];

tmp1_center = YfitSQR1_90(1,1) + (YfitSQR1_90(1,2)-YfitSQR1_90(1,1))/2;
tmp2_center = YfitSQR2_90(1,1) + (YfitSQR2_90(1,2)-YfitSQR2_90(1,1))/2;

figure
set(gcf, 'Color', 'w');
hold on
a1 = scatter(Yfit(1,1),Yfit(1,2),80,'ko','MarkerFaceColor','k');
a2 = scatter(Yfit_boot(:,1),Yfit_boot(:,2),'bo');
scatter(Yfit(1,1),Yfit(1,2),80,'ko','MarkerFaceColor','k');
a4 = plot([tmp1_center,tmp1_center],[YfitSQR2_90(1,1),YfitSQR2_90(1,2)],'r','LineWidth',1.5);
plot([YfitSQR1_90(1,1),YfitSQR1_90(1,2)],[tmp2_center,tmp2_center],'r','LineWidth',1.5)
a3 = scatter(tmp1_center,YfitSQR2,80,'ro','MarkerFaceColor','r');
scatter(YfitSQR1,tmp2_center,80,'ro','MarkerFaceColor','r');
ax = gca;
ax.XAxisLocation = 'origin';
ax.YAxisLocation = 'origin';

xlabel('House price growth (%)');
ylabel('Household credit growth (%)');

title('[0.5,0.5]-th geometric quantile');
legend([a1,a2,a3,a4],'Geometric quantile','Bootstrapped values',...
    '0.75-th quantile of household credit growth/house price growth','90% confidence interval');