# practicalGarch.ssc Examples for Practical Garch paper # # author: Eric Zivot # created: 9/1/2006 # last updated: 4/25/2008 # load FinMetrics module module(finmetrics) # # simulating ARCH models # arch8.mod = list(a.value=1, arch=0.8, garch=0) arch8.sim = simulate.garch(model=arch8.mod, n=500, rseed=124) names(arch8.sim) par(mfrow=c(4,1)) tsplot(arch8.sim$et, main="Simulated ARCH(1) with a1 = 0.8") abline(h=0) tsplot(arch8.sim$et^2, main="Squared values") tsplot(abs(arch8.sim$et), main="Absolute values") tsplot(arch8.sim$sigma.t, main="Conditional Volatility") par(mfrow=c(1,1)) par(mfrow=c(3,1)) tmp = acf(arch8.sim$et, lag=20) tmp = acf(arch8.sim$et^2, lag=20) tmp = acf(abs(arch8.sim$et), lag=20) par(mfrow=c(1,1)) summaryStats(arch8.sim$et) qqPlot(arch8.sim$et, strip="Simulated ARCH(1)") title("Normal qq-plot") # garch(1,1) garch8.mod = list(a.value=1, arch=0.1, garch=0.7) garch8.sim = simulate.garch(model=garch8.mod, n=500, rseed=124) names(garch8.sim) par(mfrow=c(4,1)) tsplot(garch8.sim$et, main="Simulated GARCH(1,1) with a1 = 0.1 and b1=0.7") abline(h=0) tsplot(garch8.sim$et^2, main="Squared values") tsplot(abs(garch8.sim$et), main="Absolute values") tsplot(garch8.sim$sigma.t, main="Conditional Volatility") par(mfrow=c(1,1)) par(mfrow=c(3,1)) tmp = acf(garch8.sim$et, lag=20) tmp = acf(garch8.sim$et^2, lag=20) tmp = acf(abs(garch8.sim$et), lag=20) par(mfrow=c(1,1)) summaryStats(garch8.sim$et) qqPlot(garch8.sim$et, strip="Simulated GARCH(1,1)") title("Normal qq-plot") # # Detecting ARCH-GARCH behavior in raw data # # analysis of Microsoft stock start(msft.ts) end(msft.ts) # plot returns par(mfrow=c(3,1)) plot(msft.ts, reference.grid=F, main="Returns") abline(h=0) plot(msft.ts^2, reference.grid=F, main="Squared Returns") plot(abs(msft.ts), reference.grid=F, main="Absolute Returns") par(mfrow=c(1,1)) # plot SACF # note: no serial correlation in returns => set conditional mean = const in # garch fit par(mfrow=c(3,1)) tmp = acf(msft.ts, lag=20) tmp = acf(msft.ts^2, lag=20) tmp = acf(abs(msft.ts), lag=20) par(mfrow=c(1,1)) # compute summary statistics and normality test summaryStats(msft.ts) normalTest(msft.ts, method="jb") # plot normal qq-plot qqPlot(msft.ts) # smoothed density superimposed with normal distribution msft.density = density(msft.ts) plot(msft.density, type="l") msft.mean = mean(msft.ts) msft.sd = colStdevs(msft.ts) sort.s = sort(as.matrix(seriesData(msft.ts))) points(sort.s, dnorm(sort.s, mean=msft.mean, sd=msft.sd), type="l") # # testing for ARCH/GARCH effects # # modified q-stats autocorTest(seriesMerge(msft.ts,msft.ts^2,abs(msft.ts)), lag.n=1) autocorTest(seriesMerge(msft.ts,msft.ts^2,abs(msft.ts)), lag.n=5) autocorTest(seriesMerge(msft.ts,msft.ts^2,abs(msft.ts)), lag.n=10) autocorTest(seriesMerge(msft.ts,msft.ts^2,abs(msft.ts)), lag.n=20) # Engle's LM test for ARCH archTest(msft.ts, lag.n=1) archTest(msft.ts, lag.n=5) archTest(msft.ts, lag.n=10) archTest(msft.ts, lag.n=20) # BDS test for non iid behavior args(BDSTest) BDSTest(msft.ts, m=5) # # GARCH estimation # # estimate garch(1,1) model msft.g11 = garch(msft.ts~1, ~garch(1,1)) # compute OP standard errors summary(msft.g11) # compute QML standard errors summary(msft.g11, method="QMLE") # compute persistence parameter coef(msft.g11)[3] + coef(msft.g11)[4] # half-life of volatility shock log(0.5)/log(coef(msft.g11)[3] + coef(msft.g11)[4]) # asymptotic standard deviation msft.g11$asymp.sd # show graphical diagnostics # SACF of squared residuals plot(msft.g11, which.plot=10) tmp = acf(residuals(msft.g11, standardized=T)^2, main="Squared Residuals") # normal qqplot of residuals plot(msft.g11, which.plot=12) # series with conditional standard deviations plot(msft.g11, which.plot=1) # BDS test on garch residuals and log squared residuals BDSTest(residuals(msft.g11, standard=T), m=5, eps=c(0.5, 1, 1.5)) BDSTest(log(residuals(msft.g11, standard=T)^2), m=5, eps=c(0.5, 1, 1.5)) # # fit variety of ARCH/GARCH models and compare using model selection criteria # GARCH(1,1) wins according to BIC # msft.a1 = garch(msft.ts~1, ~garch(1,0)) msft.a2 = garch(msft.ts~1, ~garch(2,0)) msft.a3 = garch(msft.ts~1, ~garch(3,0)) msft.a4 = garch(msft.ts~1, ~garch(4,0)) msft.a5 = garch(msft.ts~1, ~garch(5,0)) msft.g21 = garch(msft.ts~1, ~garch(2,1)) msft.g22 = garch(msft.ts~1, ~garch(2,2)) msft.g12 = garch(msft.ts~1, ~garch(1,2)) compare.mgarch(msft.a1, msft.a2, msft.a3, msft.a4, msft.a5, msft.g11, msft.g12, msft.g21, msft.g22) # compare volatilities between ARCH(5) and GARCH(1,1) # note: garch volatilities are smoother than ARCH(5) volatilities par(mfrow=c(3,1)) plot(msft.ts, reference.grid=F, main="") title("MSFT Daily Returns") abline(h=0) plot(sigma.t(msft.g11), reference.grid=F, main="") title("Conditional Volatility from GARCH(1,1)") plot(sigma.t(msft.a5), reference.grid=F, main="") title("Conditional Volatility from ARCH(5)") par(mfrow=c(1,1)) # # asymmetric GARCH models # # diagnostics for leverage effects plot(msft.g11, which.plot=5) # note: Microsoft does not have any debt outstanding # correlations between r(t)^2 and r(t-1) # correlations between r(t)^2 and r(t-1) tmp.ts = seriesMerge(msft.ts^2, tslag(msft.ts, trim=T)) tmp = acf(tmp.ts, lag.max=1) tmp$acf[1,1,2] # compute Engle-Ng regression tests for leverage # Sign Bias test s.minus = as.numeric(seriesData(msft.ts) < 0) s.minus.ts = timeSeries(pos=positions(msft.ts), data = s.minus) colIds(s.minus.ts) = "S.minus" ehat.ts = residuals(OLS(MSFT~1, data=msft.ts)) tmp.ts = seriesMerge(ehat.ts^2, ehat.ts, s.minus.ts, 1 - s.minus.ts) colIds(tmp.ts) = c("ehat2", "ehat", "S.minus", "S.plus") SignBias.fit = OLS(ehat2~tslag(S.minus), data=tmp.ts) summary(SignBias.fit) NegativeSizeBias.fit = OLS(ehat2~tslag(S.minus*ehat), data=tmp.ts) summary(NegativeSizeBias.fit) PositiveSizeBias.fit = OLS(ehat2~tslag(S.plus*ehat), data=tmp.ts) summary(PositiveSizeBias.fit) All.fit = OLS(ehat2~S.minus+tslag(S.minus*ehat)+tslag(S.plus*ehat), data=tmp.ts) summary(All.fit) # estimate asymmetric models for MSFT # EGARCH(1,1): expect gamma to be negative msft.egarch = garch(msft.ts~1, ~egarch(1,1), leverage=T, trace=F) summary(msft.egarch, method="QMLE") msft.egarch$likelihood msft.egarch$asymp.sd # TGARCH(1,1): expect gamma to be positive msft.tgarch = garch(msft.ts~1, ~tgarch(1,1), trace=F) summary(msft.tgarch, method="QMLE") msft.tgarch$likelihood msft.tgarch$asymp.sd # PGARCH(1,1) with p=2; expect LEV to be negative msft.pgarch.2 = garch(msft.ts~1, ~garch(1,1), leverage=T, trace=F) summary(msft.pgarch.2, method="QMLE") msft.pgarch.2$likelihood msft.pgarch.2$asymp.sd # PGARCH(1,1) with p=1; expect LEV to be negative msft.pgarch.1 = garch(msft.ts~1,~pgarch(1,1,1),leverage=T,trace=F) summary(msft.pgarch.1, method="QMLE") msft.pgarch.1$likelihood msft.pgarch.1$asymp.sd # model comparison # TGARCH(1,1) or GARCH(1,1) with leverage seems to win compare.mgarch(msft.g11, msft.egarch, msft.tgarch, msft.pgarch.2, msft.pgarch.1) # plot news impact curves # TGARCH a0 = msft.tgarch$coef[2] a1 = msft.tgarch$coef[3] b1 = msft.tgarch$coef[4] g1 = msft.tgarch$coef[5] A = a0 + b1 * msft.tgarch$asymp.sd^2 epsilon = seq(-0.21, 0.21, length=100) sigma2.t.TGARCH = A + (a1+g1*(epsilon < 0))*(epsilon^2) # PGARCH 1 a0 = msft.pgarch.1$coef[2] a1 = msft.pgarch.1$coef[3] b1 = msft.pgarch.1$coef[4] g1 = msft.pgarch.1$coef[5] A = (a0 + b1 * msft.pgarch.1$asymp.sd)^2 error = abs(epsilon) + g1*epsilon sigma2.t.PGARCH = A + 2*sqrt(A)*a1*error + (a1*error)^2 # PGARCH 2 a0 = msft.pgarch.2$coef[2] a1 = msft.pgarch.2$coef[3] b1 = msft.pgarch.2$coef[4] g1 = msft.pgarch.2$coef[5] A = a0 + b1 * msft.pgarch.2$asymp.sd^2 error = abs(epsilon) + g1*epsilon sigma2.t.GARCH = A + a1*(error^2) # EGARCH a0 = msft.egarch$coef[2] a1 = msft.egarch$coef[3] b1 = msft.egarch$coef[4] g1 = msft.egarch$coef[5] A = ( (msft.egarch$asymp.sd^2)^b1 )*exp(a0) error = abs(epsilon) + g1*epsilon sigma2.egarch = A*exp(a1*error/msft.egarch$asymp.sd) matplot(cbind(epsilon, epsilon, epsilon, epsilon), cbind( sigma2.t.TGARCH, sigma2.t.PGARCH, sigma2.t.GARCH, sigma2.egarch), type="l") title("Asymmetric GARCH(1,1) Models for Microsoft") key(-0.05, 0.0045, lines=list(type="l", lty=1:4), text= list(c("TGARCH", "PGARCH 1", "PGARCH 2", "EGARCH")), border=1, adj=1) # # GARCH models with non Gaussian error distributions # # garch(1,1) with student-t errors # note: qq-plot looks pretty linear! msft.g11t = garch(msft.ts~1, ~garch(1,1), cond.dist="t") summary(msft.g11t, method="QMLE") plot(msft.g11t, which.plot=12) # garch(1,1) with GED errors # qq-plot for GED distribution is not available! This is a bug msft.g11ged = garch(msft.ts~1, ~garch(1,1), cond.dist="ged") summary(msft.g11ged, method="QMLE") plot(msft.g11ged, which.plot=12) # garch(1,1) with Double exponential errors msft.g11de = garch(msft.ts~1, ~garch(1,1), cond.dist="double.exp") summary(msft.g11de, method="QMLE") plot(msft.g11de, which.plot=12) # model comparison across alternative distributions # t-distribution wins! compare.mgarch(msft.g11, msft.g11t, msft.g11ged, msft.g11de) # # asymmetric models with t errors # # pgarch 2 with t errors msft.pgarch2.t = garch(msft.ts~1, ~garch(1,1), leverage=T, cond.dist="t", trace=F) summary(msft.pgarch2.t, method="QMLE") plot(msft.pgarch2.t, which.plot=12) # tgarch with t errors msft.tgarch.t = garch(msft.ts~1, ~tgarch(1,1), cond.dist="t", trace=F) summary(msft.tgarch.t, method="QMLE") plot(msft.tgarch.t, which.plot=12) # compare asymmetric models compare.mgarch(msft.g11, msft.g11t, msft.pgarch.2, msft.pgarch2.t, msft.tgarch, msft.tgarch.t) # # Fractionally integrated GARCH models # # # Compute diagnostics for long memory behavior # n.obs = seriesLength(msft.ts) sacf.msft.sq = acf(msft.ts^2, lag.max=1000) sum(sacf.msft.sq$acf[1:500]) sum(sacf.msft.sq$acf) sacf.msft.abs = acf(abs(msft.ts), lag.max=1000) sum(sacf.msft.abs$acf[1:500]) sum(sacf.msft.abs$acf) # # estimate 2 components model # msft.2comp2 = garch(msft.ts~1, ~pgarch.2comp(2)) summary(msft.2comp2, method="QMLE") msft.2comp2.lev = garch(msft.ts~1, ~pgarch.2comp(2), leverage=T) summary(msft.2comp2.lev, method="QMLE") msft.2comp2.lev.t = garch(msft.ts~1, ~pgarch.2comp(2), leverage=T, cond.dist="t") summary(msft.2comp2.lev.t, method="QMLE") plot(msft.2comp2.lev.t, which.plot=12) # # estimate FIGARCH with and without leverage # note: leverage only available with FIEGARCH # # scale the data for better convergence msft.fia1 = fgarch(msft.ts*100~1, ~figarch(1,0)) summary(msft.fia1) msft.fia2 = fgarch(msft.ts*100~1, ~figarch(2,0)) summary(msft.fia2) msft.fig11 = fgarch(msft.ts*100~1, ~figarch(1,1)) summary(msft.fig11, method="QMLE") msft.fieg11 = fgarch(msft.ts*100~1, ~fiegarch(1,1)) summary(msft.fieg11, method="QMLE") msft.fieg11.lev = fgarch(msft.ts*100~1, ~fiegarch(1,1), leverage=T) summary(msft.fieg11.lev, method="QMLE") plot(msft.fieg11.lev, which.plot=12) # compare volatilities between FIGARCH and FIEGARCH # note: par(mfrow=c(3,1)) plot(msft.ts*100, reference.grid=F, main="") title("Microsoft Daily Returns") abline(h=0) plot(sigma.t(msft.fig11), reference.grid=F, main="") title("Conditional Volatility from FIGARCH(1,1)") plot(sigma.t(msft.fieg11.lev), reference.grid=F, main="") title("Conditional Volatility from FIEGARCH(1,1)+LEV") par(mfrow=c(1,1)) # # GARCH model prediction # # compare long-run standard deviations n.obs = seriesLength(msft.ts) msft.sd = colStdevs(msft.ts) sd.vals = c(msft.sd, msft.a5$asymp.sd, msft.g11$asymp.sd, msft.g11t$asymp.sd, msft.pgarch.2$asymp.sd, msft.pgarch2.t$asymp.sd, msft.tgarch$asymp.sd, msft.tgarch.t$asymp.sd) names(sd.vals) = c("Sample","ARCH(5)","GARCH(1,1)","GARCH(1,1) t", "GARCH(1,1) LEV", "GARCH(1,1) t LEV", "TGARCH(1,1)", "TGARCH(1,1) t") # out of sample predictions using entire data set # predicted volatility of last sample observation sigma.t.last = c(seriesData(sigma.t(msft.a5)[n.obs]), seriesData(sigma.t(msft.g11)[n.obs])) names(sigma.t.last) = c("ARCH(5)","GARCH(1,1)") td = timeSeq(from="07/01/2003", by="bizdays", length.out=25) msft.a5.pred = predict(msft.a5, n.predict=25) msft.g11.pred = predict(msft.g11, n.predict=25) msft.g11t.pred = predict(msft.g11t, n.predict=25) msft.pgarch.2.pred = predict(msft.pgarch.2, n.predict=25) msft.pgarch2.t.pred = predict(msft.pgarch2.t, n.predict=25) msft.tgarch.pred = predict(msft.tgarch, n.predict=25) msft.tgarch.t.pred = predict(msft.tgarch.t, n.predict=25) msft.vol.fcst = timeSeries(pos=td, data=cbind(msft.a5.pred$sigma.pred, msft.g11.pred$sigma.pred, msft.g11t.pred$sigma.pred, msft.pgarch.2.pred$sigma.pred, msft.pgarch2.t.pred$sigma.pred, msft.tgarch.pred$sigma.pred, msft.tgarch.t.pred$sigma.pred)) colIds(msft.vol.fcst) = names(sd.vals)[-1] # show in-sample and out-of-sample volatilities par(mfrow=c(2,1)) tmp1 = timeSeries(pos=positions(msft.ts), data=as.matrix(rep(NA,n.obs))) tmp2 = timeSeries(pos=td, data=as.matrix(rep(NA,25))) msft.vol.fcst.a5 = timeSeries(pos=td, data=as.matrix(msft.a5.pred$sigma.pred)) msft.vol.fcst.g11 = timeSeries(pos=td, data=as.matrix(msft.g11.pred$sigma.pred)) tmp3.a5 = concat(tmp1,msft.vol.fcst.a5) tmp3.a5[timeEvent("06/30/2003")] = sigma.t(msft.a5)[timeEvent("06/30/2003")] tmp3.g11 = concat(tmp1,msft.vol.fcst.g11) tmp3.g11[timeEvent("06/30/2003")] = sigma.t(msft.g11)[timeEvent("06/30/2003")] tmp4.a5 = concat(sigma.t(msft.a5), tmp2) tmp4.g11 = concat(sigma.t(msft.g11), tmp2) smpl = timeEvent("02/01/2003","08/04/2003") plot(tmp4.a5[smpl], tmp3.a5[smpl], reference.grid=F, main="Predicted Volatility from ARCH(5) for Microsoft", plot.args=list(lty=c(1,2), lwd=2, col=c(2,5))) abline(h=sd.vals[2]) plot(tmp4.g11[smpl], tmp3.g11[smpl], reference.grid=F, main="Predicted Volatility from GARCH(1,1) for Microsoft", plot.args=list(lty=c(1,2), lwd=2, col=c(2,5))) abline(h=sd.vals[3]) par(mfrow=c(1,1)) # comparison of predictions. Only the ARCH(5) is different plot(msft.vol.fcst, reference.grid=F, main="Volatility Forecasts for Microsoft from Competing GARCH Models", plot.args=list(lty=1:7, col=1:7, lwd=2)) legend(0.6, 0.021, legend=colIds(msft.vol.fcst), lty=1:7, col=1:7, lwd=2) # predictions starting from June 23, 2003 are likely to be different between models # that allow and don't allow leverage effects because # there is a large negative residual at that date plot(residuals(msft.g11)[4360:4365], residuals(msft.pgarch.2)[4360:4365], reference.grid=F) abline(h=0) # TODO: forecasts from long memory garch models # riskMetrics # h-step ahead forecast is extrapolation of last observation msft.EWMA = EWMA(msft.ts^2, lambda = 0.94) msft.EWMA.last = sqrt(as.numeric(seriesData(msft.EWMA[n.obs]))) td = timeSeq(from="07/01/2003", by="bizdays", length.out=25) msft.EWMA.pred = timeSeries(pos=td, data=as.matrix(rep(msft.EWMA.last,25))) # 2 component model msft.2comp2.pred = predict(msft.2comp2, n.predict=25) msft.2comp2.lev.pred = predict(msft.2comp2.lev, n.predict=25) msft.2comp2.lev.t.pred = predict(msft.2comp2.lev.t, n.predict=25) # figarch and fiegarch msft.fig11.pred = predict(msft.fig11, n.predict=25) msft.fieg11.pred = predict(msft.fieg11, n.predict=25) msft.fieg11.lev.pred = predict(msft.fieg11.lev, n.predict=25) td = timeSeq(from="07/01/2003", by="bizdays", length.out=25) msft.vol.fcst.lm = timeSeries(pos=td, data=cbind(seriesData(msft.EWMA.pred), msft.2comp2.pred$sigma.pred, msft.fig11.pred$sigma.pred/100, msft.fieg11.lev.pred$sigma.pred/100)) colIds(msft.vol.fcst.lm) = c("EWMA", "2 COMP", "FIGARCH", "FIEGARCH") plot(msft.vol.fcst.lm, reference.grid=F) # comparison of predictions. Asymmetric models are very different plot(msft.vol.fcst[,c(2,6)], msft.vol.fcst.lm[,c(2,3,4)], reference.grid=F, main="Volatility Forecasts for Microsoft from Competing GARCH Models", plot.args=list(lty=1:5, col=1:5, lwd=2)) legend(0.6, 0.0103, legend=colIds(msft.vol.fcst), lty=1:7, col=1:7, lwd=2) # simulation-based forecasts # ARCH(5) # these seem to consistently underpredict and start at a lower value # than the analytic predictions sigma.start = as.numeric(msft.a5$sigma.t[n.obs]) eps.start = as.numeric(msft.a5$residuals[n.obs]) eps.start = matrix(eps.start, 1, 1000) set.seed(111) error = rbind(eps.start, matrix(rnorm(25*1000), 25)) # different simulations are the columns of the following matrix set.seed(10) msft.a5.sim.pred = simulate(msft.a5, n=25, n.rep=1000, n.start=0, etat=error, sigma.start=sigma.start)$sigma.t msft.a5.sim.sigma.pred = rowMeans(msft.a5.sim.pred) cbind(msft.a5.pred$sigma.pred, msft.a5.sim.sigma.pred) tsplot(msft.a5.pred$sigma.pred, msft.a5.sim.sigma.pred) # GARCH(1,1) sigma.start = as.numeric(msft.g11$sigma.t[n.obs]) eps.start = as.numeric(msft.g11$residuals[n.obs]) eps.start = matrix(eps.start, 1, 1000) set.seed(111) error = rbind(eps.start, matrix(rnorm(25*1000), 25)) set.seed(10) msft.g11.sim.pred = simulate(msft.g11, n=25, n.rep=1000, n.start=0, etat=error, sigma.start=sigma.start)$sigma.t # incorrect forecasts msft.g11.sim.sigma.pred = rowMeans(msft.g11.sim.pred) # correct forecasts msft.g11.sim.sigma.pred2 = sqrt(rowMeans(msft.g11.sim.pred^2)) cbind(msft.g11.pred$sigma.pred, msft.g11.sim.sigma.pred, msft.g11.sim.sigma.pred2) tsplot(msft.g11.pred$sigma.pred, msft.g11.sim.sigma.pred, msft.g11.sim.sigma.pred2) # GARCH(1,1) with leverage sigma.start = as.numeric(msft.pgarch.2$sigma.t[n.obs]) eps.start = as.numeric(msft.pgarch.2$residuals[n.obs]) eps.start = matrix(eps.start, 1, 1000) set.seed(111) error = rbind(eps.start, matrix(rnorm(25*1000), 25)) set.seed(10) msft.pgarch.2.sim.pred = simulate(msft.pgarch.2, n=25, n.rep=1000, n.start=0, etat=error, sigma.start=sigma.start)$sigma.t msft.pgarch.2.sim.sigma.pred2 = sqrt(rowMeans(msft.pgarch.2.sim.pred^2)) msft.pgarch.2.sim.sigma.pred = rowMeans(msft.pgarch.2.sim.pred) cbind(msft.pgarch.2.pred$sigma.pred, msft.pgarch.2.sim.sigma.pred) tsplot(msft.pgarch.2.pred$sigma.pred, msft.pgarch.2.sim.sigma.pred, msft.pgarch.2.sim.sigma.pred2) # TODO: simulation based forecasts from egarch, two component and garch(2,2) # # Volatility term structures # term.dat = 1:200 msft.volTermStruc = cumsum(predict(msft.g11, n.predict=200)$sigma.pred^2)/term.dat tsplot(msft.volTermStruc) abline(h=msft.g11$asymp.sd^2) # TODO: volatility term structures # # GARCH models and VaR calculations # # TODO: GARCH VaR calculations # # estimation of GARCH models adjusted for outliers # # # Structural stability of GARCH estimates - rolling analysis # # TODO: rolling analysis of garch ############################################################################### # analysis of S&P 500 returns ############################################################################### start(sp500.ts) end(sp500.ts) # plot returns par(mfrow=c(3,1)) plot(sp500.ts, reference.grid=F, main="Returns") abline(h=0) plot(sp500.ts^2, reference.grid=F, main="Squared Returns") plot(abs(sp500.ts), reference.grid=F, main="Absolute Returns") par(mfrow=c(1,1)) # plot SACF par(mfrow=c(3,1)) tmp = acf(sp500.ts, lag=20) tmp = acf(sp500.ts^2, lag=20) tmp = acf(abs(sp500.ts), lag=20) par(mfrow=c(1,1)) # compute summary statistics and normality test summaryStats(sp500.ts) normalTest(sp500.ts, method="jb") # plot normal qq-plot qqPlot(sp500.ts) qqnorm(as.matrix(seriesData(sp500.ts)), ylab="sp500.ts") qqline(sp500.ts) # smoothed density superimposed with normal distribution # need to adjust scale on x-axis sp500.density = density(sp500.ts) plot(sp500.density, type="l") sp500.mean = mean(sp500.ts) sp500.sd = colStdevs(sp500.ts) sort.s = sort(as.matrix(seriesData(sp500.ts))) points(sort.s, dnorm(sort.s, mean=sp500.mean, sd=sp500.sd), type="l") # # testing for ARCH/GARCH effects # # modified q-stats autocorTest(seriesMerge(sp500.ts,sp500.ts^2,abs(sp500.ts)), lag.n=1) autocorTest(seriesMerge(sp500.ts,sp500.ts^2,abs(sp500.ts)), lag.n=5) autocorTest(seriesMerge(sp500.ts,sp500.ts^2,abs(sp500.ts)), lag.n=10) autocorTest(seriesMerge(sp500.ts,sp500.ts^2,abs(sp500.ts)), lag.n=20) args(BDSTest) BDSTest(sp500.ts, m=5) # Engle's LM test for ARCH archTest(sp500.ts, lag.n=1) archTest(sp500.ts, lag.n=5) archTest(sp500.ts, lag.n=10) archTest(sp500.ts, lag.n=20) # # GARCH estimation # # estimate garch(1,1) model sp500.g11 = garch(sp500.ts~1, ~garch(1,1)) # compute OP standard errors summary(sp500.g11) # compute QML standard errors summary(sp500.g11, method="QMLE") # compute persistence parameter coef(sp500.g11)[3] + coef(sp500.g11)[4] # half-life of volatility shock log(0.5)/log(coef(sp500.g11)[3] + coef(sp500.g11)[4]) # asymptotic standard deviation sp500.g11$asymp.sd # show graphical diagnostics # SACF of squared residuals plot(sp500.g11, which.plot=10) tmp = acf(residuals(sp500.g11, standardized=T)^2) # normal qqplot of residuals plot(sp500.g11, which.plot=12) # BDS test on garch residuals and log squared residuals BDSTest(residuals(sp500.g11, standard=T), m=5, eps=c(0.5, 1, 1.5)) BDSTest(log(residuals(sp500.g11, standard=T)^2), m=5, eps=c(0.5, 1, 1.5)) # # fit variety of ARCH/GARCH models and compare using model selection criteria # sp500.a1 = garch(sp500.ts~1, ~garch(1,0)) sp500.a2 = garch(sp500.ts~1, ~garch(2,0)) sp500.a3 = garch(sp500.ts~1, ~garch(3,0)) sp500.a4 = garch(sp500.ts~1, ~garch(4,0)) sp500.a5 = garch(sp500.ts~1, ~garch(5,0)) sp500.g21 = garch(sp500.ts~1, ~garch(2,1)) sp500.g22 = garch(sp500.ts~1, ~garch(2,2)) sp500.g12 = garch(sp500.ts~1, ~garch(1,2)) compare.mgarch(sp500.a1, sp500.a2, sp500.a3, sp500.a4, sp500.a5, sp500.g11, sp500.g12, sp500.g21, sp500.g22) # compare volatilities between ARCH(5) and GARCH(1,1) # note: garch volatilities are smoother than ARCH(5) volatilities par(mfrow=c(3,1)) plot(sp500.ts, reference.grid=F, main="") title("S&P 500 Daily Returns") abline(h=0) plot(sigma.t(sp500.g11), reference.grid=F, main="") title("Conditional Volatility from GARCH(1,1)") plot(sigma.t(sp500.a5), reference.grid=F, main="") title("Conditional Volatility from ARCH(5)") par(mfrow=c(1,1)) # # spurious inference in GARCH models # See Ma, Nelson, Startz (2007) SNDE # # create monthly sp500 returns sp500m.ts = aggregateSeries(sp500.ts, FUN=sum, by="months") par(mfrow=c(3,1)) plot(sp500m.ts, reference.grid=F, main="Returns") abline(h=0) plot(sp500m.ts^2, reference.grid=F, main="Squared Returns") plot(abs(sp500m.ts), reference.grid=F, main="Absolute Returns") par(mfrow=c(1,1)) # plot SACF par(mfrow=c(3,1)) tmp = acf(sp500m.ts, lag=20) tmp = acf(sp500m.ts^2, lag=20) tmp = acf(abs(sp500m.ts), lag=20) par(mfrow=c(1,1)) # # testing for ARCH/GARCH effects # # modified q-stats autocorTest(seriesMerge(sp500m.ts,sp500m.ts^2,abs(sp500m.ts)), lag.n=1) autocorTest(seriesMerge(sp500m.ts,sp500m.ts^2,abs(sp500m.ts)), lag.n=5) autocorTest(seriesMerge(sp500m.ts,sp500m.ts^2,abs(sp500m.ts)), lag.n=10) autocorTest(seriesMerge(sp500m.ts,sp500m.ts^2,abs(sp500m.ts)), lag.n=20) args(BDSTest) BDSTest(sp500m.ts, m=5) # Engle's LM test for ARCH archTest(sp500m.ts, lag.n=1) archTest(sp500m.ts, lag.n=5) archTest(sp500m.ts, lag.n=10) archTest(sp500m.ts, lag.n=20) # fit garch(1,1) to monthly returns sp500m.g11 = garch(sp500m.ts~1, ~garch(1,1)) summary(sp500m.g11) # fit arch(5,0) to monthly returns sp500m.a1 = garch(sp500m.ts~1, ~garch(1,0)) summary(sp500m.a1) # compute profile loglikelihood # set starting value for arch = 0.04 # set starting value for garch = values bebween 0 and 0.95 # set starting value for a equal to sample var * (1 - (arch + garch)) sp500m.var = colVars(sp500m.ts) start.garch = sp500m.g11$model garch.values = seq(from=0, to=0.92, by=0.01) likelihood.values = matrix(0,length(garch.values),1) start.garch$garch$which = F start.garch$arch$value = 0.04 start.garch$c.value = colMeans(sp500m.ts) numGarch = length(garch.values) for(i in 1:numGarch) { start.garch$garch$value=garch.values[i] start.garch$a.value = sp500m.var*(1 - (0.04 + start.garch$garch$value)) tmp.fit = garch(series=sp500m.ts, model=start.garch, trace=F) likelihood.values[i] = tmp.fit$likelihood } plot(garch.values, likelihood.values, type="l", main="Profile likelihood of GARCH(1,1)") # # asymmetric GARCH models # # diagnostics for leverage effects plot(sp500.g11, which.plot=5) # correlations between r(t)^2 and r(t-1) tmp.ts = seriesMerge(sp500.ts^2, tslag(sp500.ts, trim=T)) tmp = acf(tmp.ts, lag.max=1) tmp$acf[1,1,2] # compute Engle-Ng regression tests for leverage # Sign Bias test s.minus = as.numeric(seriesData(sp500.ts) < 0) s.minus.ts = timeSeries(pos=positions(sp500.ts), data = s.minus) colIds(s.minus.ts) = "S.minus" ehat.ts = residuals(OLS(SP500~1, data=sp500.ts)) tmp.ts = seriesMerge(ehat.ts^2, ehat.ts, s.minus.ts, 1 - s.minus.ts) colIds(tmp.ts) = c("ehat2", "ehat", "S.minus", "S.plus") SignBias.fit = OLS(ehat2~tslag(S.minus), data=tmp.ts) summary(SignBias.fit) NegativeSizeBias.fit = OLS(ehat2~tslag(S.minus*ehat), data=tmp.ts) summary(NegativeSizeBias.fit) PositiveSizeBias.fit = OLS(ehat2~tslag(S.plus*ehat), data=tmp.ts) summary(PositiveSizeBias.fit) All.fit = OLS(ehat2~S.minus+tslag(S.minus*ehat)+tslag(S.plus*ehat), data=tmp.ts) summary(All.fit) # estimate asymmetric models for sp500 # EGARCH(1,1): expect gamma to be negative sp500.egarch = garch(sp500.ts~1, ~egarch(1,1), leverage=T, trace=F) summary(sp500.egarch, method="QMLE") sp500.egarch$likelihood sp500.egarch$asymp.sd # TGARCH(1,1): expect gamma to be positive sp500.tgarch = garch(sp500.ts~1, ~tgarch(1,1), trace=F) summary(sp500.tgarch, method="QMLE") sp500.tgarch$likelihood sp500.tgarch$asymp.sd # PGARCH(1,1) with p=2; expect LEV to be negative sp500.pgarch.2 = garch(sp500.ts~1, ~garch(1,1), leverage=T, trace=F) summary(sp500.pgarch.2, method="QMLE") sp500.pgarch.2$likelihood sp500.pgarch.2$asymp.sd # PGARCH(1,1) with p=1; expect LEV to be negative sp500.pgarch.1 = garch(sp500.ts~1,~pgarch(1,1,1),leverage=T,trace=F) summary(sp500.pgarch.1, method="QMLE") sp500.pgarch.1$likelihood sp500.pgarch.1$asymp.sd # model comparison # TGARCH(1,1) or GARCH(1,1) with leverage seems to win compare.mgarch(sp500.g11, sp500.egarch, sp500.tgarch, sp500.pgarch.2, sp500.pgarch.1) # plot news impact curves # TARCH a0 = sp500.tgarch$coef[2] a1 = sp500.tgarch$coef[3] b1 = sp500.tgarch$coef[4] g1 = sp500.tgarch$coef[5] A = a0 + b1 * sp500.tgarch$asymp.sd^2 epsilon = seq(-0.21, 0.21, length=100) sigma2.t.TGARCH = A + (a1+g1*(epsilon < 0))*(epsilon^2) # PGARCH 1 a0 = sp500.pgarch.1$coef[2] a1 = sp500.pgarch.1$coef[3] b1 = sp500.pgarch.1$coef[4] g1 = sp500.pgarch.1$coef[5] A = (a0 + b1 * sp500.pgarch.1$asymp.sd)^2 error = abs(epsilon) + g1*epsilon sigma2.t.PGARCH = A + 2*sqrt(A)*a1*error + (a1*error)^2 # PGARCH 2 a0 = sp500.pgarch.2$coef[2] a1 = sp500.pgarch.2$coef[3] b1 = sp500.pgarch.2$coef[4] g1 = sp500.pgarch.2$coef[5] A = a0 + b1 * sp500.pgarch.2$asymp.sd^2 error = abs(epsilon) + g1*epsilon sigma2.t.GARCH = A + a1*(error^2) # EGARCH a0 = sp500.egarch$coef[2] a1 = sp500.egarch$coef[3] b1 = sp500.egarch$coef[4] g1 = sp500.egarch$coef[5] A = ( (sp500.egarch$asymp.sd^2)^b1 )*exp(a0) error = abs(epsilon) + g1*epsilon sigma2.egarch = A*exp(a1*error/sp500.egarch$asymp.sd) matplot(cbind(epsilon, epsilon, epsilon, epsilon), cbind( sigma2.t.TGARCH, sigma2.t.PGARCH, sigma2.t.GARCH, sigma2.egarch), type="l") title("Asymmetric GARCH(1,1) Models for S&P 500") key(-0.05, 0.0045, lines=list(type="l", lty=1:4), text= list(c("TGARCH", "PGARCH 1", "PGARCH 2", "EGARCH")), border=1, adj=1) # # GARCH models with non Gaussian error distributions # # garch(1,1) with student-t errors # note: qq-plot looks pretty linear! sp500.g11t = garch(sp500.ts~1, ~garch(1,1), cond.dist="t") summary(sp500.g11t, method="QMLE") plot(sp500.g11t, which.plot=12) # garch(1,1) with GED errors # qq-plot for GED distribution is not available! This is a bug sp500.g11ged = garch(sp500.ts~1, ~garch(1,1), cond.dist="ged") summary(sp500.g11ged, method="QMLE") plot(sp500.g11ged, which.plot=12) # garch(1,1) with Double exponential errors sp500.g11de = garch(sp500.ts~1, ~garch(1,1), cond.dist="double.exp") summary(sp500.g11de, method="QMLE") plot(sp500.g11de, which.plot=12) # model comparison across alternative distributions # t-distribution wins! compare.mgarch(sp500.g11, sp500.g11t, sp500.g11ged, sp500.g11de) # # asymmetric models with t errors # # pgarch 2 with t errors sp500.pgarch2.t = garch(sp500.ts~1, ~garch(1,1), leverage=T, cond.dist="t", trace=F) summary(sp500.pgarch2.t, method="QMLE") plot(sp500.pgarch2.t, which.plot=12) # pgarch 1 with t errors sp500.pgarch1.t = garch(sp500.ts~1, ~pgarch(1,1,1), leverage=T, cond.dist="t", trace=F) summary(sp500.pgarch1.t, method="QMLE") plot(sp500.pgarch1.t, which.plot=12) # tgarch with t errors sp500.tgarch.t = garch(sp500.ts~1, ~tgarch(1,1), cond.dist="t", trace=F) summary(sp500.tgarch.t, method="QMLE") plot(sp500.tgarch.t, which.plot=12) # compare asymmetric models compare.mgarch(sp500.g11, sp500.g11t, sp500.pgarch.2, sp500.pgarch2.t, sp500.pgarch1.t,sp500.tgarch, sp500.tgarch.t) # # Fractionally integrated GARCH models # # # Compute diagnostics for long memory behavior # # # estimate 2 components model # sp500.2comp2 = garch(sp500.ts~1, ~pgarch.2comp(2)) summary(sp500.2comp2, method="QMLE") sp500.2comp2.lev = garch(sp500.ts~1, ~pgarch.2comp(2), leverage=T) summary(sp500.2comp2.lev, method="QMLE") sp500.2comp2.lev.t = garch(sp500.ts~1, ~pgarch.2comp(2), leverage=T, cond.dist="t") summary(sp500.2comp2.lev.t, method="QMLE") plot(sp500.2comp2.lev.t, which.plot=12) # # estimate FIGARCH with and without leverage # # scale the data for better convergence sp500.fia1 = fgarch(sp500.ts*100~1, ~figarch(1,0)) summary(sp500.fia1) sp500.fia2 = fgarch(sp500.ts*100~1, ~figarch(2,0)) summary(sp500.fia2) sp500.fig11 = fgarch(sp500.ts*100~1, ~figarch(1,1)) summary(sp500.fig11, method="QMLE") sp500.fig11.lev = fgarch(sp500.ts*100~1, ~figarch(1,1), leverage=T) summary(sp500.fig11.lev, method="QMLE") sp500.fieg11 = fgarch(sp500.ts*100~1, ~fiegarch(1,1)) summary(sp500.fieg11, method="QMLE") sp500.fieg11.lev = fgarch(sp500.ts*100~1, ~fiegarch(1,1), leverage=T) summary(sp500.fieg11.lev, method="QMLE") plot(sp500.fieg11.lev, which.plot=12) # compare volatilities between FIGARCH and FIEGARCH # note: garch volatilities are smoother than ARCH(5) volatilities par(mfrow=c(3,1)) plot(sp500.ts*100, reference.grid=F, main="") title("S&P 500 Daily Returns") abline(h=0) plot(sigma.t(sp500.fig11.lev), reference.grid=F, main="") title("Conditional Volatility from FIGARCH(1,1)+LEV") plot(sigma.t(sp500.fieg11.lev), reference.grid=F, main="") title("Conditional Volatility from FIEGARCH(1,1)+LEV") par(mfrow=c(1,1)) # # Fractionally integrated GARCH models # # # Compute diagnostics for long memory behavior # n.obs = seriesLength(sp500.ts) sacf.sp500.sq = acf(sp500.ts^2, lag.max=1000) sum(sacf.sp500.sq$acf[1:500]) sum(sacf.sp500.sq$acf) sacf.sp500.abs = acf(abs(sp500.ts), lag.max=1000) sum(sacf.sp500.abs$acf[1:500]) sum(sacf.sp500.abs$acf) # # estimate 2 components model # sp500.2comp2 = garch(sp500.ts~1, ~pgarch.2comp(2)) summary(sp500.2comp2, method="QMLE") sp500.2comp2.lev = garch(sp500.ts~1, ~pgarch.2comp(2), leverage=T) summary(sp500.2comp2.lev, method="QMLE") sp500.2comp2.t = garch(sp500.ts~1, ~pgarch.2comp(2), leverage=F, cond.dist="t") summary(sp500.2comp2.t, method="QMLE") sp500.2comp2.lev.t = garch(sp500.ts~1, ~pgarch.2comp(2), leverage=T, cond.dist="t") summary(sp500.2comp2.lev.t, method="QMLE") plot(sp500.2comp2.t, which.plot=12) # # estimate FIGARCH with and without leverage # note: leverage effects only available with FIEGARCH # # scale the data for better convergence sp500.fia1 = fgarch(sp500.ts*100~1, ~figarch(1,0)) summary(sp500.fia1) sp500.fia2 = fgarch(sp500.ts*100~1, ~figarch(2,0)) summary(sp500.fia2) sp500.fig11 = fgarch(sp500.ts*100~1, ~figarch(1,1)) summary(sp500.fig11, method="QMLE") sp500.fieg11 = fgarch(sp500.ts*100~1, ~fiegarch(1,1)) summary(sp500.fieg11, method="QMLE") sp500.fieg11.lev = fgarch(sp500.ts*100~1, ~fiegarch(1,1), leverage=T) summary(sp500.fieg11.lev, method="QMLE") plot(sp500.fieg11.lev, which.plot=12) # compare volatilities between FIGARCH and FIEGARCH # note: par(mfrow=c(3,1)) plot(sp500.ts*100, reference.grid=F, main="") title("S&P 500 Daily Returns") abline(h=0) plot(sigma.t(sp500.fig11), reference.grid=F, main="") title("Conditional Volatility from FIGARCH(1,1)") plot(sigma.t(sp500.fieg11.lev), reference.grid=F, main="") title("Conditional Volatility from FIEGARCH(1,1)+LEV") par(mfrow=c(1,1)) # # GARCH model prediction # # compare long-run standard deviations n.obs = seriesLength(sp500.ts) sp500.sd = colStdevs(sp500.ts) sd.vals = c(sp500.sd, sp500.a5$asymp.sd, sp500.g11$asymp.sd, sp500.g11t$asymp.sd, sp500.pgarch.2$asymp.sd, sp500.pgarch2.t$asymp.sd, sp500.tgarch$asymp.sd, sp500.tgarch.t$asymp.sd) names(sd.vals) = c("Sample","ARCH(5)","GARCH(1,1)","GARCH(1,1) t", "GARCH(1,1) LEV", "GARCH(1,1) t LEV", "TGARCH(1,1)", "TGARCH(1,1) t") # out of sample predictions using entire data set # predicted volatility of last sample observation sigma.t.last = c(seriesData(sigma.t(sp500.a5)[n.obs]), seriesData(sigma.t(sp500.g11)[n.obs])) names(sigma.t.last) = c("ARCH(5)","GARCH(1,1)") td = timeSeq(from="07/01/2003", by="bizdays", length.out=75) sp500.a5.pred = predict(sp500.a5, n.predict=75) sp500.g11.pred = predict(sp500.g11, n.predict=75) sp500.g11t.pred = predict(sp500.g11t, n.predict=75) sp500.pgarch.2.pred = predict(sp500.pgarch.2, n.predict=75) sp500.pgarch2.t.pred = predict(sp500.pgarch2.t, n.predict=75) sp500.tgarch.pred = predict(sp500.tgarch, n.predict=75) sp500.tgarch.t.pred = predict(sp500.tgarch.t, n.predict=75) sp500.vol.fcst = timeSeries(pos=td, data=cbind(sp500.a5.pred$sigma.pred, sp500.g11.pred$sigma.pred, sp500.g11t.pred$sigma.pred, sp500.pgarch.2.pred$sigma.pred, sp500.pgarch2.t.pred$sigma.pred, sp500.tgarch.pred$sigma.pred, sp500.tgarch.t.pred$sigma.pred)) colIds(sp500.vol.fcst) = names(sd.vals)[-1] # show in-sample and out-of-sample volatilities par(mfrow=c(2,1)) tmp1 = timeSeries(pos=positions(sp500.ts), data=as.matrix(rep(NA,n.obs))) tmp2 = timeSeries(pos=td, data=as.matrix(rep(NA,75))) sp500.vol.fcst.a5 = timeSeries(pos=td, data=as.matrix(sp500.a5.pred$sigma.pred)) sp500.vol.fcst.g11 = timeSeries(pos=td, data=as.matrix(sp500.g11.pred$sigma.pred)) tmp3.a5 = concat(tmp1,sp500.vol.fcst.a5) tmp3.a5[timeEvent("06/30/2003")] = sigma.t(sp500.a5)[timeEvent("06/30/2003")] tmp3.g11 = concat(tmp1,sp500.vol.fcst.g11) tmp3.g11[timeEvent("06/30/2003")] = sigma.t(sp500.g11)[timeEvent("06/30/2003")] tmp4.a5 = concat(sigma.t(sp500.a5), tmp2) tmp4.g11 = concat(sigma.t(sp500.g11), tmp2) smpl = timeEvent("02/01/2003","08/04/2003") plot(tmp4.a5[smpl], tmp3.a5[smpl], reference.grid=F, main="Predicted Volatility from ARCH(5) for S&P 500", plot.args=list(lty=c(1,2), lwd=2, col=c(2,5))) abline(h=sd.vals[2]) plot(tmp4.g11[smpl], tmp3.g11[smpl], reference.grid=F, main="Predicted Volatility from GARCH(1,1) for S&P 500", plot.args=list(lty=c(1,2), lwd=2, col=c(2,5))) abline(h=sd.vals[3]) par(mfrow=c(1,1)) # comparison of predictions. Asymmetric models are very different plot(sp500.vol.fcst, reference.grid=F, main="Volatility Forecasts for S&P 500 from Competing GARCH Models", plot.args=list(lty=1:7, col=1:7, lwd=2)) legend(0.6, 0.0103, legend=colIds(sp500.vol.fcst), lty=1:7, col=1:7, lwd=2) # predictions starting from June 23, 2003 are likely to be different between models # that allow and don't allow leverage effects because # there is a large negative residual at that date plot(residuals(sp500.g11)[4360:4365], residuals(sp500.pgarch.2)[4360:4365], reference.grid=F) abline(h=0) # TODO: forecasts from long memory garch models # riskMetrics # h-step ahead forecast is extrapolation of last observation sp500.EWMA = EWMA(sp500.ts^2, lambda = 0.94) sp500.EWMA.last = sqrt(as.numeric(seriesData(sp500.EWMA[n.obs]))) td = timeSeq(from="07/01/2003", by="bizdays", length.out=75) sp500.EWMA.pred = timeSeries(pos=td, data=as.matrix(rep(sp500.EWMA.last,75))) # 2 component model sp500.2comp2.pred = predict(sp500.2comp2, n.predict=75) sp500.2comp2.lev.pred = predict(sp500.2comp2.lev, n.predict=75) sp500.2comp2.lev.t.pred = predict(sp500.2comp2.lev.t, n.predict=75) # figarch and fiegarch sp500.fig11.pred = predict(sp500.fig11, n.predict=75) sp500.fieg11.pred = predict(sp500.fieg11, n.predict=75) sp500.fieg11.lev.pred = predict(sp500.fieg11.lev, n.predict=75) td = timeSeq(from="07/01/2003", by="bizdays", length.out=75) sp500.vol.fcst.lm = timeSeries(pos=td, data=cbind(seriesData(sp500.EWMA.pred), sp500.2comp2.pred$sigma.pred, sp500.fig11.pred$sigma.pred/100, sp500.fieg11.lev.pred$sigma.pred/100)) colIds(sp500.vol.fcst.lm) = c("EWMA", "2 COMP", "FIGARCH", "FIEGARCH") plot(sp500.vol.fcst.lm, reference.grid=F) # comparison of predictions. Asymmetric models are very different plot(sp500.vol.fcst[,c(2,6)], sp500.vol.fcst.lm[,c(2,3,4)], reference.grid=F, main="Volatility Forecasts for S&P 500 from Competing GARCH Models", plot.args=list(lty=1:5, col=1:5, lwd=2)) legend(0.6, 0.0103, legend=colIds(sp500.vol.fcst), lty=1:7, col=1:7, lwd=2) # simulation-based forecasts # ARCH(5) # these seem to consistently underpredict and start at a lower value # than the analytic predictions sigma.start = as.numeric(sp500.a5$sigma.t[n.obs]) eps.start = as.numeric(sp500.a5$residuals[n.obs]) eps.start = matrix(eps.start, 1, 1000) set.seed(111) error = rbind(eps.start, matrix(rnorm(25*1000), 25)) # different simulations are the columns of the following matrix set.seed(10) sp500.a5.sim.pred = simulate(sp500.a5, n=25, n.rep=1000, n.start=0, etat=error, sigma.start=sigma.start)$sigma.t sp500.a5.sim.sigma.pred = rowMeans(sp500.a5.sim.pred) cbind(sp500.a5.pred$sigma.pred, sp500.a5.sim.sigma.pred) tsplot(sp500.a5.pred$sigma.pred, sp500.a5.sim.sigma.pred) # GARCH(1,1) sigma.start = as.numeric(sp500.g11$sigma.t[n.obs]) eps.start = as.numeric(sp500.g11$residuals[n.obs]) eps.start = matrix(eps.start, 1, 1000) set.seed(111) error = rbind(eps.start, matrix(rnorm(25*1000), 25)) set.seed(10) sp500.g11.sim.pred = simulate(sp500.g11, n=25, n.rep=1000, n.start=0, etat=error, sigma.start=sigma.start)$sigma.t # incorrect forecasts sp500.g11.sim.sigma.pred = rowMeans(sp500.g11.sim.pred) # correct forecasts sp500.g11.sim.sigma.pred2 = sqrt(rowMeans(sp500.g11.sim.pred^2)) cbind(sp500.g11.pred$sigma.pred, sp500.g11.sim.sigma.pred, sp500.g11.sim.sigma.pred2) tsplot(sp500.g11.pred$sigma.pred, sp500.g11.sim.sigma.pred, sp500.g11.sim.sigma.pred2) # GARCH(1,1) with leverage sigma.start = as.numeric(sp500.pgarch.2$sigma.t[n.obs]) eps.start = as.numeric(sp500.pgarch.2$residuals[n.obs]) eps.start = matrix(eps.start, 1, 1000) set.seed(111) error = rbind(eps.start, matrix(rnorm(25*1000), 25)) set.seed(10) sp500.pgarch.2.sim.pred = simulate(sp500.pgarch.2, n=25, n.rep=1000, n.start=0, etat=error, sigma.start=sigma.start)$sigma.t sp500.pgarch.2.sim.sigma.pred2 = sqrt(rowMeans(sp500.pgarch.2.sim.pred^2)) sp500.pgarch.2.sim.sigma.pred = rowMeans(sp500.pgarch.2.sim.pred) cbind(sp500.pgarch.2.pred$sigma.pred, sp500.pgarch.2.sim.sigma.pred) tsplot(sp500.pgarch.2.pred$sigma.pred, sp500.pgarch.2.sim.sigma.pred, sp500.pgarch.2.sim.sigma.pred2) # TODO: simulation based forecasts from egarch, two component and garch(2,2) # # Volatility term structures # # TODO: volatility term structures # # GARCH models and VaR calculations # # TODO: GARCH VaR calculations # # estimation of GARCH models adjusted for outliers # # # Structural stability of GARCH estimates - rolling analysis # # TODO: rolling analysis of garch