数据处理一共可以分为三个方面，一是数据的回归分类，而是时间序列数据，三是网络型数据处理。本文将要来讨论一下时间序列的应用。

一．ARMA 模型

Arma是用来讨论时间序列里面回报率的情况，假设t时刻的回报率与t时刻之前的回报率有关。同时，也与之前的误差有关。

这模型就是AR模型和MA模型的结合，非常好理解。我们在matlab中画出序列的ACF图和PACF图来找出具有相关性的时间节点。

（以ANZ公司为例）

超出两条线的部分就是具有相关性的部分，我们在这里就取ARMA（1,1）然后我们做一个LB test，为的是确定之前的return确实与t时刻具有相关性。

求出ARMA的表达式即可。

Matlab 实现：

% Plot  series and ACF

figure;subplot(3,1,1);plot(ret_is(:,1));title('ANZ returns');xlim([0,length(ret_is)]);

subplot(3,1,2);autocorr(ret_is(:,1), 25);subplot(3,1,3);parcorr(ret_is(:,1), 20)

% LB test for AR effects on ANZ

[H5, pval5, Qs5, CV5] = lbqtest(ret_is(:,1), 5, 0.05);

[H10, pval10, Qs10, CV10] = lbqtest(ret_is(:,1), 10, 0.05);

[pval5 pval10]

% I choose ARMA(1,1) for ANZ, now fit it and forecast

Mdl=arima(1,0,1);

[EstMdl,EstParamCov,logL,info] = estimate(Mdl,ret_is(:,1));

[ret_f7(1), FMSE] = forecast(EstMdl,1,'Y0',ret_is(:,1));% 1-period forecasts

二．ARCH 模型

Arch模型与之前其他模型最大的不同是他是用来预测方差的。

关于阿尔法的取值，是由SIC和AIC来决定的。

SIC或AIC计算出最小的误差值的点，因此这里可以取阿尔法等于7。

同理，以ANZ为例计算出表达式出来:

算出计算结果后，我们要做两个验证。因为这个模型建立在一个假设条件上，即残差是iid和服从正态分布的（若不是正态分布就要服从t分布）。因此，我们要做两个test，一个是LB test.

LB test 用来检验残差是否服从iid  (H0:p1=p2=p3=^=Pn=0   H1:至少有一个不等于0)， 当p-value 越大，说明原假设成立，这说明是iid的。

JB test 用来检测残差是否服从正态性（H0：是正态分布的 H1：不是）当P-value越大越好，说明原假设成立，是正态分布的。

而当我们选取残差是t分布的时候，最后做检验的时候也是将残差的分布转化为正态分布来做检验。

Matlab的实现：

LLF=0;aic1=0;sic1=0;b=0.05

for p=1:20

Mdl = garch(0,p);Mdl.Offset=NaN;

[EstMdl,EstParamCov,logL,info] = estimate(Mdl,ANZ,'display','off');

% Adding 'display','off' to the previous command suppresses the estimation

% results from appearing in the output

aic1(p)=-2*logL+2*p;sic1(p)=-2*logL+log(length(ANZ))*p;

end

figure;plot(aic1,'b+-');title('AIC & SIC for ARCH models');hold on;plot(sic1,'r+-');

legend('AIC','SIC');

% Fit ARCH(7) model

Mdl = garch(0,7);Mdl.Offset=NaN;

EstMdl = estimate(Mdl,ANZ);

v7=infer(EstMdl,ANZ);s7=sqrt(v7); %infer the conditional variance and calculate standard deviations

a7 = ANZ-EstMdl.Offset;             %calculate innovations

% Plot the innovations, conditional standard deviations and log returns

% above one another

figure;

subplot(3,1,1);

plot(ANZdates(2:end,1),a7);

%xlim([ANZdates(2) ANZdates(end)]);              % set range of x-axis

title('ARCH(7) Innovations')

subplot(3,1,2);

plot(ANZdates(2:end,1),s7);

ylim([0,0.05]);

%xlim([BHPdates(2) BHPdates(end)]);              % set range of x-axis

title('ARCH(7) Conditional Standard Deviations')

subplot(3,1,3);

plot(ANZdates(2:end,1),ANZ);

%xlim([BHPdates(2) BHPdates(end)]);              % set range of x-axis

