%% Produce Summary Statistics & Density Plots 

clear all
set(0,'defaultAxesFontName', 'Times New Roman')                 % Set font for axes
set(0,'defaultTextFontName', 'Times New Roman')                 % Set font for text
set(0,'defaultAxesFontSize', 16)                                % Set font size for axes
set(0,'defaultTextFontSize', 16)                                % Set font size for text

delete('SummaryStats.out')
diary('SummaryStats.out')

%% 1953:2 -- 2023:4

load('YSRep.mat','Spread')

% Summary Stats
Spread1953 = [ mean(Spread) ; std(Spread) ; skewness(Spread) ; kurtosis(Spread) ] ;

% SW Test
[~, SW_S_pv, SW_S_ts ] = swtest(Spread) ;
Spread1953 = [ Spread1953 ; SW_S_ts ; SW_S_pv ] ;

clear SW_S_pv SW_S_ts Spread

%% 1971:1 - 2023:4

load('YSRep1971.mat','Spread','NFCI')

% Summary Stats
Spread1971 = [ mean(Spread) ; std(Spread) ; skewness(Spread) ; kurtosis(Spread) ] ;
NFCI1971 = [ mean(NFCI) ; std(NFCI) ; skewness(NFCI) ; kurtosis(NFCI) ] ;

% SW Tests
[~, SW_S_pv, SW_S_ts ] = swtest(Spread) ;
Spread1971 = [ Spread1971 ; SW_S_ts ; SW_S_pv ] ;
[~, SW_F_pv, SW_F_ts ] = swtest(NFCI) ;
NFCI1971 = [ NFCI1971 ; SW_F_ts ; SW_F_pv ] ;

clear SW_S_pv SW_S_ts SW_F_pv SW_F_ts Spread NFCI

%% Print Summary Statistics

stat = strvcat("Mean ","StdDev ","Skewness ","Kurtosis ","SW Stat ","SW pval") ;

disp('1953:2 - 2023:4')
disp('===============') 
disp(blanks(1))
disp('Spread')
disp('------')
disp([ stat num2str(Spread1953)])
disp(blanks(2))

disp('1971:1 - 2023:4')
disp('===============') 
disp(blanks(1))
disp('Spread')
disp('------')
disp([ stat num2str(Spread1971)])
disp(blanks(1))
disp('NFCI')
disp('------')
disp([ stat num2str(NFCI1971)])
disp(blanks(2))

clear stat Spread1953 Spread1971 NFCI1971

%% Produce Density Plots 1071:1 - 2023:4

load('YSRep1971.mat','Spread','NFCI')

% Spread 
figure('position', [100, 100, 1000, 600])  
SpreadPts = linspace(-3,6,1000) ;
SpreadEmp = ksdensity(Spread,SpreadPts) ;
SpreadNorm = mean(Spread) + std(Spread)*randn(1000000,1) ;
SpreadNor = ksdensity(SpreadNorm,SpreadPts) ;
plot(SpreadPts,SpreadEmp,'LineWidth',2,'Color','r') ;
hold 
plot(SpreadPts,SpreadNor,'LineWidth',2,'Color','b') ;
legend('Empirical','Normal')
clear SpreadPts SpreadEmp SpreadNorm SpreadNor

% NFCI
figure('position', [100, 100, 1000, 600])  
NFCIPts = linspace(-4,6,1000) ;
NFCIEmp = ksdensity(NFCI,NFCIPts) ;
NFCINorm = mean(NFCI) + std(NFCI)*randn(1000000,1) ;
NFCINor = ksdensity(NFCINorm,NFCIPts) ;
plot(NFCIPts,NFCIEmp,'LineWidth',2,'Color','r') ;
hold 
plot(NFCIPts,NFCINor,'LineWidth',2,'Color','b') ;
legend('Empirical','Normal')
clear NFCIPts NFCIEmp NFCINorm NFCINor
