Variance Dynamics

Code
library(tsissm)
library(xts)
library(data.table)
library(tsaux)
library(knitr)
library(zoo)

Introduction

In Section 19.2.2 of Hyndman et al. (2008), the authors illustrate an extension of the innovations state space model incorporating exponential GARCH dynamics. The appeal of the exponential GARCH model lies in its ability to ensure positivity of the variance without requiring parameter constraints.

In the tsissm package, we instead opt for the standard (vanilla) GARCH model, applying bounds on both the GARCH parameters and the overall persistence of the variance process. While this may initially appear more complex, it is actually simpler in practice—especially when accommodating non-Gaussian distributions. For example, using Johnson’s SU distribution complicates the calculation of the expected value of the absolute standardized innovations, which lacks a closed form and is required in the exponential GARCH model. With a standard GARCH structure, this complexity is avoided.

That said, the decision to include variance dynamics should be made with care. Many nominal economic time series that exhibit apparent heteroscedasticity become much more homoscedastic when deflated by the Consumer Price Index (CPI) or a similar price index. Moreover, structural breaks and transitory changes can mimic heteroscedastic behavior, misleading standard statistical tests into detecting heteroscedasticity when none is truly present.

SPY Closing Price

Heteroscedasticity in financial returns typically stems from volatility clustering, asymmetric responses to news, and shifting market conditions. These characteristics make GARCH-type models a natural choice in financial econometrics.

In this example, we model the adjusted closing prices of the S&P 500 ETF (SPY) using the innovations state space model. The specification includes:

We also account for structural level shifts, such as those occurring around the COVID-19 pandemic, by including them as regressors.

Code
data("spy", package = "tsissm")
y <- as.xts(spy)
xreg <- auto_regressors(y["2014/"], frequency = 1, lambda = 0, sampling = "days", method = "full",
                         check.rank = TRUE, discard.cval = 3.5, maxit.iloop = 10, maxit.oloop = 10, types = "LS")

exc <- which(xreg$xreg["2020-02-03"] == 1)
xreg$xreg <- xreg$xreg[,-exc]
xreg$init <- xreg$init[-exc]

We begin by estimating two models—one with constant variance, and one with time-varying (GARCH) variance:

Code
spec_constant <- issm_modelspec(y["2014/"], slope = FALSE, seasonal = FALSE, ar = 2, ma = 0, xreg = xreg$xreg,
                       lambda = 0, variance = "constant", distribution = "jsu")
spec_constant$parmatrix[grepl("^kappa", parameters), initial := xreg$init]
mod_constant <- estimate(spec_constant, scores = FALSE)

spec_dynamic <- issm_modelspec(y["2014/"], slope = TRUE, seasonal = FALSE, ar = 1, ma = 0, xreg = xreg$xreg,
                       lambda = 0, variance = "dynamic", distribution = "jsu")
spec_dynamic$parmatrix[grepl("^kappa",parameters), initial := xreg$init]
mod_dynamic <- estimate(spec_dynamic, scores = FALSE)
print(paste0("AIC(Dynamic): ", round(AIC(mod_dynamic),1)," | AIC(Constant): ", round(AIC(mod_constant),1)))
#> [1] "AIC(Dynamic): 5059.5 | AIC(Constant): 5343.7"

Based on the AIC values, the model with dynamic variance is preferred—even though it involves more parameters—demonstrating the usefulness of capturing volatility dynamics.

Model Summary

You can obtain a clean, professional summary of the estimated model using the as_flextable() method:

Code
mod_dynamic |> summary() |> as_flextable()
ISSM Model: AAN+ARMA(1,0)+X(13)+GARCH(1,1)

Estimate

Std. Error

t value

Pr(>|t|)

α\alpha

0.1949

0.0634

3.0736

0.0021

**

β\beta

0.0013

0.0005

2.5449

0.0109

*

θ1\theta_{1}

0.7285

0.0742

9.8228

0.0000

***

κ1\kappa_{1}

-0.0398

0.0071

-5.5796

0.0000

***

κ2\kappa_{2}

-0.0603

0.0140

-4.2993

0.0000

***

κ3\kappa_{3}

-0.1387

0.0336

-4.1331

0.0000

***

κ4\kappa_{4}

0.0608

0.0300

2.0285

0.0425

*

κ5\kappa_{5}

0.0374

0.0269

1.3873

0.1654

κ6\kappa_{6}

-0.0587

0.0229

-2.5633

0.0104

*

κ7\kappa_{7}

0.0549

0.0202

2.7171

0.0066

**

κ8\kappa_{8}

0.0244

0.0188

1.2956

0.1951

κ9\kappa_{9}

-0.0612

0.0112

-5.4834

0.0000

***

κ10\kappa_{10}

-0.0457

0.0155

-2.9600

0.0031

**

κ11\kappa_{11}

-0.0489

0.0148

-3.3034

0.0010

***

κ12\kappa_{12}

-0.0476

0.0097

-4.9156

0.0000

***

κ13\kappa_{13}

0.0488

0.0108

4.5240

0.0000

***

η\eta

0.1491

0.0165

9.0197

0.0000

***

δ\delta

0.8337

0.0193

43.3078

0.0000

***

ζ\zeta

-0.8668

0.1687

-5.1394

0.0000

***

ν\nu

2.0636

0.1570

13.1431

0.0000

***

Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

sigma^2: 0.0001102

LogLik: -5011.492

AIC: 5059 | BIC: 1.021e+04

DAC : 48 | MAPE : 0.71

Model Equation

ytω=wxt1+j=113κjξj,t+ε,εJSU(0,σt,ζ,ν)y^{\omega}_t=w'x_{t-1} + \sum_{j=1}^{13} \kappa_j \xi_{j,t} + \varepsilon, \varepsilon \sim \text{JSU}\left(0,\sigma_t,\zeta, \nu\right)

xt=Fxt1+gεtx_t = Fx_{t-1}+g\varepsilon_{t}

σt2=σ^2(1P)+η1εt12+δ1σt12\sigma^2_{t} = \hat\sigma^2 (1 - P)+\eta_1\varepsilon^2_{t-1}+\delta_1\sigma^2_{t-1}

Visuals

The plot of the model components clearly displays the time-varying volatility typically associated with GARCH models:

Code
plot(mod_dynamic)

Conclusion

This vignette demonstrated how to model time-varying volatility using a financial time series example with the tsissm package. We showed how standard innovations state space models can be extended with GARCH dynamics and non-Gaussian distributions such as Johnson’s SU, allowing for more flexible modeling of conditional variance and fat-tailed behavior. This provides a robust framework for capturing some of the complexities of real-world time series, especially when volatility is an essential component of the dynamics.

References

Hyndman, Rob, Anne B Koehler, J Keith Ord, and Ralph D Snyder. 2008. Forecasting with Exponential Smoothing: The State Space Approach.