title('ANZ Log Returns')

figure;plot(ANZdates(2:end,1),ANZ,'c');hold on;

plot(ANZdates(2:end,1),s7)

%xlim([BHPdates(2) BHPdates(end)]);              % set range of x-axis

title('Log Returns and ARCH(7) Conditional Standard Deviations')

legend('ANZ returns', 'conditional standard deviations','location','South','Orientation','horizontal' );

a1A=EstMdl.ARCH; % stores the 15 ARCH coefficients for later use

%% Q1(b) Assess the fit of the ARCH model

e7=a7./s7;

figure;subplot(2,1,1);plot(ANZdates(2:end,1),e7);

title('ARCH(7) Standardised Residuals');

subplot(2,1,2);autocorr(e7);

title('ACF of ARCH(7) Standardised Residuals');

figure;subplot(2,1,1);hist(e7,25)

title('Histogram of ARCH(7) Standardised Residuals');

subplot(2,1,2);qqplot(e7);

title('QQ plot ARCH(7) Standardised Residuals');

figure;autocorr(e7.^2);

title('ACF of  ARCH(7) Squared Standardised Residuals');

% LB test on standardised residuals

[H, pValue, Qstat, CriticalValue] = lbqtest(e7, [20 25], 0.05, [13 18])

%LB test on squared standardised residuals

[H, pValue, Qstat, CriticalValue] = lbqtest(e7.^2, [20 25], 0.05, [13 18])

% JB test

[skewness(e7) kurtosis(e7)]

[h,p] = jbtest(e7)

三． GARCH 模型

类似于ARCH 模型，但是方差多加了一项，即t时刻之前的方差影响。

同理，先由sic和aic来确定P，q的值，然后开始建模和检验。

Matlab的实现：

aic0=0;sic0=0;aic1=0;sic1=0;aict0=0;sict0=0;aict1=0;sict1=0;

LLF=0;

for p=1:5

for q=1:5

% GARCH(p,q) with Gaussian errors

Mdl = garch(p,q);Mdl.Offset=NaN;Mdl.Distribution='Gaussian'; % GARCH(p,q) with Gaussian errors

[EstMdl,EstParamCov,LLF,info] = estimate(Mdl,ANZris,'display','off');

aic0(p,q)=-2*LLF+2*(p+q);sic0(p,q)=-2*LLF+log(nis)*(p+q);

% AR(1)-GARCH(p,q) with Gaussian errors

Mdl = arima('ARLags',1,'Variance',garch(p,q));

Mdl.Distribution='Gaussian';

[EstMdl,EstParamCov,LLF,info] = estimate(Mdl,ANZris,'display','off');

aic1(p,q)=-2*LLF+2*(p+q+1);sic1(p,q)=-2*LLF+log(nis)*(p+q+1);

% GARCH(p,q) with Student t errors

Mdl = garch(p,q);Mdl.Offset=NaN;Mdl.Distribution='t';

[EstMdl,EstParamCov,LLF,info] = estimate(Mdl,ANZris,'display','off');

aic0t(p,q)=-2*LLF+2*(p+q+1);sic0t(p,q)=-2*LLF+log(nis)*(p+q+1);

% AR(1)-GARCH(p,q) with Student t errors

Mdl = arima('ARLags',1,'Variance',garch(p,q));Mdl.Distribution='t';

[EstMdl,EstParamCov,LLF,info] = estimate(Mdl,ANZris,'display','off');

aic1t(p,q)=-2*LLF+2*(p+q+2);sic1t(p,q)=-2*LLF+log(nis)*(p+q+2);

end

end

% Fit GARCH(1,1)-t model, as chosen by SIC

Mdl = garch(1,1);Mdl.Offset=NaN;Mdl.Distribution='t';

[EstMdl,EstParamCov,LLF,info] = estimate(Mdl,ANZris);

% infer the conditional variance and calculate standard deviations

v=infer(EstMdl,ANZris);s=sqrt(v);

% innovations

a = ANZris-EstMdl.Offset;

% standardised residuals

e=a./s;

% Plot the standardised residuals, conditional standard deviations and log

% returns above one another

figure;subplot(3,1,1);plot(ANZdates(2:nis+1,1),e);

%xlim([BHPdates(2) BHPdates(nis+1)]);

