This vignette is intended to show the current ability to parse objects generated by forecast::Arima() by {equatiomatic}, thanks to contributions from Jay Sumners. The output uses notation from Hyndman with the exception of \(a\) and \(b\) for intercept and drift, respectively, being replace by \(\mu\) and \(delta\). To improve readability and reduce formula size, Arima functions are presented in Backshift (sometimes called lag) notation. forecast::Arima() will automatically generate either an Arima model or a linear regression model with Arima errors. We’ll address both in the examples below.
First we need to generate a model. In this example, we won’t worry too much about whether or not the model is appropriate. Rather, we aim to illustrate that the corresponding equations are parsed correctly.
library(equatiomatic)
library(forecast)
# Build Arima (no regression)
simple_ts_mod <- Arima(simple_ts,
order = c(1, 1, 1),
seasonal = c(1, 0, 1),
include.constant = TRUE
)To extract the equation, we just call equatiomatic::extract_eq() on the resulting model object, equivalent to other model types. The model above outputs the following.
extract_eq(simple_ts_mod)\[ (1 -\phi_{1}\operatorname{B} )\ (1 -\Phi_{1}\operatorname{B}^{\operatorname{4}} )\ (1 - \operatorname{B}) (y_{t} -\delta\operatorname{t}) = (1 +\theta_{1}\operatorname{B} )\ (1 +\Theta_{1}\operatorname{B}^{\operatorname{4}} )\ \varepsilon_{t} \]
Next, we’ll illustrate a slightly more complicated example that includes a linear regression. We’ll use the ts_reg_list object that is exported from {equatiomatic} to build up our data.
# Build Exogenous Regressors
xregs <- as.matrix(data.frame(
x1 = ts_reg_list$x1 + 5,
x2 = ts_reg_list$x2 * 5
))
# Build Regression Model
ts_reg_mod <- Arima(ts(ts_reg_list$ts_rnorm, freq = 4),
order = c(1, 1, 1),
seasonal = c(1, 0, 1),
xreg = xregs,
include.constant = TRUE
)Despite the extra complexity of the model, the code to pull the equation remains equivalent.
extract_eq(ts_reg_mod)\[ \begin{alignat}{2} &\text{let}\quad &&y_{t} = \operatorname{y}_{\operatorname{0}} +\delta\operatorname{t} +\beta_{1}\operatorname{x1}_{\operatorname{t}} +\beta_{2}\operatorname{x2}_{\operatorname{t}} +\eta_{t} \\ &\text{where}\quad &&(1 -\phi_{1}\operatorname{B} )\ (1 -\Phi_{1}\operatorname{B}^{\operatorname{4}} )\ (1 - \operatorname{B}) \eta_{t} \\ & &&= (1 +\theta_{1}\operatorname{B} )\ (1 +\Theta_{1}\operatorname{B}^{\operatorname{4}} )\ \varepsilon_{t} \\ &\text{where}\quad &&\varepsilon_{t} \sim{WN(0, \sigma^{2})} \end{alignat} \]
Although currently still in development, the other functionalities of extract_eq() for lm objects should generally work for forecast::Arima objects as well (e.g. interaction terms). Please let us know (preferably with a reproducible example) if you run into any issues.
Development items on the docket include:
Testing additional functionality from extract_eq() for lm objects.
Implementation of arguments that will allow the user to rename variables and/or specify a set of greek letters.
If there is interest, the ability to switch to an expanded version of the equation. However, this could potentially generate extremely long and difficult to read equations.
This is the first of (hopefully) most of the models in forecast being implemented in {equatiomatic}. Arima is one of the most used modeling techniques for time-series and, as such, got first treatment. There’s a lot we can do from here and I’m happy for any help or feedback from the community.