

clear; clc; 
optionsIVAN = optimset('Display', 'iter-detailed','LargeScale', 'off','MaxFunEvals',10000);

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

LogGDP = xlsread([CurrentFolder '\Data\GDPC1.xls'],'GDP','c12:c289'); 
DATES = [1947+1/4:1/4:2020]';
DATES=DATES(1:size(LogGDP,1),1);

%% FIND RECESSIONS 
% Recessions	
% July 1953(II)	May 1954 (II)
% August 1957(III)	April 1958 (II)
% April 1960(II)	February 1961 (I)
% December 1969(IV)	November 1970 (IV)
% November 1973(IV)	March 1975 (I)
% January 1980(I)	July 1980 (III)
% July 1981(III)	November 1982 (IV)
% July 1990(III)	March 1991(I)
% March 2001(I)	November 2001 (IV)
% December 2007 (IV)	June 2009 (II)

 
NBER_Recessions = [ 37 41;...
                    54 57;...
                    65 68;...
                    103 107;...  
                    119 124;...
                    144 146;...
                    150 155;...
%                     144 155;...
                    186 188;...
                    228 231;...
                    255 261];
NBER_Recessions = NBER_Recessions - 11*ones(size(NBER_Recessions));                

ContractionsDurations = NBER_Recessions(:,2)-NBER_Recessions(:,1);
ExpansionsDurations = NBER_Recessions(2:end,1)-NBER_Recessions(1:end-1,2);

%% Calculate Violence 

InitalDate = 1947+2/4;%
SplittingDate = 1984+3/4;


OptionLocalVariance = 1; % (1) 5 years std prior to recession; (2) pre-post great moderation; (3) Use IGARCH; (4) TVP-AR+SV;
CaseAR = 'AR2'; % 'AR1'; % 

TABELLA_AllStdViolence =[];

for OptionLocalVariance=OptionLocalVariance; %1:4
for jjj=1:2
    if jjj==2
        %% QUARTER ON QUARTER GROWTH
        D_GDP = 4*(LogGDP(2:end,1)-LogGDP(2-1:end-1,1)); 
        DATES_gr = DATES(2:end,1);
        AnnualizeConstant = 1; 
    elseif jjj==1
        %% Year on year GROWTH
        D_GDP = LogGDP(5:end,1)-LogGDP(1:end-4,1);
        DATES_gr = DATES(5:end,1); 
        AnnualizeConstant = 1; 
    end
%     FindInitalDate = find(DATES_gr==InitalDate);
    FindInitalDate = 1;
    FindSplittingDate = find(DATES_gr==SplittingDate);
    
    %% CALCULATE VARIANCE WITH IGARCH 
    if OptionLocalVariance ==3; 
    if CaseAR=='AR1';
    %% AR(1)
    y=D_GDP(1+1:end,1);
    X=D_GDP(1:end-1,1);
    X = [ones(size(y)) X];
    Rho = X\y;
    %% ESTIMATE
    vpars_init = [-1; Rho];
    lossFun = @(vpar) AR_plusIGARCH(vpar,y,X);
    [vpars,fval,exitflag,output,grad,hessianInitN] = fminunc(lossFun,vpars_init,optionsIVAN);
    [LL, Res] = AR_plusIGARCH(vpars,y,X);
    VolatilityIGARCH = [NaN; Res.VarY(1:end-1,1).^.5];
    elseif CaseAR=='AR2';
    %% PROVA AR2
    y=D_GDP(1+2:end,1);
    X=[D_GDP(3-1:end-1,1) D_GDP(3-2:end-2,1)];
    X = [ones(size(y)) X];
    Rho = X\y;
    %% ESTIMATE
    vpars_init = [-1; Rho];
    lossFun = @(vpar) AR_plusIGARCH(vpar,y,X);
    [vpars,fval,exitflag,output,grad,hessianInitN] = fminunc(lossFun,vpars_init,optionsIVAN);
    [LL, Res] = AR_plusIGARCH(vpars,y,X);
    VolatilityIGARCH = [NaN; NaN; Res.VarY(1:end-1,1).^.5];
    end
    end
    
    if OptionLocalVariance ==4;
    %% CALCULATE VARIANCE WITH IGARCH 
    if CaseAR=='AR1';
    %% AR(1)
    y=D_GDP(1+1:end,1);
    X=D_GDP(1:end-1,1);
    X = [ones(size(y)) X];
    Rho = X\y;
    %% ESTIMATE
    vpars_init = [-1; -1; -1];
%     lossFun = @(vpar) AR_plusIGARCH(vpar,y,X);
    lossFun = @(vpar) ScoreFilterGauss2_noINITpar(vpar,...
    y,X,0,0,0,5,3,Rho,std(y(1:48,1)-X(1:48,:)*Rho));
    [vpars,fval,exitflag,output,grad,hessianInitN] = fminunc(lossFun,vpars_init,optionsIVAN);