title('GARCH(1,1)-t Standardised Residuals');

subplot(3,1,2);plot(ANZdates(2:nis+1,1),s);

%xlim([ANZdates(2) BHPdates(nis+1)]);

title('GARCH(1,1)-t Conditional Standard Deviations');

subplot(3,1,3);plot(ANZdates(2:nis+1,1),ANZris);

%xlim([BHPdates(2) BHPdates(nis+1)]);title('Log Returns');

figure;plot(ANZdates(2:nis+1,1),ANZris,'c');

%xlim([BHPdates(2) BHPdates(nis+1)]);

hold on;plot(ANZdates(2:nis+1,1),s);

title('Log Returns and GARCH(1,1)-t Conditional Standard Deviations');

legend('ANZ returns','Conditional standard deviations',...

'Location','South','Orientation','horizontal');

df=EstMdl.Distribution.DoF;

% Assess fit to data

e1=sqrt(df)/sqrt(df-2)*e;   % this should have a Student-t with df degrees of freedom

eg=norminv(tcdf(e1,df));    % transform to normal errors for diagnostic analysis

figure;subplot(2,1,1);plot(ANZdates(2:nis+1,1),eg);

%xlim([BHPdates(2) BHPdates(nis+1)]);

title('GARCH(1,1)-t Transformed Standardised Residuals');

subplot(2,1,2);autocorr(eg);

figure;subplot(2,1,1);hist(eg,25);

title('Histogram of GARCH(1,1)-t Transformed Standardised Residuals');

subplot(2,1,2);qqplot(eg);

title('QQ plot of GARCH(1,1)-t Transformed Standardised Residuals');

figure;autocorr(eg.^2);

title('ACF of GARCH(1,1)-t Squared Transformed Standardised Residuals');

%LB test on standardised residuals

[H, pValue, Qstat, CriticalValue] = lbqtest(eg, [8 13], 0.05, [5 10])

%LB test on squared standardised residuals

[H, pValue, Qstat, CriticalValue] = lbqtest(eg.^2, [8 13], 0.05, [5 10])

% JB test

[skewness(eg) kurtosis(eg)]

[h,p] = jbtest(eg)

sigGt=s; % store conditional standard deviations as sigGt

四：GJR-GARCH 模型

GJR-GARCH模型的加强就在于他对于之前的mean correct return的正负会有不同的反馈。也就是前一天的回报率的盈亏情况也会对后面的方差产生很大的影响。我个人认为这是对GARCH 模型一个很大的提升，他讲考虑的层面提高了一个层级，虽然实用性还是有待商榷。当然建模之后还是要做LB和JB检验。

Matlab的实现：

Mdl = gjr(1,1);Mdl.Distribution='t';Mdl.Offset=NaN; % specify model

[EstMdl,EstParamCov,LLF,info] = estimate(Mdl,ANZris);

v = infer(EstMdl,ANZris);s=sqrt(v); %infer the conditional variance and calculate standard deviations innovations

a = ANZris-EstMdl.Offset; % innovations

e=a./s; % standardised residuals

% Plot the innovations, conditional standard deviations and log returns above one another

figure;subplot(3,1,1);plot(ANZdates(2:nis+1,1),e);

%xlim([BHPdates(2) BHPdates(nis+1)]);

title('GJR-GARCH-t Standardised Residuals');

subplot(3,1,2);plot(ANZdates(2:nis+1,1),s);

%xlim([BHPdates(2) BHPdates(nis+1)]);

title('GJR-GARCH-t Conditional Standard Deviations');

subplot(3,1,3);plot(ANZdates(2:nis+1,1),ANZris);

%xlim([BHPdates(2) BHPdates(nis+1)]);

title('Log Returns');

figure;plot(ANZdates(2:nis+1,1),ANZris,'c');

hold on;plot(ANZdates(2:nis+1,1),s);

%xlim([BHPdates(2) BHPdates(nis+1)]);

title('GJR-GARCH-t Conditional Standard Deviations and Log Returns');

legend('ANZ returns','Conditional standard deviations',...

'Location','South','Orientation','horizontal');

dfGJt=EstMdl.Distribution.DoF;

%Assess fit to data

% First transform to normal errors

teGJt=sqrt(dfGJt)/sqrt(dfGJt-2)*e;

