Does anything NOT beat the GARCH(1,1)?

In their paper on GARCH model comparison, Hansen and Lunde (2005) present evidence that among 330 different models, and using daily data on the DM/$ rate and IBM stock returns, no model does significantly better at predicting volatility (based on a realized measure) than the GARCH(1,1) model, for an out of sample period of about 250 days (between 1992/1993 for the currency, and 1999/2000 for the stock). In this demonstration, I instead ask if anything does NOT beat the GARCH(1,1) using a range of operational tests such as VaR exceedances, and model comparison using the model confidence set of Hansen, Lunde and Nason (2011) on a loss function based on the conditional quantiles. Using the S&P 500 index, and the latest 1500 days as the out of sample forecast period, I find that there are few models that do not beat the GARCH(1,1)-Normal.

Data Setup

Given that it is the most widely cited and followed equity index in the world, the S&P 500 is a natural choice for testing a hypothesis, and for the purposes of this demonstration the log returns of the SPY index tracker were used. Starting on 23-01-2007, a rolling forecast and re-estimation scheme was adopted, with a moving window size of 1500 periods, and a re-estimation period of 50. The data was obtained from Yahoo finance and adjusted for dividends and splits.

library(rugarch)
library(shape)
library(xts)
library(quantmod)
getSymbols('SPY', from = '1997-01-01')
SPY = adjustOHLC(SPY, use.Adjusted = TRUE)
R = ROC(Cl(SPY), type = 'continuous', na.pad = FALSE)

The following plot highlights the importance of the forecast period under consideration, which includes the credit crisis and an overall exceptional testing ground for risk model backtesting.

plot(Cl(SPY), col = 'black', screens = 1, blocks = list(start.time = paste(head(tail(index(SPY),1500), 1)), end.time = paste(index(tail(SPY, 1))), col = 'whitesmoke'), main = 'SPY', minor.ticks = FALSE)

beat_garch_1

Model Comparison

For the backtest the following GARCH models were chosen:

  • GARCH
  • GJR
  • EGARCH
  • APARCH
  • component GARCH
  • AVGARCH
  • NGARCH
  • NAGARCH

with four distributions, the Normal, Skew-Student (Fernandez and Steel version), Normal Inverse Gaussian (NIG) and Johnson’s SU (JSU). The conditional mean was based on an ARMA(2,1) model, while the GARCH order was set at (1,1) and (2,1), giving a total combination of 64 models. To create the rolling forecasts, the ugarchroll function was used, which has been completely re-written in the latest version of rugarch and now also includes a resume method for resuming from non-converged estimation windows.

model = c('sGARCH', 'gjrGARCH', 'eGARCH', 'apARCH', 'csGARCH', 'fGARCH', 'fGARCH', 'fGARCH')
submodel = c(NA, NA, NA, NA, NA, 'AVGARCH', 'NGARCH', 'NAGARCH')
spec1 = vector(mode = 'list', length = 16)
for (i in 1:8) spec1[[i]] = ugarchspec(mean.model = list(armaOrder = c(2, 1)), variance.model = list(model = model[i], submodel = if (i > 5) submodel[i] else NULL))
for (i in 9:16) spec1[[i]] = ugarchspec(mean.model = list(armaOrder = c(2, 1)), variance.model = list(garchOrder = c(2, 1), model = model[i - 8], submodel = if ((i - 8) > 5) submodel[i - 8] else NULL))
spec2 = vector(mode = 'list', length = 16)
for (i in 1:8) spec2[[i]] = ugarchspec(mean.model = list(armaOrder = c(2, 1)), variance.model = list(model = model[i], submodel = if (i > 5) submodel[i] else NULL), distribution = 'sstd')
for (i in 9:16) spec2[[i]] = ugarchspec(mean.model = list(armaOrder = c(2, 1)), variance.model = list(garchOrder = c(2, 1), model = model[i - 8], submodel = if ((i - 8) > 5) submodel[i - 8] else NULL), distribution = 'sstd')
spec3 = vector(mode = 'list', length = 16)
for (i in 1:8) spec3[[i]] = ugarchspec(mean.model = list(armaOrder = c(2, 1)), variance.model = list(model = model[i], submodel = if (i > 5) submodel[i] else NULL), distribution = 'nig')
for (i in 9:16) spec3[[i]] = ugarchspec(mean.model = list(armaOrder = c(2, 1)), variance.model = list(garchOrder = c(2, 1), model = model[i - 8], submodel = if ((i - 8) > 5) submodel[i - 8] else NULL), distribution = 'nig')
spec4 = vector(mode = 'list', length = 16)
for (i in 1:8) spec4[[i]] = ugarchspec(mean.model = list(armaOrder = c(2, 1)), variance.model = list(model = model[i], submodel = if (i > 5) submodel[i] else NULL), distribution = 'jsu')
for (i in 9:16) spec4[[i]] = ugarchspec(mean.model = list(armaOrder = c(2, 1)), variance.model = list(garchOrder = c(2, 1), model = model[i - 8], submodel = if ((i - 8) > 5) submodel[i - 8] else NULL), distribution = 'jsu')
spec = c(spec1, spec2, spec3, spec4)
cluster = makePSOCKcluster(15)
clusterExport(cluster, c('spec', 'R'))
clusterEvalQ(cluster, library(rugarch))
# Out of sample estimation
n = length(spec)
fitlist = vector(mode = 'list', length = n)
for (i in 1:n) {
tmp = ugarchroll(spec[[i]], R, n.ahead = 1, forecast.length = 1500, refit.every = 50, refit.window = 'moving', windows.size = 1500, solver = 'hybrid', calculate.VaR = FALSE, cluster = cluster, keep.coef = FALSE)
if (!is.null(tmp@model$noncidx)) {
tmp = resume(tmp, solver = 'solnp', fit.control = list(scale = 1), solver.control = list(tol = 1e-07, delta = 1e-06), cluster = cluster)
if (!is.null(tmp@model$noncidx))
fitlist[[i]] = NA
} else {
fitlist[[i]] = as.data.frame(tmp, which = 'density')
}
}

