Stationary Time Series

The following examples cover the main stationary time series processes presented in the lecture. Model generation, identification, estimation, forecasts and residual analysis examples are provided via R.

Generating a Simple Stationary Process

We can generate the stationary process:

  • \(Y_t = \epsilon_t\), where \(\epsilon_t \sim \mathcal{N}(0,1)\)

with the following code:

set.seed(1)
par(mfrow = c(1, 1))
x1 <- rnorm(200, mean = 0, sd = 1)
plot.ts(x1, main = "WN, mean = 0")
abline(0, 0, col = "red")

By drawing the red horizontal line we can visually inspect whether the time series fluctuates around the mean.

Generating Non-Stationary Processes

We will generate the three examples of non-stationary time series:

  • \(Y_t = t + \epsilon_t\), where \(\epsilon_t \sim \mathcal{N}(0,1)\);
  • \(Y_t = \epsilon_t \cdot t\), where \(\epsilon_t \sim \mathcal{N}(0,1)\);
  • \(Y_t = \sum_{j = 1}^t Z_j\), where each independent variable \(Z_j\) is either \(1\) or \(-1\), with a 50% probability for either value.
n = 150
set.seed(123)
#Set custom plot layout
m <- rbind(c(1, 1), c(2, 3))
layout(m)
#layout.show(3) #See how would three figures look in this layout
ns1 = rnorm(n) + 1:n
plot.ts(ns1, main = "non stationary in mean")
ns2 = rnorm(n)*(1:n)
plot.ts(ns2, main = "non stationary in variance")
ns3 = cumsum(sample(c(-1,1), n, replace = TRUE))
plot.ts(ns3, main = "no clear tendency")

The three main cases of a non-stationary process: a non-constant mean, a changing variance. The last example is also non-stationary though it is not immediately apparent from a small sample size.

We can check how this kind of process would look like with different sample sizes:

set.seed(123)
ns31 = cumsum(sample(c(-1,1), 1000, replace = TRUE))
par(mfrow = c(2, 2))
plot.ts(ns31[1:100])
plot.ts(ns31[1:300])
plot.ts(ns31[1:600])
plot.ts(ns31[1:1000])

We can see that, as we increase the sample size, the variance of the data increases.

Because we know that for this type of process \(Var(Y_t) = t\), we can check this using a Monte Carlo simulation:

set.seed(123)
Y <- matrix(NA, nrow = 1000, ncol = 11)
colnames(Y) <- paste0("X", 1:ncol(Y))
for(j in 1:nrow(Y)){
  Y[j, ] <- cumsum(sample(c(-1,1), ncol(Y), replace = TRUE))
}
#Our variance:
apply(Y, 2, var)
##        X1        X2        X3        X4        X5        X6        X7 
##  1.000324  1.968368  3.058735  3.939684  4.809894  5.835836  7.098102 
##        X8        X9       X10       X11 
##  8.151816  8.989389  9.808384 10.675792

We can see that as \(t\) increases and we look at \(X_1, X_2, ..., X_{11}\), the variance also increases, which is in line with the theoretical calculations.

Here are some of the sample series plotted with different colors:

plot.ts(t(Y[c(1, 3, 5, 8),]), plot.type = "single", col = 1:4)

And here is how all of the series look like on the same plot:

plot.ts(t(Y), plot.type = "single") # , col = 1:ncol(Y)

We see that we get a plot that resembles a board. If we increase our simulation number from 1000 to more (5000, 10000, etc.) then our plot would have some additional directories, which were not present in the simulation.

ACF and the PACF Plots

Plot the sample ACF and PACF of the non-stationary time series:

par(mfrow = c(1, 3))
acf(ns1, main = "non stationary in mean")
acf(ns2, main = "non stationary in variance")
acf(ns3, main = "no clear tendency")

pacf(ns1, main = "non stationary in mean")
pacf(ns2, main = "non stationary in variance")
pacf(ns3, main = "no clear tendency")

We can also draw the confidence intervals by hand:

par(mfrow = c(1, 2))
n = length(x1)
acf(x1, main = "WN, mean = 0")
abline(h = qnorm(c(0.025, 0.975))/sqrt(n), col = "red")
pacf(x1, main = "WN, mean = 0")
abline(h = qnorm(c(0.025, 0.975))/sqrt(n), col = "red")

Testing if a Time Series is White Noise

As an example, read the quarterly Canadian employment index data. We can quickly view the series plot, its ACF and PACF by using the tsdisplay function from the forecast package:

suppressPackageStartupMessages({require("forecast")})
txt1 <- "http://uosis.mif.vu.lt/~rlapinskas/(data%20R&GRETL/caemp.txt"
caemp <- read.csv(url(txt1), header = TRUE, as.is = TRUE)
caemp <- ts(caemp, start = c(1960, 1), freq = 4)
tsdisplay(caemp, lwd = 2)

We can test the WN hypothesis by running the Ljung-Box test statistic. Our null hypothesis:

\[H_0: \rho(1) = 0, ..., \rho(k) = 0\] Let us test it for ACF up to lag 5, one by one:

for(k in 1:5){
  print(Box.test(caemp, lag = k, type = "Ljung-Box"))
}
## 
##  Box-Ljung test
## 
## data:  caemp
## X-squared = 127.73, df = 1, p-value < 2.2e-16
## 
## 
##  Box-Ljung test
## 
## data:  caemp
## X-squared = 240.45, df = 2, p-value < 2.2e-16
## 
## 
##  Box-Ljung test
## 
## data:  caemp
## X-squared = 337.32, df = 3, p-value < 2.2e-16
## 
## 
##  Box-Ljung test
## 
## data:  caemp
## X-squared = 418.26, df = 4, p-value < 2.2e-16
## 
## 
##  Box-Ljung test
## 
## data:  caemp
## X-squared = 484.16, df = 5, p-value < 2.2e-16

In all cases, we reject the null hypothesis that all of the autocorrelation functions up to k are zero.

Note that, if we reject the null hypothesis \(H_0: \rho(1) = 0, ..., \rho(k) = 0\) - this means that at least one ACF is non-zero - this does not necessarily mean that all are non-zero.

An example of performing the Ljung-Box test on \(\epsilon_t \sim \mathcal{N}(0,1)\):