geGJt=norminv(tcdf(teGJt,dfGJt));

figure;subplot(2,1,1);plot(ANZdates(2:nis+1,1),geGJt);

%xlim([BHPdates(2) BHPdates(nis+1)]);

title('GJR-GARCH-t Transformed Standardised Residuals');

subplot(2,1,2);autocorr(geGJt)

figure;subplot(2,1,1);hist(geGJt,25);

title('Histogram of GJR-GARCH-t Transformed Standardised Residuals');

subplot(2,1,2);qqplot(geGJt);

title('QQ plot for GJR-GARCH-t Transformed Standardised Residuals');

figure;autocorr(geGJt.^2);

title('ACF of GJR-GARCH-t Squared Transformed Standardised Residuals');

sigGJt=s;aGJt=a;

%LB test on standardised residuals

[H, pValue, Qstat, CriticalValue] = lbqtest(geGJt, [10 15], 0.05, [5 10])

%LB test on squared standardised residuals

[H, pValue, Qstat, CriticalValue] = lbqtest(geGJt.^2,[10 15], 0.05, [5 10])

% JB test

[skewness(geGJt) kurtosis(geGJt)]

[h,p] = jbtest(geGJt)

figure;plot(ANZdates(2:nis+1,1),sigGt);

hold on;plot(ANZdates(2:nis+1,1),sigGJt,'r');

title('Conditional Standard Deviations for GARCH-t and GJR-t');

%xlim([BHPdates(2) BHPdates(nis+1)]);legend('GARCH-t','GJR-t')

四．最后一部分讲一下预测的部分

不多解释，就是用proxy1,2,3来检验方差的bias。

Matlab的实现：

%% Lab Sheet 12: Forecasting with ARCH & GARCH Model

% Import IBM daily price index data

% Import columns 2:5 of "IBM80-97.xls" as numeric matrix called 'IBMdata'

% Import column 1 of "IBM80-97.xls" as a column vector named 'IBMdates'

% load lab12.mat;

%% Q1(a) Use SIC (and here AIC) to determine ARCH(p)

IBMp=IBMdata(:,4); % IBMdata is columns 2:5 of the "IBM80-97.xls" file imported as a numeric matrix

IBMr = 100*price2ret(IBMp); % calculate returns series

% Plot IBM prices and returns

figure;subplot(2,1,1);plot(IBMdates,IBMp);title('ANZ Prices');

%xlim([min(IBMdates) max(IBMdates)]);

subplot(2,1,2);plot(IBMdates(2:end,:),IBMr);title('ANZ Log Returns');

%xlim([min(IBMdates)+1 max(IBMdates)]);

ylim([-30 20]);

% Leave out last 550 for forecasting

n=length(IBMr);

nis=1769;

IBMris=IBMr(1:nis); % 'learning' sample

IBMrf=IBMr(nis+1:n); % forecast sample

% Fit model to in-sample data (first 4000 obs)

LLF=0;aic1=0;sic1=0;

for p=1:20

Mdl = garch(0,p);Mdl.Offset=NaN;Mdl.Distribution='t';

[EstMdl,EstParamCov,LLF,info] = estimate(Mdl,IBMris,'display','off');

aic1(p)=-2*LLF+2*p;sic1(p)=-2*LLF+log(length(IBMris))*p;

end

figure;subplot(2,1,1);plot(aic1,'r+-');title('AIC')

subplot(2,1,2);plot(sic1,'o-');title('SIC')

%SIC chooses p=7, (AIC chooses p>20)

%Define models

Mdla = garch(0,7);Mdla.Offset=NaN;Mdla.Distribution='t'; %ARCH(7)

Mdlg = garch(1,1);Mdlg.Offset=NaN;Mdlg.Distribution='t'; %GARCH(1,1)

Mdlgj = gjr(1,1);Mdlgj.Offset=NaN;Mdlgj.Distribution='t'; %GJRGARCH(1,1)

Mdleg = egarch(1,1);Mdleg.Offset=NaN;Mdleg.Distribution='t'; %EGARCH(1,0)

% Fit ARCH(7) model

[EstMdla,EstParamCov,LLF,info] = estimate(Mdla,IBMris);

va=infer(EstMdla,IBMris);sa=sqrt(va);

