ARMA-GARCH Estimation

The Model

The observed series $y_t$ is decomposed into a conditional mean $\mu_t$ and an innovation $\varepsilon_t$:

\[ y_t = \mu_t + \varepsilon_t \]

By default (arma = c(0,0)) $\mu_t$ reduces to either a constant $\mu$ or zero, controlled independently by the constant argument to garch_modelspec. Setting arma = c(p,q) instead activates a jointly estimated ARMA($p,q$) mean equation, in which $\mu_t$ is time-varying and follows the mean-adjusted recursion

\[ \mu_t = \mu + \sum_{i=1}^{p}\phi_i\left(y_{t-i}-\mu\right) + \sum_{j=1}^{q}\theta_j\varepsilon_{t-j} \]

so that $\mu$ retains its interpretation as the unconditional mean of $y$, and the AR feedback term is centered on $y_{t-i}-\mu$ rather than the raw $y_{t-i}$. This is the same mean-adjusted form used in rugarch’s ARFIMA-X mean equation.

Combining the two equations above, the innovation $\varepsilon_t$ is obtained recursively as

\[ \varepsilon_t = y_t - \mu_t = \left(y_t-\mu\right) - \sum_{i=1}^{p}\phi_i\left(y_{t-i}-\mu\right) - \sum_{j=1}^{q}\theta_j\varepsilon_{t-j} \]

which only requires $\varepsilon_{t-1},\ldots,\varepsilon_{t-q}$ from previous time steps, so the whole recursion can be evaluated forward in a single pass once $\mu$, $\phi_i$ and $\theta_j$ are known. Setting arma = c(0,0) collapses both equations exactly to the constant-mean case, $\varepsilon_t = y_t - \mu$, matching prior (pre 1.0.5) behavior with no change in results.

It is $\varepsilon_t$ (rather than $y_t-\mu$) that is passed on to the conditional variance equation of whichever GARCH flavor is chosen (see the GARCH Models vignette for the flavor-specific variance equations), decomposed as

\[ \varepsilon_t = \sigma_t z_t, \qquad z_t \overset{\text{iid}}{\sim} D(0,1;\ldots) \]

with $\sigma_t$ the conditional standard deviation implied by the chosen GARCH variance equation, and $z_t$ an iid draw from the chosen conditional distribution $D$ (normal, Student-t, skew distributions, etc.), with any shape/skew parameters estimated jointly with everything else.

The combined pre-sample/burn-in length used to initialize this recursion is $\max\left(\text{GARCH order}, \text{ARMA order}\right)$, so that e.g. an ARMA(2,2) mean equation together with a GARCH(1,1) variance equation initializes correctly; over that pre-sample the process is assumed to start at its unconditional mean with zero shocks, i.e. $y_t = \mu$ and $\varepsilon_t = 0$.

Regressors in the Mean Equation (ARMAX)

Regressors can also enter the conditional mean via the xreg argument to garch_modelspec, with coefficients named tau1, tau2, … in the parmatrix ($\tau_k$; the symbol $\xi$ is already used in this package for the variance equation regressors). Two conventions are available through xreg_type. The default, "arma_errors", runs the ARMA recursion on the regressor-adjusted deviations

\[ y_t - \mu - \sum_{k=1}^{s}\tau_k x_{k,t} = \sum_{i=1}^{p}\phi_i\left(y_{t-i} - \mu - \sum_{k=1}^{s}\tau_k x_{k,t-i}\right) + \varepsilon_t + \sum_{j=1}^{q}\theta_j\varepsilon_{t-j} \]

matching stats::arima(xreg = ) and the sibling tsarma package, so that $\tau_k$ is the long-run marginal effect of $x_{k,t}$. The alternative, "armax" (the rugarch convention), adds the regressor contribution to the conditional mean at time $t$ only,

\[ y_t = \mu + \sum_{i=1}^{p}\phi_i\left(y_{t-i}-\mu\right) + \sum_{k=1}^{s}\tau_k x_{k,t} + \varepsilon_t + \sum_{j=1}^{q}\theta_j\varepsilon_{t-j} \]

so that $\tau_k$ is the impact effect and the long-run effect of a sustained change is $\tau_k/\left(1-\sum_i\phi_i\right)$. The two conventions differ only by the lagged-regressor term inside the AR feedback and are algebraically identical whenever $p = 0$.