set.seed(123)
Box.test(rnorm(100, mean = 0, sd = 1), lag = 1, type = "Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  rnorm(100, mean = 0, sd = 1)
## X-squared = 0.067503, df = 1, p-value = 0.795

p-value = 0.795 > 0.05 - we do not reject the null hypothesis \(H_0: \rho(1) = 0\).

Moving-Average MA(q) Models

We will generate three MA models with \(\epsilon_t \sim \mathcal{N}(0,1)\):

  • \(Y_t = 3 + \epsilon_t + 0.5 \epsilon_{t-1}\);
  • \(Y_t = \epsilon_t + 1.3 \epsilon_{t-1}\);
  • \(Y_t = \epsilon_t - 1.3 \epsilon_{t-1} + 0.4 \epsilon_{t-2}\);

We can generate the processes in three ways: manually or using arima.sim from the forecast package. Generating manually involves writing the code equivalent of the previously provided mathematical formulas:

set.seed(100)
n = 5000
epsilon <- rnorm(n, mean = 0, sd = 1)
MA.m1 <- 3 + epsilon + 0.5 * c(0, epsilon[-n])
MA.m2 <- epsilon + 1.3 * c(0, epsilon[-n])
MA.m3 <- epsilon - 1.3 * c(0, epsilon[-n]) + 0.4 * c(0, 0, epsilon[-((n-1):n)])

The mean of each sample is:

mean(MA.m1)
## [1] 3.002095
mean(MA.m2)
## [1] 0.003159988
mean(MA.m3)
## [1] 0.0001168011

which is what we expect (Note: you can check this by calculating \(\mathbb{E}Y_t\) for each case).

The variance of the second process:

var(MA.m2)
## [1] 2.685624

is very close to the theoretical value \(Var(Y_t) = (1 + \theta^2)\sigma^2 = 1 + 1.3^2 = 2.69\).

The ACF and PACF or the three different models:

par(mfrow = c(2,3))
acf(MA.m1, main = expression(paste(Y[t], " = 3 + ", epsilon[t], " + 0.5", epsilon[t-1])))
acf(MA.m2, main = expression(paste(Y[t], " = ", epsilon[t], " + 1.3", epsilon[t-1])))
acf(MA.m3, main = expression(paste(Y[t], " = ",
                                   epsilon[t], " - 1.3", epsilon[t-1], " + 0.4", epsilon[t-2])))

pacf(MA.m1, main = expression(paste(Y[t], " = 3 + ", epsilon[t], " + 0.5", epsilon[t-1])))
pacf(MA.m2, main = expression(paste(Y[t], " = ", epsilon[t], " + 1.3", epsilon[t-1])))
pacf(MA.m3, main = expression(paste(Y[t], " = ", 
                                    epsilon[t], " - 1.3", epsilon[t-1], " + 0.4", epsilon[t-2])))

For MA(1) the ACF has a sharp cutoff after lag = 1, for MA(2) - the sharp cutoff is after lag = 2. The PACF of the MA(q) process is slowly decaying to zero. Note that the second equation does not meet the requirement that \(|\theta| < 1\). Therefore, it is not invertible, i.e. it cannot be expressed as an infinite-order AR process.

Autoregressive AR(p) Models

We will generate three AR models with \(\epsilon_t \sim \mathcal{N}(0,1)\):

  • \(Y_t = 3 + 0.5 Y_{t-1} + \epsilon_t\);
  • \(Y_t = 1.2 Y_{t-1} + \epsilon_t\);
  • \(Y_t = 1.2 Y_{t-1} - 0.5 Y_{t-2} + \epsilon_t\);

We will also give an example of model generation using arima.sim function:

set.seed(100)
n = 5000
epsilon <- rnorm(n, mean = 0, sd = 1)
AR.m1 <- AR.m2 <- NULL
AR.m1[1] <- 3 + epsilon[1]
AR.m2[1] <- epsilon[1]
for(j in 2:n){
  AR.m1[j] <- 3 + 0.5 * AR.m1[j-1] + epsilon[j]
  AR.m2[j] <- 1.2 * AR.m2[j-1] + epsilon[j]
}
#We can simulate an AR, MA or ARIMA model using arima.sim:
AR.m3 <- arima.sim(list(order = c(2,0,0), ar = c(1.2, -0.5)), n = n)

Note: The ARMA model is checked for stationarity when using arima.sim !.

To make it easier to analyse the second model, we will reduce the sample of the generated data:

AR.m2 <- AR.m2[1:50]
plot.ts(AR.m2)

Because the coefficient \(\phi = 1.2 \geq 1\) the time series grows exponentially.

The mean of each sample is:

mean(AR.m1)
## [1] 6.001387
mean(AR.m2)
## [1] 106.433
mean(AR.m3)
## [1] 0.02753754

Remember that for the generalized AR(1) case (with constant component) the mean is \(m = \alpha/(1-\phi)\). For the first model, we have that \(\alpha = 3\) and \(\phi = 0.5\) \(\Rightarrow\) \(m = 6\) (which is similar to what the sample mean of the first model is).

Note that the second model violates the requirement that \(|\phi| < 1\). This kind of AR process is not stationary.

Trying to generate this kind of model using arima.sim would result in an error:

arima.sim(list(order = c(1,0,0), ar = c(1.2)), n = n)
## Error in arima.sim(list(order = c(1, 0, 0), ar = c(1.2)), n = n): 'ar' part of model is not stationary

This is on of the differences when generating data manually and using arima.sim - sometimes we want to generate non-stationary data and see how our R script would handle this situation.

The ACF and PACF:

par(mfrow = c(2,3))
acf(AR.m1, main = expression(paste(Y[t], " = 3 + 0.5", Y[t-1], " + ", epsilon[t])))
acf(AR.m2, main = expression(paste(Y[t], " = 1.2", Y[t-1], " + ", epsilon[t])))
acf(AR.m3, main = expression(paste(Y[t], " = 1.2", Y[t-1], " - 0.5", Y[t-2], " + ", epsilon[t])))

pacf(AR.m1, main = expression(paste(Y[t], " = 3 + 0.5", Y[t-1], " + ", epsilon[t])))
pacf(AR.m2, main = expression(paste(Y[t], " = 1.2", Y[t-1], " + ", epsilon[t])))
pacf(AR.m3, main = expression(paste(Y[t], " = 1.2", Y[t-1], " - 0.5", Y[t-2], " + ", epsilon[t])))

For the AR(p) models, the ACF decays gradually to zero. For the AR(1), the PACF has a sharp cutoff after lag = 1. For the AR(2), the PACF has a sharp cutoff after lag = 2. Note that even though the second model is not stationary, the PACF is cut-off after the first lag and the ACF is slowly decaying (because the coefficient \(\phi > 1\) the decay is slower compared to the first model with \(\phi = 0.5\))

Autoregressive Moving-Average ARMA(p,q) Models

We will generate an ARMA model:

  • \(Y_t = 0.4 Y_{t-1} + \epsilon_t + 0.6 \epsilon_{t-1}\)
set.seed(100)
n = 5000
ARMA.m1 <- arima.sim(list(order = c(1,0,1), ar = 0.4, ma = 0.6), n = n)

The ACF and PACF:

acf(ARMA.m1, main = expression(paste(Y[t], " = 0.4", 
                                     Y[t-1], " + ", epsilon[t], " + 0.6", epsilon[t-1])))

pacf(ARMA.m1, main = expression(paste(Y[t], " = 0.4", 
                                      Y[t-1], " + ", epsilon[t], " + 0.6", epsilon[t-1])))

Both ACF and PACF decline for the ARMA process.

Model Estimation with Arima

Let’s say we do not know the data generating process that generated out previously mentioned AR, MA and ARMA models. If we examined the ACF and PACF graphs and we have an idea of what our underlying model might be, we can use the Arima function from the forecast package (NOT arima from stats package) to estimate our model.

MA.m1.est <- Arima(MA.m1, order = c(0, 0, 1), include.mean = TRUE)
MA.m2.est <- Arima(MA.m2, order = c(0, 0, 1), include.mean = FALSE)
MA.m3.est <- Arima(MA.m3, order = c(0, 0, 2), include.mean = FALSE)
AR.m1.est <- Arima(AR.m1, order = c(1, 0, 0), include.mean = TRUE)
AR.m2.est <- Arima(AR.m2, order = c(1, 0, 0), include.mean = FALSE, method = "ML")
AR.m3.est <- Arima(AR.m3, order = c(2, 0, 0), include.mean = FALSE)
ARMA.m1.est<-Arima(ARMA.m1, order = c(1, 0, 1), include.mean = FALSE)

The summary statistics of our models:

summary(MA.m1.est)
## Series: MA.m1 
## ARIMA(0,0,1) with non-zero mean 
## 
## Coefficients:
##          ma1    mean
##       0.5083  3.0021
## s.e.  0.0122  0.0213
## 
## sigma^2 estimated as 0.9999:  log likelihood=-7093.67
## AIC=14193.33   AICc=14193.34   BIC=14212.88
## 
## Training set error measures:
##                        ME      RMSE       MAE       MPE     MAPE      MASE
## Training set 4.945214e-05 0.9997648 0.8003394 -31.77517 66.31426 0.8106853
##                      ACF1
## Training set -0.009559869

To only leave the relevant info, we only take the coefficients and their standard errors:

print(data.frame(coef = coef(MA.m2.est), "s.e" = sqrt(diag(vcov(MA.m2.est)))))
##          coef         s.e
## ma1 0.7712797 0.008735223
print(data.frame(coef = coef(MA.m3.est), `s.e` = sqrt(diag(vcov(MA.m3.est)))))
##           coef        s.e
## ma1 -1.3100158 0.01311805
## ma2  0.4017394 0.01315534

We can see that our estimated coefficients are close to the actual values, except for the second process, which is not invertible.

For the AR models:

print(data.frame(coef = coef(AR.m1.est), "s.e" = sqrt(diag(vcov(AR.m1.est)))))
##                coef        s.e
## ar1       0.4854857 0.01237243
## intercept 6.0010377 0.02749717
print(data.frame(coef = coef(AR.m2.est), "s.e" = sqrt(diag(vcov(AR.m2.est)))))
##          coef         s.e
## ar1 0.9979733 0.002871269
print(data.frame(coef = coef(AR.m3.est), "s.e" = sqrt(diag(vcov(AR.m3.est)))))
##           coef        s.e
## ar1  1.1678721 0.01241928
## ar2 -0.4781098 0.01242592

Note that the second process is not stationary. As such the estimated value of \(\phi\) is close to 1. We also use method = "ML" when estimating the non-stationary series to get a stationary model. If we do not use this option, then we get an error that our time series is not stationary:

Arima(AR.m2, order = c(1, 0, 0), include.mean = FALSE)
## Error in stats::arima(x = x, order = order, seasonal = seasonal, include.mean = include.mean, : non-stationary AR part from CSS

When estimating AR, MA or ARMA processes we assume that they are stationary (if they are not, then we need to create such models or transform the data such that it would be stationary).

Note that the intercept coefficient of the first \(AR(1)\) model is close to the process mean m = 6.

summary(ARMA.m1.est)
## Series: ARMA.m1 
## ARIMA(1,0,1) with zero mean 
## 
## Coefficients:
##          ar1     ma1
##       0.3735  0.6199
## s.e.  0.0160  0.0131
## 
## sigma^2 estimated as 0.9996:  log likelihood=-7093.15
## AIC=14192.3   AICc=14192.31   BIC=14211.85
## 
## Training set error measures:
##                       ME      RMSE       MAE      MPE     MAPE      MASE
## Training set 0.001018066 0.9995867 0.8001791 -79.8131 424.9294 0.8285089
##                     ACF1
## Training set 0.005094493

The ARIMA model coefficients are also close to the actual values used when generating the model.

Model Estimation with auto.arima

If we are not sure what order of AR, MA or ARMA model would best fit our data, we can use the auto.arima function.

MA.m1.auto <- auto.arima(MA.m1)
MA.m2.auto <- auto.arima(MA.m2)
MA.m3.auto <- auto.arima(MA.m3)
AR.m1.auto <- auto.arima(AR.m1)
AR.m2.auto <- auto.arima(AR.m2)
AR.m3.auto <- auto.arima(AR.m3)
ARMA.m1.auto<-auto.arima(ARMA.m1)

We can see, that for simple AR or MA processes auto.arima works fine, however, in some cases the best model differs from the actual model used when generating the data:

summary(MA.m1.auto)
## Series: MA.m1 
## ARIMA(0,0,1) with non-zero mean 
## 
## Coefficients:
##          ma1    mean
##       0.5083  3.0021
## s.e.  0.0122  0.0213
## 
## sigma^2 estimated as 0.9999:  log likelihood=-7093.67
## AIC=14193.33   AICc=14193.34   BIC=14212.88
## 
## Training set error measures:
##                        ME      RMSE       MAE       MPE     MAPE      MASE
## Training set 4.945214e-05 0.9997648 0.8003394 -31.77517 66.31426 0.8106853
##                      ACF1
## Training set -0.009559869
print(data.frame(coef = coef(MA.m2.auto), "s.e" = sqrt(diag(vcov(MA.m2.auto)))))
##          coef         s.e
## ma1 0.7712797 0.008735223
print(data.frame(coef = coef(MA.m3.auto), "s.e" = sqrt(diag(vcov(MA.m3.auto)))))
##             coef        s.e
## ar1  0.002381094 0.06571518
## ar2 -0.024920403 0.03035926
## ma1 -1.304041029 0.06444188
## ma2  0.398838343 0.05764883

For the MA process, the first two models are identical to the Arima specified models. However, the last model is different.

Similar results can be seen for the AR(2) models:

print(data.frame(coef = coef(AR.m1.auto), "s.e" = sqrt(diag(vcov(AR.m1.auto)))))
##                  coef        s.e
## ar1        0.49872222 0.01414085
## ar2       -0.02732634 0.01415136
## intercept  6.00100454 0.02675482
print(data.frame(coef = coef(AR.m2.auto), "s.e" = sqrt(diag(vcov(AR.m2.auto)))))
##          coef       s.e
## ma1 0.8915756 0.1164140
## ma2 1.0606306 0.1011557
## ma3 0.4085473 0.1054677
print(data.frame(coef = coef(AR.m3.auto), "s.e" = sqrt(diag(vcov(AR.m3.auto)))))
##           coef        s.e
## ar1  1.1678721 0.01241928
## ar2 -0.4781098 0.01242592

Because the second model is not stationary, auto.arima tries to fit the best stationary model.

summary(ARMA.m1.auto)
## Series: ARMA.m1 
## ARIMA(2,0,3) with non-zero mean 
## 
## Coefficients:
##           ar1     ar2     ma1     ma2     ma3    mean
##       -0.6637  0.3137  1.6620  0.7175  0.0481  0.0058
## s.e.   0.0469  0.0450  0.0484  0.0836  0.0378  0.0359
## 
## sigma^2 estimated as 0.9997:  log likelihood=-7091.58
## AIC=14197.16   AICc=14197.19   BIC=14242.78
## 
## Training set error measures:
##                        ME      RMSE       MAE       MPE     MAPE     MASE
## Training set -0.001232387 0.9992723 0.7999059 -92.63542 440.1049 0.828226
##                      ACF1
## Training set 0.0004421623

The ARMA model differs the most from the actual model even though we used arima.sim to generate the model.

Model Forecasts

Having estimated the models, we would like to create our time series forecasts:

Lets say we want to forecast 10 periods into the future. Since we have 5000 historic data points we would also want to reduce the historic data in the plot to, say 100.

Our code would look like:

#?plot.forecast
plot(forecast(MA.m1.auto, h = 10), include = 100)

Some other forecast plots:

par(mfrow = c(2, 2))
plot(forecast(MA.m3.auto, h = 10), include = 100)
plot(forecast(AR.m2.auto, h = 10), include = 100)
plot(forecast(AR.m3.auto, h = 10), include = 100)
plot(forecast(ARMA.m1.auto, h = 10), include = 100)

AR and MA model forecasts perform as we would expect: MA forecasts quickly tend to the mean, whereas AR forecasts tend to the average but never reach it.

Model Residuals

After estimating our models, it is also important to check if the residuals are WN. To get the model residuals use the residuals function:

res1.MA <- residuals(MA.m1.auto)
res2.MA <- residuals(MA.m2.auto)
res3.MA <- residuals(MA.m3.auto)
res1.AR <- residuals(AR.m1.auto)
res2.AR <- residuals(AR.m2.auto)
res3.AR <- residuals(AR.m3.auto)
res1.ARMA <- residuals(ARMA.m1.auto)

We can plot the residuals (we will plot the last 100):

par(mfrow = c(1,2))
plot.ts(tail(res2.AR, 100))
plot.ts(tail(res1.ARMA, 100))

We can check their ACF and PACF:

par(mfrow = c(2,2))
acf(res2.AR, 100)
acf(res1.ARMA, 100)
pacf(res2.AR, 100)
pacf(res1.ARMA, 100)

And perform the Ljung-Box test:

Box.test(res2.AR, lag = 20, type = "Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  res2.AR
## X-squared = 50.262, df = 20, p-value = 0.0002031

Because p-value < 0.05, we reject the null hypothesis of no serial correlation, i.e. the residuals are not white noise. This is to be expected, because the second model is estimated on non-stationary data.

Box.test(res1.ARMA, lag = 20, type = "Ljung-Box")$p.value
## [1] 0.7141348

We do not reject the null hypothesis that the first 20 lag autocorrelations among the residuals are zero (p-value > 0.05) for the ARMA model.

We can also perform the Box-Pierce test:

Box.test(res2.AR, lag = 20, type = "Box-Pierce")
## 
##  Box-Pierce test
## 
## data:  res2.AR
## X-squared = 41.936, df = 20, p-value = 0.002819
Box.test(res1.ARMA, lag = 20, type = "Box-Pierce")
## 
##  Box-Pierce test
## 
## data:  res1.ARMA
## X-squared = 15.993, df = 20, p-value = 0.7171

Which lets us draw the same conclusions as the Ljung-Box test results.

Remember, the null hypothesis is:

\[H_0: \rho(1) = 0, ..., \rho(k) = 0\]

The test statistics:

  • Ljung-Box test statistic: \[ Q = T(T+2)\sum_{\tau = 1}^k \dfrac{\hat{\rho}^2(\tau)}{T - \tau}\]
  • Box-Pierce test statistic: \[ Q = T\sum_{\tau = 1}^k \hat{\rho}^2(\tau)\]

are calculated differently but have the same critical region - if \(Q > \chi_{K}^2\), we reject the null hypothesis of independently distributed residuals.

Model Misspecification

Another example, this time by fitting an incorrect model.

Remember our MA(1) model which is not invertible \(Y_t = \epsilon_t + 1.3 \epsilon_{t-1}\)?

Let us try to fit an AR(1) model and check its residuals:

par(mfrow = c(1,2))
bad.model <- Arima(MA.m2, order = c(1,0,0))
acf(residuals(bad.model))
pacf(residuals(bad.model))

Box.test(residuals(bad.model), lag = 10, type = "Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  residuals(bad.model)
## X-squared = 633.87, df = 10, p-value < 2.2e-16

The ACF and PACF graphs show that the residuals exhibit autocorrelation. The Ljung-Box test rejects the null hypothesis (because p-value < 0.05) that the first 10 lag autocorrelations among residuals are zero. We can conclude that the residuals are not WN, which indicates that the specified model is not correct.

To correct this, with our current knowledge of time series, we would need to specify a different model or try to leverage the results of auto.arima.

An Example of how to Calculate the Mean of an ARMA Process

Lets say our \(ARMA(1,1)\) model is of the following general form: \[Y_t = \alpha + \phi_1 Y_{t-1} + \epsilon_t + \theta_1 \epsilon_{t-1}\]

Let \(\mathbb{E}Y_t = \mu\), then we have:

\[\begin{align} \mathbb{E}(Y_t) &= \alpha + \phi_1 \mathbb{E} Y_{t-1} + \mathbb{E}\epsilon_t + \theta_1 \mathbb{E}\epsilon_{t-1}\\ \mu &= \alpha + \phi_1 \mu \\ \mu &= \dfrac{\alpha}{1 - \phi_1} \end{align}\]

Which is the same as the generalized \(AR(1)\) case.

Let’s take an example with \(\alpha = 2\), \(\phi_1 = 0.2\) and \(\theta_1 = 0.3\). then our \(\mu\) is:

2/(1 - 0.2)
## [1] 2.5

Why would we need tho separately calculate the mean? Because arima.sim generates a model which has a zero mean. To get around this, we either generate the model manually:

set.seed(1)
general.arima.v1 <- NULL
epsilon <- rnorm(mean = 0, sd = 1, n = 2500)

general.arima.v1[1] <- 2 + epsilon[1]
for(j in 2:length(epsilon)){
  general.arima.v1[j] <- 2 + 0.2 * general.arima.v1[j-1] + epsilon[j] + 0.3 * epsilon[j-1]
}
mean(general.arima.v1)
## [1] 2.484964

or we can generate the data using arima.sim:

set.seed(1)
general.arima.v2 <- arima.sim(list(order = c(1,0,1), ar = 0.2, ma = 0.3), n = 2500)
mean(general.arima.v2)
## [1] -0.01640406

and add the mean to the series:

general.arima.v2 <- 2.5 + general.arima.v2
mean(general.arima.v2)
## [1] 2.483596

We can plot the first 100 observation the two series side by side:

plot.ts(cbind(general.arima.v1, general.arima.v2)[1:100, ])

Next, we can estimate the two different datasets:

mdl.arma.v1 <- Arima(general.arima.v1, order = c(1, 0, 1), include.mean = TRUE)
mdl.arma.v2 <- Arima(general.arima.v2, order = c(1, 0, 1), include.mean = TRUE)

coef(mdl.arma.v1)
##       ar1       ma1 intercept 
## 0.1860932 0.3085092 2.4844429
coef(mdl.arma.v2)
##       ar1       ma1 intercept 
## 0.1911828 0.3046751 2.4837918

We note that the intercept in both cases is close to the process mean. The difference is in the naming - intercept is not the \(\alpha\) coefficient, it is the \(\mu\)!

Interestingly, if we print the summary of the model:

summary(mdl.arma.v2)
## Series: general.arima.v2 
## ARIMA(1,0,1) with non-zero mean 
## 
## Coefficients:
##          ar1     ma1    mean
##       0.1912  0.3047  2.4838
## s.e.  0.0419  0.0407  0.0335
## 
## sigma^2 estimated as 1.077:  log likelihood=-3638.9
## AIC=7285.8   AICc=7285.81   BIC=7309.09
## 
## Training set error measures:
##                        ME     RMSE       MAE       MPE    MAPE      MASE
## Training set 3.903516e-05 1.037248 0.8270324 -90.77095 179.704 0.8361583
##                      ACF1
## Training set 0.0007538847

The intercept is now labeled as the mean! This one of the problems when using very generalized functions for model estimation - since they are written in such a way that they cover a variety of models and cases, their output is not always as intuitive as it may seem. This is further discussed here and here

Process Invertibility

If an \(MA(q)\) process is invertible, then:

\[Y_t = \epsilon_t + \theta \epsilon_{t-1}\]

Can be written as an \(AR(\infty )\) process:

\[Y_t = \sum_{j = 1}^\infty \tilde{\phi}_j Y_{t-j} + \epsilon_t\]

Since the \(MA\) process is a finite approximation of the Wold representation, we could also approximate its infinite \(AR\) representation reasonably well. This could also be the case for the infinite \(MA\) process approximation.

This means that, in practice, we may come across cases where more than one model could be used, e.g. a lower order \(ARMA\) process, which is stationary and invertible, could be approximated by a higher order \(MA\) or \(AR\) process, especially if the sample size is small.

We will provide an example using data generated from an \(ARMA(2, 3)\) process:

set.seed(100)
# We will use a smaller data sample:
n = 100
# We will generate an additional 5 periods ahead, which we will use to compare the predictions to:
ARMA_data <- arima.sim(list(order = c(2,0,3), ar = c(0.3, -0.2), ma = c(0.4, 0.2, -0.1)), n = n + 5)
# Save the last five periods:
actual_data = ts(ARMA_data[(n+1):(n+5)], start = n + 1)
# Exclude the last 5 values from the data:
ARMA_data <- ARMA_data[1:n]

We will estimate 4 models:

  • The true model using Arma;
  • The best model using auto.arima;
  • The best \(AR\) model using auto.arima;
  • The best \(MA\) model using auto.arima;
ARMA_model = Arima(ARMA_data, order = c(2,0,3), include.mean = FALSE)
ARMA_auto_arma = auto.arima(ARMA_data)
ARMA_auto_ar   = auto.arima(ARMA_data, max.p = 10, max.q = 0)
ARMA_auto_ma   = auto.arima(ARMA_data, max.p = 0, max.q = 10)

The resulting models:

print(data.frame(ARMA = ARMA_model$coef))
##            ARMA
## ar1  0.53388006
## ar2 -0.05686963
## ma1 -0.11185555
## ma2 -0.09679085
## ma3 -0.34581676
print(data.frame(ARMA = ARMA_auto_arma$coef))
##           ARMA
## ar1  0.8294162
## ar2 -0.6277636
## ma1 -0.4047657
## ma2  0.3237735
print(data.frame(AR = ARMA_auto_ar$coef))
##             AR
## ar1  0.4009224
## ar2 -0.1016042
## ar3 -0.2739891
print(data.frame(MA = ARMA_auto_ma$coef))
##              MA
## ma1  0.42249089
## ma2  0.05592508
## ma3 -0.33983869
## ma4 -0.14555316

Gives us different models:

  • The best model using auto.arima is an \(ARMA(2, 2)\);
  • The best \(AR\) model using auto.arima is an \(AR(3)\);
  • The best \(MA\) model using auto.arima is an \(MA(4)\);

The residuals ACF and PACF plots:

par(mfrow = c(2, 2))
Acf(ARMA_model$residuals);Pacf(ARMA_model$residuals)
Acf(ARMA_auto_arma$residuals);Pacf(ARMA_auto_arma$residuals)

par(mfrow = c(2, 2))
Acf(ARMA_auto_ar$residuals);Pacf(ARMA_auto_ar$residuals)
Acf(ARMA_auto_ma$residuals);Pacf(ARMA_auto_ma$residuals)

We can see that the residuals are not autocorrelated in either of these models, so each of them is a valid specification for our data sample.

And from the model AIC values:

print(paste0("True ARMA AIC = ", ARMA_model$aic))
## [1] "True ARMA AIC = 294.427693105794"
print(paste0("AUTO ARMA AIC = ", ARMA_auto_arma$aic))
## [1] "AUTO ARMA AIC = 293.348716616519"
print(paste0("AUTO AR   AIC = ", ARMA_auto_ar$aic))
## [1] "AUTO AR   AIC = 290.187669396479"
print(paste0("AUTO MA   AIC = ", ARMA_auto_ma$aic))
## [1] "AUTO MA   AIC = 292.598059624327"

we can see that:

  • Our true \(ARMA(2,3)\) model has the largest AIC, therefore it is the worst one of the 4 possible models.

  • The \(ARMA(2, 2)\) model estimated via auto.arima is the second-worst as it has the second largest AIC from the estimated models.

  • The \(AR(3)\) model has the lowest \(AIC\) value from the estimated models.

Note: we can also examine their BIC values with $bic but we will get the same results.

So, the best model for our data sample is the \(AR(3)\) model. Let’s compare the forecasts for 5 periods ahead:

forc_1 = forecast(ARMA_model, h = 5)$mean
forc_2 = forecast(ARMA_auto_arma, h = 5)$mean
forc_3 = forecast(ARMA_auto_ar, h = 5)$mean
forc_4 = forecast(ARMA_auto_ma, h = 5)$mean

We can check how the fitted data looks like for each model:

plot.ts(ARMA_data)
lines(fitted(ARMA_model), col = "red", lty = 2)
lines(fitted(ARMA_auto_arma), col = "blue", lty = 2)
lines(fitted(ARMA_auto_ar), col = "darkgreen", lty = 2)
lines(fitted(ARMA_auto_ma), col = "darkorange", lty = 2)

legend("bottomleft", #x = 0, y = -1, 
       legend = c("Actual", "ARMA_true", "ARMA_auto", "AR_auto", "MA_auto"),
       col = c("black", "red", "blue", "darkgreen", "darkorange"), 
       lty = c(1, 2, 2, 2, 2), cex = 0.8)

Graphically, there is little difference between the models. This can also be said for their forecasts:

plot.ts(ARMA_data, xlim = c(98, 105))
lines(actual_data, col = "black", lty = 2)
lines(forc_1, col = "red", lty = 2)
lines(forc_2, col = "blue", lty = 2)
lines(forc_3, col = "darkgreen", lty = 2)
lines(forc_4, col = "darkorange", lty = 2)

legend("bottomleft", #x = 0, y = -1, 
       legend = c("Historical", "Forecast_actual", "ARMA_true", "ARMA_auto", "AR_auto", "MA_auto"),
       col = c("black", "black", "red", "blue", "darkgreen", "darkorange"), 
       lty = c(1, 2, 2, 2, 2, 2), cex = 0.8)

Financial Volatility

Example Data Introduction

The following data contains:

  • S&P Composite index returns Ret over the period 1954:1-2001:9;
  • the consumer price index (CPI) inflation rate Inf.
  • the change in the three-month Treasury-Bill (T-bill) rate DT.bill.
txt1 <- "http://uosis.mif.vu.lt/~rlapinskas/(data%20R&GRETL/"
txt2 <- "S&P.txt"
#Read it as an excel file
SP <- read.table(paste0(txt1, txt2), header = TRUE)
SP <- read.table(paste0(txt1, txt2), header = TRUE)
SP <- ts(SP, start = 1954, freq = 12)

We begin by modelling the return series Ret as a function of a constant, one lag of returns Ret, one lag of the inflation rate Inf and one lag of the first-difference of the three-month T-bill rate DT.bill. We expect our model to be of the form:

\[ \begin{cases} Y_t &= \gamma + \phi_1 Y_{t-1} + ... \beta_1 X_{1,t} + ... + \epsilon_t \\ \epsilon_t &= \sigma_t w_t \\ \sigma_t^2 &= \omega + \sum_{j = 1}^q \alpha_j \epsilon_{t-j}^2 + \sum_{i = 1}^p \beta_i \sigma_{t-i}^2 \end{cases} \]

where:

  • \(w_t\) are i.i.d. with conditional mean 0 and conditional variance 1.
  • \(X_t\) are exogenous variables;

We will begin by adding the first lags of the aforementioned variables and create the usual OLS model. Then we will test the null hypothesis:

\[H_0: \text{residuals of the model are normal}\]

In conventional regression models, rejection of the null hypothesis would lead us to change the model specification. However, non-normality is an inherent feature of the errors from the regression models for financial data. As such, robust standard errors are needed.

Mean equation specification

suppressPackageStartupMessages({library(dynlm)})
mean_mdl <- dynlm(Ret ~ L(Ret) + L(Inf.) + L(DT.bill), data = SP)
summary(mean_mdl)
## 
## Time series regression with "ts" data:
## Start = 1954(2), End = 2001(8)
## 
## Call:
## dynlm(formula = Ret ~ L(Ret) + L(Inf.) + L(DT.bill), data = SP)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -12.1199  -1.8907   0.1141   2.0607  10.3608 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.17715    0.21004   5.604 3.26e-08 ***
## L(Ret)       0.20799    0.04102   5.070 5.39e-07 ***
## L(Inf.)     -1.17945    0.44805  -2.632  0.00871 ** 
## L(DT.bill)  -1.25007    0.29135  -4.291 2.10e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.259 on 567 degrees of freedom
## Multiple R-squared:  0.1011, Adjusted R-squared:  0.0963 
## F-statistic: 21.25 on 3 and 567 DF,  p-value: 4.688e-13

As mentioned before, we need to correct the standard errors of our model:

suppressPackageStartupMessages({
  library(sandwich)
  library(lmtest)
})
coeftest(mean_mdl, df = Inf, vcov = vcovHC(mean_mdl, type = "HC"))
## 
## z test of coefficients:
## 
##              Estimate Std. Error z value  Pr(>|z|)    
## (Intercept)  1.177149   0.228283  5.1565 2.516e-07 ***
## L(Ret)       0.207991   0.048629  4.2771 1.893e-05 ***
## L(Inf.)     -1.179451   0.521536 -2.2615 0.0237285 *  
## L(DT.bill)  -1.250071   0.328941 -3.8003 0.0001445 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

where:

  • Variance-covariance matrix with type = "HC" means that the standard errors are corrected to be Heteroskedasticity-Consistent.
  • coeftest - a z-test (using a normal approximation (df = Inf)).

As we can see, all the coefficients are statistically significant because their p-values are less than 0.05.

We can test the residual normality with the Shapiro-Wilk test:

shapiro.test(mean_mdl$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  mean_mdl$residuals
## W = 0.97873, p-value = 2.233e-07
# tseries::jarque.bera.test(mdl.1$residuals)

Since p-value < 0.05 we reject the null hypothesis that our residuals are normal.

Let us look at the residual correlograms:

suppressPackageStartupMessages({library(forecast)})
tsdisplay(mean_mdl$residuals)

Box.test(mean_mdl$residuals, lag = 15)
## 
##  Box-Pierce test
## 
## data:  mean_mdl$residuals
## X-squared = 23.766, df = 15, p-value = 0.06917

We do not reject the null hypothesis that the autocorrelations are not statistically significantly different from zero.

Let us test the same with the squared residuals:

tsdisplay(mean_mdl$residuals^2)

Box.test(mean_mdl$residuals^2, lag = 6)
## 
##  Box-Pierce test
## 
## data:  mean_mdl$residuals^2
## X-squared = 17.567, df = 6, p-value = 0.007412

We reject the null hypothesis, which means that our squared residuals are autocorrelated.

From the PACF of the squared residuals we see that the an \(ARCH(6)\) model might be appropriate. So, we will now test the residuals for \(ARCH(6)\):

suppressPackageStartupMessages({library(fGarch)})
mdl.res <- mean_mdl$residuals
mdl.arch <- garchFit(formula = ~ garch(6, 0), trace = FALSE, data = mdl.res)
summary(mdl.arch)
## 
## Title:
##  GARCH Modelling 
## 
## Call:
##  garchFit(formula = ~garch(6, 0), data = mdl.res, trace = FALSE) 
## 
## Mean and Variance Equation:
##  data ~ garch(6, 0)
## <environment: 0x000000003bbf2f48>
##  [data = mdl.res]
## 
## Conditional Distribution:
##  norm 
## 
## Coefficient(s):
##          mu        omega       alpha1       alpha2       alpha3  
## -1.5103e-16   5.6728e+00   1.2585e-01   1.2245e-02   1.2945e-01  
##      alpha4       alpha5       alpha6  
##  1.1232e-02   9.3974e-02   1.2805e-01  
## 
## Std. Errors:
##  based on Hessian 
## 
## Error Analysis:
##          Estimate  Std. Error  t value Pr(>|t|)    
## mu     -1.510e-16   1.286e-01    0.000   1.0000    
## omega   5.673e+00   1.072e+00    5.294  1.2e-07 ***
## alpha1  1.259e-01   5.986e-02    2.102   0.0355 *  
## alpha2  1.225e-02   3.999e-02    0.306   0.7595    
## alpha3  1.294e-01   6.952e-02    1.862   0.0626 .  
## alpha4  1.123e-02   3.852e-02    0.292   0.7706    
## alpha5  9.397e-02   6.653e-02    1.412   0.1578    
## alpha6  1.280e-01   5.793e-02    2.210   0.0271 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log Likelihood:
##  -1468.767    normalized:  -2.572272 
## 
## Description:
##  Mon Feb 19 23:28:51 2018 by user: Andrius 
## 
## 
## Standardised Residuals Tests:
##                                 Statistic p-Value     
##  Jarque-Bera Test   R    Chi^2  28.08916  7.952716e-07
##  Shapiro-Wilk Test  R    W      0.9883611 0.0001652405
##  Ljung-Box Test     R    Q(10)  10.2739   0.416801    
##  Ljung-Box Test     R    Q(15)  22.56042  0.09392936  
##  Ljung-Box Test     R    Q(20)  27.36542  0.1252926   
##  Ljung-Box Test     R^2  Q(10)  1.268283  0.9994941   
##  Ljung-Box Test     R^2  Q(15)  6.411728  0.9719713   
##  Ljung-Box Test     R^2  Q(20)  9.804126  0.9715862   
##  LM Arch Test       R    TR^2   2.103424  0.9992301   
## 
## Information Criterion Statistics:
##      AIC      BIC      SIC     HQIC 
## 5.172564 5.233473 5.172179 5.196327

We see that alpha6 coefficient is statistically significant and we reject the null hypothesis \(H_0: \alpha_1 =...=\alpha_6=0\) in:

\[ \begin{cases} \epsilon_t = \sigma_t w_t \\ \sigma_t^2 = \omega + \alpha_1 \epsilon_{t-1}^2 + ... + \alpha_6 \epsilon_{t-6}^2 \end{cases} \] The LM Arch Test on the standardized residuals from this ARCH model says that we do not reject the null hypothesis of uncorrelated residuals (i.e. the standardized residuals from the ARCH model do not have any ARCH effects).

Thus, the residuals of our model make \(ARCH(6)\). The \(ARCH(6)\) model has many parameters and one common advice is to replace the model by \(GARCH(1, 1)\):

mdl.garch <- garchFit(formula = ~ garch(1, 1), trace = FALSE, data = mdl.res)
summary(mdl.garch)
## 
## Title:
##  GARCH Modelling 
## 
## Call:
##  garchFit(formula = ~garch(1, 1), data = mdl.res, trace = FALSE) 
## 
## Mean and Variance Equation:
##  data ~ garch(1, 1)
## <environment: 0x000000001db25b60>
##  [data = mdl.res]
## 
## Conditional Distribution:
##  norm 
## 
## Coefficient(s):
##          mu        omega       alpha1        beta1  
## -1.5103e-16   1.1511e+00   1.1392e-01   7.8510e-01  
## 
## Std. Errors:
##  based on Hessian 
## 
## Error Analysis:
##          Estimate  Std. Error  t value Pr(>|t|)    
## mu     -1.510e-16   1.292e-01    0.000  1.00000    
## omega   1.151e+00   5.183e-01    2.221  0.02636 *  
## alpha1  1.139e-01   4.002e-02    2.846  0.00442 ** 
## beta1   7.851e-01   6.584e-02   11.925  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log Likelihood:
##  -1472.623    normalized:  -2.579025 
## 
## Description:
##  Mon Feb 19 23:28:51 2018 by user: Andrius 
## 
## 
## Standardised Residuals Tests:
##                                 Statistic p-Value     
##  Jarque-Bera Test   R    Chi^2  36.679    1.084565e-08
##  Shapiro-Wilk Test  R    W      0.9851709 1.483964e-05
##  Ljung-Box Test     R    Q(10)  10.01125  0.4395069   
##  Ljung-Box Test     R    Q(15)  22.05404  0.1063948   
##  Ljung-Box Test     R    Q(20)  27.39449  0.1245256   
##  Ljung-Box Test     R^2  Q(10)  7.909245  0.6377014   
##  Ljung-Box Test     R^2  Q(15)  13.07198  0.5967378   
##  Ljung-Box Test     R^2  Q(20)  17.00733  0.6524979   
##  LM Arch Test       R    TR^2   8.202932  0.7690776   
## 
## Information Criterion Statistics:
##      AIC      BIC      SIC     HQIC 
## 5.172060 5.202514 5.171962 5.183941

  • From the summary we see that Shapiro-Wilk Test and Jarque-Bera Test reject the null hypothesis that residuals are normal.

  • The Ljung-Box Test does not reject the null hypothesis that residuals form \(WN\). They also do not reject the null hypothesis that the squared residuals R^2 are WN.

Note on output formatting

We should also note that the output from \(GARCH\) models is quite large. If we want to only get the coefficients, the fit statistics, we can do the following:

mdl.garch@fit$matcoef
##             Estimate  Std. Error       t value   Pr(>|t|)
## mu     -1.510320e-16  0.12923004 -1.168707e-15 1.00000000
## omega   1.151120e+00  0.51831280  2.220897e+00 0.02635790
## alpha1  1.139233e-01  0.04002464  2.846328e+00 0.00442266
## beta1   7.851049e-01  0.06583855  1.192470e+01 0.00000000
mdl.garch@fit$ics
##      AIC      BIC      SIC     HQIC 
## 5.172060 5.202514 5.171962 5.183941

The value of the garchFit function is an not a list (S3 object), but an S4 object; its components are accessed not with the $ symbol but with @. Each accessed component can be a list, vector, data.frame, etc. which is why we also used the $ symbol to extract the coefficient matrix.

If need be, we can also use capture.output to save each output from the summary function into a vector of strings. Then we can print any summary lines that we need to reduce the output clutter:

capture.output(summary(mdl.garch))[39:49]
##  [1] "Standardised Residuals Tests:"                         
##  [2] "                                Statistic p-Value     "
##  [3] " Jarque-Bera Test   R    Chi^2  36.679    1.084565e-08"
##  [4] " Shapiro-Wilk Test  R    W      0.9851709 1.483964e-05"
##  [5] " Ljung-Box Test     R    Q(10)  10.01125  0.4395069   "
##  [6] " Ljung-Box Test     R    Q(15)  22.05404  0.1063948   "
##  [7] " Ljung-Box Test     R    Q(20)  27.39449  0.1245256   "
##  [8] " Ljung-Box Test     R^2  Q(10)  7.909245  0.6377014   "
##  [9] " Ljung-Box Test     R^2  Q(15)  13.07198  0.5967378   "
## [10] " Ljung-Box Test     R^2  Q(20)  17.00733  0.6524979   "
## [11] " LM Arch Test       R    TR^2   8.202932  0.7690776   "

Combining Regression and GARCH equations

As a final step, we combine the regression model for conditional mean and the GARCH model for conditional variance:

\[ \begin{cases} \text{Ret}_t &= \gamma + \gamma_1\cdot \text{Ret}_{t-1} + \gamma_2\cdot \text{Inf}_{t-1} + \gamma_3\cdot \text{DT.bill}_{t-1} + \epsilon_t \\ \epsilon_t &= \sigma_t \cdot w_t \\ \sigma_t^2 &= \omega + \alpha_1 \epsilon_{t-1}^2 + \beta_1 \sigma_{t-1}^2 \end{cases} \]

Save the data subset of exogenous variables. Because we need to include their lags in the right-hand side - we need to exclude the last row:

spec.dt <- as.matrix(SP)[-nrow(SP), 1:3]
head(spec.dt)
##           Ret       Inf. DT.bill
## [1,] 2.678492  0.3710579   -0.21
## [2,] 2.584550 -0.3710579    0.06
## [3,] 4.448416  0.0000000   -0.06
## [4,] 4.420518  0.0000000   -0.21
## [5,] 1.221139  0.0000000   -0.12
## [6,] 4.459217  0.0000000    0.08

Specify the model in R:

suppressPackageStartupMessages({library(rugarch)})
mdl.spec.garch <- ugarchspec(
  variance.model = list(garchOrder = c(1, 1)),
  mean.model = list(armaOrder = c(0,0),
                    include.mean = TRUE,
                    external.regressors = spec.dt)
  )
mdl.spec.garch
## 
## *---------------------------------*
## *       GARCH Model Spec          *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## ------------------------------------
## GARCH Model      : sGARCH(1,1)
## Variance Targeting   : FALSE 
## 
## Conditional Mean Dynamics
## ------------------------------------
## Mean Model       : ARFIMA(0,0,0)
## Include Mean     : TRUE 
## GARCH-in-Mean        : FALSE 
## Exogenous Regressor Dimension: 3
## 
## Conditional Distribution
## ------------------------------------
## Distribution :  norm 
## Includes Skew    :  FALSE 
## Includes Shape   :  FALSE 
## Includes Lambda  :  FALSE

where:

  • variance.model is the \(GARCH(1,1)\) model specification for our residual volatility;
  • mean.model is the mean equation for our returns with:
    • armaOrder - an \(ARMA\) model in case we need to specify one;
    • include.mean - if we would like to include a constant/intercept in our model;
    • external.regressors - we specify the exogenous variables. Note that we are including the lag of the endogenous regressor as well, so it also needs to be included in the model.

Because we are including lags, our endogenous data will start from \(t = 2,...,T\), so that the lagged exogenous data would be from \(t = 1,..., T-1\):

mdl.spec.garch.fit <- ugarchfit(data = SP[-1, "Ret"], spec = mdl.spec.garch)
mdl.spec.garch.fit
## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : sGARCH(1,1)
## Mean Model   : ARFIMA(0,0,0)
## Distribution : norm 
## 
## Optimal Parameters
## ------------------------------------
##         Estimate  Std. Error  t value Pr(>|t|)
## mu       1.31263    0.204030   6.4335 0.000000
## mxreg1   0.18528    0.045487   4.0731 0.000046
## mxreg2  -1.31787    0.456021  -2.8899 0.003853
## mxreg3  -1.19690    0.309909  -3.8621 0.000112
## omega    1.19923    0.536503   2.2353 0.025400
## alpha1   0.11996    0.041406   2.8970 0.003767
## beta1    0.77559    0.068169  11.3775 0.000000
## 
## Robust Standard Errors:
##         Estimate  Std. Error  t value Pr(>|t|)
## mu       1.31263    0.230400   5.6972 0.000000
## mxreg1   0.18528    0.043918   4.2187 0.000025
## mxreg2  -1.31787    0.547165  -2.4085 0.016016
## mxreg3  -1.19690    0.361031  -3.3152 0.000916
## omega    1.19923    0.598210   2.0047 0.044996
## alpha1   0.11996    0.047766   2.5113 0.012028
## beta1    0.77559    0.069713  11.1254 0.000000
## 
## LogLikelihood : -1472.298 
## 
## Information Criteria
## ------------------------------------
##                    
## Akaike       5.1814
## Bayes        5.2347
## Shibata      5.1811
## Hannan-Quinn 5.2022
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                    0.08685  0.7682
## Lag[2*(p+q)+(p+q)-1][2]   0.26510  0.8155
## Lag[4*(p+q)+(p+q)-1][5]   1.49164  0.7420
## d.o.f=0
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                   0.002882  0.9572
## Lag[2*(p+q)+(p+q)-1][5]  1.322472  0.7834
## Lag[4*(p+q)+(p+q)-1][9]  4.093256  0.5731
## d.o.f=2
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[3] 0.0002154 0.500 2.000  0.9883
## ARCH Lag[5] 0.8432633 1.440 1.667  0.7800
## ARCH Lag[7] 3.7926024 2.315 1.543  0.3775
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  0.8629
## Individual Statistics:              
## mu     0.13023
## mxreg1 0.06168
## mxreg2 0.13919
## mxreg3 0.02595
## omega  0.11633
## alpha1 0.09033
## beta1  0.15680
## 
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic:          1.69 1.9 2.35
## Individual Statistic:     0.35 0.47 0.75
## 
## Sign Bias Test
## ------------------------------------
##                    t-value    prob sig
## Sign Bias           0.8170 0.41426    
## Negative Sign Bias  0.8844 0.37684    
## Positive Sign Bias  0.7499 0.45362    
## Joint Effect        8.1065 0.04386  **
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     27.77      0.08787
## 2    30     47.65      0.01597
## 3    40     55.76      0.03990
## 4    50     59.56      0.14348
## 
## 
## Elapsed time : 0.7686911

We get a large output. Let’s extract only the relevant information.

The coefficients:

round(mdl.spec.garch.fit@fit$robust.matcoef, 5)
##         Estimate  Std. Error  t value Pr(>|t|)
## mu       1.31263     0.23040  5.69719  0.00000
## mxreg1   0.18528     0.04392  4.21873  0.00002
## mxreg2  -1.31787     0.54716 -2.40855  0.01602
## mxreg3  -1.19690     0.36103 -3.31523  0.00092
## omega    1.19923     0.59821  2.00469  0.04500
## alpha1   0.11996     0.04777  2.51132  0.01203
## beta1    0.77559     0.06971 11.12540  0.00000

Here we extracted the coefficients with Robust Standard Errors - unbiased standard errors of coefficients under the hypothesis of residual heteroscedasticity (i.e. when the residual variance is not the same across all observation points). We see that all the p-values are less than 0.05, so all of the coefficients are statistically significant.

We also get Ljung-Box tests on the standardized (squared) residuals we conclude that the residuals are \(WN\) p-value > 0.05 for all cases).

The Weighted ARCH LM Tests for ARCH effects testing conclude that the model residuals do not have any more ARCH effects.

So, we have finalized our model. Now, we can also create its forecasts:

forc.garch.mdl <- ugarchforecast(mdl.spec.garch.fit, n.ahead = 5)
forc.garch.mdl
## 
## *------------------------------------*
## *       GARCH Model Forecast         *
## *------------------------------------*
## Model: sGARCH
## Horizon: 5
## Roll Steps: 0
## Out of Sample: 0
## 
## 0-roll forecast [T0=1971-07-26 02:00:00]:
##     Series Sigma
## T+1  1.313 5.580
## T+2  1.313 5.393
## T+3  1.313 5.220
## T+4  1.313 5.059
## T+5  1.313 4.912

Where we forecast both the return Ret and the volatility Sigma. Finally, we can plot the data:

par(mfrow = c(2, 1))
plot(forc.garch.mdl, which = 1)
plot(forc.garch.mdl, which = 3)

Note - if we do not specify which = ..., we will be greeted with the following prompt:

Make a plot selection (or 0 to exit): 

1:   Time Series Prediction (unconditional)
2:   Time Series Prediction (rolling)
3:   Sigma Prediction (unconditional)
4:   Sigma Prediction (rolling)

where we would have to input the number 1-4 to specify which plot we would like to see.

Regarding the rolling forecasts - they are used for cross-validation. From the documentation:

The forecast is based on the expected value of the innovations and hence the density chosen. One step ahead forecasts are based on the value of the previous data, while n-step ahead (n>1) are based on the unconditional expectation of the models.
The ability to roll the forecast 1 step at a time is implemented with the n.roll argument which controls how many times to roll the n.ahead forecast. The default argument of n.roll = 0 denotes no rolling and returns the standard n.ahead forecast. Critically, since n.roll depends on data being available from which to base the rolling forecast, the ugarchfit function needs to be called with the argument out.sample being at least as large as the n.roll argument, or in the case of a specification being used instead of a fit object, the out.sample argument directly in the forecast function.

Some info on rolling forecasts can be found here. The general idea is:

  1. Select a sub-sample size of your data \(m < T\).
  2. Select an increment between successive rolling windows. Let’s say we select it equal to 1.
  3. The total number of sub-samples generated by this is \(T - m + 1\):
    • the first sub-sample is for periods \(t = 1,...,m\)
    • the second sub-sample is for periods \(t = 2, ..., m + 1\), etc.
  4. Select the forecast period \(h\).
  5. For each sub-sample \(n = 1, ..., T - m + 1\):
    • Exclude the last \(h\) observations from your sub-sample and estimate the model.
    • Calculate the forecasts for \(h\) periods.
    • Calculate the error between the actual values and the forecasted values : \(e_{n,j} = y_{n + m - h + j} - \hat{y}_{n + m - h + j}\), \(j = 1,..., h\)
  6. Calculate the forecast accuracy. For example, for each forecast period \(j = 1, ..., h\) calculate the forecast \(RMSE_j = \sqrt{\dfrac{1}{T-m+1} \sum_{j = 1}^{T-m+1} e_{n,j}^2}\).

This allows us to examine the forecast accuracy of the model.

A similar test for the coefficients and their confidence intervals of the model lets us test the model stability over time.