%% Q1(b) Fit GARCH(1,1), GJR-GARCH(1,1) and EGARCH(1,0) models, all with Student-t errors, to the in-sample data

%Fit GARCH(1,1) model

[EstMdlg,EstParamCov,LLF,info] = estimate(Mdlg,IBMris);

vg=infer(EstMdlg,IBMris);sg=sqrt(vg);

%Fit GJR-GARCH model

[EstMdlgj,EstParamCov,LLF,info] = estimate(Mdlgj,IBMris);

vgj=infer(EstMdlgj,IBMris);sgj=sqrt(vgj);

%Fit E-GARCH model

[EstMdleg,EstParamCov,LLF,info] = estimate(Mdleg,IBMris);

veg=infer(EstMdleg,IBMris);seg=sqrt(veg);

%% Q1(c) IGARCH model - Risk Metrics (no estimation of model required, just conditional vols)

% Estimate sample variance of returns to initialise variance

lr_var=var(IBMris);

% initialise the volatility series using long-run sample variance

% for time 0 and assuming a(0)=0

vrm=0.06*0+0.94*lr_var;

for t=2:nis+1

vrm(t)=0.06*(IBMris(t-1)^2)+0.94*vrm(t-1);

end

srm=sqrt(vrm)';

% Plot conditional volatility estimates for all models

allmod_is= [sa sg sgj seg];

figure;plot(IBMdates(2:nis+1),allmod_is);

title('Conditional Standard Deviations');

legend('ARCH(7)','GARCH','GJR','EGARCH');

%xlim([IBMdates(2) IBMdates(nis+1)]);

% Plot focusing on 1987 crash period

figure;plot(IBMdates(2:nis+1),allmod_is);

title('Conditional Standard Deviations - 1987 Crash');

legend('ARCH(7)','GARCH','GJR','EGARCH','RiskMetrics');

%xlim([IBMdates(1750) IBMdates(2250)]);

% Volatility proxies

prox1C=abs(IBMr-mean(IBMr));

prox1Ci=prox1C(1:nis);prox1Cf=prox1C(nis+1:2273);

highC=IBMdata(:,2);lowC=IBMdata(:,3);openC=IBMdata(:,1);

rangC=100*log(highC./lowC);rangC(rangC <= 0)=mean([0 min(rangC(rangC>0))]);

prox2C=sqrt(0.3607*(rangC.^2));prox2C=prox2C(2:end);

prox3C=1.107*prox2C.^2 + 0.68*((100*log(openC(2:end)./Close(1:end-1))).^2);

prox3C=sqrt(prox3C);

prox4C=exp(2*log(rangC)-0.86+2*0.29^2);

prox4C=sqrt(prox4C(2:end));

[length(prox1C) length(prox2C) length(prox3C) length(prox4C)]

% first 1515 proxies

prox1Ci=prox1C(1:nis);prox2Ci=prox2C(1:nis);

prox3Ci=prox3C(1:nis);prox4Ci=prox4C(1:nis);

% accuracy of volatility series estimates (in-sample).

[sqrt(mean((sa-prox1Ci).^2)) sqrt(mean((sa-prox2Ci).^2)) ...

sqrt(mean((sa-prox3Ci).^2)) sqrt(mean((sa-prox4Ci).^2));

sqrt(mean((sg-prox1Ci).^2)) sqrt(mean((sg-prox2Ci).^2)) ...

sqrt(mean((sg-prox3Ci).^2)) sqrt(mean((sg-prox4Ci).^2));

sqrt(mean((sgj-prox1Ci).^2)) sqrt(mean((sgj-prox2Ci).^2)) ...

sqrt(mean((sgj-prox3Ci).^2)) sqrt(mean((sgj-prox4Ci).^2));

sqrt(mean((seg-prox1Ci).^2)) sqrt(mean((seg-prox2Ci).^2)) ...

sqrt(mean((seg-prox3Ci).^2)) sqrt(mean((seg-prox4Ci).^2));

sqrt(mean((srm(1:nis,:)-prox1Ci).^2)) sqrt(mean((srm(1:nis,:)-prox2Ci).^2)) ...

sqrt(mean((srm(1:nis,:)-prox3Ci).^2)) sqrt(mean((srm(1:nis,:)-prox4Ci).^2))]

