% Construct GARCH Filtered Standardized Daily Returns
% Save Daily Returns and Standardized Returns

clear all;
small = 1.0e-10;
pinv_tol = 1.0e-05;
big = 1.0e+8;
global datadir;

% -- File Directories   
outdir = 'out/';
figdir = 'fig/';
matdir = 'mat/';
datadir = '../Data/Returns/';

% Read in Data
str_data = [datadir 'sp_500_Daily.xlsx'];
ret_date = readmatrix(str_data,'Range','A10:A25557');
ret = readmatrix(str_data,'Range','B10:B25557');

ret_save = ret;
ret_date_save = ret_date;

% Find years
yr = floor(ret_date/10000);
fy = 1926;
ii = yr >= fy;
yr = yr(ii==1);
ret = ret(ii==1);
ret_date = ret_date(ii==1);

% Estimate GARCH(1,1) model
T = length(ret);
cal = linspace(fy,2023,T)';
Mdl = garch('GARCHLags',1,'ARCHLags',1,'Offset',NaN);
EstMdl = estimate(Mdl,ret);
V=infer(EstMdl,ret);
std_resid = sqrt(V);
ret_std = (ret-mean(ret))./std_resid;  

% Plot Raw Returns and Standardized Returns
figure;
subplot(3,1,1)
 ax_fs = 10;
 plot(cal,ret);
 xlim([1920,2030]);
 xlabel('Year');
 ax = gca;
 ax.FontSize = ax_fs;
 title('Daily Returns');
subplot(3,1,2)
 plot(cal,std_resid);
 xlim([1920,2030]);
 %% ylim([-12 10]);
 xlabel('Year');
 ax = gca;
 ax.FontSize = ax_fs;
 title('GARCH(1,1) Standard Deviation');
subplot(3,1,3)
 plot(cal,ret_std);
 xlim([1920,2030]);
 ylim([-12 10]);
 xlabel('Year');
 ax = gca;
 ax.FontSize = ax_fs;
 title('Daily Standardized Returns');

 % Save Year, Month, return and standardized return
 yr = floor(ret_date/10000);
 mth = floor((ret_date-yr*10000)/100);
 fname = [matdir 'SP500_VW_Daily_StandardizedReturns_1926_2022.mat'];
 save(fname,'yr','mth','ret','ret_std','ret_date','cal');
 
 