# Time Series ARIMA GARCH Models in R
# Copyright 2020 Agust by Tolga Omay
#install.packages(c("rugarch","rmgarch","MASS")) # only needed in case you have not yet installed these packages
library(rugarch)
library(rmgarch)
library(MASS)
library(tseries)
plot(SP500, type = 'l')
# we took SP500 from MASS library
data(SP500)
plot(SP500, type = 'l')
par(mfrow=c(1,2))
acf(SP500)
pacf(SP500)
#View(SP500)
y <- SP500
adf.test(y, alternative="stationary")
arima(y, order = c(1,0,0), include.mean = FALSE)
fit1 <- arima(y, order = c(1,0,0))
jarque.bera.test(residuals(fit1))
res <- fit1$residuals
Box.test(res, lag=15)
Box.test(res, lag=1, type = "Ljung")
#install.packages("aTSA")
library(aTSA)
arch.test(fit1)
par(mfrow=c(1,1))
hist(y)
summary(y)
res2 <- res*res
par(mfrow=c(1,2))
acf(res2)
pacf(res2)
#When we are estimating volatility models we work with returns.
#There is a function that transforms the data to returns.
#Here we are using the functionality provided by the rugarch package written by Alexios Galanos
#The first thing you need to do is to ensure you know what type of GARCH model
#you want to estimate and then let R know about this. It is the
#ugarchspec( ) function which is used to let R know about the model type.
#There is in fact a default specification and the way to invoke this is as follows
ug_spec = ugarchspec()
#ug_spec is now a list which contains all the relevant model specifications.
#Let's look at them:
ug_spec
#The key issues here are the spec for the Mean Model (here an ARMA(1,1) model) and the specification for the GARCH Model, here an sGARCH(1,1) which is basically a GARCH(1,1). To get details on all the possible specifications and how to change them it is best to consult the documentation of the rugarch package.
#Let's say you want to change the mean model from an ARMA(1,1) to an ARMA(1,0), i.e. an AR(1) model.
ug_spec <- ugarchspec(mean.model=list(armaOrder=c(2,1)))
ug_spec
ug_spec <- ugarchspec(variance.model=list(model="sGARCH", garchOrder=c(1,1)))
ug_spec
ug_spec <- ugarchspec(mean.model=list(armaOrder=c(1,0)), variance.model=list(model="sGARCH", garchOrder=c(1,1)))
ug_spec
#ewma_spec = ugarchspec(variance.model=list(model="iGARCH", garchOrder=c(1,1)),
# mean.model=list(armaOrder=c(0,0), include.mean=TRUE),
#
#Odistribution.model="norm", fixed.pars=list(omega=0))
#ewma_spec
ugfit = ugarchfit(spec = ug_spec, data = y)
#fit is now a list that contains a range of results from the estimation.
#Let's have a look at the results
ugfit
#GARCH models you will recognise some of the parameters.
#ar1 is the AR1 coefficient of the mean model
#(here very small and basically insignificant),
#alpha1 is the coefficient to the squared residuals in the GARCH equation
#and beta1 is the coefficient to the lagged variance.
ugfit@fit$coef
ug_var <- ugfit@fit$var # save the estimated conditional variances
ug_res2 <- (ugfit@fit$residuals)^2 # save the estimated squared residuals
#Let's plot the squared residuals and the estimated conditional variance:
plot(ug_res2, type = "l")
lines(ug_var, col = "green")
View(ug_var)
#chartSeries(y)
#Explanation of the test
#GARCH models provide a way of modelling conditional volatility. I.e. They are useful in situations where the volatility of a time series is a function of previous levels of volatility AKA volatility clustering.
#A GARCH model is typically of the following form:
# \sigma_t^2 = \omega+ \sum_{i=1}^{p} \alpha_i\epsilon_{t-i}^2+\sum_{i=1}^{q} \beta_i\sigma_{t-i}^2
#which means that the variance ( \sigma_t^2) of the time series today is equal to a constant (\omega), plus some amount ( \alpha) of the previous residual (\epsilon_{t-1}), plus some amount ( \beta) of the previous variance (\sigma^2_{t-1}).
#2. FITTING
#Fitting GARCH models is usually trivial using modern software such as the rugarch package for R. However checking whether the fitted model is any good is less trivial since a range of diagnostic tests must be applied and most importantly understood to ensure that the model captures the intended behaviour.
#3. GARCH DIAGNOSTICS
#LJUNG-BOX TEST
#The Ljung-Box test provides a means of testing for auto-correlation within the GARCH model’s standardized residuals. If the GARCH model has done its job there should be NO auto-correlation within the residuals. The null-hypothesis of the Ljung-Box test is that the auto-correlation between the residuals for a set of lags k = 0. If at least one auto-correlation for a set of lags k > 0 then the test statistic indicates that the null-hypothesis may be rejected.
#The rugarch package for R applies a weighted Ljung-Box Test on the standardized Residuals and the standardized Squared Residuals. Should the p-value be <= 0.05 (your significance level \alpha) then the null hypothesis should be rejected meaning that the GARCH model has not captured the auto-correlation.
#ARCH LM TEST
#Similar to the Ljung-Box Test, the ARCH LM test provides a means of testing for serial dependence (auto-correlation) due to a conditional variance process by testing for auto-correlation within the squared residuals.The null hypothesis is that the auto-correlation between the residuals for a set of lags k = 0.
#NYBLOM STABILITY TEST
#The Nyblom stability test provides a means of testing for structural change within a time series. A structural change implies that the relationship between variables changes overtime e.g. for the regression y= \beta x beta changes over time. The null hypothesis is that the parameter values are constant i.e. zero variance, the alternative hypothesis is that their variance > 0.
#SIGN BIAS TEST
#Engle & Ng sign bias tests provide a means of testing for mispecification of conditional volatility models. Specifically they examine whether the standardized squared residual is predictable using (dummy) variables indicative of certain information.
#The sign bias test has a dummy variable S^- that is 1 when \epsilon_{t-1}<0. It tests for the impact of positive & negative shocks on volatility not predicted by the model.
#The negative sign bias test uses a dummy variable S^- \epsilon_{t-1} – it focuses on the effect of large and small negative shocks.
#The positive sign bias test uses a dummy variable S^+ \epsilon_{t-1} where S^+ = 1-S^- – it focuses on the effect of large and small positive shocks.
#The null hypothesis for these tests is that additional parameters corresponding to the additional (dummy) variables = 0. The alternative hypothesis is that the addition parameters are non-zero indicating mispecification of the model.
#ADJUSTED PEARSON GOODNESS-OF-FIT TEST
#The adjusted pearson goodness-of-fit test compares the empirical distribution of the standardized residuals with the selected theoretical distribution. The null hypothesis is that the empirical and theoretical distribution is identical.