[mean(abs(sa-prox1Ci)) mean(abs(sa-prox2Ci)) ...

mean(abs(sa-prox3Ci)) mean(abs(sa-prox4Ci));

mean(abs(sg-prox1Ci)) mean(abs(sg-prox2Ci)) ...

mean(abs(sg-prox3Ci)) mean(abs(sg-prox4Ci));

mean(abs(sgj-prox1Ci)) mean(abs(sgj-prox2Ci)) ...

mean(abs(sgj-prox3Ci)) mean(abs(sgj-prox4Ci));

mean(abs(seg-prox1Ci)) mean(abs(seg-prox2Ci)) ...

mean(abs(seg-prox3Ci)) mean(abs(seg-prox4Ci));

mean(abs(srm(1:nis,:)-prox1Ci)) mean(abs(srm(1:nis,:)-prox2Ci)) ...

mean(abs(srm(1:nis,:)-prox3Ci)) mean(abs(srm(1:nis,:)-prox4Ci))]

% Plots for proxies 1, 2, and 3 separately with all models

% Proxy 1

figure;plot(IBMdates(2:nis+1),prox1Ci,'c+');hold on;

plot(IBMdates(2:nis+1),sa,'m');plot(IBMdates(2:nis+1),sg,'r');

plot(IBMdates(2:nis+1),sgj);plot(IBMdates(2:nis+1),seg,'k');

%plot(IBMdates(2:nis+1),srm(1:nis,:),'g');

title('Proxy 1 and conditional volatility from all models');

legend('Proxy 1','ARCH','GARCH','GJR','EGARCH');

% Proxy 2

figure;plot(IBMdates(2:nis+1),prox2Ci,'c+');hold on;

plot(IBMdates(2:nis+1),sa,'m');plot(IBMdates(2:nis+1),sg,'r');

plot(IBMdates(2:nis+1),sgj);plot(IBMdates(2:nis+1),seg,'k');

%plot(IBMdates(2:nis+1),srm(1:nis,:),'g');

title('Proxy 2 and conditional volatility from all models');

legend('Proxy 2','ARCH','GARCH','GJR','EGARCH');

% Proxy 3

figure;plot(IBMdates(2:nis+1),prox3Ci,'c+');hold on;

plot(IBMdates(2:nis+1),sa,'m');plot(IBMdates(2:nis+1),sg,'r');

plot(IBMdates(2:nis+1),sgj);plot(IBMdates(2:nis+1),seg,'k');

%plot(IBMdates(2:nis+1),srm(1:nis,:),'g');

title('Proxy 3 and conditional volatility from all models');

legend('Proxy 3','ARCH','GARCH','GJR','EGARCH');

%% Q1(d)&(e) One step ahead volatility forecasting for last 550 days

% Proxies for last 550 returns

prox1Cf=prox1C(end-503:end);prox2Cf=prox2C(end-503:end);

prox3Cf=prox3C(end-503:end);prox4Cf=prox4C(end-503:end);

% Calculate forecasts for last 550 days

% Initialise vectors for keeping forecast results

sigAf=0;sigGf=0;sigEGf=0;sigLTf=0;sig100f=0;sig25f=0;sigIGf=0;

for t=1:504

IBMris=IBMr(1:n-504+t-1);ncis=length(IBMris);

if mod(t,10)==0     % below code to reestimate models will run if the

% remainder when divided t is divided by 25 is 0, i.e. every 25 days

[EstMdla,EstParamCov,LLF,info]=estimate(Mdla,IBMris,'display','off');

va=infer(EstMdla,IBMris);sa=sqrt(va);

[EstMdlg,EstParamCov,LLF,info]=estimate(Mdlg,IBMris,'display','off');

vg=infer(EstMdlg,IBMris);sg=sqrt(vg);

[EstMdlgj,EstParamCov,LLF,info]=estimate(Mdlgj,IBMris,'display','off');

vgj=infer(EstMdlgj,IBMris);sgj=sqrt(vgj);

[EstMdleg,EstParamCov,LLF,info]=estimate(Mdleg,IBMris,'display','off');

veg=infer(EstMdleg,IBMris);seg=sqrt(veg);

end

sigAf(t)=sqrt(forecast(EstMdla,1,'Y0',IBMris));

sigGf(t)=sqrt(forecast(EstMdlg,1,'Y0',IBMris));