xreg must be an xts matrix aligned with y: same number of rows, matching time index, no NA/NaN/Inf values, and full column rank (including against the constant when constant = TRUE) - violations are rejected at specification time. Because the recursion consumes a regressor value at every step, predict() and simulate() require future values (newxreg with h rows, xreg with h + burn rows respectively), and tsfilter() requires newxreg covering the appended observations; if the model was specified with xreg but no future values are supplied, a zero matrix is substituted with a warning rather than an error.

A short demo with a day-of-week effect in the Nikkei returns:

Code

library(tsgarch)
#> Loading required package: tsmethods
suppressMessages(library(xts))
#> Warning: package 'xts' was built under R version 4.6.1
#> Warning: package 'zoo' was built under R version 4.6.1
data(nikkei)
nikkei <- xts(nikkei$value, as.Date(nikkei$index))
monday <- xts(as.numeric(weekdays(index(nikkei)) == "Monday"), index(nikkei))
colnames(monday) <- "monday"
spec_x <- garch_modelspec(nikkei, arma = c(1,1), xreg = monday,
                          xreg_type = "arma_errors", model = 'garch',
                          constant = TRUE, init = 'unconditional',
                          distribution = 'jsu')
mod_x <- estimate(spec_x)
coef(mod_x)["tau1"]
#>       tau1 
#> -0.0384228

Forecasting requires the future regressor values:

Code

newx <- xts(matrix(as.numeric(weekdays(index(nikkei)[NROW(nikkei)] + 1:10) == "Monday"),
                   ncol = 1), index(nikkei)[NROW(nikkei)] + 1:10)
colnames(newx) <- "monday"
predict(mod_x, h = 10, newxreg = newx)$mean
#>                    [,1]
#> 2000-12-22 -0.002029333
#> 2000-12-23  0.112332341
#> 2000-12-24  0.028983951
#> 2000-12-25  0.051306625
#> 2000-12-26  0.045457273
#> 2000-12-27  0.077723438
#> 2000-12-28  0.054207405
#> 2000-12-29  0.071346219
#> 2000-12-30  0.058855211
#> 2000-12-31  0.067958832

Which convention should I use?

Written out in levels for an AR(1) mean equation, arma_errors is

\[ y_t = \mu\left(1-\phi_1\right) + \phi_1 y_{t-1} + \tau_1 x_{1,t} - \phi_1\tau_1 x_{1,t-1} + \varepsilon_t \]

and armax is the same expression without the final regressor lag. Both are restrictions of the autoregressive distributed lag model $y_t = c + \phi_1 y_{t-1} + \beta_0 x_{1,t} + \beta_1 x_{1,t-1} + \varepsilon_t$: arma_errors imposes $\beta_1 = -\phi_1\beta_0$ and armax imposes $\beta_1 = 0$. Neither is nested in the other.

The difference that matters in practice is propagation. Under arma_errors a one period blip in a regressor shifts $y_t$ by $\tau_k$ and nothing afterwards, because the ARMA dynamics apply only to the error; $\tau_k$ is both the impact and the total effect. Under armax the same blip is carried forward by the autoregressive term as $\tau_k, \phi_1\tau_k, \phi_1^2\tau_k,\ldots$, so $\tau_k$ is the impact effect and $\tau_k/\left(1-\sum_i\phi_i\right)$ the cumulative one. Event and calendar dummies usually suit the former, covariates whose influence builds and decays with the series the latter. With no AR terms the two are identical, so the choice only matters when the mean equation has autoregressive dynamics and the regressor is itself persistent.

Since armax with both $x_{k,t}$ and $x_{k,t-1}$ as separate xreg columns spans the unrestricted model, the restriction each convention imposes can be tested directly by fitting that larger model and inspecting $\hat\beta_1$. Simpler still, the likelihood discriminates sharply between the two: in a Monte Carlo experiment reported in the GARCH Models vignette the correctly specified convention had the higher likelihood in all 1000 comparisons. That is worth doing rather than guessing, because the parameter most reliably damaged is the autoregressive one: the bias in $\phi_1$ is roughly $\pm 0.36$ on a true value of 0.6 in three of the four misspecified cells considered there. How the error divides between $\tau$ and $\phi$ turns on how persistent the regressor is, so a plausible looking regressor coefficient is no evidence that the convention is right.

