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} + \epsilon_t + \theta_1 \epsilon_{t-1} \end{align*}\] with GARCH error \[\begin{align*} \epsilon_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 \epsilon_t \hspace{10mm} \epsilon_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) ) )


2a. Fit ARMA-GARCH

  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)
## Warning in arima(.series$x, order = c(u, 0, v), include.mean = include.mean):
## possible convergence problem: optim gave code = 1
## Warning in sqrt(diag(fit$cvar)): NaNs produced
  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: 0x160dc3d58>
##  [data = Y]
## 
## Conditional Distribution:
##  norm 
## 
## Coefficient(s):
##         ar1          ar2          ar3          ar4          ar5          ma1  
## -7.4255e-01  -1.5549e-01  -8.7181e-01  -5.6246e-01   3.7109e-01   7.1967e-01  
##         ma2          ma3          ma4          ma5        omega       alpha1  
##  1.3684e-01   8.6172e-01   5.3251e-01  -3.9581e-01   3.7362e-06   1.4904e-01  
##       beta1  
##  8.1578e-01  
## 
## Std. Errors:
##  based on Hessian 
## 
## Error Analysis:
##          Estimate  Std. Error  t value Pr(>|t|)    
## ar1    -7.425e-01         NaN      NaN      NaN    
## ar2    -1.555e-01         NaN      NaN      NaN    
## ar3    -8.718e-01         NaN      NaN      NaN    
## ar4    -5.625e-01         NaN      NaN      NaN    
## ar5     3.711e-01         NaN      NaN      NaN    
## ma1     7.197e-01         NaN      NaN      NaN    
## ma2     1.368e-01         NaN      NaN      NaN    
## ma3     8.617e-01         NaN      NaN      NaN    
## ma4     5.325e-01         NaN      NaN      NaN    
## ma5    -3.958e-01         NaN      NaN      NaN    
## omega   3.736e-06   4.809e-07    7.769 7.99e-15 ***
## alpha1  1.490e-01   1.292e-02   11.531  < 2e-16 ***
## beta1   8.158e-01   1.414e-02   57.710  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log Likelihood:
##  13926.6    normalized:  3.318227 
## 
## Description:
##  Sun Sep 13 15:24:57 2026 by user:  
## 
## 
## 
## Standardised Residuals Tests:
##                                   Statistic   p-Value
##  Jarque-Bera Test   R    Chi^2  1015.142128 0.0000000
##  Shapiro-Wilk Test  R    W         0.972619 0.0000000
##  Ljung-Box Test     R    Q(10)     7.188338 0.7075531
##  Ljung-Box Test     R    Q(15)    11.849116 0.6904120
##  Ljung-Box Test     R    Q(20)    18.354837 0.5640458
##  Ljung-Box Test     R^2  Q(10)    13.867852 0.1791051
##  Ljung-Box Test     R^2  Q(15)    16.191197 0.3694597
##  Ljung-Box Test     R^2  Q(20)    16.790092 0.6665588
##  LM Arch Test       R    TR^2     15.163040 0.2326412
## 
## Information Criterion Statistics:
##       AIC       BIC       SIC      HQIC 
## -6.630260 -6.610615 -6.630279 -6.623313
  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.69 0.564 0.072 0.369 0.667  0 0.998
  ##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: 0x13c81c9b0>
##  [data = Y]
## 
## Conditional Distribution:
##  norm 
## 
## Coefficient(s):
##      omega      alpha1       beta1  
## 3.7798e-06  1.5098e-01  8.1387e-01  
## 
## Std. Errors:
##  based on Hessian 
## 
## Error Analysis:
##         Estimate  Std. Error  t value Pr(>|t|)    
## omega  3.780e-06   4.755e-07    7.949 1.78e-15 ***
## alpha1 1.510e-01   1.290e-02   11.703  < 2e-16 ***
## beta1  8.139e-01   1.396e-02   58.291  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log Likelihood:
##  13917.88    normalized:  3.316149 
## 
## Description:
##  Sun Sep 13 15:24:57 2026 by user:  
## 
## 
## 
## Standardised Residuals Tests:
##                                    Statistic   p-Value
##  Jarque-Bera Test   R    Chi^2  1020.8847125 0.0000000
##  Shapiro-Wilk Test  R    W         0.9727087 0.0000000
##  Ljung-Box Test     R    Q(10)    10.3053323 0.4141262
##  Ljung-Box Test     R    Q(15)    16.2668287 0.3645432
##  Ljung-Box Test     R    Q(20)    24.4955019 0.2214173
##  Ljung-Box Test     R^2  Q(10)    13.9808894 0.1738651
##  Ljung-Box Test     R^2  Q(15)    16.1999190 0.3688908
##  Ljung-Box Test     R^2  Q(20)    16.8389205 0.6634071
##  LM Arch Test       R    TR^2     14.8604609 0.2491563
## 
## Information Criterion Statistics:
##       AIC       BIC       SIC      HQIC 
## -6.630869 -6.626336 -6.630870 -6.629266
  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.365 0.221 0.03 0.369 0.663  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
## 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.4901347
  mean( sign(Y[index(Y.pred)])==1 )             #- num of days SPY went up
## [1] 0.5502662
## [1] 0.2466391
## [1] 0.4779874
## [1] 0.5702306


3. Cheking Conditional Distribution

3a. 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] -2.1015322  0.8825195 -1.2327865  1.5668140  0.5881110  1.0515181
##  [7] -0.7412360 -1.1190486 -0.1183287 -0.3939993
  plot(ecdf(X22))


3b. 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.

3c. 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 \(\epsilon_t\) are i.i.d. noise, if we can observe \(\epsilon_t\), K-S test can be directly applicable.

However, in Time Series analysis \(\epsilon_t\) are not observable, and only residuals \(\hat \epsilon_t\) are available.

Can we still use K-S test on \(\hat \epsilon_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.default(as.numeric(res01), "pnorm"): ties should not be
## present for the one-sample Kolmogorov-Smirnov test
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  as.numeric(res01)
## D = 0.081575, p-value < 2.2e-16
## alternative hypothesis: two-sided