Forecasting Time Series

Bear Stearns Data Analysis

Preamble

library("forecast")
library("rugarch")

Loading required package: methods

Loading required package: parallel

Attaching package: 'rugarch'

The following object is masked from 'package:stats':

    sigma

library("tseries")

Bear Stearns Daily Returns

January 2, 2008 - March 28, 2008

data <- read.csv("http://ptrckprry.com/course/forecasting/data/bearstearns.csv")
date <- as.Date(data$date, format="%m/%d/%Y")
time <- 1:length(date)
BSCret <- data$BSCret

Time Series Plot

plot(date, BSCret, type="l", col=2)

ACF, PACF

Returns

par(mfrow=c(1,2))
Acf(BSCret)
Pacf(BSCret)

Diff. Returns

par(mfrow=c(1,2))
Acf(diff(BSCret))
Pacf(diff(BSCret))

ARIMA Model for Mean

AICC Values

d <- 0
for (include.constant in c(FALSE, TRUE)) {
    for (p in 0:4) {
        for (q in 0:4) {
            fit <- Arima(BSCret, c(p,d,q), method="ML",
                         include.constant=include.constant)
             cat("ARIMA",
                 "(", p, ",", d, ",", q, ")",
                 "(constant=", include.constant, ")",
                 " : ", fit$aicc, "\n", sep="")
        }
    }
}

ARIMA(0,0,0)(constant=FALSE) : -35.16031
ARIMA(0,0,1)(constant=FALSE) : -34.62698
ARIMA(0,0,2)(constant=FALSE) : -35.99712
ARIMA(0,0,3)(constant=FALSE) : -42.84813
ARIMA(0,0,4)(constant=FALSE) : -47.21773
ARIMA(1,0,0)(constant=FALSE) : -34.53045
ARIMA(1,0,1)(constant=FALSE) : -40.18723
ARIMA(1,0,2)(constant=FALSE) : -38.01908
ARIMA(1,0,3)(constant=FALSE) : -42.31904
ARIMA(1,0,4)(constant=FALSE) : -49.44943
ARIMA(2,0,0)(constant=FALSE) : -32.32843
ARIMA(2,0,1)(constant=FALSE) : -38.12093
ARIMA(2,0,2)(constant=FALSE) : -37.41074
ARIMA(2,0,3)(constant=FALSE) : -44.39582
ARIMA(2,0,4)(constant=FALSE) : -47.69503
ARIMA(3,0,0)(constant=FALSE) : -30.47742
ARIMA(3,0,1)(constant=FALSE) : -39.23293
ARIMA(3,0,2)(constant=FALSE) : -42.51048
ARIMA(3,0,3)(constant=FALSE) : -42.79488
ARIMA(3,0,4)(constant=FALSE) : -45.05117
ARIMA(4,0,0)(constant=FALSE) : -41.74055
ARIMA(4,0,1)(constant=FALSE) : -39.2671
ARIMA(4,0,2)(constant=FALSE) : -41.3395
ARIMA(4,0,3)(constant=FALSE) : -44.81976
ARIMA(4,0,4)(constant=FALSE) : -42.341
ARIMA(0,0,0)(constant=TRUE) : -33.21241
ARIMA(0,0,1)(constant=TRUE) : -32.55884
ARIMA(0,0,2)(constant=TRUE) : -33.89324
ARIMA(0,0,3)(constant=TRUE) : -40.47506
ARIMA(0,0,4)(constant=TRUE) : -45.0509
ARIMA(1,0,0)(constant=TRUE) : -32.45962
ARIMA(1,0,1)(constant=TRUE) : -38.02689
ARIMA(1,0,2)(constant=TRUE) : -35.78304
ARIMA(1,0,3)(constant=TRUE) : -39.87188
ARIMA(1,0,4)(constant=TRUE) : -48.92238
ARIMA(2,0,0)(constant=TRUE) : -30.18139
ARIMA(2,0,1)(constant=TRUE) : -35.89054
ARIMA(2,0,2)(constant=TRUE) : -36.90732
ARIMA(2,0,3)(constant=TRUE) : -41.89337
ARIMA(2,0,4)(constant=TRUE) : -46.60697
ARIMA(3,0,0)(constant=TRUE) : -28.22613
ARIMA(3,0,1)(constant=TRUE) : -36.84817
ARIMA(3,0,2)(constant=TRUE) : -40.06975
ARIMA(3,0,3)(constant=TRUE) : -40.20443
ARIMA(3,0,4)(constant=TRUE) : -44.18157
ARIMA(4,0,0)(constant=TRUE) : -39.55506
ARIMA(4,0,1)(constant=TRUE) : -37.01267
ARIMA(4,0,2)(constant=TRUE) : -39.44428
ARIMA(4,0,3)(constant=TRUE) : -42.1212
ARIMA(4,0,4)(constant=TRUE) : -41.30252