The method of using as.data.frame with option “density”, returns a data.frame with the forecast conditional mean, sigma, any skew or shape parameters for the estimation/forecast window, and the realized data for the period. Because I have made use of a time based data series, the returned data.frame has rownames the forecast dates.

DateMuSigmaSkewShapeShape (GIG)Realized
23/01/20070.000660.004960000.00295
24/01/20070.000580.004870000.008
25/01/20070.000250.00517000-0.01182
26/01/20070.000740.00604000-0.00088
29/01/20070.000870.00587000-0.00056
30/01/20070.000860.005710000.00518
31/01/20070.000630.005670000.00666
01/02/20070.000330.005750000.00599
02/02/20070.00010.00580000.00141
05/02/20070.00010.005640000.00024

In order to compare the models, the conditional coverage test for VaR exceedances of Christoffersen (1998) and the tail test of Berkowitz (2001) were used. The following code calculates the quantiles and applies the relevant tests described above.

vmodels = c('sGARCH(1,1)', 'gjrGARCH(1,1)', 'eGARCH(1,1)', 'apARCH(1,1)', 'csGARCH(1,1)','AVGARCH(1,1)', 'NGARCH(1,1)', 'NAGARCH(1,1)', 'sGARCH(2,1)', 'gjrGARCH(2,1)','eGARCH(2,1)', 'apARCH(2,1)', 'csGARCH(2,1)', 'AVGARCH(2,1)', 'NGARCH(2,1)','NAGARCH(2,1)')
modelnames = c(paste(vmodels, '-N', sep = ''), paste(vmodels, '-sstd', sep = ''), paste(vmodels, '-nig', sep = ''), paste(vmodels, '-jsu', sep = ''))
q1 = q5 = px = matrix(NA, ncol = 64, nrow = 1500)
dist = c(rep('norm', 16), rep('sstd', 16), rep('nig', 16), rep('jsu', 16))
# use apply since nig and gh distributions are not yet vectorized
for (i in 1:64) {
q1[, i] = as.numeric(apply(fitlist[[i]], 1, function(x) qdist(dist[i], 0.01, mu = x['Mu'], sigma = x['Sigma'], skew = x['Skew'], shape = x['Shape'])))
q5[, i] = as.numeric(apply(fitlist[[i]], 1, function(x) qdist(dist[i], 0.05, mu = x['Mu'], sigma = x['Sigma'], skew = x['Skew'], shape = x['Shape'])))
px[, i] = as.numeric(apply(fitlist[[i]], 1, function(x) pdist(dist[i], x['Realized'], mu = x['Mu'], sigma = x['Sigma'], skew = x['Skew'], shape = x['Shape'])))
}
VaR1cc = apply(q1, 2, function(x) VaRTest(0.01, actual = fitlist[[1]][, 'Realized'], VaR = x)$cc.LRp)
VaR5cc = apply(q5, 2, function(x) VaRTest(0.05, actual = fitlist[[1]][, 'Realized'], VaR = x)$cc.LRp)
BT1 = apply(px, 2, function(x) BerkowitzTest(qnorm(x), tail.test = TRUE, alpha = 0.01)$LRp)
BT5 = apply(px, 2, function(x) BerkowitzTest(qnorm(x), tail.test = TRUE, alpha = 0.05)$LRp)
VTable = cbind(VaR1cc, VaR5cc, BT1, BT5)
rownames(VTable) = modelnames
colnames(VTable) = c('VaR.CC(1%)', 'VaR.CC(5%)', 'BT(1%)', 'BT(5%)')

From a summary look at the head and tails of the tables, it is quite obvious that the GARCH(1,1)-Normal model is NOT hard to beat. More generally, the (1,1)-Normal models fare worst, and this is not surprising since financial markets, particularly for this crisis prone period under study, cannot possibly be well modelled without fat tailed and asymmetric distributions. A more interesting result is that the (2,1)-Non Normal models beat (1,1)-Non Normal models on average. Given the widely adopted (1,1) approach almost everywhere, this is a little surprising. It should however be noted that 1500 points where used for each estimation window…with less data, parameter uncertainty is likely to induce noise in a higher order GARCH model.

Top 10

ModelVaR.CC(1%)VaR.CC(5%)BT(1%)BT(5%)
AVGARCH(2,1)-sstd0.010.270.010.01
AVGARCH(2,1)-nig0.170.180.120.03
NAGARCH(2,1)-sstd0.110.180.10.02
gjrGARCH(2,1)-sstd0.60.140.410.07
NAGARCH(2,1)-nig0.250.140.130.04
NAGARCH(2,1)-jsu0.360.140.230.05
gjrGARCH(2,1)-nig0.60.140.450.13
gjrGARCH(1,1)-jsu0.60.140.480.11
gjrGARCH(2,1)-jsu0.60.140.550.12
gjrGARCH(1,1)-nig0.60.140.390.13

