We already learned about the basics of ARIMA
For seasonal data, like quarterly and monthly time series, the ARIMA model can be extended to Seasonal ARIMA (SARIMA) models. See https://rpubs.com/SunnyBingoMe/sarima for some sample work.
Another popular extension to ARIMA models is called ARIMAX, implemented by incorporating additional exogenous variables (regressors) that are external to and different from the forecast variable
\[ ARMAX(p,q,r): r_t = \sum_{i=1}^p r_{t-i} + \sum_{j=1}^q \varepsilon_{t-j}+\sum_{k=1}^r \gamma_k X_{tk} + \varepsilon_t \] - where: \(\gamma_k\) are the parameters for the exogenous variables \(X_{tk}\).
# Hypothesis: I believe that AAPL daily volume can explain AAPL future prices
psych::describe(s1v)
psych::describe(log(s1v))
plot(s1v)
plot(log(s1v))
hist(log(s1v))
cor(s1p, log(s1v))
cor.test(s1p, log(s1v))
# Model Estimation
arima <- auto.arima(s1p)
arima
arima.x <- auto.arima(s1p, xreg = log(s1v) )
arima.x# Setting Time Parameters
tt <- dim(mk) # Total Time Period
hh <- 200 # Out-sample Period for Prediction
# ARIMA
mk.arima.1 <- auto.arima(mk[1:(tt-hh),])
mk.arima.1
mk.arima.f <- forecast(mk.arima.1, h = hh)
# ARIMAX with VIX
mk.arima.2 <- auto.arima(mk[1:(tt-hh),], xreg = vix[1:(tt-hh),])
mk.arima.2
mk.arima.f2 <- forecast(mk.arima.2, h = hh, xreg = vix[(tt-hh+1):tt])
# Forecast Plot
par(mfrow=c(1,3))
plot(mk.arima.f)
plot(mk.arima.f2)
plot(mk)# Setting Time Parameters
tt <- dim(s1p) # Total Time Period
hh <- 100 # Out-sample Period for Prediction
# ARIMA
s1p.arima.1 <- auto.arima(s1p[1:(tt-hh),])
s1p.arima.1
s1p.arima.f <- forecast(s1p.arima.1, h = hh)
# ARIMAX with VIX
s1p.arima.2 <- auto.arima(s1p[1:(tt-hh),], xreg = log(s1v)[1:(tt-hh),])
s1p.arima.2
s1p.arima.f2 <- forecast(s1p.arima.2, h = hh, xreg = log(s1v)[(tt-hh+1):tt])
# Forecast Plot
par(mfrow=c(1,2))
plot(s1p.arima.f)
plot(s1p.arima.f2)\[ r_t = ARMA/ARMAX + \sigma_t Z_t \\ GARCH(p, q, r) : \sigma_t^2 = \omega + \sum_{i=1}^p \alpha_i \epsilon_{t-i}^2 + \sum_{j=1}^q \beta_j \sigma_{t-j}^2 + \sum_{k = 1}^r \gamma_k X{tk} \]
plot(vix)
dim(vix)
dim(vix)
?ugarchspec
?ugarchfit
garchs <- ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(1, 1)),
mean.model = list(armaOrder = c(0, 0), include.mean = TRUE))
garchs.res <- ugarchfit(garchs, mkr)
garchx <- ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(1, 1),
external.regressors = vix),
mean.model = list(armaOrder = c(0, 0), include.mean = TRUE))
garchx.res <- ugarchfit(garchx, mkr, xreg = vix)
print(garchs.res)
print(garchx.res)err <- garchx.res@fit[["residuals"]]
# Residual Visualization
par(mfrow = c(2,3))
plot(err)
hist(err, breaks = 50)
qqnorm(err)
acf(err)
pacf(err)
# Ljung-Box Test -> H1: Error is not stationary
Box.test(err, lag = 16, type = "Ljung")
# ADF Test -> H1: Error has no unit root
adf.test(err)
# Jarque-Bera Test -> H1: Error is not normal
jarque.bera.test(err)\[ CME(x)=\mathbb{E}[X-x|X \geq x] \] - Examples of the heavy tail distribution according to this definition include Pareto distribution, Weibull distribution, Cauchy distribution etc.
There’s a whole documentations on probability distributions used in r: https://cran.r-project.org/web/views/Distributions.html
Which model is the best? Why?
Look at the parameters, are they significant?
Are the information criteria (AIC, BIC… ) the lowest?
garchs.n <- ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(1, 1)),
mean.model = list(armaOrder = c(0, 0), include.mean = TRUE),
distribution.model = "norm")
garchs.sn <- ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(1, 1)),
mean.model = list(armaOrder = c(0, 0), include.mean = TRUE),
distribution.model = "snorm")
garchs.t <- ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(1, 1)),
mean.model = list(armaOrder = c(0, 0), include.mean = TRUE),
distribution.model = "std")
garchs.st <- ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(1, 1)),
mean.model = list(armaOrder = c(0, 0), include.mean = TRUE),
distribution.model = "sstd")
garchs.g <- ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(1, 1)),
mean.model = list(armaOrder = c(0, 0), include.mean = TRUE),
distribution.model = "ged")
garchs.sg <- ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(1, 1)),
mean.model = list(armaOrder = c(0, 0), include.mean = TRUE),
distribution.model = "sged")
garchs.h <- ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(1, 1)),
mean.model = list(armaOrder = c(0, 0), include.mean = TRUE),
distribution.model = "ghyp")
garchs.n <- ugarchfit(garchs.n, s1r)
garchs.sn <- ugarchfit(garchs.sn, s1r)
garchs.t <- ugarchfit(garchs.t, s1r)
garchs.st <- ugarchfit(garchs.st, s1r)
garchs.g <- ugarchfit(garchs.g, s1r)
garchs.sg <- ugarchfit(garchs.sg, s1r)
garchs.h <- ugarchfit(garchs.h, s1r)
garchs.n@fit[["coef"]]
garchs.sn@fit[["coef"]]
garchs.t@fit[["coef"]]
garchs.st@fit[["coef"]]
garchs.g@fit[["coef"]]
garchs.sg@fit[["coef"]]
garchs.h@fit[["coef"]]
### Remarks ###
# This is how they compute the information criteria in R
# Check this out -> rugarch:::.information.test
info.test <- function (LLH, nObs, nPars)
{
AIC = (-2 * LLH)/nObs + 2 * nPars/nObs
BIC = (-2 * LLH)/nObs + nPars * log(nObs)/nObs
SIC = (-2 * LLH)/nObs + log((nObs + 2 * nPars)/nObs)
HQIC = (-2 * LLH)/nObs + (2 * nPars * log(log(nObs)))/nObs
informationTests = list(AIC = AIC, BIC = BIC, SIC = SIC, HQIC = HQIC)
return(informationTests)
}\[ r_t = \mu + \sigma_t z_t \]
\[ ln(\sigma^2_t) = \omega + \Sigma_{i=1}^p \alpha_i (|z_{t-i}|-\mathbb{E}[|z_{t-i}|] + \gamma z_{t-i}) + \Sigma_{j=1}^q \beta_j ln (\sigma^2_{t-j}) \]
eg1 <- ugarchspec(variance.model = list(model = "eGARCH", garchOrder = c(1, 1)),
mean.model = list(armaOrder = c(0,0)),
distribution.model = "norm")
eg1fit <- ugarchfit(eg1,s1r)
eg1fiterr <- eg1fit@fit[["residuals"]]
par(mfrow = c(2,3))
plot(err)
hist(err, breaks = 50)
qqnorm(err)
acf(err)
pacf(err)
Box.test(err, lag = 16, type = "Ljung")
adf.test(err) # make sure you have the tseries package to use this function
jarque.bera.test(err)Another volatility model commonly used to handle leverage effects is the threshold GARCH (gjrGARCH or TGARCH) model; see Glosten, Jagannathan, and Runkle (1993) and Zakoian (1994).9 10
If the shock is negative, there will be a leverage effect on the future volatility (negative shock \(\rightarrow\) more volatile than positive shock).
This model looks `cleaner’ than EGARCH but may incur negative \(\sigma\) in some cases.
\[ r_t = \mu + \sigma_t z_t \]
\[ \sigma^2_t = \omega + \Sigma_{i=1}^p (\alpha_i + \gamma_i N_{t-i})z_{t-i}^2 + \Sigma_{j=1}^q \beta_j \sigma_{t-j}^2 \]- Where \(N_{t-i} = 1; z_{t-i}< 0\) and 0 otherwise.
- Again, you can use normal distribution or other fatter-tailed distributions for \(z_t\).
tg1 <- ugarchspec(variance.model = list(model = "gjrGARCH", garchOrder = c(1, 1)),
mean.model = list(armaOrder = c(0,0)),
distribution.model = "norm")
tg1fit <- ugarchfit(tg1,s1r)
tg1fiterr <- tg1fit@fit[["residuals"]]
par(mfrow = c(2,3))
plot(err)
hist(err, breaks = 50)
qqnorm(err)
acf(err)
pacf(err)
Box.test(err, lag = 16, type = "Ljung")
adf.test(err) # make sure you have the tseries package to use this function
jarque.bera.test(err)\[ \sigma_t^\delta = \omega + \Sigma_{i=1}^p \alpha_i (|\varepsilon_{t-i}| - \gamma_i \varepsilon_{t-i} )^\delta + \Sigma_{j=1}^q \beta_j \sigma_{t-j}^\delta \]
ag1 <- ugarchspec(variance.model = list(model = "apARCH", garchOrder = c(1, 1)),
mean.model = list(armaOrder = c(0,0)),
distribution.model = "norm")
ag1fit <- ugarchfit(ag1,s1r)
ag1fiterr <- ag1fit@fit[["residuals"]]
par(mfrow = c(2,3))
plot(err)
hist(err, breaks = 50)
qqnorm(err)
acf(err)
pacf(err)
Box.test(err, lag = 16, type = "Ljung")
adf.test(err) # make sure you have the tseries package to use this function
jarque.bera.test(err)Use “ugarchforecast” to get the future predicted “expected values” of \(r_{t+h}\) and future predicted \(\sigma_{t+h}\) for h steps ahead
Use “ugarchboot” to get the simulated future returns and future volatility
?ugarchforecast
garchs.n.f <- ugarchforecast(garchs.n, n.ahead = 100)
plot(garchs.n.f)
garchs.sn.f <- ugarchforecast(garchs.sn, n.ahead = 100)
plot(garchs.sn.f)“ugarchforecast” will only give you the expected values of \(r_t\) and \(\sigma_t\)
I will use “ugarchboot” function in the rugarch package to make a simulated forecast using empirical distributions.
“ugarchsim” can do similar stuff but slower and buggy (imo)
print(eg1fit)
print(tg1fit)
print(ag1fit)
?ugarchboot
# EGARCH
eg.boot <- ugarchboot(eg1fit, method = "Partial", n.ahead = 10, n.bootpred = 2000)
par(mfrow = c(1,2))
boxplot(eg.boot@fseries)
boxplot(eg.boot@fsigma)
step <- 10
par(mfrow = c(1,2))
hist(eg.boot@fseries[,step])
hist(eg.boot@fsigma[,step])
# gjrGARCH
tg.boot <- ugarchboot(tg1fit, method = "Partial", n.ahead = 10, n.bootpred = 2000)
par(mfrow = c(1,2))
boxplot(tg.boot@fseries)
boxplot(tg.boot@fsigma)
step <- 2
hist(tg.boot@fseries[,step])
hist(ag.boot@fsigma[,step])
# APARCH
ag.boot <- ugarchboot(ag1fit, method = "Partial", n.ahead = 10, n.bootpred = 2000)
par(mfrow = c(1,2))
boxplot(ag.boot@fseries)
boxplot(ag.boot@fsigma)
step <- 2
par(mfrow = c(1,2))
hist(ag.boot@fseries[,step])
hist(ag.boot@fsigma[,step])Let’s look at the stock market index. The sign of regime switch could arise from (but not limited to):
A ‘switch’ in trends
A ‘switch’ in volatility
(Not shown here) A ‘switch’ in correlation among stocks
And many others!
par(mfrow = c(1,2))
plot(mk)
plot(mkr)library(MSGARCH)
dat <- mkr
spec <- CreateSpec(variance.spec = list(model = c("eGARCH","eGARCH")),
distribution.spec = list(distribution = c("std","std")),
switch.spec = list(do.mix = FALSE))
msg.res <- FitML(spec, dat)
print(msg.res)
msg.sprob <- State(msg.res)
msg.state <- msg.sprob[["Viterbi"]]
plot(msg.state, type = "l")Szakmary, A., Ors, E., Kim, J. K., & Davidson III, W. N. (2003). The predictive power of implied volatility: Evidence from 35 futures markets. Journal of Banking & Finance, 27(11), 2151-2175.↩︎
Bekaert, G., & Hoerova, M. (2014). The VIX, the variance premium and stock market volatility. Journal of econometrics, 183(2), 181-192.↩︎
Giot, P. (2003). The information content of implied volatility in agricultural commodity markets. Journal of Futures Markets: Futures, Options, and Other Derivative Products, 23(5), 441-454.↩︎
Fernández, C., & Steel, M. F. (1998). On Bayesian modeling of fat tails and skewness. Journal of the american statistical association, 93(441), 359-371. https://www.jstor.org/stable/2669632↩︎
See Nelson, D. B. (1991). Conditional heteroskedasticity in asset returns: A new approach. Econometrica: Journal of the econometric society, 347-370.↩︎
Glosten, L. R., R. Jagannathan, and D. E. Runkle, 1993. On The Relation between The Expected Value and The Volatility of Nominal Excess Return on stocks. Journal of Finance 48: 1779-1801. https://www.jstor.org/stable/2329067↩︎
Zakoian, J. M., 1994. Threshold Heteroscedastic Models. Journal of Economic Dynamics and Control 18: 931-955. https://doi.org/10.1016/0165-1889(94)90039-6↩︎
Ding, Z., Granger C.W., Engle R.F., 1991. A Long Memory Property of Stock Market Returns and a New Model. Journal of Empirical Finance 1 (1993) 83-106. https://doi.org/10.1016/S0304-4076(97)00072-9↩︎
Haas, M., Mittnik, S., & Paolella, M. S. (2004). A new approach to Markov-switching GARCH models. Journal of financial Econometrics, 2(4), 493-530.↩︎
https://www.weforum.org/agenda/2023/03/charted-the-link-between-unemployment-and-recessions/↩︎