(9-3) ARMA-GARCH

Nao Mimoto - Dept. of Statistics : The University of Akron

1. ARMA-GARCH model

Let $Y_t$ be the observatoin. (Log-difference of the financial time series) Then, \[\begin{align*} Y_t & \sim \mbox{ARMA}(p,q)\\\\ Y_t &= \phi_1 Y_{t-1} + e_t + \theta_1 e_{t-1} \end{align*}\] with GARCH error \[\begin{align*} e_t &= \hspace{3mm} \sigma_t \epsilon_t \hspace{5mm} \epsilon_t \sim N(0,1)\\\\ \sigma_t^2 \hspace{3mm} = \hspace{3mm} \omega + \alpha Y_t^2 + \beta \sigma_t^2 \end{align*}\]

$Y_t$ is ARMA with error $\epsilon_t$. $\epsilon_t$ is GARCH. \[\begin{align*} \Phi(B) Y_t &= \Theta(B) \epsilon_t \\ \\ \epsilon_t &= \sigma_t e_t \hspace{10mm} e_t \sim_{iid} N(0,1) \\\\ \sigma_t^2 &= \omega + \alpha \epsilon_{t-1}^2 + \beta \sigma^2_{t-1} \end{align*}\]

2. Example: SPY

Daily Price of SP500 ETF (SPY) from Jan 02 2000 to Dec 31 2014

  library(quantmod)
  source('https://nmimoto.github.io/R/TS-00.txt')


  getSymbols("SPY")    #- SP500 download from Yahoo!

## [1] "SPY"

  SPY <- Ad(SPY)['2010::']

  is.ts(SPY)      # not ts object

## [1] FALSE

  is.xts(SPY)     # its xts object

## [1] TRUE

  plot( SPY )

  plot( log(SPY) )

  plot( diff( log(SPY) ) )

  library(fGarch)

  Y <- diff( log(SPY) )[-1]     # remove the first diff for NA
  
  #- Estimae parameters of ARMA and GARCH at the same time 
  Fit01 <-  garchFit(~ arma(5,5) + garch(1,1), data=Y, cond.dist="norm",  include.mean = FALSE, trace = FALSE)
  summary(Fit01)

## 
## Title:
##  GARCH Modelling 
## 
## Call:
##  garchFit(formula = ~arma(5, 5) + garch(1, 1), data = Y, cond.dist = "norm", 
##     include.mean = FALSE, trace = FALSE) 
## 
## Mean and Variance Equation:
##  data ~ arma(5, 5) + garch(1, 1)
## <environment: 0x7ff513b0c208>
##  [data = Y]
## 
## Conditional Distribution:
##  norm 
## 
## Coefficient(s):
##         ar1          ar2          ar3          ar4          ar5  
##  5.7637e-02   2.5862e-01  -2.9950e-01   2.6887e-02   9.5685e-01  
##         ma1          ma2          ma3          ma4          ma5  
## -7.1717e-02  -2.3916e-01   3.0287e-01  -3.2324e-02  -9.5271e-01  
##       omega       alpha1        beta1  
##  3.8302e-06   1.8922e-01   7.7602e-01  
## 
## Std. Errors:
##  based on Hessian 
## 
## Error Analysis:
##          Estimate  Std. Error    t value Pr(>|t|)    
## ar1     5.764e-02   1.587e-05   3631.373  < 2e-16 ***
## ar2     2.586e-01   1.586e-05  16307.742  < 2e-16 ***
## ar3    -2.995e-01   1.546e-05 -19374.512  < 2e-16 ***
## ar4     2.689e-02   1.581e-05   1700.958  < 2e-16 ***
## ar5     9.569e-01   1.588e-05  60245.048  < 2e-16 ***
## ma1    -7.172e-02   1.684e-05  -4258.421  < 2e-16 ***
## ma2    -2.392e-01   1.685e-05 -14194.731  < 2e-16 ***
## ma3     3.029e-01   1.644e-05  18418.088  < 2e-16 ***
## ma4    -3.232e-02   1.684e-05  -1919.603  < 2e-16 ***
## ma5    -9.527e-01   1.685e-05 -56538.900  < 2e-16 ***
## omega   3.830e-06   5.436e-07      7.046 1.84e-12 ***
## alpha1  1.892e-01   1.887e-02     10.030  < 2e-16 ***
## beta1   7.760e-01   1.843e-02     42.100  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log Likelihood:
##  8779.44    normalized:  3.391055 
## 
## Description:
##  Mon Apr 20 20:44:11 2020 by user:  
## 
## 
## Standardised Residuals Tests:
##                                 Statistic p-Value  
##  Jarque-Bera Test   R    Chi^2  597.138   0        
##  Shapiro-Wilk Test  R    W      0.9728211 0        
##  Ljung-Box Test     R    Q(10)  11.01984  0.3559744
##  Ljung-Box Test     R    Q(15)  17.83918  0.2712147
##  Ljung-Box Test     R    Q(20)  26.22266  0.1585834
##  Ljung-Box Test     R^2  Q(10)  7.343764  0.6926416
##  Ljung-Box Test     R^2  Q(15)  11.54014  0.7134504
##  Ljung-Box Test     R^2  Q(20)  12.36018  0.9031213
##  LM Arch Test       R    TR^2   9.118898  0.6927428
## 
## Information Criterion Statistics:
##       AIC       BIC       SIC      HQIC 
## -6.772067 -6.742647 -6.772117 -6.761405

  sig.t  <- xts(Fit01@sigma.t,                 order.by=index(Y))
                               #-- estimated sig_t fo GARCH
  res01  <- xts(Fit01@residuals/Fit01@sigma.t, order.by=index(Y))
                               #-- this is the (standardized) ARMA-GARCH residuals

  Randomness.tests(res01)

