clear
clc

WhereAreWeNow = cd;
WhereAreSimulations = 'C:\Users\dlx832\Desktop\Step 2 Simulation\Computestats';

WhereToSavePlots = [WhereAreWeNow '\FIGURES\'];
Sbar_Values = [ 67 70 75 80 85 90]'; 


LocateGDP = 1;
LocateCons = 2; 
LocateInv =3; 

CalibrateGDPgrowth = 0.031267; %Average GDP growth rate 1948-2016

for iii=1:size(Sbar_Values,1)
    NameFile = ['SimOutputLevels_yci_sbar=' num2str(Sbar_Values(iii)) '_pct.mat'];
    NameFileWithLoc = [WhereAreSimulations '\' NameFile];
    load(NameFileWithLoc);
    
    for sim = 1:size(SimOutputLevels_yci,2);
        GDPSim = squeeze(SimOutputLevels_yci(LocateGDP,sim,:));
        ConsSim = squeeze(SimOutputLevels_yci(LocateCons,sim,:));
        InvSim = squeeze(SimOutputLevels_yci(LocateInv,sim,:));
        
        %% RECESSIONS STATISTICS
        TT = size(GDPSim,1);

        TrendGrowth = (((1+CalibrateGDPgrowth)^(1/4))-1)*[1:TT]';
        GDPlev = 100*(GDPSim+TrendGrowth);
        [PT,RECESSIONS_PT,DURATIONS,VIOLENCE_TABLE,MeanViolence,EXCESS_v2]=BusinessCycleDating_v2(GDPlev);
        SAVE_RESULTS_Recessions_AllSim(sim,:)=[DURATIONS MeanViolence(2,:)];

        
        ReoraganizeExpansions = [RECESSIONS_PT(1:end-1,2) RECESSIONS_PT(2:end,1)]; 
        ReoraganizeExpansions = [ReoraganizeExpansions ceil(mean(ReoraganizeExpansions,2))];
        for jjvv = 1:size(ReoraganizeExpansions,1)
            SelectTP = ReoraganizeExpansions(jjvv,:);
            AVERAGEGrowth_EXP(jjvv,:)= 4*((GDPlev(SelectTP(2))-GDPlev(SelectTP(1)))/(SelectTP(2)-SelectTP(1)));
        end
        SAVE_AVERAGEGrowthEXP_AllSim(sim,1)=[mean(AVERAGEGrowth_EXP)];        

        
        for jjvv = 1:size(RECESSIONS_PT,1)
        GDPlevPeak = GDPlev(RECESSIONS_PT(jjvv,1));
        GDPlevPostPeak = GDPlev(RECESSIONS_PT(jjvv,2)+1:end);
        [~,b] = max((GDPlevPostPeak>=GDPlevPeak));
        TIMETOPEAK(jjvv,1) = b; 
        end
        
        
        %% SKEWNESS STATISTICS
        Variables = [GDPSim ConsSim InvSim]; 
        VariablesYoYGr = Variables(4+1:end,:)-Variables(1:end-4,:);
        SAVE_RESULTS_Skewness_AllSim(sim,:) = skewness(VariablesYoYGr);
        
        conta=1; for jj=1:3; correl(conta)=corr(VariablesYoYGr(:,1),VariablesYoYGr(:,jj)); conta=conta+1; end 
        
        %% Volatility STATISTICS
        GDP_YoYGr = [GDPSim(4+1:end,:)-GDPSim(1:end-4,:)]; 
        GDP_level = [GDPSim(4+1:end,:)];
        VariablesForSTD = NaN*zeros(size(GDP_YoYGr,1),3);
        VariablesForSTD(:,1) = GDP_YoYGr;     
        VariablesForSTD((GDP_level>=0),2) = GDP_YoYGr((GDP_level>=0));     
        VariablesForSTD((GDP_level<0),3) = GDP_YoYGr((GDP_level<0));     
        
        SAVE_RESULTS_STD_AllSim(sim,:) = nanstd(VariablesForSTD);

        
        %% SemiVolatility STATISTICS
        GDP_YoYGr = [GDPSim(4+1:end,:)-GDPSim(1:end-4,:)]; 
        SemiStd_All = sqrt(sum((GDP_YoYGr-mean(GDP_YoYGr)).^2)/(size(GDP_YoYGr,1)-1));     
        SemiStd_Pos = sqrt(sum((GDP_YoYGr(GDP_YoYGr>=mean(GDP_YoYGr))-mean(GDP_YoYGr)).^2)/(size(GDP_YoYGr,1)-1));     
        SemiStd_Neg = sqrt(sum((GDP_YoYGr(GDP_YoYGr<mean(GDP_YoYGr))-mean(GDP_YoYGr)).^2)/(size(GDP_YoYGr,1)-1));     
        
        SAVE_RESULTS_SemiSTD_AllSim(sim,:) = [SemiStd_All SemiStd_Pos SemiStd_Neg];
               
        
    end
    SAVE_RESULTS_Recessions_AllSim_AllSbar(:,:,iii) = SAVE_RESULTS_Recessions_AllSim;
    SAVE_RESULTS_Skewness_AllSim_AllSbar(:,:,iii) = SAVE_RESULTS_Skewness_AllSim;
    SAVE_RESULTS_STD_AllSim_AllSbar(:,:,iii) = SAVE_RESULTS_STD_AllSim;
    SAVE_AVERAGEGrowthEXP_AllSim_AllSim(:,iii) = SAVE_AVERAGEGrowthEXP_AllSim;
    SAVE_RESULTS_SemiSTD_AllSim_AllSim(:,:,iii) = SAVE_RESULTS_SemiSTD_AllSim; 
    
    SAVE_RESULTS_Recessions_Median(:,iii) = median(SAVE_RESULTS_Recessions_AllSim,1)';
    SAVE_RESULTS_Skewness_Median(:,iii) = median(SAVE_RESULTS_Skewness_AllSim,1)';
    SAVE_RESULTS_STD_Median(:,iii) = median(SAVE_RESULTS_STD_AllSim,1)';


    SAVE_AVERAGEGrowthEXP(iii,1) = median(SAVE_AVERAGEGrowthEXP_AllSim,1)';
    SAVE_RESULTS_SemiSTD_Median(:,iii) = median(SAVE_RESULTS_SemiSTD_AllSim,1)'; 

    clear SAVE_RESULTS_Recessions_AllSim SAVE_RESULTS_Skewness_AllSim SAVE_RESULTS_STD_AllSim...
        SAVE_EXCESS_AllSim SAVE_AVERAGEGrowthEXP_AllSim
end


%% PLOT FIGURS

Sbar_Values = Sbar_Values/100;


SPESSORE =2.5; 
FontSIZE = 14;



figure(71)
plot1 = plot(Sbar_Values,SAVE_RESULTS_Skewness_Median(1:3,:)',...
    'linewidth',SPESSORE,'color','b'); 
%set(plot1(2),'color','g','linestyle','--')
set(plot1(2),'color',[0.54 0.81 0.94],'linestyle','--')
%set(plot1(3),'color','m','linestyle',':')
set(plot1(3),'color',[0.67 0.88 0.69],'linestyle',':')
xlim([min(Sbar_Values) max(Sbar_Values)]);
xlabel('${s}^{I}$','FontSize',12,'interpreter','Latex');
grid on;
legend1 = legend('Output','Consumption','Investment');
set(legend1,'EdgeColor',[1 1 1],'Location','southwest','FontSize',10);
title('Skewness','fontsize',FontSIZE);
% tightfig;
saveas(gcf,[WhereToSavePlots 'AllVariablesSkewness.pdf']);
saveas(gcf,[WhereToSavePlots 'AllVariablesSkewness.tif']);

figure(83)
plot1 = plot(Sbar_Values,SAVE_RESULTS_Recessions_Median(4,:),...
    'linewidth',SPESSORE,'color','b'); 
xlim([min(Sbar_Values) max(Sbar_Values)]);
xlabel('${s}^{I}$','FontSize',12,'interpreter','Latex');
grid on;
title('(Standardized) Violence','fontsize',FontSIZE);
% tightfig;
saveas(gcf,[WhereToSavePlots 'SimRec_Violence.pdf']);
saveas(gcf,[WhereToSavePlots 'SimRec_Violence.tif']);

figure(81)
plot1 = plot(Sbar_Values,SAVE_RESULTS_Recessions_Median(2,:),...
    'linewidth',SPESSORE,'color','b'); 
xlim([min(Sbar_Values) max(Sbar_Values)]);
xlabel('${s}^{I}$','FontSize',12,'interpreter','Latex');
grid on;
title('Expansions','fontsize',FontSIZE);
% tightfig;
saveas(gcf,[WhereToSavePlots 'SimExpDuration.pdf']);
saveas(gcf,[WhereToSavePlots 'SimExpDuration.tif']);

figure(82)
plot1 = plot(Sbar_Values,SAVE_RESULTS_Recessions_Median(1,:),...
    'linewidth',SPESSORE,'color','b'); 
xlim([min(Sbar_Values) max(Sbar_Values)]);
xlabel('${s}^{I}$','FontSize',12,'interpreter','Latex');
ylim([2 8]); 
%set(gca,'YTick',2:0.5:4);

grid on;
title('Contractions','fontsize',FontSIZE);
% tightfig;
saveas(gcf,[WhereToSavePlots 'SimRecDuration.pdf']);
saveas(gcf,[WhereToSavePlots 'SimRecDuration.tif']);

figure(91)
plot1 = plot(Sbar_Values,SAVE_RESULTS_STD_Median(1,:),...
    'linewidth',SPESSORE,'color','b'); 
xlim([min(Sbar_Values) max(Sbar_Values)]);
xlabel('${s}^{I}$','FontSize',12,'interpreter','Latex');
grid on;
title('Business Cycle Volatility','fontsize',FontSIZE);
% tightfig;
saveas(gcf,[WhereToSavePlots 'SimGDP_Volatility.pdf']);
saveas(gcf,[WhereToSavePlots 'SimGDP_Volatility.tif']);



figure(92)
plot1 = plot(Sbar_Values,SAVE_RESULTS_STD_Median(2:3,:)',...
    'linewidth',SPESSORE,'color','b'); 
set(plot1(2),'color','b','linestyle','--')
xlim([min(Sbar_Values) max(Sbar_Values)]);
xlabel('${s}^{I}$','FontSize',12,'interpreter','Latex');
grid on;
legend1 = legend('Booms','Busts');
set(legend1,'EdgeColor',[1 1 1],'Location','southwest','FontSize',10);
title('Volatility of Booms and Busts','fontsize',FontSIZE);
% tightfig;
saveas(gcf,[WhereToSavePlots 'SimGDP_VolatilityExpCont.pdf']);
saveas(gcf,[WhereToSavePlots 'SimGDP_VolatilityExpCont.tif']);



SAVE_RESULTS_SemiSTD_Median(4,:)=SAVE_RESULTS_SemiSTD_Median(3,:)./SAVE_RESULTS_SemiSTD_Median(2,:);
figure(72)
plot1 = plot(Sbar_Values,SAVE_RESULTS_SemiSTD_Median(4,:)',...
    'linewidth',SPESSORE,'color','b'); 
%set(plot1(2),'color','b','linestyle','--')
xlim([min(Sbar_Values) max(Sbar_Values)]);
grid on;
legend1 = legend('Ratio');
set(legend1,'EdgeColor',[1 1 1],'Location','southwest');
title('Ratio of Negative to Positive Semivariance','fontsize',FontSIZE);
% tightfig;
saveas(gcf,[WhereToSavePlots 'SimGDP_VolatilityPosNegSemiRatio.pdf']);
saveas(gcf,[WhereToSavePlots 'SimGDP_VolatilityPosNegSemiRatio.tif']);



%save('simulresults_all','Sbar_Values','SAVE_RESULTS_Skewness_Median','SAVE_RESULTS_Recessions_Median','SAVE_RESULTS_STD_Median','SAVE_RESULTS_SemiSTD_Median');

