%% Application Plot File
% This file plots the main graphs for the application.
% File is divided into one section per 
%% Plot IP growth
fig=figure('Name','IP growth');
fig.PaperPositionMode = 'manual';
fig.PaperOrientation = 'landscape';
fig.PaperSize = [ 10  6];
fig.PaperPosition = [ 0 0 10 6]; 

L0table = readtable('IP98updated_formulas.xlsx','Sheet','L0','Range',[4 2]);
L0table = table2array(L0table);
IPgrowth = 1200*log(L0table(2:end)./L0table(1:end-1));
date = datetime(1972,1+(1:length(IPgrowth)),1);

plot(date(1:ts(end)-DataOffset),IPgrowth(1:ts(end)-DataOffset),'k','LineWidth',1.75)
fig.CurrentAxes.XLim = [ date(1) - years(3) date(ts(end)-DataOffset) + years(3)];
for t = 1:length(ts)-1
    xline(date(ts(t)-DataOffset),'k--','LineWidth',1)
    text(date(ts(t) -DataOffset -round(RW/2)),-40,['\sigma = ' num2str(std(IPgrowth(ts(t)-DataOffset -RW:ts(t)-DataOffset )),'%.2f')])
end
text(date(ts(end) -DataOffset -round(RW/2)),-40,['\sigma = ' num2str(std(IPgrowth(ts(end)-DataOffset -RW:ts(end)-DataOffset )),'%.2f')])
print('-dpdf', './Graphs/IPgrowth.pdf','-fillpage')

%% MSE coefficient
PaperSizeMSE=[ 10 7];
titleFontSize= 14;
FontSize = 13;
fig=figure('Name','MSE Coef-CV');
fig.PaperPositionMode = 'manual';
fig.PaperOrientation = 'landscape';
fig.PaperSize = PaperSizeMSE;
fig.PaperPosition = [ 0 0 1 1];
plotMSE = 100*MSE(:,:,ts);
TITLE={'1972 to 1983','1984 to 1995','1996 to 2007','2008 to 2019'};
tiledlayout(2,2, 'Padding', 'none', 'TileSpacing', 'compact');

k=0;
for t=ts
    k=k+1;
    
    ax = nexttile;
    h=plot(squeeze(plotMSE(:,:,k)),'LineWidth',2);
    h(2).LineStyle = ':';
    h(3).LineStyle = '--';
    h(3).Color = 'k';
    hold on
    
    [y,x]=min(squeeze(plotMSE(:,:,k)));
    scatter(x,y,50,[h(1).Color ;h(2).Color; h(3).Color],'d','fill')
    
    xticks(linspace(0,1, size(plotMSE,1)))
    xlim([1 size(plotMSE,1)])
    xticklabels([0,1])
    if k == 1 || k == 3
        ylabel('Minimal MSE', 'FontSize',FontSize )
        if k == 3
                xlabel('Hyperparameter grid', 'FontSize',FontSize )
        end
    elseif k == 4
        XLABEL = xlabel('alpha in adapt. ENET', 'FontSize',FontSize );
        %         XLABEL.Position(3)=-1;
        xlabel('Hyperparameter grid', 'FontSize',FontSize )
        ax.YTickLabel={''};
    else
        ax.YTickLabel={''};
    end
    ylim([0.99*min(min(min(plotMSE))) 1.01*max(max(max(plotMSE)))])
    title(TITLE{k}, 'FontSize',titleFontSize )
    ytickformat('%.1g%%')
    ax.YGrid = 'on';
    
end

legend('AENET: OLS initial guess','AENET: e.-net initial guess','Elastic-net: no inital guess', 'FontSize',FontSize,'Location',  'east')
% legend boxoff
print('-dpdf', './Graphs/A_MSE.pdf','-fillpage')

%% MSE sigma
fig = figure('Name','MSE Sigma CV');
fig.PaperPositionMode = 'manual';
fig.PaperOrientation = 'landscape';
fig.PaperSize = PaperSizeMSE;
fig.PaperPosition = [ 0 0 1 1];
tiledlayout(2,2, 'Padding', 'none', 'TileSpacing', 'compact');