Bottom 10

ModelVaR.CC(1%)VaR.CC(5%)BT(1%)BT(5%)
sGARCH(2,1)-N0000
NGARCH(2,1)-N0000
NGARCH(1,1)-N0000
csGARCH(1,1)-N0000
csGARCH(2,1)-N0000
apARCH(2,1)-N0000
sGARCH(1,1)-N0000
apARCH(1,1)-N0000
eGARCH(1,1)-N0000
eGARCH(2,1)-N0000

 

The following plots provide for a clearer (hopefully!) visual representation of the results.
beat_garch_2

The models with conditional distribution the Normal are excluded since none of them are above the 5% significance cutoff. All other barplots start at 5% so that any bar which shows up indicates acceptance of the null of the hypothesis (for both tests this is phrased so that it indicates a correctly specified model). A careful examination of all the models indicates that the GJR (1,1) and (2,1) model from any of the 3 skewed and shaped conditional distributions passes all 4 tests. The SSTD distribution fares worst among the three, with the NIG and JSU appearing to provide an equally good fit. The final plot shows the VaR performance of the GARCH(1,1)-Normal and the GJR(2,1)-NIG for the period.
beat_garch_3

MCS Test on Loss Function

It is often said that the VaR exceedance tests are rather crude and require a lot of data in order to offer any significant insight into model differences. An alternative method for evaluating competing models is based on the comparison of some relevant loss function using a test such as the Model Confidence Set (MCS) of Hansen, Lunde and Nason (2011). In keeping with the focus on the tail of the distribution, a VaR loss function at the 1% and 5% coverage levels was used, described in González-Rivera et al.(2004) and defined as:
$${Q_{loss}} \equiv {N^{ – 1}}\sum\limits_{t = R}^T {\left( {\alpha  – 1\left( {{r_{t + 1}} < VaR_{t + 1}^\alpha } \right)} \right)} \left( {{r_{t + 1}} – VaR_{t + 1}^\alpha } \right)$$
where \( P = T – R \) is the out-of-sample forecast horizon, \( T \) the total horizon to include in estimation, and \( R \) the start of the out-of-sample forecast. This is an asymmetric loss function, linearly penalizing exceedances more heavily by \( (1-\alpha) \). This is not yet exported in rugarch but can still be called as follows:

LossV5 = apply(q5, 2, function(x) rugarch:::.varloss(0.05, fitlist[[1]][, 'Realized'], x))
LossV1 = apply(q1, 2, function(x) rugarch:::.varloss(0.01, fitlist[[1]][, 'Realized'], x))

The MCS function is also not yet included, but may be in a future version. The following barplot shows the result of running the MCS test using a stationary bootstrap with 5000 replications.
beat_garch_4

 

beat_garch_5

 

While no models are excluded from the model confidence set, at the 95% significance level, the ranking is still somewhat preserved for the top and bottom models, with the (2,1) NIG and JSU models being at the top. Based on this test, it would appear that the eGARCH(2,1)-NIG is the superior model at both coverage rates, with the GARCH(1,1)-Normal being at the very bottom with a 7% and 20% probability of belonging to the set at the 1% and 5% coverage rates respectively.

Conclusion

The normality assumption does not realistically capture the observed market dynamics, and neither does the GARCH(1,1)-N model. There were few models which did not beat it in this application. In fact, higher order GARCH models were shown to provide significant out performance on a range of tail related measures, and distributions such as the NIG and JSU appeared to provide the most realistic representation of the observed dynamics.

References

Berkowitz, J. (2001). Testing density forecasts, with applications to risk management. Journal of Business & Economic Statistics, 19(4), 465-474.
Christoffersen, P. F. (1998). Evaluating Interval Forecasts. International Economic Review, 39(4), 841-62.
González-Rivera, G., Lee, T. H., & Mishra, S. (2004). Forecasting volatility: A reality check based on option pricing, utility function, value-at-risk, and predictive likelihood. International Journal of Forecasting, 20(4), 629-645.
Hansen, P. R., & Lunde, A. (2005). A forecast comparison of volatility models: does anything beat a GARCH (1, 1)?. Journal of applied econometrics, 20(7), 873-889.
Hansen, P. R., Lunde, A., & Nason, J. M. (2011). The model confidence set. Econometrica, 79(2), 453-497.

The GARCH-DCC Model and 2-stage DCC(MVT) estimation.

This short demonstration illustrates the use of the DCC model and its methods using the rmgarch package, and in particular an alternative method for 2-stage DCC estimation in the presence of the MVT distribution shape (nuisance) parameter. The theoretical background and representation of the model is detailed in the package’s vignette. The dataset and period used are purely for illustration purposes.

