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

%"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
 
%First compute statistics for the LTV distributions. 
%Then plot annual time series
tic
clear; clc; close all;
 
whichplot = 2; %1 if compute stats for LTV distributions; 2 if plot macro series

%For the part on LTV distributions, all variables labeled original_ltv_... 
%represent distributions of loan-level LTVs for the period indicated in the 
%variable name. Because these data are confidential, this variables are
%empty here.
 
stats = 1; %1: we compute the share of observations which we treat as 
%exceptions to the limit (see Section 3.2 of the paper); in addition, we 
%compute statistics for "sprinkling", i.e. for the alternative case where 
%we impose an LTV limit such that we move some LTVs to a level close to but 
%below the limit (see Section 3.3 of the paper for motivation and detail)

dist = 1; 
%This variable marks which distributions we use for computing the share of exceptions
%1: pre-cap 2010M8-2011M7, 106-cap 2011M11-2012M12. Others: 2013M4-M12 through 2017M4M12. 
%2: pre-cap 2010M4-2010M12, 106-cap 2012M4-2012M12. Others: same.
%3: pre-cap 2010, 106-cap 2011M11-2012M12. Others: same.
 

if whichplot == 1  
     
    if stats == 1
              
        %FIRST, COMPUTE STATISTICS RELEVANT FOR EXCEPTIONS
        
        %First compute fractions of observations above various marks
        himat = nan(7,8); 
        %(1,1) entry: nan. Rest of first row: dates for distributions
        %Rest of first column: marks for which we compute percentages above
        %the mark
        %The rest of the matrix fills in the share of observations (on a
        %0-1 scale) in that column's distribution for which the LTV exceeds
        %that row's mark
    
        %Fill in codes for empirical distributions, e.g. original_ltv_1011M8M7 original_ltv_1112M11M12 original_ltv_13M4M12 original_ltv_14M4M12 original_ltv_15M4M12 original_ltv_16M4M12 original_ltv_17M4M12 
        if dist == 1
            himat(1,2:8) = [101187 11121112 2013412 2014412 2015412 2016412 2017412];
        elseif dist == 2
            himat(1,2:8) = [2010412 2012412 2013412 2014412 2015412 2016412 2017412];
        elseif dist == 3
            himat(1,2:8) = [2010 11121112 2013412 2014412 2015412 2016412 2017412];
        end
            
        %Fill in the marks
        himat(2:7,1) = [106:-1:101]';
        
        %Compute fraction of observations exceeding the various marks in
        %the pre-cap distribution
        if dist == 1
            for j = 2:7 %for each of the six marks
                clear hi
                hi = nan(size(original_ltv_1011M8M7)); %vector to be filled in with hi/not marks for that loop
                for i = 1 : size(original_ltv_1011M8M7,1)
                    if original_ltv_1011M8M7(i) > himat(j,1)
                        hi(i) = 1; %For every LTV above the mark, set the high marker to 1
                    else
                        hi(i) = 0; %High marker 0: not above the mark.
                    end
                end
                himat(j,2) = mean(hi); %Share of observations above the mark (in units from 0 to 1)
            end
        elseif dist == 2
            for j = 2:7 %for each of the six marks
                clear hi
                hi = nan(size(original_ltv_10M4M12)); %vector to be filled in with hi/not marks for that loop
                for i = 1 : size(original_ltv_10M4M12,1)
                    if original_ltv_10M4M12(i) > himat(j,1)
                        hi(i) = 1; %For every LTV above the mark, set the high marker to 1
                    else
                        hi(i) = 0; %High marker 0: not above the mark.
                    end
                end
                himat(j,2) = mean(hi); %Share of observations above the mark (in units from 0 to 1)
            end            
        elseif dist == 3
            for j = 2:7 %for each of the six marks
                clear hi
                hi = nan(size(original_ltv_2010)); %vector to be filled in with hi/not marks for that loop
                for i = 1 : size(original_ltv_2010,1)
                    if original_ltv_2010(i) > himat(j,1)
                        hi(i) = 1; %For every LTV above the mark, set the high marker to 1
                    else
                        hi(i) = 0; %High marker 0: not above the mark.
                    end
                end
                himat(j,2) = mean(hi); %Share of observations above the mark (in units from 0 to 1)
            end
        end
        
        %Compute fraction of observations exceeding the various marks in empirical 106-cap distribution        
        if dist == 1 | dist == 3
            for j = 2:7 %for each of the six marks
                clear hi
                hi = nan(size(original_ltv_1112M11M12)); %vector to be filled in with hi/not marks for that loop
                for i = 1 : size(original_ltv_1112M11M12,1)
                    if original_ltv_1112M11M12(i) > himat(j,1)
                        hi(i) = 1; %For every LTV above the mark, set the high marker to 1
                    else
                        hi(i) = 0; %High marker 0: not above the mark.
                    end
                end
                himat(j,3) = mean(hi); %Share of observations above the mark (in units from 0 to 1)
            end
        elseif dist == 2
            for j = 2:7 %for each of the six marks
                clear hi
                hi = nan(size(original_ltv_12M4M12)); %vector to be filled in with hi/not marks for that loop
                for i = 1 : size(original_ltv_12M4M12,1)
                    if original_ltv_12M4M12(i) > himat(j,1)
                        hi(i) = 1; %For every LTV above the mark, set the high marker to 1
                    else
                        hi(i) = 0; %High marker 0: not above the mark.
                    end
                end
                himat(j,3) = mean(hi); %Share of observations above the mark (in units from 0 to 1)
            end           
        end
        
        %Compute fraction of observations exceeding the various marks in empirical 105-cap distribution        
        for j = 2:7 %for each of the six marks
            clear hi
            hi = nan(size(original_ltv_13M4M12)); %vector to be filled in with hi/not marks for that loop
            for i = 1 : size(original_ltv_13M4M12,1)
                if original_ltv_13M4M12(i) > himat(j,1)
                    hi(i) = 1; %For every LTV above the mark, set the high marker to 1
                else
                    hi(i) = 0; %High marker 0: not above the mark.
                end
            end
            himat(j,4) = mean(hi); %Share of observations above the mark (in units from 0 to 1)
        end
        
        %Compute fraction of observations exceeding the various marks in empirical 104-cap distribution        
        for j = 2:7 %for each of the six marks
            clear hi
            hi = nan(size(original_ltv_14M4M12)); %vector to be filled in with hi/not marks for that loop
            for i = 1 : size(original_ltv_14M4M12,1)
                if original_ltv_14M4M12(i) > himat(j,1)
                    hi(i) = 1; %For every LTV above the mark, set the high marker to 1
                else
                    hi(i) = 0; %High marker 0: not above the mark.
                end
            end
            himat(j,5) = mean(hi); %Share of observations above the mark (in units from 0 to 1)
        end
        
        %Compute fraction of observations exceeding the various marks in empirical 103-cap distribution        
        for j = 2:7 %for each of the six marks
            clear hi
            hi = nan(size(original_ltv_15M4M12)); %vector to be filled in with hi/not marks for that loop
            for i = 1 : size(original_ltv_15M4M12,1)
                if original_ltv_15M4M12(i) > himat(j,1)
                    hi(i) = 1; %For every LTV above the mark, set the high marker to 1
                else
                    hi(i) = 0; %High marker 0: not above the mark.
                end
            end
            himat(j,6) = mean(hi); %Share of observations above the mark (in units from 0 to 1)
        end
        
        %Compute fraction of observations exceeding the various marks in empirical 102-cap distribution        
        for j = 2:7 %for each of the six marks
            clear hi
            hi = nan(size(original_ltv_16M4M12)); %vector to be filled in with hi/not marks for that loop
            for i = 1 : size(original_ltv_16M4M12,1)
                if original_ltv_16M4M12(i) > himat(j,1)
                    hi(i) = 1; %For every LTV above the mark, set the high marker to 1
                else
                    hi(i) = 0; %High marker 0: not above the mark.
                end
            end
            himat(j,7) = mean(hi); %Share of observations above the mark (in units from 0 to 1)
        end
        
        %Compute fraction of observations exceeding the various marks in empirical 101-cap distribution        
        for j = 2:7 %for each of the six marks
            clear hi
            hi = nan(size(original_ltv_17M4M12)); %vector to be filled in with hi/not marks for that loop
            for i = 1 : size(original_ltv_17M4M12,1)
                if original_ltv_17M4M12(i) > himat(j,1)
                    hi(i) = 1; %For every LTV above the mark, set the high marker to 1
                else
                    hi(i) = 0; %High marker 0: not above the mark.
                end
            end
            himat(j,8) = mean(hi); %Share of observations above the mark (in units from 0 to 1)
        end
        clear hi
        
        %Next, compute ratios of fractions => share of observations that
        %stays behind as exceptions
        %106/precap, 105/precap, 104/precap, 103/precap, 102/precap, 101/precap,
        
        exbigmat = nan(8,22);
        %exceptions. First row: coded distributions in numerator of ratio, e.g. 106
        %Second row: coded distrubtion in nominator of ratio, e.g. 107
        %First column: mark.
        %Body: fraction of exceptions with those distributions for that
        %mark
        exbigmat(1,2:end) = [106 105 104 103 102 101 105 104 103 102 101 104 103 102 101 103 102 101 102 101 101];
        exbigmat(2,2:end) = [107 107 107 107 107 107 106 106 106 106 106 105 105 105 105 104 104 104 103 103 102];
        exbigmat(3:end,1) = [106:-1:101]';
        
        %First fill in all the ratios
        for j = 2:7
            i = j+1;
            exbigmat(i,2) = himat(j,3)/himat(j,2);
            exbigmat(i,3) = himat(j,4)/himat(j,2);
            exbigmat(i,4) = himat(j,5)/himat(j,2);
            exbigmat(i,5) = himat(j,6)/himat(j,2);
            exbigmat(i,6) = himat(j,7)/himat(j,2);
            exbigmat(i,7) = himat(j,8)/himat(j,2);
            
            exbigmat(i,8) = himat(j,4)/himat(j,3);
            exbigmat(i,9) = himat(j,5)/himat(j,3);
            exbigmat(i,10) = himat(j,6)/himat(j,3);
            exbigmat(i,11) = himat(j,7)/himat(j,3);
            exbigmat(i,12) = himat(j,8)/himat(j,3);
            
            exbigmat(i,13) = himat(j,5)/himat(j,4);
            exbigmat(i,14) = himat(j,6)/himat(j,4);
            exbigmat(i,15) = himat(j,7)/himat(j,4);
            exbigmat(i,16) = himat(j,8)/himat(j,4);
            
            exbigmat(i,17) = himat(j,6)/himat(j,5);
            exbigmat(i,18) = himat(j,7)/himat(j,5);
            exbigmat(i,19) = himat(j,8)/himat(j,5);
            
            exbigmat(i,20) = himat(j,7)/himat(j,6);
            exbigmat(i,21) = himat(j,8)/himat(j,6);
            
            exbigmat(i,22) = himat(j,8)/himat(j,7);           
        end
               
        %Then set ratios to nan when the distributions are on the same 
        %side of the mark so they don't reflect exceptions
        exbigmatnan = exbigmat; %Give it a different name so I can still check the other ratios when I need to
        for i = 3:8 %rows
            for j = 2:22 %columns
                if (exbigmatnan(1,j) > exbigmatnan(i,1) & exbigmatnan(2,j) > exbigmatnan(i,1)) | (exbigmatnan(1,j) <= exbigmatnan(i,1) & exbigmatnan(2,j) <= exbigmatnan(i,1))
                   exbigmatnan(i,j) = nan; 
                end
            end
        end
        %Add a column at the left with means for every mark
        col = [nan(2,1); mean(rmmissing(exbigmatnan(3,2:end))) ; mean(rmmissing(exbigmatnan(4,2:end))) ; mean(rmmissing(exbigmatnan(5,2:end))); mean(rmmissing(exbigmatnan(6,2:end))); mean(rmmissing(exbigmatnan(7,2:end))); mean(rmmissing(exbigmatnan(8,2:end)))];
        exbigmatnan = [col exbigmatnan];
        
        %Add a row at the bottom with the mean of these means in the first
        %entry
        row = [mean(exbigmatnan(3:8,1)) nan(1,22)];
        exbigmatnan = [exbigmatnan; row];
        
        %Add another column at the beginning with the share of exceptions
        %for only the occasion were the mark just became binding as a
        %limit, relative to the distribution from the spell right before that        
        %The last entry is the mean
        col2 = [nan(2,1); exbigmatnan(3,3); exbigmatnan(4,9); exbigmatnan(5,14); exbigmatnan(6,18); exbigmatnan(7,21); exbigmatnan(8,23);nan];
        col2(end) = mean(col2(3:end-1));
        exbigmatnan = [col2 exbigmatnan];
        
        %And add a column which computes the share of exceptions above 106
        %relative to the precap distribution
        col3 = [nan(2,1); [exbigmatnan(3,4:9)]'; nan];
        col3(end) = mean(col3(3:end-1));
        exbigmatnan = [col3 exbigmatnan];
        
        %And add a column with the share of exceptions
        %for only the occasion where the mark just became binding
        %relative to the pre-cap distribution
        col4 = [nan(2,1); exbigmatnan(3,5); exbigmatnan(4,6); exbigmatnan(5,7); exbigmatnan(6,8); exbigmatnan(7,9); exbigmatnan(8,10); nan];
        col4(end) = mean(col4(3:end-1)); 
        exbigmatnan = [col4 exbigmatnan];
        %col4, i.e. the first column of exbigmatnan, IS THE MAIN VERSION.
        %exbigmatnan(9,1) is the output: share of exceptions
        
      
      
        %NEXT, COMPUTE SHARE IN SIX BINS 100-106
        %AND DO FURTHER WORK WITH THEM, SO AS TO INFORM SPRINKLING.
        
        %Only do with baseline distributions

        if dist == 1
        
            sprinkmat = nan(7,7); %Matrix to fill in the shares and derived states for sprinkling
        
            %Columns: 
            %lower edge bin | upper edge bin | share
            %original_ltv_1011M8M7 | share original_ltv_1112M11M12 | difference
            %| "share of the differences"

            %Rows: 
            %names | bin by bin

            %Fill in the distribution names
            sprinkmat (1,3:4) = [101187 11121112];
             
            %Fill in the lower and upper edges of the binds
            sprinkmat(2:7,1) = [100:105]';
            sprinkmat(2:7,2) = [101:106]';


            for j = 2:7 %for each of the six marks
                %Shares in pre-cap distribution
                xpre = nan(size(original_ltv_1011M8M7)); %vector to be filled in with yes/not marks for that loop
                for i = 1 : size(original_ltv_1011M8M7,1)
                    if original_ltv_1011M8M7(i) > sprinkmat(j,1) & original_ltv_1011M8M7(i) <= sprinkmat(j,2)
                        xpre(i) = 1; %For every LTV above the mark, set the high marker to 1
                    else
                        xpre(i) = 0; %High marker 0: not above the mark.
                    end
                end
                sprinkmat(j,3) = mean(xpre); %Share of observations above the mark (in units from 0 to 1)
                clear xpre

                %Shares in 106-cap distribution
                xpost = nan(size(original_ltv_1112M11M12)); %vector to be filled in with yes/not marks for that loop
                for i = 1 : size(original_ltv_1112M11M12,1)
                    if original_ltv_1112M11M12(i) > sprinkmat(j,1) & original_ltv_1112M11M12(i) <= sprinkmat(j,2)
                        xpost(i) = 1; %For every LTV above the mark, set the high marker to 1
                    else
                        xpost(i) = 0; %High marker 0: not above the mark.
                    end
                end
                sprinkmat(j,4) = mean(xpost); %Share of observations above the mark (in units from 0 to 1)
                clear xpost
            
                %Now compute difference of bin shares
                %106-cap - pre-cap
                sprinkmat(j,5) = sprinkmat(j,4)- sprinkmat(j,3);   
            
            end
            
            %Now compute "difference shares"
            sprinksum = sum(sprinkmat(2:7,5));
            for j = 2:7
                sprinkmat(j,6) = sprinkmat(j,5) / sprinksum;
            end
            %Column 6 is the outcome! 
            %On imposition, for LTVs on which we impose the limit, 
            %we set some to 106, ...., 101, as indicated by the shares in
            %column 6.

            %But column 7 is what we will reall use: cumulate the outcome in column6, starting from
            %105-106 bin
            sprinkmat(7,7) = sprinkmat(7,6);
            for j = 6:-1:2
                sprinkmat(j,7) = sprinkmat(j+1,7) + sprinkmat(j,6);
            end

        end        
  
    end 
elseif whichplot == 2
    %Plot annual time series
        
    %Go to VECM directory because that's where these data are
    cd('JMCB_replic_deveirman-jong/VECM estimation')
    
    %Import the time series
    E=importdata('for_annual_plot.csv'); %Matlab produces a structure  
    colheaders=E.colheaders'; %Take colheaders from structure
    E=E.data; %Rename D as the numerical part of the original structure
    year=E(:,1);ltv_comb_nz=E(:,5);hpcq_r=E(:,2);ybbpm15=E(:,3);rh_ea1=E(:,4);
    lim=nan(size(year)); lim(34:40)=[106:-1:100]';
    
    figure
    haha = tight_subplot(2,2,[.16 .09],[.10 .09],[.08 .02])
    
    %Loan-to-Value ratio
    %figure
    %subplot(2,2,1)
    axes(haha(1));
    hey = plot(year,lim,'k')
    hold on
    ha = plot(year,ltv_comb_nz,'k')
    hold off
    legend('LTV limit','average LTV','Fontsize',20,'NumColumns',2)
    lgd.NumColumns = 2;
    set(ha(1),'LineWidth',1)
    set(hey(1),'LineWidth',0.5,'Marker','*')
    axis([year(1) year(end) 80 110])
    set(gca,'FontSize',19) %This affects the font size of everything (except legends if there were legends)
    title('Loan-to-Value Ratio: Limit and Average','FontSize',25), ylabel('%','FontSize',22)%, xlabel('Year','FontSize',20)
    yticks([80:10:100])
    xticks([1980:5:2015])
    grid on; %ax=gca; ax.GridAlpha = 0.30;
    %subplot(2,2,2) %Real mortgage rate
    axes(haha(2));
    hb=plot(year,rh_ea1,'k')
    set(hb(1),'LineWidth',1)
    axis([year(1) year(end) 0 10])
    set(gca,'FontSize',19)
    title('Real Mortgage Rate','FontSize',25), ylabel('%','FontSize',22)%, xlabel('Year','FontSize',20)
    yticks([0:4:10])
    xticks([1980:5:2015])
    grid on; ax=gca; ax.GridAlpha = 0.30;
    %subplot(2,2,3) %GDP
    axes(haha(3));
    hc=plot(year,ybbpm15,'k')
    set(hc(1),'LineWidth',1)
    axis([year(1) year(end) 200 800])
    set(gca,'FontSize',19)
    title('Real Gross Domestic Product','FontSize',25), ylabel('Billions 2015 euro','FontSize',22)%, xlabel('Year','FontSize',20)
    yticks([200:200:800])
    xticks([1980:5:2015])
    grid on; ax=gca; ax.GridAlpha = 0.30;
    %subplot(2,2,4) %Relative price of housing
    axes(haha(4));
    hd=plot(year,hpcq_r,'k') %Note this is the rescaled relative price of housing so 2010=100
    set(hd(1),'LineWidth',1)
    axis([year(1) year(end) 0 150])
    set(gca,'FontSize',19)
    title('Real House Prices','FontSize',25), ylabel('2015=100','FontSize',22)%, xlabel('Year','FontSize',20)
    yticks([0:50:150])
    xticks([1980:5:2015])
    grid on; ax=gca; ax.GridAlpha = 0.30;
    
    
       
end
toc