plotLL = 100*squeeze(LL(:,:,ts));
k=0;
for t=ts
    k = k+1;
    ax = nexttile;
    h = plot(squeeze(plotLL(:,:,k)),'LineWidth',2);
    h(2).LineStyle = ':';     h(3).LineStyle = '--';     h(3).Color = 'k';
    hold on
    
    [y,x] = min(squeeze(plotLL(:,:,k)));
    hColor = [];
    for hh= 1:length(h)
        hColor = [hColor;h(hh).Color];
    end
    scatter(x,y,50,hColor,'d','fill')
    if k==1 || k == 3
        ylabel('Mean squared CV-loss', 'FontSize',FontSize )
        if k ==3
            xlabel('Hyperparameter grid', 'FontSize',FontSize )
        end
    elseif k == 4
        xlabel('Hyperparameter grid', 'FontSize',FontSize )
    end
    
    [minY,minX] = min(y);
    scatter(x(minX),minY,75,'k','LineWidth',2)
    text(x(1)+4,y(1),['eta = ' num2str(eta(t))])
    
    ylim([min(min(min(plotLL(:,:,k))))-0.05 max(max(max(plotLL(1,:,k))))])
    xlim([1 size(plotLL,1)])
    ytickformat('%.2f')

    xticks('')
    title(TITLE{k}, 'FontSize',titleFontSize )
    ax.YGrid = 'on';
end
legend('Adap.-Threshold','GLASSO','Ledoit-Wolf', 'FontSize',FontSize)
print('-dpdf', './Graphs/Sig_MSE.pdf','-fillpage')

%% Network plots Reg
networkPlot(DfinalReg,meanWeight,names,ts)
print('-dpng', './Graphs/Spillovers.png')

%% Network plots OLS
networkPlot(DfinalOLS,meanWeight,names,ts)
print('-dpng', './Graphs/SpilloversOLS.png')

%% Histogram ingoing connectedness
PaperSizePlots = [10 6];
fig = figure('Name','Histogram Ingoing Connectedness');
fig.PaperPositionMode = 'manual';
fig.PaperOrientation = 'landscape';
fig.PaperSize = PaperSizePlots;
fig.PaperPosition = [ 0 0 1 1];
tiledlayout(2,2, 'Padding', 'none', 'TileSpacing', 'compact');

k=0;

for t=ts
    k=k+1;
    weightComponent = N*meanWeight(:,k)/sum(meanWeight(:,k));
    Table=DfinalReg(:,:,t)-diag(diag(DfinalReg(:,:,t))); %
    
    deg2=sum(Table,2)';
    [totest, rankings]=sort(meanWeight(:,k),'descend');
    deg2ranked=deg2(rankings);
    deg2stairs=deg2ranked(1:floor(N/2));
    
    [~, maxI]=max(meanWeight(:,k));
    H1 = nexttile;
    if k == 1
        bins=linspace(0,max(deg2),N/2);
    end
    H = histogram(deg2,bins);
    ylim([0 15])

    hold on
    
    H2 = histogram(deg2stairs,bins);
    title(TITLE{k})
    hold off
    if k>1
        H1.XLim = limsX;
        H1.YLim = limsY;
    else
        limsX = H1.XLim;
        limsY = H1.YLim;
    end
    if k > 2
        xlabel('In-Connectedness')
    end
    
    if k ==2 || k ==4
        yticklabels({})
    end
    
end
legend('All sectors','50% sectors with highest weight')
legend boxoff
print('-dpdf', './Graphs/HistInCon.pdf','-fillpage')
%% Scatter Outgoing Connectedness
fig = figure('Name','Scatter Outgoing Connectedness');
fig.PaperPositionMode = 'manual';
fig.PaperOrientation = 'landscape';
fig.PaperSize = [ 12.5  4];
fig.PaperPosition = [ 0 0 1 1];
tiledlayout(1,length(ts)-1, 'Padding', 'none', 'TileSpacing', 'compact');