library(rmgarch)
library(parallel)
data(dji30retw)
Dat = dji30retw[, 1:10, drop = FALSE]
# define a DCCspec object: 2 stage estimation should usually always use
# Normal for 1-stage (see below for
xspec = ugarchspec(mean.model = list(armaOrder = c(1, 1)), variance.model = list(garchOrder = c(1,1), model = 'eGARCH'), distribution.model = 'norm')
uspec = multispec(replicate(10, xspec))
spec1 = dccspec(uspec = uspec, dccOrder = c(1, 1), distribution = 'mvnorm')
spec1a = dccspec(uspec = uspec, dccOrder = c(1, 1), model='aDCC', distribution = 'mvnorm')
spec2 = dccspec(uspec = uspec, dccOrder = c(1, 1), distribution = 'mvlaplace')
spec2a = dccspec(uspec = uspec, dccOrder = c(1, 1), model='aDCC', distribution = 'mvlaplace')

Since multiple DCC models are being estimated on the same dataset with the same first stage dynamics, we can estimate the first stage once and pass it to the dccfit routine (rather than re-estimating it every time dccfit is called):

cl = makePSOCKcluster(10)
multf = multifit(uspec, Dat, cluster = cl)

Next, the second stage of the DCC model is estimated.

fit1 = dccfit(spec1, data = Dat, fit.control = list(eval.se = TRUE), fit = multf, cluster = cl)
fit1a = dccfit(spec1a, data = Dat, fit.control = list(eval.se = TRUE), fit = multf, cluster = cl)
fit2 = dccfit(spec2, data = Dat, fit.control = list(eval.se = TRUE), fit = multf, cluster = cl)
fit2a = dccfit(spec2a, data = Dat, fit.control = list(eval.se = TRUE), fit = multf, cluster = cl)

To fit the DCC (MVT) model in practice, one either assumes a first stage QML, else must jointly estimate in 1 stage the common shape parameter. In the example that follows below, an alternative approach is used to approximate the common shape parameter.

# First Estimate a QML first stage model (multf already estimated). Then
# estimate the second stage shape parameter.
spec3 = dccspec(uspec = uspec, dccOrder = c(1, 1), distribution = 'mvt')
fit3 = dccfit(spec3, data = Dat, fit.control = list(eval.se = FALSE), fit = multf)
# obtain the multivariate shape parameter:
mvt.shape = rshape(fit3)
# Plug that into a fixed first stage model and iterate :
mvt.l = rep(0, 6)
mvt.s = rep(0, 6)
mvt.l[1] = likelihood(fit3)
mvt.s[1] = mvt.shape
for (i in 1:5) {
xspec = ugarchspec(mean.model = list(armaOrder = c(1, 1)), variance.model = list(garchOrder = c(1,1), model = 'eGARCH'), distribution.model = 'std', fixed.pars = list(shape = mvt.shape))
spec3 = dccspec(uspec = multispec(replicate(10, xspec)), dccOrder = c(1,1), distribution = 'mvt')
fit3 = dccfit(spec3, data = Dat, solver = 'solnp', fit.control = list(eval.se = FALSE))
mvt.shape = rshape(fit3)
mvt.l[i + 1] = likelihood(fit3)
mvt.s[i + 1] = mvt.shape
}
# Finally, once more, fixing the second stage shape parameter, and
# evaluating the standard errors
xspec = ugarchspec(mean.model = list(armaOrder = c(1, 1)), variance.model = list(garchOrder = c(1,1), model = 'eGARCH'), distribution.model = 'std', fixed.pars = list(shape = mvt.shape))
spec3 = dccspec(uspec = multispec(replicate(10, xspec)), dccOrder = c(1, 1), distribution = 'mvt', fixed.pars = list(shape = mvt.shape))
fit3 = dccfit(spec3, data = Dat, solver = 'solnp', fit.control = list(eval.se = TRUE), cluster = cl)

The plot in the change of the likelihood and shape shows that only a few iterations are necessary to converge to a stable value.
dccmvt_1

The value of the shape parameter implies an excess kurtosis of 1.06. The exercise is repeated for the asymmetric DCC (MVT) model.

xspec = ugarchspec(mean.model = list(armaOrder = c(1,1)), variance.model = list(garchOrder = c(1,1), model = "eGARCH"),  distribution.model = "norm")
spec3a  = dccspec(uspec = multispec( replicate(10, xspec) ), dccOrder = c(1,1), distribution = "mvt", model="aDCC")
fit3a = dccfit(spec3a, data = Dat, fit.control = list(eval.se=FALSE), fit = multf)
# obtain the multivariate shape parameter:
mvtx.shape = rshape(fit3a)
# Plug that into a fixed first stage model and iterate :
mvtx.l = rep(0, 6)
mvtx.s = rep(0, 6)
mvtx.l[1] = likelihood(fit3)
mvtx.s[1] = mvtx.shape
for(i in 1:5){
xspec = ugarchspec(mean.model = list(armaOrder = c(1,1)), variance.model = list(garchOrder = c(1,1), model = "eGARCH"),  distribution.model = "std", fixed.pars = list(shape=mvtx.shape))
spec3a = dccspec(uspec = multispec( replicate(10, xspec) ), dccOrder = c(1,1), model="aDCC", distribution = "mvt")
fit3a = dccfit(spec3a, data = Dat, solver = "solnp", fit.control = list(eval.se=FALSE))
mvtx.shape = rshape(fit3a)
mvtx.l[i+1] = likelihood(fit3a)
mvtx.s[i+1] = mvtx.shape
}
# Finally, once more, fixing the second stage shaoe parameter, and evaluating the standard errors
xspec = ugarchspec(mean.model = list(armaOrder = c(1,1)), variance.model = list(garchOrder = c(1,1), model = "eGARCH"),  distribution.model = "std", fixed.pars = list(shape=mvtx.shape))
spec3a = dccspec(uspec = multispec( replicate(10, xspec) ), dccOrder = c(1,1), model="aDCC", distribution = "mvt", fixed.pars=list(shape=mvtx.shape))
fit3a = dccfit(spec3a, data = Dat, solver = "solnp", fit.control = list(eval.se=TRUE), cluster = cl)