%     [LL, Res] = AR_plusIGARCH(vpars,y,X);
    [LL, Res] = ScoreFilterGauss2_noINITpar(vpars,...
    y,X,0,0,0,5,3,Rho,std(y(1:48,1)-X(1:48,:)*Rho));
    VolatilityTVPSV = [NaN; Res.VarY(1:end-1,1).^.5];
    elseif CaseAR=='AR2';
    %% PROVA AR2
    y=D_GDP(1+2:end,1);
    X=[D_GDP(3-1:end-1,1) D_GDP(3-2:end-2,1)];
    X = [ones(size(y)) X];
    Rho = X\y;
    %% ESTIMATE
    vpars_init = [-1; -1; -1];
%     lossFun = @(vpar) AR_plusIGARCH(vpar,y,X);
    lossFun = @(vpar) ScoreFilterGauss2_noINITpar(vpar,...
    y,X,0,0,0,5,3,Rho,std(y(1:48,1)-X(1:48,:)*Rho));
    [vpars,fval,exitflag,output,grad,hessianInitN] = fminunc(lossFun,vpars_init,optionsIVAN);
%     [LL, Res] = AR_plusIGARCH(vpars,y,X);
    [LL, Res] = ScoreFilterGauss2_noINITpar(vpars,...
    y,X,0,0,0,5,3,Rho,std(y(1:48,1)-X(1:48,:)*Rho));
    VolatilityTVPSV = [NaN; NaN; Res.VarY(1:end-1,1).^.5];
    end
        
    end
    
%%

    

    for ii=1:size(NBER_Recessions,1) 
        % (b) SELECT Peak and Trough from the data 
        Peak_Date = NBER_Recessions(ii,1); 
        Trough_Date = NBER_Recessions(ii,2); 
        [Peak_Select(ii,1), bb] = max(LogGDP(Peak_Date-1:Peak_Date+1,1)); 
        PeakTrough_DateAuto(ii,1) = bb + Peak_Date - 1 - 1; 
        [Trough_Select(ii,1), cc] = min(LogGDP(Trough_Date-1:Trough_Date+1,1)); 
        PeakTrough_DateAuto(ii,2) = cc + Trough_Date - 1 - 1; 
        RecessionLENGHT(ii,1) = (PeakTrough_DateAuto(ii,2)-PeakTrough_DateAuto(ii,1)); % RECESSION LENGHT
        
        RecessionDept_2(ii,1) = Peak_Select(ii,1) - Trough_Select(ii,1);
        RecessionDept_3(ii,1) = 4*RecessionDept_2(ii,1)/RecessionLENGHT(ii,1); % ADJUSTED FOR RECESSION LENGHT
        
        if OptionLocalVariance==1;
%         LocalVariance(ii,1) = std(D_GDP(max(1,Peak_Date-5*4-1):Peak_Date-1),1); 
        LocalVariance(ii,1) = std(D_GDP(max(1,Peak_Date-4-5*4):Peak_Date-4),1); 
        elseif OptionLocalVariance==2;
            if Peak_Date<FindSplittingDate
            LocalVariance(ii,1) = std(D_GDP(FindInitalDate:FindSplittingDate-1),1); 
            else
            LocalVariance(ii,1) = std(D_GDP(FindSplittingDate:end),1); 
            end    
        elseif OptionLocalVariance==3;
        LocalVariance(ii,1) = VolatilityIGARCH(Peak_Date-1); 
%         LocalVariance(ii,1) = mean(VolatilityIGARCH(max(1,Peak_Date-4-5*4):Peak_Date-4)); 
        elseif OptionLocalVariance==4;
        LocalVariance(ii,1) = VolatilityTVPSV(Peak_Date-1); 
%         LocalVariance(ii,1) = mean(VolatilityTVPSV(max(1,Peak_Date-4-5*4):Peak_Date-4)); 
        end
        
        RecessionDept_5(ii,1) =  RecessionDept_3(ii,1)/LocalVariance(ii,1);
        
    end 

%% SELECT VIOLENCE    
Violence = RecessionDept_3;    
StandardizedViolence = RecessionDept_5;    

TABELLA_AllStdViolence = [ TABELLA_AllStdViolence...
            StandardizedViolence];

        
        

end
end

TABELLA = [ContractionsDurations [NaN; ExpansionsDurations] Violence TABELLA_AllStdViolence]; 
Pre84_Mean = nanmean(TABELLA(1:end-3,:));
Post84_Mean = nanmean(TABELLA(end-2:end,:));
TABELLA = [TABELLA; NaN*TABELLA(1,:); Pre84_Mean; Post84_Mean];
         
MainTable_Recessions = TABELLA(:,[3,4,1,2]); 
RobustnessTable_Recessions = TABELLA(:,5:end); 

%% Save and Print Results
swp2(MainTable_Recessions,[CurrentFolder '\Tables_And_Figures\TABELLA_Recessions_Main.tex']);
swp2(RobustnessTable_Recessions,[CurrentFolder '\Tables_And_Figures\TABELLA_Recessions_Robustness.tex']);