##   B-L test H0: the series is uncorrelated
##   M-L test H0: the square of the series is uncorrelated
##   J-B test H0: the series came from Normal distribution
##   SD         : Standard Deviation of the series

##       BL15  BL20  BL25  ML15  ML20 JB    SD
## [1,] 0.271 0.159 0.037 0.713 0.903  0 0.997

  ##plot(Fit1)            #-- choose option 13 for residual qq plot

    ##  Fit01@fit$par                  # estimated parameters
    ##  Fit01@residuals                # this is not GARCH residuals! This is same as Y.
    ##  Fit01@sigma.t                  # estimated sig_t
    ##  Fit01@fit$ics                   # AIC and BIC are here
    ##  Fit01@residuals/Fit1@sigma.t   # this is the (standardized) GARCH residuals
    


  #- Estimae parameters of ARMA and GARCH at the same time 
  Fit02 <-  garchFit(~ garch(1,1), data=Y, cond.dist="norm",  include.mean = FALSE, trace = FALSE)
  summary(Fit02)

## 
## Title:
##  GARCH Modelling 
## 
## Call:
##  garchFit(formula = ~garch(1, 1), data = Y, cond.dist = "norm", 
##     include.mean = FALSE, trace = FALSE) 
## 
## Mean and Variance Equation:
##  data ~ garch(1, 1)
## <environment: 0x7ff51770c2a0>
##  [data = Y]
## 
## Conditional Distribution:
##  norm 
## 
## Coefficient(s):
##      omega      alpha1       beta1  
## 3.8225e-06  1.7627e-01  7.8681e-01  
## 
## Std. Errors:
##  based on Hessian 
## 
## Error Analysis:
##         Estimate  Std. Error  t value Pr(>|t|)    
## omega  3.823e-06   5.542e-07    6.897 5.32e-12 ***
## alpha1 1.763e-01   1.763e-02    9.997  < 2e-16 ***
## beta1  7.868e-01   1.808e-02   43.512  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log Likelihood:
##  8748.313    normalized:  3.379032 
## 
## Description:
##  Mon Apr 20 20:44:12 2020 by user:  
## 
## 
## Standardised Residuals Tests:
##                                 Statistic p-Value  
##  Jarque-Bera Test   R    Chi^2  632.927   0        
##  Shapiro-Wilk Test  R    W      0.9711714 0        
##  Ljung-Box Test     R    Q(10)  10.48633  0.3989101
##  Ljung-Box Test     R    Q(15)  15.75681  0.3983959
##  Ljung-Box Test     R    Q(20)  23.56569  0.2618783
##  Ljung-Box Test     R^2  Q(10)  10.6118   0.3885507
##  Ljung-Box Test     R^2  Q(15)  14.87802  0.4602388
##  Ljung-Box Test     R^2  Q(20)  15.63522  0.7389809
##  LM Arch Test       R    TR^2   13.2704   0.3496974
## 
## Information Criterion Statistics:
##       AIC       BIC       SIC      HQIC 
## -6.755746 -6.748957 -6.755749 -6.753285

  res02 <- Y/Fit02@sigma.t
  Randomness.tests(res02)