The table below displays the summary of the estimated models, where the stars next to the coefficients indicate the level of significance (*** 1%, ** 5%, * 10%). The asymmetry parameter is everywhere insignificant and, not surprisingly, the MVT distribution provides for the best overall fit (even when accounting for one extra parameter estimated).

##           DCC-MVN   aDCC-MVN    DCC-MVL   aDCC-MVL      DCC-MVT     aDCC-MVT
## a      0.00784*** 0.00639*** 0.00618***  0.0055***   0.00665***    0.00623***
## b      0.97119*** 0.96956*** 0.97624*** 0.97468***   0.97841***    0.97784***
## g                    0.00439               0.00237                 0.00134
## shape                                                9.63947***    9.72587***
## LogLik      22812      22814      22858      22859        23188         23188

To complete the demonstration, the plots below illustrate some of the dynamic correlations from the different models:
dccmvt_2

 

…and remember to terminate the cluster object:

stopCluster(cl)

Rolling GARCH Forecasts

The rugarch package contains a rolling volatility forecast function called ugarchroll, but in this example I will show how easy it is to create a quick custom function. Having a rolling forecast of volatility can prove an invaluable indicator for use in trading systems. an example of which is also included.

# The function makes exclusive use of xts based timeseries indexing.
garchv = function(x, start = 600, roll = 25, n.ahead = 1, spec = NULL, cluster = NULL)
{
# x : is a series of closing prices

# start : is the number of periods which are used for the first rolling
# estimation window

# roll : the number of periods to roll the 1-ahead forecast before
# re-estimating the model parameters

# n.ahead : the forecast horizon

# spec : the univariate GARCH model specification (uGARCHspec object)

# cluster : a pre-created cluster object from the parallel package
require(rugarch)
if (!is.xts(x))
stop('\nx must be an xts object')
# create continuous return series
r = na.omit(diff(log(x), 1))
n = nrow(r)
if (n &lt;= start)
stop('\nstart cannot be greater than length of x')
# the ending points of the estimation window
s = seq(start + roll, n, by = roll)
m = length(s)
# the rolling forecast length
out.sample = rep(roll, m)
# adjustment to include all the datapoints from the end
if (s[m] &lt; n) {
s = c(s, n)
m = length(s)
out.sample = c(out.sample, s[m] - s[m - 1])
}
# the default specification if none is provided is set to an eGARCH-std
# model
if (is.null(spec))
spec = ugarchspec(variance.model = list(model = 'eGARCH'), distribution = 'std')
if (!is.null(cluster)) {
clusterEvalQ(cluster, library(rugarch))
clusterEvalQ(cluster, library(xts))
clusterExport(cluster, varlist = c('r', 's', 'roll', 'spec', 'out.sample'),
envir = environment())
tmp = clusterApply(cl = cluster, 1:m, fun = function(i) {
fit = ugarchfit(spec, r[1:s[i], ], out.sample = out.sample[i], solver = 'hybrid')
f = ugarchforecast(fit, n.ahead = 1, n.roll = out.sample[i] - 1)
y = sigma(f)
# use xts for indexing the forecasts
y = xts(as.numeric(y), tail(time(r[1:s[i], ]), out.sample[i]))
return(y)
})
} else{
tmp = lapply(1:m, FUN = function(i) {
fit = ugarchfit(spec, r[1:s[i], ], out.sample = out.sample[i], solver = 'hybrid')
f = ugarchforecast(fit, n.ahead = 1, n.roll = out.sample[i] - 1)
y = sigma(f)
# use xts for indexing the forecasts
y = xts(as.numeric(y), tail(time(r[1:s[i], ]), out.sample[i]))
return(y)
})
}
# re-assemble the forecasts
forc = tmp[[1]]
for (i in 2:m) {
forc = rbind(forc, tmp[[i]])
}
colnames(forc) = 'V'
return(forc)
}

As an application of a rolling volatility forecast, I will consider the variable moving average (VMA) indicator of Chande (see Technical Analysis of Stocks & Commodities, March 1992), implemented in the TTR package.

getSymbols('SPY', src = 'yahoo', from = '1990-01-01')
SPY = adjustOHLC(SPY, use.Adjusted = TRUE)
library(parallel)
cluster = makePSOCKcluster(10)
garch_V = garchv(Cl(SPY), start = 600, roll = 25, n.ahead = 1, spec = NULL,
cluster = cluster)
stopCluster(cluster)

The VMA function uses a weighting vector, scaled between 0 and 1, with higher values (volatility) leading to faster adjustment of the moving average. Formally, it is defined as $$VM{A_t} = \left( {r{w_t}} \right){P_t} + \left( {1 – r{w_t}} \right)VM{A_{t – 1}}$$ where r is the decay ratio which adjusts the smoothness of the moving average, and w the weighting vector. In order to use the forecast volatility with the VMA indicator we need to rescale it to lie between the required bounds. A simple transformation which uses the min and max values can be used, such that $\bar h=\left(max(h)-h\right)/\left(\max(h)-min(h)\right)$. A slight problem with this approach is that the max and min of the forecast for the period is not known (not without looking ahead), but we ignore this problem at present.

