AIB Stock Price, EGARCH-M, and rgarch

Posted on May 17, 2011 by timeseriesireland in R bloggers, Uncategorized | 0 Comments

[This article was first published on TimeSeriesIreland » R, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

This post examines conditional heteroskedasticity models in the context of daily stock price data for Allied Irish Banks (AIB), specifically how to test for conditional heteroskedasticity in a series, how to approach model specification and estimation when time-varying volatility is present, and how to forecast with these models; all of this is done in R, with relevant R code provided at the bottom of the post.

Now follows a necessary disclaimer: This is not intended to be financial advice, and should not in any way be interpreted as such; this is purely an academic post, and the subtle irony of analysing price data on a bank that no longer trades on the ISEQ index should not go unnoticed by the reader. I am not responsible for any losses, financial or otherwise, incurred by anyone who dives head-first into stock/options trading with no understanding and/or appreciation of the financial risks involved, especially with any leveraged investments/products. By continuing to read this post you agree to release me of any legal liability and responsibility for your actions based on the content of this post.

Vanilla Black-Scholes-Merton option pricing formulae are based on underlying asset prices being driven by geometric Brownian motion (GBM), which implies log-normally distributed prices and normally distributed returns. The strict assumptions imposed by GBM are often a point of amusement by those looking in from outside of finance and economics, but it’s worth noting that the field does not end with this simple model. Forecasting implied volatility from pricing data is a point of interest in applied financial econometrics, and a class of models utilised for this purpose are conditional heteroskedasticity models. For those not familiar with this term, ‘heteroskedasticity’ simply means non-constant variance, where the case of constant variance called ‘homoskedasticity’.

If a stock, $S$ , follows a random walk (without drift), we can denote this process by $s_t = s_{t-1} + \epsilon_t$ , where lower case letters denote logs and $\epsilon_t \sim IID(0, \sigma_{\epsilon}^2)$ . This implications of this model are ubiquitious: our best guess of tomorrow’s stock price, $E_t(s_{t+1})$ , is today’s price $s_t$ , as $E_t (\epsilon_{t+1}) = 0$ . If we take first-differences, $\Delta s_t = \epsilon_t$ , where the log change approximates the growth rate in the stock price, or the return, we can re-write as $r_t = \epsilon_t$ , where $r_t$ is the return at time $t$ . This implies that the return of an asset is an IID process with mean zero and variance $\sigma_{\epsilon}^2$ . Tests of the random walk hypothesis come in various form, and also depend on how strict our assumptions are for the model, for example $\epsilon_t \sim IID(0, \sigma_{\epsilon}^2)$ . Independent increments in the error term can be easily tested by means of Box-Ljung Q-statistics.

Some other features of the random walk model without drift are worth noting: $s_{t-1} = s_{t-2} + \epsilon_{t-1}$ , so $s_t = (s_{t-2} + \epsilon_{t-1}) + \epsilon_t$ . Continuing with this process of backward-iteration, we see that $s_t = s_{0} + \sum_{j=1}^{t} \epsilon_{j}$ . The mean is given by $E(s_t) = s_0$ , as $\epsilon_t \sim (0,\sigma_{\epsilon}^2)$ by definition, and the variance is $V(s_t) = E[ s_t - E(s_t) ]^2$ , or $V(s_t) = E [\sum_{j=1}^{t} \epsilon_{j}]^2 = t \sigma_{\epsilon}^2$ .

Moving to an empirical example, below is a graph of Allied Irish Banks’ (AIB) daily stock price, from the beginning of 2003 to the end of 2010:

For those unfamiliar with this bank, it is essentially a government-owned entity of the Irish State, where government ownership occurred due to large losses on loans extended to property developers. Below is a graph of the log daily returns of the stock:

The marked drop in the series occurred on January 19, 2009; an astonishing 88% fall over one trading day. Volatility clustering is evident in the series, with relatively higher volatility from the middle of 2008 than the comparatively tranquil 2003 to 2008 period.

Our first task is to identify the correct (mean) model specification, here from the autoregressive moving average (ARMA) suite of models. An ARMA(p,q) model is given by

$r_t = \delta + \sum_{i = 1}^p \phi_i r_{t - i} + \epsilon + \sum_{j = 1}^q \theta_{j} \epsilon_{t - j}$ ,

where I’ve included a drift term, $\delta$ , for maximum generality. When dealing with time-varying volatility, Engle’s autoregressive conditional heteroskedasticity (ARCH) model, and its generalised cousin, Bollerslev’s generalised ARCH (GARCH), are staples of model financial econometrics and are frequently deployed to forecast volatility in financial time-series. If we take the simplest form of the ARMA class, ARMA(0,0),

$r_t = \epsilon_t$ ,

we specify the error term as

$\epsilon_t = \nu_t h_t$ , $\nu_t \sim IID(0,1)$ ,

where

$h_t^2 = \omega + \alpha_1 \epsilon_{t-1}^2$ .

$h_t^2$ is the conditional variance of $\epsilon_t$ . This is an ARCH(1) model, which easily generalises to ARCH(q),

$h_t^2 = \omega + \sum_{k=1}^q \alpha_k \epsilon_{t-k}^2$

The GARCH(p,q) model is given by,

$h_t^2 = \omega + \sum_{k=1}^q \alpha_k \epsilon_{t-k}^2 + \sum_{l=1}^p \beta_l h_{t-l}^2$

where we have simply added autoregressive terms to the conditional variance process. Low-order GARCH(p,q) models have been shown to fit data better than high-order ARCH(q) models. The GARCH(1,1) specification, for example, is frequently used in finance, and one can find this model discussed in popular options textbooks, such as John Hull’s Options, Futures, and Other Derivatives.

If ARCH effects are not present in a series, we would expect past squared values of the residuals to be insignificant covariates of $\epsilon_t^2$ , today’s squared residual. Thus, a simple and popular test for ARCH effects is to, first, estimate an ARMA model, then take the residuals, $\widehat{\epsilon_{t}}$ , where the circumflex denotes an estimate, square them and estimate an autoregressive model of order q for the squared residual series, including a constant.

To identify the appropriate ARMA model, the autocorrelation function (ACF) and partial autocorrelation (PACF) of the log-returns series are given below.

We see significant sample lags out to 50 lags in both the ACF and PACF, a rather unfortunate fact that will hinder estimation of a parsimonious model. (The erudite reader might point out that this is a candidate for fractional integration, with the series exhibiting so-called long-memory, even after integer-order differencing.) It’s arbitrary, really, what model is chosen here, due to the indefinite form of the ACF and PACF graph above. However, the oscillatory nature of the ACF and positive to negative move from the first to second sample lag in the PACF are indicative of a pure AR model; though, again, any model can be justified given the ACF and PACF above. I have selected an AR(9) model, based on minimising information criteria and analysing the ACF of the resulting residuals.

With the mean model estimated to control for linear dependence, we now look to the squared residuals produced by fitting our AR(9) model. Two tests are commonly used to check for ARCH effects. The first is to estimate an autoregressive model of the squared residuals, and employ a test of joint insignificance for all the coefficients on lagged values in the model, i.e., $\widehat{\epsilon_t}^2 = \alpha_0 + \alpha_1 \widehat{\epsilon}_{t-1}^2 + \cdots + \alpha_m \widehat{\epsilon}_{t-m}^2 + error$ . The null hypothesis is such that $\alpha_1 = \alpha_2 = \cdots = \alpha_m = 0$ . An ‘F-test’ is automatically printed along with a linear regression in most (good) statistical packages, and R being an awesome statistical package simple use of the ‘lm’ function will suffice. A regression of the form $\widehat{\epsilon_t}^2 = \alpha_0 + \sum_{i = 1}^8 \alpha_i \widehat{\epsilon}_{t-i}^2 + error$ gives an F-statistic of 15.48, with 8 and 2048 degrees of freedom (a p-value of roughly zero), so we reject the null of no ARCH effects. The second test involves the Box-Ljung Q-statistics to test for independence in the squared residuals: we reject the null of independence at 8 lags, with a $\chi_{8}^2 = 170.2$ , a p-value of roughly zero.

The PACF of the squared residuals can be used to identify the order of an ARCH model, which is given below:

This implies a high-order ARCH(q) model, with significant sample lags at 19 and beyond. Hence, the GARCH model is preferred to estimating an AR(9)-ARCH(19) model. Computationally, this may be a non-trivial exercise and convergence may not occur depending on the suite of algorithms used. However, identifying the order of a GARCH model is essentially a guess-and-go process, with GARCH(1,1), GARCH(1,2), GARCH (2,2) (and higher) being plausible specifications. One could use information criteria here to determine the correct model specification, though some authors do caution on the exact meaning of these for GARCH processes.

Now that we have identified the presence of ARCH effects, and determined that GARCH is a preferable approach than pure ARCH, we proceed to estimate our ARMA-GARCH model. This can be easily achieved in the ‘fGarch’ package, which is part of the wider Rmetrics project. However, I will use ‘rgarch’, a relatively more flexible (beta) package available on R-Forge, that also allows for estimation of GARCH in mean models (GARCH-M) and asymmetric GARCH specifications. Whatever model one estimates, the standardised residuals produced by the estimated model should be white noise, if correctly specified.

Before I estimate the model, it’s interesting to note that the relatively more volatile period in the log-return series occurred when the stock price was falling. This can be incorporated by using a GARCH in mean, or GARCH-M, model, where the condicational variance appears in the mean equation; e.g.,

$r_t = \delta + \mu h_t^2 + \sum_{i = 1}^p \phi_i r_{t - i} + \epsilon + \sum_{j = 1}^q \theta_{j} \epsilon_{t - j}$ ,

$\epsilon_t = \nu_t h_t$ , $\nu_t \sim IID(0,1)$ ,

and

$h_t^2 = \omega + \sum_{k=1}^q \alpha_k \epsilon_{t-k}^2 + \sum_{l=1}^p \beta_l h_{t-l}^2$ .

The $\mu$ coefficient is interpreted as a risk premium parameter. Our a priori expectation is that this will be negative, i.e., as the return increases, volatility decreases, and vice versa. I’ve also chosen to use an exponential GARCH model, or EGARCH, to account for asymmetric effects in the return series by means of a parameter $\gamma$ that allows for a leverage effect in the lagged residuals in $h_t^2$ , which can now be rewritten as:

$\ln(h_t^2) = \omega + \sum_{k=1}^q \alpha_k \frac{ | \epsilon_{t-k} | + \gamma_k \epsilon_{t-k}}{h_{t - k}} + \sum_{l=1}^p \beta_l \ln (h_{t-l}^2)$ .

I estimate an AR(9)-EGARCH(4,2)-M model based on several factors. First, the stability of the model is undermined with higher-order speficications, noticeably so for any EGARCH(5,q). The coefficients and their respective standard errors are given in the tables below, the first table being the mean equation and the second being the variance equation. Nyblom’s parameter stability test yield a joint test statistic of 3.3267, so we cannot reject the null hypothesis of parameter stability at the 10, 5 or 1% significance levels.

Mean Equation:

	Estimate	Standard Error
$\delta$	0.001166	0.003648
$\phi_1$	0.043452	0.038624
$\phi_2$	-0.063144	0.023096
$\phi_3$	-0.012754	0.025271
$\phi_4$	-0.026951	0.024755
$\phi_5$	-0.064749	0.024209
$\phi_6$	-0.024585	0.048999
$\phi_7$	-0.028065	0.021314
$\phi_8$	0.020022	0.037065
$\phi_9$	0.024686	0.020654
$\mu$	-0.101332	0.176526

Variance equation:

$\omega$	-0.097284	0.633358
$\alpha_1$	-0.136455	0.035170
$\alpha_2$	-0.087314	0.060027
$\gamma_1$	0.373658	0.131616
$\gamma_2$	0.407474	0.165395
$\beta_1$	-0.541935	0.038899
$\beta_2$	0.772521	0.003665
$\beta_3$	0.585102	0.062753
$\beta_4$	0.164900	0.020075

To analyse if the model is correctly specified, we first construct the standardised residuals, $\widetilde{\epsilon}_t$ , where $\widetilde{\epsilon}_t = \epsilon_t/h_t$ . If the standardised residuals are white noise, the Q-statistics should be statistically insignificant, which is the case for our model (see the first table below), so we continue under the assumption that the mean model is correctly specified. Similarly, a test for remaining ARCH effects can be conducted by testing the significance of Q-statistics on the squared standardised residuals. ARCH LM tests are significant for two lags at the 5% significance level, though only marginally so.

Q-Statistics on Standardized Residuals:

	Statistic	p-value
Q(10)	8.361	0.5936
Q(15)	11.348	0.7275
Q(20)	20.644	0.4184

Q-Statistics on Standardized Squared Residuals:

	Statistic	p-value
Q(10)	13.42	0.2013
Q(15)	19.26	0.2024
Q(20)	22.78	0.2996

Below is an empirical density plot of the standardised residuals, which isn’t too far from what we would expect for correct model specification.

Below is the daily return series with 2 conditional standard deviations imposed on either side.

Forecasting is easily done with the rgarch package, using the ‘ugarchforecast’ function. If one re-estimates the model leaving out the last ten observations, evaluation of volatility forecasts are based on four mean loss functions: the mean square error (MSE); mean absolute deviation (MAD); QLIKE; and R2Log. I invite the reader to estimate another model of their choice and compare 10-day ahead forecasts with:

MSE	MAD	QLIKE	R2Log
0.00021026	0.009247	-4.0351	4.0524

One might prefer to use, say, weekly or monthly data rather than daily data, which is usually easier to work with, and higher-order models are rarely necessary with weekly and monthly frequencies.

R Code:

require(tseries)
require(rgarch)
require(urca)
require(ggplot2)
#
a<-read.csv(file="AIB.csv")
#
AIB<-a$Close
#
AIBts<-ts(a$Close)
#
dates <- as.Date(a$Date, "%d/%m/%Y")
#
LAIBts<-log(AIBts)
#
retAIB<-diff(LAIBts)
#
df1<-data.frame(dates,AIB)
#
### Plot AIB stock price:
gg1.1<-ggplot(df1,aes(dates,AIB)) + xlab(NULL) + ylab("AIB Stock Price in Euro")
gg1.2<-gg1.1+geom_line(colour="darkblue") + opts(title="Daily AIB Stock Price")
#
png(filename = "aib1.png", width = 580, height = 400, units = "px", pointsize = 12, bg = "white")
gg1.2
dev.off()
#
### Plot AIB stock price returns:
dates2<-dates[2:length(dates)]
ReturnsAIB<-as.numeric(retAIB)
df2<-data.frame(dates2,ReturnsAIB)
#
gg2.1<-ggplot(df2,aes(dates2,ReturnsAIB)) + xlab(NULL) + ylab("Log Changes")
gg2.2<-gg2.1+geom_line(colour="darkred") + opts(title="Daily AIB Stock Price Return")
#
png(filename = "aib2.png", width = 580, height = 400, units = "px", pointsize = 12, bg = "white")
gg2.2
dev.off()
#
### ACFs and PACFs
png(filename = "acfpacf.png", width = 580, height = 700, units = "px", pointsize = 12, bg = "white")
par(mfrow=c(2,1))
acf(retAIB, main="ACF of AIB Log Returns", lag = 50)
pacf(retAIB, main="PACF of AIB Log Returns", lag = 50)
dev.off()
#
ar9<-arima(retAIB, order=c(9,0,0))
acf(ar9$residuals)
#
ressq<-(ar9$residuals)^2
Box.test(ressq, lag = 8, type = "Ljung-Box")
#
png(filename = "pacfressq.png", width = 580, height = 350, units = "px", pointsize = 12, bg = "white")
pacf(ressq, main="PACF of Squared Residuals", lag = 30)
dev.off()
#
#
retAIB2<-as.numeric(retAIB)
# Note that the GARCH order is revered from what I have discussed above
specm1 <- ugarchspec(variance.model=list(model="eGARCH", garchOrder=c(2,4), submodel = NULL),
		mean.model=list(armaOrder=c(9,0), include.mean=TRUE, garchInMean = TRUE))
fitm1 <- ugarchfit(data = retAIB2, spec = specm1)
fitm1
#
fittedmodel <- fitm1@fit
sigma1<-fittedmodel$sigma
#
df2<-data.frame(dates2,ReturnsAIB,sigma1)
#
gg3.1<-ggplot(df2,aes(dates2)) + xlab(NULL) + ylab("Log Changes")
gg3.2<-gg3.1+geom_line(aes(y = ReturnsAIB, colour="Log Returns")) + opts(title="Daily Log Return with 2 Conditional Standard Deviations")
gg3.3<-gg3.2 + geom_line(aes(y = sigma1*2, colour="2 S.D.")) + geom_line(aes(y = sigma1*-2, colour="2 S.D.")) + scale_colour_hue("Series:") + opts(legend.position=c(.18,0.8))
#
png(filename = "aib3.png", width = 580, height = 400, units = "px", pointsize = 12, bg = "white")
gg3.3
dev.off()
#
fitm2 <- ugarchfit(data = retAIB2,out.sample = 10, spec = specm1)
fitm2
pred <- ugarchforecast(fitm2, n.ahead = 10,n.roll = 1)
pred.fpm <- fpm(pred)
#