##   B-L test H0: the series is uncorrelated
##   M-L test H0: the square of the series is uncorrelated
##   J-B test H0: the series came from Normal distribution
##   SD         : Standard Deviation of the series

##       BL15  BL20  BL25 ML15  ML20 JB    SD
## [1,] 0.398 0.262 0.047 0.46 0.739  0 0.998

  #- plot log-return and sigma_t -
    sigma.upper = xts(1.96*Fit01@sigma.t,  order.by=index(Y))
    sigma.lower = xts(-1.96*Fit01@sigma.t, order.by=index(Y))
    
    plot( cbind(Y, sigma.upper, sigma.lower),
        col=c("black", "red", "red"), lwd=c(2,1,1) )

    plot(  as.numeric(          Y["2018::"]), type="h", lwd=2)
    lines( as.numeric(sigma.upper["2018::"]), col="red")
    lines( as.numeric(sigma.lower["2018::"]), col="red")

  #--- h-step ahead prediction from ARMA-GARCH ---
  h = 15    #- predict 15 days ahead

  X.pred  <- xts(predict(Fit2, n.ahead=h)[,1], order.by=index(X1)[length(X1)]+(1:h))
  SD.pred <- xts(predict(Fit2, n.ahead=h)[,3], order.by=index(X1)[length(X1)]+(1:h))

  #- plot log-return overlay with sigma_t -
  plot(rbind(X1["2014"], X.pred))
  lines(X.pred, col="red")
  lines( 1.96*SD.pred, type="l", col="red")
  lines(-1.96*SD.pred, type="l", col="red")

  #- Rolling 1-step prediction of sig.t and ARMA
  Y <- diff( log(SPY["2013::"]) )[-1]    # Original data
  window.size <- 250                         # Window size for estimation

    Y.pred1   <- numeric(0)
    sig.pred1 <- numeric(0)
    date.stp  <- numeric(0)
    for (i in 1:(length(Y)-window.size)){
      T <- Y[(1:window.size)+(i-1)]     
      n <- length(T)
      out1  <-  garchFit(~ arma(0,3) + garch(1,1), data=T, cond.dist="norm", include.mean = FALSE, trace = FALSE)
      Y.pred1[i]   <- predict(out1)[1,1]   #  ARMA prediction
      sig.pred1[i] <- predict(out1)[1,3]   #  sig.t prediction
      date.stp[i]  <- index(Y)[window.size+i]
    }

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

## Warning in sqrt(diag(fit$cvar)): NaNs produced

    Y.pred   <- xts(Y.pred1,   order.by=as.Date(date.stp))
    sig.pred <- xts(sig.pred1, order.by=as.Date(date.stp))



  #- Plot and compare estimated vs rolling predicted
  plot(cbind(Y[             '2016'],
             Y.pred[        '2016'],
             1.96*sig.pred[ '2016'],
            -1.96*sig.pred[ '2016']),
     col=c("black", "red", "blue"), lwd=c(1,1,1),
     main="ARMA-GARCH")

  time.period="2019::"
    plot(cbind(Y[          time.period],
               Y.pred[     time.period],
             1.96*sig.pred[time.period],
            -1.96*sig.pred[time.period]),
     col=c("black", "red", "blue", "blue"), lwd=c(1,1,1),
     main="ARMA-GARCH")

  time.period="2019::"  
    plot(  as.numeric(             Y[time.period]), type="h", lwd=2, ylim=c(-.12, .12),
            xlab="diff( log(Y) )")
    lines( as.numeric(        Y.pred[time.period]), col="red")        
    lines( as.numeric( 1.96*sig.pred[time.period]), col="blue")
    lines( as.numeric(-1.96*sig.pred[time.period]), col="blue")

  mean( sign(Y[index(Y.pred)])==sign(Y.pred) )  #- Can Y.pred guess up/down ?