barh = (garch_V - min(garch_V))/(max(garch_V) - min(garch_V))
tmp = cbind(Cl(SPY), barh)
tmp = na.omit(tmp)
VMA_garch1 = VMA(tmp[, 1], w = tmp[, 2]^2, ratio = 0.1)
VMA_garch2 = VMA(tmp[, 1], w = tmp[, 2]^2, ratio = 1)

For illustration, I’ll plot 2 VMA, a slow and fast version, the calibration of which is done through the decay ‘ratio’ option in the VMA function.

par(mfrow = c(2, 1))
plot(garch_V, main = 'SPY Volatility')
plot(tmp[, 1], main = 'SPY w/th VMA Indicator')
lines(VMA_garch1, col = 2)
lines(VMA_garch2, col = 4)

rolling_var_1

As can be visually inferred, the effect of the spikes in volatility on the moving average indicator are clearly illustrated, particularly around the 2008 period. A simple system would probably seek to be long the security when the fast moving average is above the slow moving average, when markets are trending.

How Good Are Your VaR Estimates?

Despite its shortcomings, Value at Risk (VaR) is still the most widely used measure for measuring the risk of a portfolio, and the preferred measure in the Basel II Accord. In this demonstration, I backtest a group of indices using a GARCH-DCC(MVT) model and test the VaR obtained from randomly weighted portfolio margins using a number of relevant tests, based on functions in the rugarch and rmgarch packages. The use of randomly generated weights is in order to avoid any bias from testing a model with only equal weighted margins which is sometimes used. In practise, the weight vector is pre-determined, but for the backtesting, this approach appears a reasonable compromise for evaluating multivariate models using only univariate measures.

library(rmgarch)
library(quantmod)
# International Equity ETFs
Symbols = c('SPY', 'EWC', 'EWW', 'EWA', 'EWH', 'EWJ', 'EWS', 'EWG',
'EWQ', 'EWP', 'EWI', 'EWU', 'EWL', 'EWD')
m = length(Symbols)
Countries = c('USA', 'Canada', 'Mexico', 'Australia', 'Hong.Kong',
'Japan', 'Singapore', 'Germany', 'France', 'Spain', 'Italy', 'UK',
'Switzerland', 'Sweden')
getSymbols(Symbols, src = 'yahoo', from = '1990-01-01')
Y = Ad(SPY)
for (i in 2:m) Y = cbind(R, Ad(get(Symbols[i])))
Y = na.omit(Y)
R = ROC(Y, na.pad = FALSE)

It’s always a good idea to check the data visually for outliers or bad data:

par(mfrow = c(3, 5))
for (i in 1:m) hist(R[, i], breaks = 100, main = Countries[i])

var_estimates_1

Not suprisingly, there is a large spike for the zeros in many countries and this is likely related to low volume when these ETFs where first introduced in the mid 1990s. I’ll exclude about 100 days from the start, since consecutive zeros in a dataset can seriously affect estimation (particularly a rolling backtest which does not use all the sample at once).

R = R[-c(1:100), ]

This is a ‘large’ scale backtest, so full use of parallel functionality is made, which is one of the features of the multivariate GARCH models included in the rmgarch package (the separability of the univariate from the multivariate dynamics).
A 2-stage QML DCC-(MVT) model (with first stage based ARMA(1,1)-gjrGARCH dynamics) is adopted, in the same spirit as in Bauwens and Laurent (2005). A separate post provides for an alternative approach using an iterated procedure for joint shape parameter estimation in the first stage conditional student density.
The model is estimated every 25 days (approx. 1 month), using all the data available up to that point, and the rolling 1-ahead forecast for the next 25 days is evaluated. The backtest procedure starts on 02-06-1999, thus providing an initial window for estimation of size 800, and ends on 10-19-2012 (the code adjusts to the last date available from the downloaded data.)

n = nrow(R)
roll = 25
start = 800
s = seq(start + roll, n, by = roll)
p = length(s)
out.sample = rep(roll, p)
if (s[p] &lt; n) {
out.sample = c(out.sample, n - s[p])
s = c(s, n)
p = length(s)
}

500 sets of random weights (no shorts) are generated with full investment constraint (and for replicability the random seed is initialized to a predefined integer).

set.seed(111111)
w = matrix(rexp(14 * 500), 500, 14)
w = t(apply(w, 1, FUN = function(x) x/sum(x)))
w[500, ] = rep(1/m, m)
# create a required function for use later.
wsig = function(s, w) apply(w, 1, FUN = function(x) sqrt(x %*% s %*% x))

The backtest returns a list of length sum(out.sample)=3288 with elements the 25-rolling 1-ahead conditional forecasts of the mean and sigma weighted by the 500 random weight sets, and a vector of the shape parameter for each forecast window. While 2-stage DCC models offer the possibility of large scale parallel estimation because of the separation of the univariate and multivariate likelihoods, the calculation of the partitioned standard error (s.e.) matrix is quite slow (though this too can take advantage of parallel methods). For this application, s.e. are not calculated.