Stationarity and Invertibility

For $\phi=\left(\phi_1,\ldots,\phi_p\right)$ and $\theta=\left(\theta_1,\ldots,\theta_q\right)$ to be well defined ARMA coefficients, the AR polynomial

\[ \Phi(z) = 1 - \phi_1 z - \cdots - \phi_p z^p \]

must have all roots strictly outside the unit circle (stationarity), and the MA polynomial

\[ \Theta(z) = 1 + \theta_1 z + \cdots + \theta_q z^q \]

must have all roots strictly outside the unit circle (invertibility). Equivalently, the inverse roots of both polynomials must lie strictly inside the unit circle - this is exactly what the first panel of plot(object, type = "arma") displays, and what arma_inverse_roots() computes programmatically.

Rather than impose these two conditions as nonlinear inequality constraints on $\phi_i$ and $\theta_j$ directly during optimization (which would require either finite-difference Jacobians, or an autodiff-aware eigendecomposition of a companion matrix), tsgarch instead reparameterizes the raw quantities being optimized over as partial-autocorrelation-like parameters $r^{ar}_i, r^{ma}_j \in \left(-1,1\right)$, and recovers $\phi$ and $\theta$ from them via the Durbin-Levinson / Jones recursion, which is guaranteed by construction to always land in the stationarity/invertibility region (Barndorff-Nielsen and Schou, 1973; Jones, 1980) - the same device used internally by stats::arima(..., transform.pars = TRUE).

Concretely, given raw parameters $r_1,\ldots,r_p\in\left(-1,1\right)$, the AR coefficients are obtained by the forward recursion

\[ \phi^{(1)}_1 = r_1, \qquad \phi^{(k)}_k = r_k, \qquad \phi^{(k)}_i = \phi^{(k-1)}_i - r_k\,\phi^{(k-1)}_{k-i} \quad \left(i=1,\ldots,k-1\right) \]

for $k=2,\ldots,p$, with the final AR coefficients given by $\phi_i = \phi^{(p)}_i$. The same recursion applied to the negated raw MA parameters, with the result negated again, yields MA coefficients $\theta_j$ that satisfy invertibility, via the well-known AR/MA duality:

\[ \theta = -\,\text{DL}\left(-r^{ma}_1,\ldots,-r^{ma}_q\right) \]

where $\text{DL}(\cdot)$ denotes the recursion above. Because this map is a smooth (indeed, polynomial) function of the raw parameters, it differentiates exactly under automatic differentiation, so TMB provides exact (not finite-difference) Jacobians of the whole ARMA mean equation throughout estimation, prediction, simulation and filtering. Only simple box constraints $r^{ar}_i, r^{ma}_j \in \left(-1,1\right)$ (tightened slightly to $\left(-0.995, 0.995\right)$ for numerical safety near the boundary) are needed during optimization - no nonlinear stationarity/invertibility constraint is ever evaluated or differentiated.

The raw parameters themselves (arpacf/mapacf in the model’s parmatrix) are not directly interpretable; the transformed $\phi$/$\theta$ coefficients are what summary() displays, and are also available via arma_coefficients().

Diagnostic helpers and the “ARMA(1,1) trap”

Four helpers are exported for inspecting the estimated mean equation:

arma_coefficients() returns the transformed AR and MA coefficients $\phi_i$ and $\theta_j$.
arma_inverse_roots() returns the inverse roots of $\Phi(z)$ and $\Theta(z)$; all of them must lie inside the unit circle.
arma_irf() computes the impulse-response function of the mean equation (how a one-unit shock to $\varepsilon_t$ propagates through $\mu_t$).
arma_near_cancellation() flags pairs of AR and MA inverse roots that are closer than a chosen tolerance in the complex plane.

The first three are mostly self-explanatory, but the last one is worth a longer warning because ARMA($1,1$) is so often used as a default without checking whether it is actually identified by the data.

Why near-cancellation matters

For an ARMA($1,1$) with our sign convention,

\[ \Phi(z) = 1 - \phi_1 z, \qquad \Theta(z) = 1 + \theta_1 z. \]

