clear; clc; 

WhereAmI = cd; 
WhereToSavePlots = [WhereAmI '\FiguresAndTables\'];

SPESSORE =2.5;

CurrentFolder = cd;
path(path,[CurrentFolder '\Data']);
path(path,[CurrentFolder '\Subroutines']);

NomDebt_PP = xlsread('NY FED STATE LEVEL DATA ON DEBT.xls','total','b2:n53');
Population = xlsread('NY FED STATE LEVEL DATA ON DEBT.xls','population','b2:n53');
ReOrder = xlsread('NY FED STATE LEVEL DATA ON DEBT.xls','Reorder','e1:e51');
load StatesNames

NomDebt_PP = NomDebt_PP(ReOrder,:);
Population = Population(ReOrder,:);

NomDebt_tot = NomDebt_PP.*Population; 

NomIncome = xlsread('BEA STATE LEVEL DATA ON INCOME_ANNUAL','Sheet0','BY8:CK58');
NomIncome_PerCapita = xlsread('BEA STATE LEVEL DATA ON INCOME_ANNUAL_PerCapita.xls','Sheet0','BY7:CK57');

Debt_to_Income = NomDebt_PP./NomIncome_PerCapita;
Mean_Debt_to_Income = nanmean(Debt_to_Income,2);

WhichYEarPreRecession =4;
PreRecession_Debt_to_Income = Debt_to_Income(:,WhichYEarPreRecession);


%% 

LogGDP = log(xlsread('BEA STATE LEVEL DATA REAL GDP.xls','Sheet0','C8:AU58')); 
LogGDP_US = log(xlsread('BEA STATE LEVEL DATA REAL GDP.xls','Sheet0','C7:AU7'))'; 
Dates_LogGDP = x2mdate(xlsread('BEA STATE LEVEL DATA REAL GDP.xls','Sheet0','C5:AU5')); 
Dyoy_GDP = LogGDP(:,5:end)-LogGDP(:,5-4:end-4);

% 1 - Calculate Skewness 
for ii=1:size(LogGDP,1) 
    SkewLogGDP(ii,:) = skewness(Dyoy_GDP(ii,:)');
    STDLogGDP(ii,:) = std(Dyoy_GDP(ii,:)');
end


% 2 - Recessions Dept 
% (a) USE NBER DATES 
Peak_Date = 12; 
Trough_Date = 18;

RecessionDept_1 = LogGDP(:,Peak_Date)-LogGDP(:,Trough_Date);
RecessionDept_1b = LogGDP(:,Trough_Date)-LogGDP(:,Peak_Date);

% (b) SELECT Peak and Trough from the data 
for ii=1:size(LogGDP,1) 
    [Peak_Select(ii,1), bb] = max(LogGDP(ii,Peak_Date-5:Peak_Date+1)); 
    PeakTrough_DateAuto(ii,1) = bb + Peak_Date - 5 - 1; 
    [Trough_Select(ii,1), cc] = min(LogGDP(ii,Trough_Date-1:Trough_Date+5)); 
    PeakTrough_DateAuto(ii,2) = cc + Trough_Date - 1 - 1; 
    RecessionDept_2(ii,1) = Peak_Select(ii,1) - Trough_Select(ii,1);
    RecessionLENGHT(ii,1) = (PeakTrough_DateAuto(ii,2)-PeakTrough_DateAuto(ii,1)); % RECESSION LENGHT
    RecessionDept_3(ii,1) = 4*(Peak_Select(ii,1) - Trough_Select(ii,1))/RecessionLENGHT(ii,1); % ADJUSTED FOR RECESSION LENGHT
    RecessionDept_4(ii,1) = RecessionDept_3(ii,1)/STDLogGDP(ii,:); % ACCOUNTING FOR DIFFERENT VOLATILITIES 
end
 

%% CORRELAZIONE CON SKEWNESS

NAME_Y_AXIS{1} = 'Skewness (2005 - 2016)';
NAME_Y_AXIS{2} = 'Great Recession Violence';

NAMEFIGURESAVE{1}='Plot_ScatterSkewness';
NAMEFIGURESAVE{2}='Plot_ScatterViolence';

for jjj=1:2
    if jjj==1
        SelectX = PreRecession_Debt_to_Income;
        SelectY = SkewLogGDP;
    elseif jjj==2
        SelectX = PreRecession_Debt_to_Income;
        SelectY = RecessionDept_3;
%         SelectY = RecessionDept_2;
    end


xdata=SelectX;
ydata=SelectY;
%% Run regression to check significance of trendline
StatsREG = regstats2(ydata,xdata,'linear','all');
RegrREG=[StatsREG.beta(2) StatsREG.hc0.se(2) StatsREG.hc0.t(2) StatsREG.hc0.pval(2)]; % [slope, standard error, t-statistic, p-value]

%% Outlier Robust Regression
CUToffPerc = .95; CUToff = norminv(CUToffPerc);
outliers=abs(StatsREG.standres)>CUToff;
StatsOUT = regstats2(ydata(outliers==0,:),xdata(outliers==0,:),'linear','all');
RegrOUT=[StatsOUT.beta(2) StatsOUT.hc0.se(2) StatsOUT.hc0.t(2) StatsOUT.hc0.pval(2)]; % [slope, standard error, t-statistic, p-value]

%%
B=linspace(0.6,2,length(SelectX)); B=B';
B(:,2)=StatsREG.beta(2)*B(:,1)+StatsREG.beta(1);


figure(10+jjj); 
% subplot(2,2,1);
plot(SelectX,SelectY,'linestyle','none','marker','.','Color','b','markersize',8); 
labelpoints (SelectX,SelectY, StatesNames,'Color','b');    
% dx = 0.01; dy = 0.01; % displacement so the text does not overlay the data points
% textfit(SelectX+dx, SelectY+dy, StatesNames);
% text(SelectX+dx, SelectY+dy, StatesNames);
ylabel(NAME_Y_AXIS{jjj});
xlabel('Debt-to-Income Ratio');
hold on
plot(B(:,1),B(:,2),'g','linewidth',SPESSORE);
ylabel(NAME_Y_AXIS{jjj},'FontSize',11);
xlabel('Debt/Income','FontSize',11);
if jjj==1; dim = [.63 .75 .25 .15]; elseif jjj==2; dim = [.6 .15 .35 .15]; end 
str = ['Pr(H_0: b  = 0) = ' num2str(round(StatsREG.hc0.pval(2)*100,1)) '%',...
    sprintf('\n'),'Pr(H_0: b_{r} = 0) = ' num2str(round(StatsOUT.hc0.pval(2)*100,1)) '%'];
annotation('textbox',dim,'String',str,'FitBoxToText','on',...
     'FontSize',11,'EdgeColor',[1 1 1]);
 grid on
tightfig
% saveas(gcf,[WhereToSavePlots [NAMEFIGURESAVE{jjj} '.pdf']]);
% saveas(gcf,[WhereToSavePlots [NAMEFIGURESAVE{jjj} '.tif']]);


end

%% 
%% INDEX SYNTHETIC GDP

NAMEFIGURESAVEINDEX{1}='Syntetic_Indexes_10States';
NAMEFIGURESAVEINDEX{2}='Syntetic_Indexes_15States';

[~,IndexIndebtness] = sort(PreRecession_Debt_to_Income);
% [~,IndexIndebtness] = sort(Mean_Debt_to_Income);
LogGDP_Reordered = LogGDP(IndexIndebtness,:)';
[T,N] = size(LogGDP_Reordered);

% 2 - Recessions Dept 
% (a) USE NBER DATES 
Peak_Date = 12; 
Trough_Date = 18;

FontSIZE = 15;
for jjj=2
    if jjj==1
        QuantiStates = 10;
    elseif jjj==2
        QuantiStates = 15;
    end
    
    LogGDP_Reordered = LogGDP_Reordered - ones(T,1)*LogGDP_Reordered(Peak_Date,:);
    Quantiles = [.5];
    LowDebtStates_GDP = median(LogGDP_Reordered(:,1:QuantiStates),2);
    HighDebtStates_GDP = median(LogGDP_Reordered(:,end-QuantiStates:end),2);
    LogGDP_US_toPlot = LogGDP_US - ones(T,1)*LogGDP_US(Peak_Date,1);

    Dates_RGDP = [2005+1/4:1/4:2016+1/4]';

    start_plot = 1; 
    end_plot = 27; 


    dates = Dates_LogGDP';

    figure
    % subplot(2,2,3:4);
    plot1 = plot(dates(start_plot:end_plot,1),100*[LowDebtStates_GDP(start_plot:end_plot,:) HighDebtStates_GDP(start_plot:end_plot,:)],'linewidth',2.5);
    set(plot1(1),'Color','g','linestyle','--');
    set(plot1(2),'Color','b','linestyle','-');
    set(gca,'Xtick',dates(start_plot:end_plot,1)); datetick('x','yyyy','keepticks'); 
    recessionplot; datetick('x','yyyy','keeplimits');
    legend1 = legend('Low Debt','High Debt');
    set(legend1,'Edgecolor',[1 1 1],'location','southeast','FontSize',FontSIZE-2);
    legend('boxoff');
    xlim([min(dates(start_plot:end_plot,1)) max(dates(start_plot:end_plot,1))]);
    %title('Time Varying Constant');
    grid('on');
    tightfig(gcf);
    saveas(gcf,[WhereToSavePlots [NAMEFIGURESAVEINDEX{jjj} '.pdf']]);
    saveas(gcf,[WhereToSavePlots [NAMEFIGURESAVEINDEX{jjj} '.tif']]);

end