# set up the parallel estimation
library(parallel)
spec = multispec(replicate(m, ugarchspec(variance.model = list(model
= 'gjrGARCH'))))
mspec = dccspec(spec, distribution = 'mvt')
cl &lt; - makePSOCKcluster(20)
clusterEvalQ(cl = cl, library(rmgarch))
clusterEvalQ(cl = cl, library(xts))
clusterExport(cl, varlist = c('R', 's', 'w', 'mspec', 'wsig', 'out.sample'),
envir = environment())
tic = Sys.time()
mod = clusterApply(cl = cl, 1:p, fun = function(i) {
fit = dccfit(mspec, R[1:s[i], ], out.sample = out.sample[i],
solver = 'solnp', fit.control = list(eval.se = FALSE))
f = dccforecast(fit, n.ahead = 1, n.roll = out.sample[i] - 1)
mu = t(apply(f@mforecast$mu, 3, FUN = function(x) x))
mu = apply(w, 1, FUN = function(x) mu %*% x)
# rcov(DCCforecast) is a list of size n.roll+1, with each
# element being an array of size m by m by n.ahead
sig = t(sapply(rcov(f), FUN = function(x) wsig(x[, , 1], w)))
shape = rep(as.numeric(rshape(fit)), out.sample[i])
dt = as.character(tail(time(R[1:s[i], ]), out.sample[i]))
ans = list(mu = mu, sig = sig, shape = shape, dates = dt)
return(ans)
})
toc = Sys.time()
stopCluster(cl)

The total time taken is about 3.62 mins, but with cloud based computing this could certainly be reduced, given that the functions themselves can take advantage of parallel estimation (not used here) to estimate the univariate models. Next, the forecasts are  re-assembled and the VaR calculated, which because of the location and scaling invariant property of this distribution, is quite fast to evaluate.

mu = sig = shape = Dates = NULL
for (i in 1:length(mod)) {
mu = rbind(mu, mod[[i]]$mu)
sig = rbind(sig, mod[[i]]$sig)
shape = c(shape, mod[[i]]$shape)
Dates = c(Dates, mod[[i]]$dates)
}
# The realized weighted returns
realized = apply(w, 1, FUN = function(x) tail(R, sum(out.sample)) %*% x)
VaR01 = VaR05 = matrix(NA, ncol = 500, nrow = sum(out.sample))
for (i in 1:500) {
VaR01[, i] = mu[, i] + sig[, i] * qdist('std', 0.01, mu = 0, sigma = 1,
shape = shape)
VaR05[, i] = mu[, i] + sig[, i] * qdist('std', 0.05, mu = 0, sigma = 1,
shape = shape)
}

Before proceeding to the tests, do some plots to see some of the results (use the equal weighting margin=500)

par(mfrow = c(2, 1))
rugarch:::.VaRplot('ETFs', 0.01, realized[, 500], Dates, VaR01[, 500])
rugarch:::.VaRplot('ETFs', 0.05, realized[, 500], Dates, VaR05[, 500])

var_estimates_2

The test of VaR exceedances and VaR Duration are used, both of which are implemented and described in detail in the rugarch vignette and documentation. Once again, parallel functionality is used to evaluate these tests. For the exceedances test only the conditional coverage test is considered.

f = function(VaR, realized, p) {
tmp1 = VaRTest(p, realized, VaR)
tmp2 = VaRDurTest(p, realized, VaR)
Exceed = tmp1$actual.exceed
Excp = tmp1$cc.LRp
Excs = tmp1$cc.LRstat
Durb = tmp2$b
Durp = tmp2$LRp
return(c(Exceed = Exceed, Excp = Excp, Excs = Excs, Durb = Durb, Durp = Durp))
}
cl &lt; - makeCluster(20)
clusterEvalQ(cl = cl, library(rugarch))
clusterExport(cl, varlist = c('VaR01', 'VaR05', 'f', 'realized'), envir = environment())
test1pct = clusterApply(cl = cl, 1:500, fun = function(i) f(VaR01[, i], realized[,
i], 0.01))
test5pct = clusterApply(cl = cl, 1:500, fun = function(i) f(VaR05[, i], realized[,
i], 0.05))
stopCluster(cl)
test1pct = t(sapply(test1pct, FUN = function(x) x))
test5pct = t(sapply(test5pct, FUN = function(x) x))
VaRresults = matrix(NA, ncol = 2, nrow = 5)
VaRDurresults = matrix(NA, ncol = 2, nrow = 3)
rownames(VaRDurresults) = c('Median_[b]', 'Median_pvalue', '%Rejections')
rownames(VaRresults) = c('Mean_Exceedances', 'E[Exceedances]', 'Median_stat',
'Median_pvalue', '%Rejections')
colnames(VaRDurresults) = colnames(VaRresults) = c('VaR[1%]', 'VaR[5%]')
VaRresults[1, ] = c(mean(test1pct[, 1]), mean(test5pct[, 1]))
VaRresults[2, ] = c(floor(sum(out.sample) * 0.01), floor(sum(out.sample) * 0.05))
VaRresults[3, ] = c(mean(test1pct[, 3]), median(test5pct[, 3]))
VaRresults[4, ] = c(mean(test1pct[, 2]), median(test5pct[, 2]))
VaRresults[5, ] = c(100 * length(which(test1pct[, 2] &lt; 0.05))/nrow(test1pct),
100 * length(which(test5pct[, 2] &lt; 0.05))/nrow(test5pct))
VaRDurresults[1, ] = c(median(test1pct[, 4]), median(test5pct[, 4]))
VaRDurresults[2, ] = c(median(test1pct[, 5]), median(test5pct[, 5]))
VaRDurresults[3, ] = c(100 * length(which(test1pct[, 5] &lt; 0.05))/nrow(test1pct),
100 * length(which(test5pct[, 5] &lt; 0.05))/nrow(test5pct))