k=0;
for t=ts
    k=k+1;
    Table=DfinalReg(:,:,t)-diag(diag(DfinalReg(:,:,t))); %
    
    if k==1
        deg2ref = sum(Table,1);
    else
        H1 = nexttile;
        deg2 = sum(Table,1);
        sz = 100*meanWeight(:,k);
        c = deg2-deg2ref;
        H2 = scatter(deg2,deg2ref,sz,c,'filled','MarkerFaceAlpha',0.85,'LineWidth',1.5);
        [totest, rankings]=sort(c);     % rank by changes
        [~, rankingsSz]=sort(sz);       % rank by weights    

        colormap(0.7*jet)
        Alignment={'left','right'};
        for i=1:length(deg2)
            if i<=4 || i >= length(deg2) -1
                top = names(rankings(i)); top=char(string(top));  
                [~,ranks2]=sort(deg2ref,'descend');
                if i<=4
                    line([deg2(rankings(i)) 0.6],[deg2ref(rankings(i)) max(deg2ref)-i*0.075] ,'LineWidth',1,'Color',[0 0 0])
                    text(0.6,max(deg2ref)-i*0.075,top); %+0.05*(-floor(rem(I,2))+1)
                elseif i >= length(deg2) -1 && deg2(rankings(i))>0.75 
                    line([deg2(rankings(i)) 1],[deg2ref(rankings(i)) min(deg2ref)-(i-length(deg2ref))*0.075+0.05] ,'LineWidth',1,'Color',[0 0 0])
                    text(1,min(deg2ref)-(i-length(deg2ref))*0.075 + 0.05,top); %+0.05*(-floor(rem(I,2))+1)
                elseif k ==4 && i >= length(deg2) -1
                    top = names(rankingsSz(i)); top=char(string(top));
                    line([deg2(rankingsSz(i)) 1],[deg2ref(rankingsSz(i)) min(deg2ref)-(i-length(deg2ref))*0.075+0.05] ,'LineWidth',1,'Color',[0 0 0])
                    text(1,min(deg2ref)-(i-length(deg2ref))*0.075 + 0.05,top); %+0.05*(-floor(rem(I,2))+1)
                end
            end
        end
        hold on
        plot(linspace(0,max(deg2ref)*1.05,10),linspace(0,max(deg2ref)*1.05,10),'k')
        H1.YLim = [0 max(deg2ref)*1.05];
        H1.XLim = [0 max(deg2ref)*1.05];
        H1.XTick = [0:0.2:1.6];
        H1.YTick = [0:0.2:1.6];
        H1.YTickMode = 'manual';
        H1.XTickMode = 'manual';

        grid on
        box on
        
    end
    ylabel(['Out-con. ' TITLE{1}])
    xlabel(['Out-con. ' TITLE{k}])
end
print('-dpdf', './Graphs/ScatterOutCon.pdf','-fillpage')
%% Testing Area Inventories
Table = DfinalReg(:,:,ts(1))-diag(diag(DfinalReg(:,:,ts(1)))); %
deg1 = sum(Table,1);
Table = DfinalReg(:,:,ts(2))-diag(diag(DfinalReg(:,:,ts(2)))); %
deg2 = sum(Table,1);

% adjust for break in data series
changeShip = Shipments(1:end-1,:)./Shipments(2:end,:);
changeShip(132,:) = 1;
for tt = 0:131
    Shipments(132-tt,:) = Shipments(133-tt,:).*changeShip(132-tt,:);
end
changeInv = Inventories(1:end-1,:)./Inventories(2:end,:);
changeInv(132,:) = 1;
for tt = 0:131
    Inventories(132-tt,:) = Inventories(133-tt,:).*changeInv(132-tt,:);
end
InvToShip = Inventories./Shipments;
InvToShip(InvToShip==0) = NaN;

% Calculate mean levels pre and post GM
InvToShipPreGM = mean(InvToShip(1:ts(1),:));  
InvToShipPostGM = mean(InvToShip(ts(1)+1:end,:));

% changes for regressions
changeInvToShip = InvToShipPostGM./InvToShipPreGM -1; log(InvToShipPostGM) - log(InvToShipPreGM );
changeDeg = deg2 - deg1;

tbl = table(changeDeg',changeInvToShip',...
    'VariableNames',{'changeDeg','dInvToShip'});

mdlAll = fitlm(tbl,'changeDeg ~ dInvToShip ')