sigGJf(t)=sqrt(forecast(EstMdlgj,1,'Y0',IBMris));

sigEGf(t)=sqrt(forecast(EstMdleg,1,'Y0',IBMris));

% Risk metrics (I-GARCH) forecasts

if t==1

sigIGf(t)=vrm(end);

else

sigIGf(t)=sqrt(0.06*IBMris(end)^2+0.94*sigIGf(t-1)^2);

end

% ad-hoc forecasts

sigLTf(t)=std(IBMris);

sig100f(t)=std(IBMris(end-99:end));

sig25f(t)=std(IBMris(end-24:end));

end

%% Q1(f) Assess accuracy against proxies 1, 2, and 3

allmod_f = [sigAf' sigGf' sigGJf' sigEGf'];

% Plots for proxies 1, 2, and 3 with forecasts from parametric models

figure;plot(IBMdates(nis+2:end,:),prox1Cf,'yd');hold on;

plot(IBMdates(nis+2:end,:),allmod_f);

title('Proxy 1 and Garch model forecasts');

legend('Proxy 1','ARCH','GARCH','GJR','EGARCH','location','NorthWest');

%xlim([IBMdates(nis+2) IBMdates(end)]);

figure;plot(IBMdates(nis+2:end,:),prox2Cf,'yd');hold on;

plot(IBMdates(nis+2:end,:),allmod_f);

title('Proxy 2 and Garch model forecasts');

legend('Proxy 2','ARCH','GARCH','GJR','EGARCH','location','NorthWest');

%xlim([IBMdates(nis+2) IBMdates(end)]);

figure;plot(IBMdates(nis+2:end,:),prox3Cf,'yd');hold on;

plot(IBMdates(nis+2:end,:),allmod_f);

title('Proxy 3 and Garch model forecasts');

legend('Proxy 3','ARCH','GARCH','GJR','EGARCH','location','NorthWest');

%xlim([IBMdates(nis+2) IBMdates(end)]);

% Plots for proxies with ad-hoc measures

figure;plot(IBMdates(nis+2:end,:),prox1Cf,'yd');hold on;

plot(IBMdates(nis+2:end,:),prox4Cf,'rx');

plot(IBMdates(nis+2:end,:),sigLTf,'b');

plot(IBMdates(nis+2:end,:),sig100f,'m');

plot(IBMdates(nis+2:end,:),sig25f,'k');

title('Ad hoc measures and proxies 1 and 4');

legend('Proxy 1','Proxy 4','s(n)','s(100)','s(25)','location','NorthWest');

%xlim([IBMdates(nis+2) IBMdates(end)]);

figure;plot(IBMdates(nis+2:end,:),prox2Cf,'go');hold on;

plot(IBMdates(nis+2:end,:),prox3Cf,'cx');

plot(IBMdates(nis+2:end,:),sigLTf,'b');

plot(IBMdates(nis+2:end,:),sig100f,'m');

plot(IBMdates(nis+2:end,:),sig25f,'k');

title('Ad hoc measures and proxies 2 and 3');

legend('Proxy 2','Proxy 3','s(n)','s(100)','s(25)','location','NorthWest');

%xlim([IBMdates(nis+2) IBMdates(end)]);

figure;plot(IBMdates(nis+2:end,:),prox1Cf,'yd');hold on;

plot(IBMdates(nis+2:end,:),prox4Cf,'rx');

plot(IBMdates(nis+2:end,:),sigAf,'b');

plot(IBMdates(nis+2:end,:),sigGf,'m');

plot(IBMdates(nis+2:end,:),sigGJf,'y');

plot(IBMdates(nis+2:end,:),sigEGf,'k');

%plot(IBMdates(nis+2:end,:),sigIGf,'r');

title('GARCH models and proxies 1 and 4');

legend('Proxy 1','Proxy 4','ARCH','GARCH','GJR','EGARCH','location','NorthWest');

%xlim([IBMdates(nis+2) IBMdates(end)]);

% Plots of proxies 2 and 3 with all parametric models

figure;plot(IBMdates(nis+2:end,:),prox2Cf,'co');hold on;

plot(IBMdates(nis+2:end,:),prox3Cf,'yd');

plot(IBMdates(nis+2:end,:),sigAf,'b');