Selected Model for Mean

(fit.mean <- Arima(BSCret, c(1, 0, 4), include.constant=FALSE))

Series: BSCret 
ARIMA(1,0,4) with zero mean     

Coefficients:
         ar1      ma1      ma2     ma3      ma4
      0.4625  -0.2547  -0.0151  0.2768  -0.6937
s.e.  0.1845   0.1467   0.1028  0.0950   0.1049

sigma^2 estimated as 0.02135:  log likelihood=31.52
AIC=-51.03   AICc=-49.45   BIC=-38.47

Residuals

resid <- as.numeric(residuals(fit.mean))
plot(date, resid, type="l", col=3)

Ljung-Box Tests

Box.test(resid, lag=12, type="Ljung-Box", fitdf=5)

    Box-Ljung test

data:  resid
X-squared = 3.4851, df = 7, p-value = 0.8368

Box.test(resid, lag=24, type="Ljung-Box", fitdf=5)

    Box-Ljung test

data:  resid
X-squared = 4.1192, df = 19, p-value = 0.9999

Box.test(resid, lag=36, type="Ljung-Box", fitdf=5)

    Box-Ljung test

data:  resid
X-squared = 6.9687, df = 31, p-value = 1

Box.test(resid, lag=48, type="Ljung-Box", fitdf=5)

    Box-Ljung test

data:  resid
X-squared = 10.745, df = 43, p-value = 1

Residual ACF, PACF

par(mfrow=c(1,2))
Acf(resid)
Pacf(resid)

Sq. Residual ACF, PACF

par(mfrow=c(1,2))
Acf(resid^2)
Pacf(resid^2)

GARCH Model for Shocks (Using tseries package)

AICC Values for ARCH Models

q <- 0:10
loglik <- rep(NA, length(q))
N <- length(resid)

for (i in 1:length(q)) {
    if (q[i] == 0) {
        loglik[i] <- -0.5 * N * (1 + log(2 * pi * mean(resid^2)))
    } else {
        fit <- garch(resid, c(0,q[i]), trace=FALSE)
        loglik[i] <- logLik(fit)
    }
}

Warning in sqrt(pred$e): NaNs produced

Warning in sqrt(pred$e): NaNs produced

k <- q + 1
aicc <- -2 * loglik  + 2 * k * N / (N - k - 1)

print(data.frame(q, loglik, aicc))

    q   loglik       aicc
1   0 32.88073  -63.69249
2   1 46.66435  -89.11817
3   2 44.76260  -83.09662
4   3 53.78042  -98.83356
5   4 59.39001 -107.66891
6   5 44.18812  -74.79133
7   6 44.87877  -73.60369
8   7 41.82517  -64.82682
9   8 40.44996  -59.29992
10  9 37.93360  -51.37741
11 10 37.58527  -47.67054

AICC Value for GARCH(1,1) Model

fit <- garch(resid, c(1,1), trace=FALSE)
loglik <- logLik(fit)
k <- 2
aicc <- -2 * loglik  + 2 * k * N / (N - k - 1)

print(data.frame(loglik, aicc))

    loglik      aicc
1 46.29492 -88.37932

Selected Model for Shocks

fit.var <- garch(resid, c(0,5), trace=FALSE)
summary(fit.var)

Call:
garch(x = resid, order = c(0, 5), trace = FALSE)

Model:
GARCH(0,5)

Residuals:
     Min       1Q   Median       3Q      Max 
-5.88399 -0.54685 -0.09173  0.39859  1.61239 

Coefficient(s):
    Estimate  Std. Error  t value Pr(>|t|)    
a0 5.419e-03   2.481e-03    2.184 0.028964 *  
a1 1.012e-01   2.614e-02    3.872 0.000108 ***
a2 7.765e-02   2.535e-01    0.306 0.759365    
a3 4.387e-02   3.264e-01    0.134 0.893080    
a4 1.523e-01   1.717e-01    0.887 0.374919    
a5 1.220e-15   1.485e-01    0.000 1.000000    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Diagnostic Tests:
    Jarque Bera Test

data:  Residuals
X-squared = 308.45, df = 2, p-value < 2.2e-16


    Box-Ljung test

data:  Squared.Residuals
X-squared = 4.6983, df = 1, p-value = 0.03019

Estimated Conditional Variances

ht <- fit.var$fit[,1]^2 # ARIMA-ARCH

par(mfrow=c(1,2))
plot(date, BSCret, type="l", col=2)
plot(date, ht, type="l", col=3)

ARCH Residuals

resid.arch <- resid / sqrt(ht)

Check Normality

boxplot(resid.arch)

hist(resid.arch)

qqnorm(resid.arch, col=2)
qqline(resid.arch, col=1, lty=2)

Check Autocorrelation