fig=figure('Name','Regression Graph');
fig.PaperPositionMode = 'manual';
fig.PaperOrientation = 'landscape';
fig.PaperSize = [ 10  4];
fig.PaperPosition = [ 0 0 1 1];
plot(mdlAll,'Color',[0.25  0.25 0.7])

ylabel('Difference Out-connectedness')
xlabel('%-difference inventory to shipments ratio')
title('')
grid on
legend('Data','Fit','95% Confidence bounds')
print('-dpdf', './Graphs/RegressionGraph.pdf','-fillpage')

%% Histogram Outgoing Connectedness
fig = figure('Name','Histogram Outgoing Connectedness');
fig.PaperPositionMode = 'manual';
fig.PaperOrientation = 'landscape';
fig.PaperSize = PaperSizePlots;
fig.PaperPosition = [ 0 0 1 1];
tiledlayout(2,2, 'Padding', 'none', 'TileSpacing', 'compact');

percentageHighlight = 25;
k=0;
for t=ts
    k=k+1;
    Table=DfinalReg(:,:,t)-diag(diag(DfinalReg(:,:,t))); %
    
    deg2 = sum(Table,1);
    if k==1
        [totest, rankings]=sort(deg2,'descend');
        [maxD, maxI]=max(deg2);
    end
    deg2ranked=deg2(rankings);
    deg2stairs=deg2ranked(1:floor(percentageHighlight*N/100));
    
    deg2stairs2=deg2ranked(end-floor(percentageHighlight*N/100):end);
    
    nexttile
    bins=linspace(0,maxD,N/2);
    H = histogram(deg2,bins);
    ylim([0 12])
     
    H2 = histogram(deg2stairs,bins);
    
    title(TITLE{k})
    hold off
    
    if k == 2 || k == 4
        yticklabels({})
    end
    if k > 2
        xlabel('Out-Connectedness')
    end
    
end
legend('All sectors',[ num2str(percentageHighlight) '% largest Out-Con. pre-1984'])
legend boxoff
print('-dpdf', './Graphs/HistOutCon.pdf','-fillpage')

%% scatter ingoing connectedness
fig=figure('Name','Ingoing vs. Sector Weight');
fig.PaperPositionMode = 'manual';
fig.PaperOrientation = 'landscape';
fig.PaperSize = PaperSizePlots;
fig.PaperPosition = [ 0 0 1 1];
tiledlayout(2,2, 'Padding', 'none', 'TileSpacing', 'compact');

k=0;
for t=ts
    k=k+1;
    weightComponent = meanWeight(:,k)/sum(meanWeight(:,k));
    Table=DfinalReg(:,:,t)-diag(diag(DfinalReg(:,:,t))); %
    deg2=sum(Table,2)';
    [totest, rankings]=sort(meanWeight(:,k),'descend');
    deg2ranked=deg2(rankings);
    deg2stairs=deg2ranked(1:floor(N/2));
    
    [~, maxI] = max(meanWeight(:,k));
    
    H1 = nexttile;
    sz = 100*meanWeight(:,k);
    H2 = scatter(deg2,weightComponent,sz,'filled','MarkerFaceAlpha',0.85,'LineWidth',1.5);
    colormap(jet)
    ls=lsline;
    ls.LineWidth=1.5; ls.Color=0.25*[1 1 1];
    title(TITLE{k})
    hold off
    
    if k == 2 || k == 4
        yticklabels({})
        H1.YLim = limsY;
        H1.XLim = limsX;
        
    else
        H1.YLim(2) = H1.YLim(2)*1.05;
        limsY = H1.YLim;
        limsX = H1.XLim;
        
        ylabel('Sector weight')
    end
    if k> 2
        xlabel('In-Connectedness')
    end
    H1.XLim = [0 1];
    box on
    grid on
end
print('-dpdf', './Graphs/ScatterInCon.pdf','-fillpage')

%% BumpChart
clear bumpRanks sortings NewSortings outDeg sortings outDegOff outDegDiag sortingsNew NewSortings
fig = figure('Name','Main Bumpchart');
fig.PaperPositionMode = 'manual';
fig.PaperOrientation = 'landscape';
fig.PaperSize = [12  7];
fig.PaperPosition = [ 0 0 1 1];
dimensionSum = 1;

