%===============================
%Load Data – Summary Statistics
%===============================
clear all
load CSFBCompEx_new.txt
data=CSFBCompEx_new;

load datafanalyse_new.txt
datafactors=datafanalyse_new;

[T,plithosmet]=size(data)
[T,plithosfactors]=size(datafactors)
%===============================

%===============================
%plot of the data
plot(data)
pause
hist(data)
pause
%calculate summary statistics
[mean(data) std(data) median(data) kurtosis(data) skewness(data)]
%===============================

%===============================
%===============================
%RUN MULTIPLE REGRESSION MODEL
%===============================
%================================================
% 1. Multiple Regression : Applying the formulas
%================================================
model=[1 2 3 4]; numfactors=length(model);
X=[]; X=[ones(T,1) datafactors(:,model)];
b=(inv(X'*X))*X'*data;
yhat=X*b;
errors=data-yhat;
s2=(1/(T-numfactors-1))*sum(errors.^2); s=sqrt(s2);
loglikelihood=-0.5*T*log(2*pi)-0.5*T*log(s2)-0.5*(1/s2)*sum((data-X*b).^2);
numberparam=1+numfactors+1; 
AIC=-2*loglikelihood+2*(numberparam);  
BIC=-2*loglikelihood+(log(T)*numberparam);

% Residual tests
% 1. Autocorrelation of Residuals
autocorr(errors,12)
[H,P,Qstat,CV]=lbqtest(errors,12)

% 2. Autocorrelation of Squared residuals
errors2=errors.^2;
autocorr(errors2,12)
[H2,P2,Qstat2,CV2]=lbqtest(errors2,12)

% 3. Normality test
alpha=0.05;
[H,P,JBstat,CVN]=jbtest(errors,alpha)

%================================================
% MODEL SELECTION : Via STEPWISEFIT Routine
%================================================
[B,SE,PVAL,INMODEL,STATS,NEXTSTEP,HISTORY]=stepwisefit(datafactors,data,'display','off');


%================================================
% 2. Multiple Regression : Via GARCHFIT Routine
%================================================
% NORMAL ERRORS
speck=garchset('P',0,'Q',0,'Regress',[1:numfactors],'display','off');
Coeff=[]; Errors=[]; LLF=[]; Innovations=[]; Sigmas=[]; Summary=[];
[Coeff,Errors,LLF,Innovations,Sigmas,Summary]=garchfit(speck,data,datafactors(:,model)); 
%Keep the parameters
param(1)=garchget(Coeff,'C');
param((1+1):(1+numfactors))=garchget(Coeff,'Regress');
param(numfactors+1+1)=garchget(Coeff,'K');
param=param';

% STUDENT-T ERRORS
speckT=garchset('P',0,'Q',0,'Regress',[1:numfactors],'Distribution','T','display','off');
CoeffT=[]; ErrorsT=[]; LLFT=[]; InnovationsT=[]; SigmasT=[]; SummaryT=[];
[CoeffT,ErrorsT,LLFT,InnovationsT,SigmasT,SummaryT]=garchfit(speckT,data,datafactors(:,model)); 


%===========================================================
% 3. Multiple Regression - GARCH Model: Via GARCHFIT Routine
%===========================================================
% NORMAL ERRORS
speckGN=garchset('P',1,'Q',1,'Regress',[1:numfactors],'display','off');
CoeffGN=[]; ErrorsGN=[]; LLFGN=[]; InnovationsGN=[]; SigmasGN=[]; SummaryGN=[];
[CoeffGN,ErrorsGN,LLFGN,InnovationsGN,SigmasGN,SummaryGN]=garchfit(speckGN,data,datafactors(:,model)); 
subplot(2,1,1)
plot(data)
subplot(2,1,2)
plot(SigmasGN)
autocorr(InnovationsGN,24)
[Hr,Pr,Qr,CVr]=lbqtest(InnovationsGN,24)
autocorr(abs(InnovationsGN),24)
[Har,Par,Qar,CVar]=lbqtest(abs(InnovationsGN),24)

% STUDENT-T ERRORS
speckGT=garchset('P',1,'Q',1,'Regress',[1:numfactors],'Distribution','T','display','off');
CoeffGT=[]; ErrorsGT=[]; LLFGT=[]; InnovationsGT=[]; SigmasGT=[]; SummaryGT=[];
[CoeffGT,ErrorsGT,LLFGT,InnovationsGT,SigmasGT,SummaryGT]=garchfit(speckGT,data,datafactors(:,model)); 
subplot(2,1,1)
plot(data)
subplot(2,1,2)
plot(SigmasGT)
autocorr(InnovationsGT,24)
[Hr,Pr,Qr,CVr]=lbqtest(InnovationsGT,24)
autocorr(abs(InnovationsGT),24)
[Har,Par,Qar,CVar]=lbqtest(abs(InnovationsGT),24)