par(mfrow=c(1,2))
Acf(resid.arch)
Pacf(resid.arch)

Check Autocorrelation in Volatility

par(mfrow=c(1,2))
Acf(resid.arch^2)
Pacf(resid.arch^2)

ARIMA Forecast

forecast(fit.mean, h=1)

   Point Forecast      Lo 80     Hi 80      Lo 95     Hi 95
61     0.03603092 -0.1512096 0.2232714 -0.2503288 0.3223906

# By hand:
f1 <- as.numeric(forecast(fit.mean, h=1)$mean)

f1 + sqrt(fit.mean$sigma2) * c(-1.96, 1.96) # 95%

[1] -0.2503305  0.3223923

ARIMA-ARCH Forecast

Get omega, alpha

coef(fit.var)

          a0           a1           a2           a3           a4           a5 
5.419432e-03 1.012155e-01 7.764522e-02 4.386777e-02 1.523309e-01 1.220438e-15 

(omega <- coef(fit.var)[[1]])

[1] 0.005419432

(alpha <- coef(fit.var)[-1])

          a1           a2           a3           a4           a5 
1.012155e-01 7.764522e-02 4.386777e-02 1.523309e-01 1.220438e-15 

Get last 5 squared residuals

resid^2 # all

 [1] 7.308631e-04 9.257799e-05 2.028942e-03 2.767860e-04 4.190183e-03
 [6] 3.659430e-03 1.047687e-04 5.073931e-04 1.525619e-03 1.284868e-03
[11] 1.031298e-03 1.049198e-03 6.299390e-04 2.307987e-03 1.535209e-02
[16] 7.597199e-04 2.222645e-03 1.410402e-03 7.163183e-03 4.962642e-04
[21] 2.159323e-04 7.583050e-04 1.170522e-03 2.149947e-03 2.401038e-03
[26] 7.034550e-04 9.431637e-05 1.931921e-04 2.471243e-03 7.717163e-04
[31] 3.230007e-04 4.443429e-03 7.076370e-03 3.242646e-03 2.012913e-03
[36] 9.781121e-03 2.261636e-03 8.951150e-04 1.657402e-03 1.029122e-03
[41] 4.774644e-03 4.888628e-05 5.552887e-04 3.360819e-04 1.317528e-02
[46] 4.344694e-04 1.690150e-02 5.438516e-03 9.134266e-03 1.367884e-03
[51] 3.143726e-01 4.710913e-01 1.426707e-01 1.334274e-04 7.157952e-04
[56] 6.032774e-02 1.371442e-02 7.127600e-04 8.305394e-03 1.934357e-02

(resid^2)[60:56]

[1] 0.019343574 0.008305394 0.000712760 0.013714422 0.060327737

Compute 1-step-ahead conditional variance

(h1 <- omega + sum(alpha * (resid^2)[60:56]))

[1] 0.01014257

Construct forecast intervals

f1 <- as.numeric(forecast(fit.mean, h=1)$mean)

# ARIMA-ARCH inteval
f1 + sqrt(h1) * c(-1.96, 1.96) # 95%

[1] -0.1613614  0.2334232

95% Forecast Intervals

ft <- fitted.values(fit.mean)
sigma2 <- fit.mean$sigma2


plot(date, BSCret, type="l", col=2, ylim=c(-1.25, 1.25))

# ARIMA intervals
lines(date, ft + sqrt(sigma2) * 1.96, col=3)
lines(date, ft - sqrt(sigma2) * 1.96, col=3)

# ARIMA-ARCH intervals
lines(date, ft + sqrt(ht) * 1.96, col=4)
lines(date, ft - sqrt(ht) * 1.96, col=4)

#legend("bottomleft", inset=0.05,
#       legend=c("Series", "ARIMA PI", "ARIMA-ARCH PI"), lty=1, col=2:4)

Zoomed In on Crash

plot(date, BSCret, type="l", col=2,
     xlim=range(date[40:60]), ylim=c(-1.25, 1.25))

# ARIMA intervals
lines(date, ft + sqrt(sigma2) * 1.96, col=3)
lines(date, ft - sqrt(sigma2) * 1.96, col=3)

# ARIMA-ARCH intervals
lines(date, ft + sqrt(ht) * 1.96, col=4)
lines(date, ft - sqrt(ht) * 1.96, col=4)

legend("bottomleft", inset=0.05,
       legend=c("Series", "ARIMA PI", "ARIMA-ARCH PI"), lty=1, col=2:4)

How Many Failures?

# Failures rate for ARIMA:
sum(abs(BSCret - ft) > sqrt(sigma2) * 1.96, na.rm=TRUE)

[1] 3

# Failures for ARIMA-ARCH:
sum(abs(BSCret - ft) > sqrt(ht) * 1.96, na.rm=TRUE)

[1] 2