BumpChartFont = 8;
k=0;
for t=ts
    k=k+1;
    weightComponent = (meanWeight(:,k)/sum(meanWeight(:,k))); weightComponent = weightComponent.^2;
    weightInnovations = ones(size(weightComponent));diag(FinalSigmaReg(:,:,t)).^(-0.5);
    
    outDeg(:,k)=sum((weightComponent)'*(DfinalRegIRF(:,:,t))*diag(weightInnovations),dimensionSum);
    outDegOff(:,k)=sum(diag(weightComponent)'*(DfinalRegIRF(:,:,t)-diag(diag(DfinalRegIRF(:,:,t))))*diag(weightInnovations),dimensionSum);
    outDegDiag(:,k)=sum(diag(weightComponent)'*(diag(diag(DfinalRegIRF(:,:,t))))*diag(weightInnovations),dimensionSum);
    
    [bumpRanks(k,:), sortings(k,:)] = sort(outDeg(:,k),'descend');
end

% % % % % Define Color for bumpchart % % % % %
ColorsScale= 0.9*[autumn(ceil(length(outDeg)/3)); flipud(summer(ceil(length(outDeg)/3))); winter(ceil(length(outDeg)/3))];
lineColors = [        0    0.4470    0.7410
    0.8500    0.3250    0.0980
    0.9290    0.6940    0.1250
    0.4940    0.1840    0.5560
    0.4660    0.6740    0.1880
    0.3010    0.7450    0.9330
    0.6350    0.0780    0.1840];
numberOfDifColors = 6;
[rline, bline, gline] = deal([]);
for kk = 1:numberOfDifColors
    rline = [rline;linspace(lineColors(kk,1), lineColors(kk+1,1),ceil(length(outDeg)/numberOfDifColors))'];
    gline = [gline;linspace(lineColors(kk,2), lineColors(kk+1,2),ceil(length(outDeg)/numberOfDifColors))'];
    bline = [bline;linspace(lineColors(kk,3), lineColors(kk+1,3),ceil(length(outDeg)/numberOfDifColors))'];
end
ColorsScale = [rline gline bline];

% jet(ceil(length(outDeg)));
colorweight = 1;
ColorsScale= colorweight*ColorsScale + (1-colorweight)*ones(size(ColorsScale));
distWhite=[0 linspace(0.75,0.125,5) zeros(1,length(outDeg)-6)];
% % % % % % % % %
for i=1:length(outDeg)
    for j=1:length(ts)
        
        sortingsNew(j,sortings(j,i),1) = sum(bumpRanks(j,:)) - sum(bumpRanks(j,i:end) )+distWhite(i);
        
        sortingsNew(j,sortings(j,i),2) = sum(bumpRanks(j,:)) - sum(bumpRanks(j,i+1:end) ); %
    end
end

k=0;
Xs=0.5:0.025:(length(ts)+0.5);
for x=Xs
    k=k+1;
    weightings = zeros(length(ts),1);
    if x<1
        firstWeight=1;
        lastWeight=ceil(x);
    elseif x>length(ts)
        firstWeight=floor(x);
        lastWeight=length(ts);
    else
        firstWeight=floor(x);
        lastWeight=ceil(x);
    end
    weightings(firstWeight:lastWeight) = pdf(makedist('Normal','mu',x,'sigma',0.2),firstWeight:lastWeight);
    weightings = weightings/sum(weightings);
    NewSortings(k,:,1) = weightings'*sortingsNew(:,:,1);
    NewSortings(k,:,2) = weightings'*sortingsNew(:,:,2);
end
hold on

for i=1:length(outDeg)-0
    fill([Xs, fliplr(Xs)],[NewSortings(:,length(outDeg)+1-i,1)' fliplr( NewSortings(:,length(outDeg)+1-i,2)' ) ],ColorsScale(i,:),'EdgeColor',0.125*[1 1 1])
end

for i=1:6
    for j=1:length(ts)
        txt = convertStringsToChars(names{sortings(j,i)});
        
        txt =[txt ' (' num2str( round( 100*outDegOff(sortings(j,i),j)/outDeg(sortings(j,i),j) ,0) ) '%)'];
        if j<=1 || i<=3
            text(j,0.5*sortingsNew(j,sortings(j,i),1)+0.5*sortingsNew(j,sortings(j,i),2),txt,'HorizontalAlignment','center','Color',[ 1 1 1],'FontSize',BumpChartFont)
        end
    end
end
xlim([0.7 length(ts)+0.25])
ylim([0 max(sum(outDeg))])
box off
% axis off

xticks(1:4)
ylabel('Contribution to IP index')
xticklabels({'1972 to 1983','1984 to 1995','1996 to 2007', '2008 to 2019'});
print('-dpdf', './Graphs/BumpChart.pdf','-fillpage')

%% BumpChart diagonal
clear bumpRanks sortings NewSortings outDeg sortings
fig = figure('Name','IRV no spillovers');
fig.PaperPositionMode = 'manual';
fig.PaperOrientation = 'landscape';
fig.PaperSize = [10  7];
fig.PaperPosition = [ 0 0 1 1];
dimensionSum = 1;
k=0;
for t=ts
    k=k+1;
    weightComponent = meanWeight(:,k)/sum(meanWeight(:,k)); weightComponent = weightComponent.^2;
    weightInnovations = ones(size(weightComponent));
    
    outDeg(:,k)=sum(diag(weightComponent)'*(diag(diag(DfinalRegIRF(:,:,t))))*diag(weightInnovations),dimensionSum);
    
    [bumpRanks(k,:), sortings(k,:)] = sort(outDeg(:,k),'descend');
end
distWhite=[0 linspace(0.1,0.05,5) zeros(1,length(outDeg)-6)];

for i=1:length(outDeg)
    for j=1:length(ts)
        
        sortingsNew(j,sortings(j,i),1) = sum(bumpRanks(j,:)) - sum(bumpRanks(j,i:end) )+distWhite(i);
        
        sortingsNew(j,sortings(j,i),2) = sum(bumpRanks(j,:)) - sum(bumpRanks(j,i+1:end) ); %
    end
end

k=0;
Xs=0.5:0.025:(length(ts)+0.5);
for x=Xs
    k=k+1;
    weightings = zeros(length(ts),1);
    if x<1
        firstWeight=1;
        lastWeight=ceil(x);
    elseif x>length(ts)
        firstWeight=floor(x);
        lastWeight=length(ts);
    else
        firstWeight=floor(x);
        lastWeight=ceil(x);
    end
    weightings(firstWeight:lastWeight) = pdf(makedist('Normal','mu',x,'sigma',0.2),firstWeight:lastWeight);
    weightings = weightings/sum(weightings);
    NewSortings(k,:,1) = weightings'*sortingsNew(:,:,1);
    NewSortings(k,:,2) = weightings'*sortingsNew(:,:,2);
end
hold on

for i=1:length(outDeg)-0
    fill([Xs, fliplr(Xs)],[NewSortings(:,length(outDeg)+1-i,1)' fliplr( NewSortings(:,length(outDeg)+1-i,2)' ) ],ColorsScale(i,:),'EdgeColor',0.125*[1 1 1])
end

for i=1:3
    for j=1:length(ts)
        txt=convertStringsToChars(names{sortings(j,i)});
        text(j,0.5*sortingsNew(j,sortings(j,i),1)+0.5*sortingsNew(j,sortings(j,i),2),txt,'HorizontalAlignment','center','Color',[ 1 1 1],'FontSize',BumpChartFont)
    end
end
xlim([0.7 length(ts)+0.25])
ylim([0 max(sum(outDeg))])
box off
% axis off
xticks(1:4)
ylabel('Contribution to IP index (no spillovers)')
xticklabels({'1972 to 1983','1984 to 1995','1996 to 2007','2008 to 2019'});
print('-dpdf', './Graphs/BumpChartDiag.pdf','-fillpage')
%% BumpChart off diag
clear bumpRanks sortings NewSortings outDeg sortings
fig = figure('Name','IRV only spillovers');
fig.PaperPositionMode = 'manual';
fig.PaperOrientation = 'landscape';
fig.PaperSize = [10  7];
fig.PaperPosition = [ 0 0 1 1];
dimensionSum = 1;
k=0;
for t=ts
    k=k+1;
    weightComponent = meanWeight(:,k)/sum(meanWeight(:,k));weightComponent = weightComponent.^2;
    weightInnovations = ones(size(weightComponent));
    
    outDeg(:,k)=sum(diag(weightComponent)'*(DfinalRegIRF(:,:,t)-diag(diag(DfinalRegIRF(:,:,t))))*diag(weightInnovations),dimensionSum);
    
    [bumpRanks(k,:), sortings(k,:)] = sort(outDeg(:,k),'descend');
    
end

distWhite=[0 linspace(0.1,0.05,5) zeros(1,length(outDeg)-6)];
for i=1:length(outDeg)
    for j=1:length(ts)
        
        sortingsNew(j,sortings(j,i),1) = sum(bumpRanks(j,:)) - sum(bumpRanks(j,i:end) )+distWhite(i);
        sortingsNew(j,sortings(j,i),2) = sum(bumpRanks(j,:)) - sum(bumpRanks(j,i+1:end) ); %
    end
end

k=0;
Xs=0.5:0.025:(length(ts)+0.5);
for x=Xs
    k=k+1;
    weightings = zeros(length(ts),1);
    if x<1
        firstWeight=1;
        lastWeight=ceil(x);
    elseif x>length(ts)
        firstWeight=floor(x);
        lastWeight=length(ts);
    else
        firstWeight=floor(x);
        lastWeight=ceil(x);
    end
    weightings(firstWeight:lastWeight) = pdf(makedist('Normal','mu',x,'sigma',0.2),firstWeight:lastWeight);
    weightings = weightings/sum(weightings);
    NewSortings(k,:,1) = weightings'*sortingsNew(:,:,1);
    NewSortings(k,:,2) = weightings'*sortingsNew(:,:,2);
end
hold on

for i=1:length(outDeg)-0
    fill([Xs, fliplr(Xs)],[NewSortings(:,length(outDeg)+1-i,1)' fliplr( NewSortings(:,length(outDeg)+1-i,2)' ) ],ColorsScale(i,:),'EdgeColor',0.125*[1 1 1])
end

for i=1:3
    for j=1:length(ts)
        txt=convertStringsToChars(names{sortings(j,i)});
        text(j,0.5*sortingsNew(j,sortings(j,i),1)+0.5*sortingsNew(j,sortings(j,i),2),txt,'HorizontalAlignment','center','Color',[ 1 1 1],'FontSize',BumpChartFont)
    end
end
xlim([0.7 length(ts)+0.25])
ylim([0 max(sum(outDeg))])
box off
% axis off
xticks(1:4)
ylabel('Contribution to IP index (spillover only)')
xticklabels({'1972 to 1983','1984 to 1995','1996 to 2007','2008 to 2019'});
print('-dpdf', './Graphs/BumpChartOffDiag.pdf','-fillpage')

%% CCDF
fig = figure('Name','CCDF');
fig.PaperPositionMode = 'manual';
fig.PaperOrientation = 'portrait';
fig.PaperSize = [ 7  5 ];
fig.PaperPosition = [ 0 0 1 1];

k=0;
lines={{'-'},{':'},{'--'},{'-.'}};

for t=ts
    k=k+1;
    weightComponent = N*meanWeight(:,k)/sum(meanWeight(:,k));
    Table=round( diag(weightComponent)*(DfinalReg(:,:,t)-diag(diag(DfinalReg(:,:,t))))*diag(weightComponent) ,3);
    Table=DfinalReg(:,:,t)-diag(diag(DfinalReg(:,:,t)));
    deg=sum(Table,2)'+eps;
    deg2=sum(Table,1)+eps;
    
    pd=fitdist(deg2','kernel','Kernel','normal','width',0.05,'support','positive');
    x = [-0.02:0.01:1.7];
    y = 1-cdf(pd,x);
    plot(x,y,char(lines{k}),'LineWidth',2)
    hold on
    
end
legend('1972 to 1983','1984 to 1995','1996 to 2007','2008 to 2019')
axis([-0.05 max(x) 0 1])
ylabel('1-CDF(Out-connectedness)')
xlabel('Out-connectedness')
hold off
% box off
legend('boxoff')
print('-dpdf', './Graphs/CDFToCon.pdf','-fillpage')