plot(IBMdates(nis+2:end,:),sigGf,'m');

plot(IBMdates(nis+2:end,:),sigGJf,'y');

plot(IBMdates(nis+2:end,:),sigEGf,'k');

%plot(IBMdates(nis+2:end,:),sigIGf,'r');

plot(IBMdates(nis+2:end,:),sigLTf,'b');

plot(IBMdates(nis+2:end,:),sig100f,'m');

plot(IBMdates(nis+2:end,:),sig25f,'k');

title('All forecasts and proxies 2 and 3');

legend('Proxy 2','Proxy 3','ARCH','GARCH','GJR','EGARCH','s(n)','s(100)','s(25)','location','NorthWest');

%xlim([IBMdates(nis+2) IBMdates(end)]);

% Calculate MSE and MAD for proxies 1, 2 and 3

% Proxy 1 measures

mseC11=sqrt(mean((sigAf'-prox1Cf).^2));madC11=mean(abs(sigAf'-prox1Cf));

mseC21=sqrt(mean((sigGf'-prox1Cf).^2));madC21=mean(abs(sigGf'-prox1Cf));

mseC31=sqrt(mean((sigEGf'-prox1Cf).^2));madC31=mean(abs(sigEGf'-prox1Cf));

mseC41=sqrt(mean((sigIGf'-prox1Cf).^2));madC41=mean(abs(sigIGf'-prox1Cf));

mseC51=sqrt(mean((sigLTf'-prox1Cf).^2));madC51=mean(abs(sigLTf'-prox1Cf));

mseC61=sqrt(mean((sig100f'-prox1Cf).^2));madC61=mean(abs(sig100f'-prox1Cf));

mseC71=sqrt(mean((sig25f'-prox1Cf).^2));madC71=mean(abs(sig25f'-prox1Cf));

% Proxy 2 measures

mseC12=sqrt(mean((sigAf'-prox2Cf).^2));madC12=mean(abs(sigAf'-prox2Cf));

mseC22=sqrt(mean((sigGf'-prox2Cf).^2));madC22=mean(abs(sigGf'-prox2Cf));

mseC32=sqrt(mean((sigEGf'-prox2Cf).^2));madC32=mean(abs(sigEGf'-prox2Cf));

mseC42=sqrt(mean((sigIGf'-prox2Cf).^2));madC42=mean(abs(sigIGf'-prox2Cf));

mseC52=sqrt(mean((sigLTf'-prox2Cf).^2));madC52=mean(abs(sigLTf'-prox2Cf));

mseC62=sqrt(mean((sig100f'-prox2Cf).^2));madC62=mean(abs(sig100f'-prox2Cf));

mseC72=sqrt(mean((sig25f'-prox2Cf).^2));madC72=mean(abs(sig25f'-prox2Cf));

% Proxy 3 measures

mseC13=sqrt(mean((sigAf'-prox3Cf).^2));madC13=mean(abs(sigAf'-prox3Cf));

mseC23=sqrt(mean((sigGf'-prox3Cf).^2));madC23=mean(abs(sigGf'-prox3Cf));

mseC33=sqrt(mean((sigEGf'-prox3Cf).^2));madC33=mean(abs(sigEGf'-prox3Cf));

mseC43=sqrt(mean((sigIGf'-prox3Cf).^2));madC43=mean(abs(sigIGf'-prox3Cf));

mseC53=sqrt(mean((sigLTf'-prox3Cf).^2));madC53=mean(abs(sigLTf'-prox3Cf));

mseC63=sqrt(mean((sig100f'-prox3Cf).^2));madC63=mean(abs(sig100f'-prox3Cf));

mseC73=sqrt(mean((sig25f'-prox3Cf).^2));madC73=mean(abs(sig25f'-prox3Cf));

[mseC11 mseC21 mseC31 mseC41 mseC51 mseC61 mseC71;

mseC12 mseC22 mseC32 mseC42 mseC52 mseC62 mseC72;

mseC13 mseC23 mseC33 mseC43 mseC53 mseC63 mseC73]

[madC11 madC21 madC31 madC41 madC51 madC61 madC71;

madC12 madC22 madC32 madC42 madC52 madC62 madC72;

madC13 madC23 madC33 madC43 madC53 madC63 madC73]