The results summarized below tell an interesting story. The VaR Exceedances test indicates that at the 1% coverage the model used is unable to capture the extremes, with a rejection of the Null (H0: Correct Exceedances & Independent) in the region of 65% across all randomly weighted portfolios. However, when it comes to the 5% coverage the model does quite well with only a 13% rejection rate across all 500 portfolios. The VaR Duration test on the other hand, which measures the independence between exceedances (H0: Duration Between Exceedances have no memory (Weibull b=1 = Exponential)) has a rejection rate of less than 1% for both coverage rates indicating that the model used (ARMA-GARCH-DCC) is adequate for the modelling of the underlying dynamics, and the quantiles in particular. Nevertheless, exceedances for the 1% coverage are high and alternative dynamics may be called for, such as an Autoregressive Conditional Density model, a feasible multivariate example of which is presented in a recent paper by myself and colleagues The IFACD model.

print(VaRresults, digit = 2)
##                VaR[1%] VaR[5%]
## Exceedances     46.930     197
## E[Exceedances]  32.000     164
## Median_stat      6.666       0
## Median_pvalue    0.063       1
## %Rejections     65.400      13
##
##
print(VaRDurresults, digit = 2)
##               VaR[1%] VaR[5%]
## Median_[b]       0.98    1.00
## Median_pvalue    0.64    0.67
## %Rejections      0.20    0.20

References

Bauwens L and Laurent S (2005). “A new class of multivariate skew densities, with application to generalized autoregressive conditional
heteroscedasticity models.” Journal of Business and Economic Statistics, 23(3), pp. 346-354.

GARCH Parameter Uncertainty and Data Size

A frequently asked question relates to the minimum size of a dataset, required to obtain ‘good’ GARCH estimates. In this demonstration, the ugarchdistribution function is used to show how this question can be addressed within the rugarch package and the relevance of $\sqrt{N}$ consistency and how this relates to this question.
First define a GARCH specification

library(rugarch)
library(parallel)
cluster = makePSOCKcluster(10)
spec = ugarchspec(mean.model = list(armaOrder = c(0, 0)))
setfixed(spec)&lt;-list(mu = 1e-04, omega = 3e-06, alpha1 = 0.05, beta1 = 0.92)
sqrt(uncvariance(spec))
## [1] 0.01

Now use this specification to simulate and estimate based on different data sizes.

mod = ugarchdistribution(spec, n.sim = 101, n.start = 1, m.sim = 100,
recursive = TRUE, recursive.length = 3000, recursive.window = 250,
cluster = cluster)
# remember to terminate the cluster
stopCluster(cluster)

The resulting object is of class uGARCHdistribution, with 3 slots:

slotNames(mod)
## [1] 'dist'     'truecoef' 'model'

The dist slot is a list of size equal to the number of estimated windows plus 1 (the last object in the list contains details of the estimation). Each list object contains the estimated parameters per window in addition to other calculated statistics for that window size. Next, we investigate the distribution of the parameters per window:

n = length(mod@dist) - 1
clr = topo.colors(n, alpha = 1)
mu = sapply(mod@dist, FUN = function(x) x$simcoef[, 1])
mu$details = NULL
omega = sapply(mod@dist, FUN = function(x) x$simcoef[, 2])
omega$details = NULL
alpha1 = sapply(mod@dist, FUN = function(x) x$simcoef[, 3])
alpha1$details = NULL
beta1 = sapply(mod@dist, FUN = function(x) x$simcoef[, 4])
beta1$details = NULL
par(mfrow = c(2, 2))
boxplot(na.omit(mu), names = paste('w[', 1:n, ']', sep = ''), col = clr)
abline(h = 1e-04, col = 2)
title('mu')
boxplot(na.omit(omega), names = paste('w[', 1:n, ']', sep = ''), col = clr)
abline(h = 3e-06, col = 2)
title('omega')
boxplot(na.omit(alpha1), names = paste('w[', 1:n, ']', sep = ''), col = clr)
abline(h = 0.05, col = 2)
title('alpha')
boxplot(na.omit(beta1), names = paste('w[', 1:n, ']', sep = ''), col = clr)
abline(h = 0.92, col = 2)
title('beta')

parameter_uncertainty_1

 

As expected, the standard deviation from the true parameter decreases as the window size increases, with the largest errors for the small window sizes. Another way to see this is via the root mean squared error (RMSE) plots:

plot(mod, which = 4)

parameter_uncertainty_2

 

The plots show the RMSE of the fitted versus true coefficients per window size, and the red line the expected RMSE under the assumption of $\sqrt{N}$ consistency. As expected, more (data) leads to less (error). Finally, it it possible to investigate some additional plots showing the distribution of other measures of interest such as persistence, half-life etc:

plot(mod, which = 3)

parameter_uncertainty_3

The ugarchdistribution function makes it easy to investigate the importance of data size on the estimated parameters and provides an initial estimate of the cost of using too little data in GARCH estimation.