If the two polynomials share a common root then the AR and MA operators cancel. After cancellation the mean equation is just white noise around $\mu$; the model has two parameters doing the job of zero and the likelihood is flat along a ridge. Even when the match is not exact, near-cancellation makes the model practically unidentifiable: the data cannot tell whether the observed persistence comes from the AR side, the MA side, or a mixture of both.

The fitted coefficients below illustrate exactly this. The AR coefficient is $\hat\phi_1 \approx -0.99$ and the MA coefficient is $\hat\theta_1 \approx 0.99$. Because $\hat\phi_1 \approx -\hat\theta_1$, the two polynomials $1 - \hat\phi_1 z$ and $1 + \hat\theta_1 z$ are almost identical, so their roots are almost the same. arma_near_cancellation() reports the distance between the corresponding inverse roots.

Demo

Code

library(tsgarch)
suppressMessages(library(data.table))
#> Warning: package 'data.table' was built under R version 4.6.1
suppressMessages(library(xts))
data(nikkei)
nikkei <- xts(nikkei$value, as.Date(nikkei$index))
spec <- garch_modelspec(nikkei, arma = c(1,1), model = 'garch', constant = TRUE, 
                        init = 'unconditional', distribution = 'jsu')
mod <- estimate(spec)
as_flextable(summary(mod))

GARCH Model Summary
	Estimate	Std. Error	t value	Pr(>\|t\|)
$\mu$	0.0560	0.0145	3.8737	0.0001	***
$\phi_1$	-0.9917	0.0082	-120.3148	0.0000	***
$\theta_1$	0.9943	0.0064	155.3514	0.0000	***
$\omega$	0.0184	0.0045	4.1035	0.0000	***
$\alpha_1$	0.1169	0.0134	8.7040	0.0000	***
$\beta_1$	0.8802	0.0125	70.3128	0.0000	***
$\zeta$	-0.1544	0.0609	-2.5338	0.0113	*
$\nu$	1.7501	0.0853	20.5168	0.0000	***
$P$	0.9971	0.0055	180.4444	0.0000	***
Signif. codes: 0 '*' 0.001 '' 0.01 '*' 0.05 '.' 0.1 ' ' 1
variance targeting: FALSE
initialization value: 1.817
LogLik: -6426.112
AIC: 1.287e+04 \| BIC: 1.293e+04
Model Equation
$\mu_t = \mu + \phi_1\left(y_{t-1} - \mu\right) + \theta_1\varepsilon_{t-1}$
$\varepsilon_t \sim JSU\left(0,\sigma_t,\zeta, \nu\right)$
$\sigma^2_t = \omega + \alpha_1\varepsilon^2_{t-1} + \beta_1\sigma^2_{t-1}$

Persistence (P) and Unconditional Variance Equations
$P = \sum_{j=1}^q \alpha_j + \sum_{j=1}^p \beta_j$
$E\left[\varepsilon^2_t\right] = \frac{\omega}{1 - P}$

The estimated AR and MA coefficients and the inverse roots:

Code

arma_coefficients(mod)
#> $ar
#>        ar1 
#> -0.9916538 
#> 
#> $ma
#>      ma1 
#> 0.994338
arma_inverse_roots(mod)
#> $ar
#> [1] -0.9916538+0i
#> 
#> $ma
#> [1] -0.994338+0i
arma_near_cancellation(arma_inverse_roots(mod), tol = 0.1)
#>    ar_index ma_index       ar_root      ma_root    distance
#>       <num>    <num>        <cplx>       <cplx>       <num>
#> 1:        1        1 -0.9916538+0i -0.994338+0i 0.002684213

The plot method now takes in an argument type which dispatches to either the GARCH (default and backwards compatible) or ARMA plots.

Code

plot(mod)

Code

plot(mod, type = "arma", which = NULL, envelope = "parametric")

The first panel of the ARMA plot is the inverse-root diagram. Because the AR and MA inverse roots are extremely close, the plot joins them with a dashed segment and the legend flags “Possible common factors”. This is the visual cue that an ARMA($1,1$) mean equation may be over-parameterized for this series.

ARMA-GARCH Estimation

Alexios Galanos

2026-10-09

The Model

Regressors in the Mean Equation (ARMAX)

Which convention should I use?

Stationarity and Invertibility

Diagnostic helpers and the “ARMA(1,1) trap”

Why near-cancellation matters

Demo