## [1] 0.4794953

  mean( sign(Y[index(Y.pred)])==1 )             #- num of days SPY went up

## [1] 0.5507886

## [1] 0.1165644

## [1] 0.4210526

## [1] 0.5526316

3. Cheking Conditional Distribution

Empirical Distribution Function

Sample version of CDF \[ \hat F_n(x) \hspace{3mm} = \hspace{3mm} \frac{1}{n} \sum_{i=1}^n I(x_i < x) \]

  X22 <- rnorm(10)
  X22

##  [1]  1.2328929  1.4598186 -0.2701371 -1.3288165 -0.1796356  1.4949569
##  [7] -0.3288075 -0.6891992 -0.2188739  3.5027069

  plot(ecdf(X22))

Kolmogolov-Smirnov Test

In iid setting, K-S test can be used to test hypothesis $X_i \sim F$, and

\[ H_0: F = F_0 \hspace{5mm} vs \hspace{5mm} H_A: F \ne F_0 \hspace{10mm} \mbox{{\scriptsize ($F_0$ is completely specified) } } \]
K-S test uses Empirical Distibution (EDF, ECDF), and calculate the test statistic \[ KS \hspace{3mm} = \hspace{3mm} \sup_x |\hat F_n(x) - F(x) | \]
Where is $(\frac{i}{n}, X_{(i)})$ and $(\frac{i-1}{n}, X_{(i)})$?

\[\begin{align*} D_n &= \sup_x |\hat F_n(x) - F_0(x)| \\\\ &= \max_i \Big( \Big|\frac{i}{n} - F_0(x_{(i)}) \Big| , \hspace{3mm} \Big|\frac{i-1}{n} - F_0(x_{(i)}) \Big| \Big) \\ \\ \end{align*}\]
If $X_i$ are indeed from $F_0$, \[ D_n \hspace{3mm} =_d \hspace{3mm} \max_i \Big( \Big|\frac{i}{n} - U_{(i)} \Big| , \hspace{3mm} \Big|\frac{i-1}{n} - U_{(i)} \Big| \Big) \] with $U_i \sim_{iid} Unif(0,1)$.
No matter what $F$ is, distribution of the test statistic KS under the null is known.
KS test with completely specified $F_0$ is a distribution-free test.
When $F_0$ contains nuisance parameter, Distribution of $D_n$ under the null depends on $F$. KS test becomes only asymptotically distribution free test. For finite $n$, null distribution must be computed by Monte Carlo simulation.

K-S test in Time Series Setting

In time sries, we want to use K-S test for checking distribution of innovations. (errors).
Since we assume $e_t$ are i.i.d. noise, if we can observe $e_t$, K-S test can be directly applicable.
However, in Time Series analysis $e_t$ are not observable, and only residuals $\hat e_t$ are available.
Can we still use K-S test on $\hat e_t$?

  Fn <- ecdf(as.numeric(res01))   #- Fn becomes an function
  plot( ecdf(as.numeric(res01)) )

  t=seq(-4,4, .01)
  plot(t,  Fn(t), type="l" )
  lines(t, pnorm(t,0,1), col="red")

  ks.test(as.numeric(res01), "pnorm")

## Warning in ks.test(as.numeric(res01), "pnorm"): ties should not be
## present for the Kolmogorov-Smirnov test

## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  as.numeric(res01)
## D = 0.048574, p-value = 9.891e-06
## alternative hypothesis: two-sided

TS Class Web Page – R resource page