# Invert SARIMA models to AR

## SARIMA

Seasonal Autoregressive Integrated Moving Average, SARIMA or Seasonal ARIMA, is an extension of ARIMA that explicitly supports univariate time series data with a seasonal component.

It adds three new hyperparameters to specify the autoregression (AR), differencing (I) and moving average (MA) for the seasonal component of the series, as well as an additional parameter for the period of the seasonality. The notation for a SARIMA model is specified as ARIMA$(p,d,q)(P,D,Q) _ {m}$

An ARIMA$(p, d, q)(P, D, Q) _ {m}$ model without intercept can be generally written as \begin{aligned} (1- \phi _ {1} B- \cdots - \phi _ {p} B^{p})(1- \Phi _ {1}B^{m}- \cdots - \Phi _ {P}B^{Pm})(1-B)^{d}(1-B^{m})^D y_{t} = (1 + \theta _ {1}B + \cdots + \theta _ {q}B^{q})(1 + \Theta _ {1}B^{m} + \cdots + \Theta _ {Q}B^{Qm})e _ {t}, \end{aligned}

where $B$ is the backshift operator. R uses the parametrization of this equation.

## Invert SARIMA to ARMA

In this part, we invert SARIAM model to ARMA model.

For autoregression, we have

• $ar = (1- \phi _ {1}B - \cdots - \phi _ {p}B^p)(1-B)^d$,

• $sar = (1- \Phi _ {1}B^m- \cdots - \Phi _ {P}B^{Pm})(1-B^m)^D$ for $m > 1$, $sar = 1$ for $m = 1$.

Then, the AR polynomial of the inverted ARMA model can be written as $AR = ar \cdot sar$.

For moving average, we have

• $ma = (1 + \theta _ {1}B + \cdots + \theta _ {q}B^q)$.

• $sma = (1 + \Theta _ {1}B^m + \cdots + \Theta _ {Q}B^{Qm})$ for $m > 1$, $sma = 1$ for $m = 1$.

Then, the MA polynomial of the inverted ARMA model can be written as $MA = ma \cdot sma$.

By this way, we can get the ARMA$(u, v)$ model inverted from the SARIMA model. It can be written as $$(1- \phi’ _ {1}B- \cdots - \phi’ _ {u}B^u)y_t = (1 + \theta’ _ {1}B + \cdots + \theta’ _ {v}B^{v})e _ {t}.$$

## Three model representations for an ARMA model

In this part, we discuss three model representations for a stationary ARMA$(u, v)$ model.

1. Using the back-shift operator. $$(1- \phi’ _ {1}B - \cdots - \phi’ _ {u}B^u)y _ t = (1 + \theta’ _ {1}B + \cdots + \theta’ _ {v}B^{v})e _ {t}.$$

This representation is compact and useful in parameter estimation. It is also useful in computing recursively multistep-ahead forecasts of $y_t$.

1. AR representation.

2. MA representation.

For the AR and MA representations, we use long division of two polynomials. Given two polynomials $\phi(B)=1-\sum _ {i=1}^{u} \phi _ {i} B^{i}$ and $\theta(B)=1-\sum _ {i=1}^{v} \theta _ {i} B^{i}$, we can obtain, by long division, that $$\frac{\theta(B)}{\phi(B)}=1+\psi _ {1} B+\psi _ {2} B^{2}+\cdots \equiv \psi(B) \tag{1}$$ and $$\frac{\phi(B)}{\theta(B)}=1-\pi _ {1} B-\pi _ {2} B^{2}-\cdots \equiv \pi(B). \tag{2}$$

### AR representation

Using the result of long division in Eq ($2$), the ARMA($u, v$) model can be written as $$\frac{\phi(B)}{\theta(B)}y _ t = e _ t,$$ So, $$y_t = \pi _ {1} y _ {t-1}+\pi _ {2}y _ {t-2}+\pi _ {3}y _ {t-3}+\cdots+e _ t.$$

### MA representation

Again, using the result of long division in Eq (1), an ARMA(u,v) model can also be written as $$y _ t = \frac{\theta(B)}{\phi(B)}e _ t,$$ So, $$y _ t = e _ t + \psi _ {1}e _ {t-1} + \psi _ {2}e _ {t-2} + \cdots .$$

## Function pi_coefficients

In the function pi_coefficients( ) in the R package gratis, the steps of inverting SARIMA model to AR representation are:

#### Step 1. Invert SARIMA model to ARMA model.

ARIMA$(p, d, q)(P, D, Q) _ {m}$ model: \begin{aligned} (1- \phi _ {1}B- \cdots - \phi _ {p}B^p)(1- \Phi _ {1}B^m- \cdots - \Phi _ {P}B^{Pm})(1-B)^d(1-B^m)^Dy_t = (1 + \theta _ {1}B + \cdots + \theta _ {q}B^q)(1 + \Theta _ {1}B^m + \cdots + \Theta _ {Q}B^{Qm})e _ {t}. \end{aligned}

ARMA($u, v$) model: $$(1- \phi’ _ {1}B- \cdots - \phi’ _ {u}B^u)y_t = (1 - \theta’ _ {1}B - \cdots - \theta’ _ {v}B^{v})e _ {t}.$$

#### Step 2. Invert ARMA model to AR model.

Given two polynomials $\phi(B)=1-\sum _ {i=1}^{u} \phi _ {i} B^{i}$ and $\theta(B)=1-\sum _ {i=1}^{v} \theta _ {i} B^{i}$, we can obtain, by long division, that $$\frac{\phi(B)}{\theta(B)}y_t = e_t.$$ Let $(1- \phi’ _ {1}B- \cdots - \phi’ _ {u}B^u) = (1 - \theta’ _ {1}B - \cdots - \theta’ _ {v}B^{v})(1 + \pi _ {1}B + \pi _ {2}B^2 + \cdots)$, then the values $\pi_{i}$ can be obtained recursively:

$\pi _ {1} = - \phi _ {1} + \theta _ {1}, \quad (-\phi _ {1}=\pi _ {1}-\theta _ {1})$

$\pi _ {2} = - \phi _ {2} + (\theta _ {2} + \pi _ {1}\theta _ {1}), \quad (-\phi _ {2}=-\pi _ {1}\theta _ {1}+ \pi _ {2} - \theta _ {2})$

$\pi _ {3} = - \phi _ {3} + (\theta _ {3} + \pi _ {1} \theta _ {2} + \pi _ {2} \theta _ {1})$

$\cdots$

$\pi _ {i} = - \phi _ {i} + (\theta _ {i} + \pi _ {1}\theta _ {i-1} + \pi _ {2}\theta _ {i-2} + \cdots \pi _ {i-1}\theta _ {1})$

$\cdots$

So, we can obtain the AR representation of ARMA($u, v$) model.

#### step 3. Obtain the AR coefficients.

From the recursion process presented above, we get the value $\pi_{i}$, then the inverted AR model can be written as $$(1 + \pi _ {1}B + \pi _ {2}B^2 + \cdots)y_t = e_t.$$ So, the coefficents $\phi_i$ of the inverted AR model is $-\pi_i$.

### Test of pi_coefficents

library(M4comp2018)
library(forecast)
library(polynom)

data(M4)
hourly_M4 <- Filter(function(l) l$period == "Hourly", M4) n_ts <- length(hourly_M4) nlength <- sapply(hourly_M4, function(lentry){ n <- lentry$n
n
})
ts <- hourly_M4[[300]]

forecast_fun <- function(ts, coef.AR){
# ts: a list contains of x, xx, n, h, period
x <- ts$x h <- ts$h
tspx <- tsp(x)
coef <- coef.AR

p <- ifelse("intercept" %in% names(coef), length(coef)-1, length(coef))
c <- ifelse("intercept" %in% names(coef), coef[p+1]*(1-sum(coef[1:p])), 0)

fits <- (c + quantmod::Lag(as.vector(x), k = 1:p) %*% coef[1:p]) %>%
as.vector() %>%
ts(frequency = tspx[3], start = tspx[1])

res <- x - fits

y <- x
for (i in 1:h){
predvalue <- c + tail(quantmod::Lag(as.vector(y), k = 0:(p-1)), 1) %*%
coef[1:p]
y_extend <- c(y, predvalue) %>%
ts(frequency = tspx[3], start = tspx[1])
y <- y_extend
}
pred <- tail(y, h)
return(
list(
x = ts$x, n = ts$n, xx= ts$xx, h = ts$h,
period = ts$period, mean = pred, fitted = fits, residuals = res ) ) } mase.fun <- function(x, xx, pred){ freq <- frequency(x) scaling <- mean(abs(diff(as.vector(x), freq))) mase <- abs(xx - pred) / scaling mase } pi_coefficients <- function(ar = 0, d = 0L, ma = 0, sar = 0, D = 0L, sma = 0, m = 1L, tol = 1e-10) { # non-seasonal AR ar <- polynomial(c(1, -ar)) * polynomial(c(1, -1))^d # seasonal AR if (m > 1) { P <- length(sar) seasonal_poly <- numeric(m * P) seasonal_poly[m * seq(P)] <- sar sar <- polynomial(c(1, -seasonal_poly)) * polynomial(c(1, rep(0, m - 1), -1))^D } else { sar <- 1 } # non-seasonal MA ma <- polynomial(c(1, ma)) # seasonal MA if (m > 1) { Q <- length(sma) seasonal_poly <- numeric(m * Q) seasonal_poly[m * seq(Q)] <- sma sma <- polynomial(c(1, seasonal_poly)) } else { sma <- 1 } n <- 500L theta <- -c(coef(ma * sma))[-1] if (length(theta) == 0L) { theta <- 0 } phi <- -c(coef(ar * sar)[-1], numeric(n)) q <- length(theta) pie <- c(numeric(q), 1, numeric(n)) for (j in seq(n)) pie[j + q + 1L] <- -phi[j] + sum(theta * pie[(q:1L) + j]) pie <- pie[(0L:n) + q + 1L] # Return up to last element greater than tol maxj <- max(which(abs(pie) > tol)) pie <- head(pie, maxj) return(-pie[-1]) } ### 1. ARMA model (include mean) # ARMA p <- 2 q <- 2 fit_arma <- Arima(ts$x, order = c(p, 0, q), seasonal = c(0, 0, 0),
method = "ML")
forec_arma <- forecast(fit_arma, h = ts$h) mase_arma <- mase.fun(x = ts$x, xx = ts$xx, pred = forec_arma$mean)
mean(mase_arma)

# invert ARMA to AR
frequency <- frequency(ts$x) phi <- fit_arma$coef[1:p]
theta <- fit_arma$coef[(p+1):(p+q)] ar <- pi_coefficients(ar = phi, d = 0L, ma = theta, sar = 0, D = 0L, sma = 0, m = frequency) ar <- 'names<-' (c(ar, tail(fit_arma$coef, 1)),
c(paste("ar", sep = "", seq(length(ar))), "intercept"))
forec_ar <- forecast_fun(ts, coef.AR = ar)
mase_ar <- mase.fun(x = ts$x, xx = ts$xx, pred = forec_ar$mean) mean(mase_ar) all.equal(forec_arma$mean, forec_ar$mean) ### 2. ARMA model (no mean) # ARMA p <- 2 q <- 2 fit_arma <- Arima(ts$x, order = c(p, 0, q), seasonal = c(0, 0, 0),
method = "ML", include.constant = FALSE)
forec_arma <- forecast(fit_arma, h = ts$h) mase_arma <- mase.fun(x = ts$x, xx = ts$xx, pred = forec_arma$mean)
mean(mase_arma)

# invert ARMA to AR
frequency <- frequency(ts$x) phi <- fit_arma$coef[1:p]
theta <- fit_arma$coef[(p+1):(p+q)] ar <- pi_coefficients(ar = phi, d = 0L, ma = theta, sar = 0, D = 0L, sma = 0, m = frequency) ar <- 'names<-' (ar, paste("ar", sep = "", seq(length(ar)))) forec_ar <- forecast_fun(ts, coef.AR = ar) mase_ar <- mase.fun(x = ts$x, xx = ts$xx, pred = forec_ar$mean)
mean(mase_ar)
all.equal(forec_arma$mean, forec_ar$mean)

### 3. ARIMA model (d=1, no mean and no drift)
# ARIMA
p <- 2
q <- 1
d <- 1
fit_arima <- Arima(ts$x, order = c(p, d, q), seasonal = c(0, 0, 0), include.drift = FALSE, method = "ML") forec_arima <- forecast(fit_arima, h = ts$h)
mase_arima <- mase.fun(x = ts$x, xx = ts$xx, pred = forec_arima$mean) mean(mase_arima) # invert ARIMA to AR frequency <- frequency(ts$x)
phi <- fit_arima$coef[1:p] theta <- fit_arima$coef[(p+1):(p+q)]
ar <- pi_coefficients(ar = phi, d = d, ma = theta, sar = 0, D = 0L,
sma = 0, m = frequency)
ar <- 'names<-' (ar, paste("ar", sep = "", seq(length(ar))))
forec_ar <- forecast_fun(ts, coef.AR = ar)
mase_ar <- mase.fun(x = ts$x, xx = ts$xx, pred = forec_ar$mean) mean(mase_ar) all.equal(forec_arima$mean, forec_ar$mean) ### 4. SARIMA model (d=D=0, include mean) # SARIMA p <- 2 q <- 1 d <- 0 P <- 1 Q <- 1 D <- 0 fit_sarima <- Arima(ts$x, order = c(p, d, q), seasonal = c(P, D, Q),
method = "ML")
forec_sarima <- forecast(fit_sarima, h = ts$h) mase_sarima <- mase.fun(x = ts$x, xx = ts$xx, pred = forec_sarima$mean)
mean(mase_sarima)

# invert SARIMA to AR
phi <- fit_sarima$coef[1:p] theta <- fit_sarima$coef[(p+1):(p+q)]
Phi <- fit_sarima$coef[(p+q+1):(p+q+P)] Theta <- fit_sarima$coef[(p+q+P+1):(p+q+P+Q)]
ar <- pi_coefficients(ar = phi, d = d, ma = theta, sar = Phi, D = D,
sma = Theta, m = frequency)
ar <- 'names<-' (c(ar, tail(fit_sarima$coef, 1)), c(paste("ar", sep = "", seq(length(ar))), "intercept")) forec_ar <- forecast_fun(ts, coef.AR = ar) mase_ar <- mase.fun(x = ts$x, xx = ts$xx, pred = forec_ar$mean)
mean(mase_ar)
all.equal(forec_sarima$mean, forec_ar$mean)

### 5. SARIMA model (d+D>0, no mean and no drift)
# SARIMA
p <- 2
q <- 1
d <- 0
P <- 1
Q <- 1
D <- 1
fit_sarima <- Arima(ts$x, order = c(p, d, q), seasonal = c(P, D, Q), method = "ML") forec_sarima <- forecast(fit_sarima, h = ts$h)
mase_sarima <- mase.fun(x = ts$x, xx = ts$xx, pred = forec_sarima$mean) mean(mase_sarima) # invert SARIMA to AR phi <- fit_sarima$coef[1:p]
theta <- fit_sarima$coef[(p+1):(p+q)] Phi <- fit_sarima$coef[(p+q+1):(p+q+P)]
Theta <- fit_sarima$coef[(p+q+P+1):(p+q+P+Q)] ar <- pi_coefficients(ar = phi, d = d, ma = theta, sar = Phi, D = D, sma = Theta, m = 24) ar <- 'names<-' (ar, paste("ar", sep = "", seq(length(ar)))) forec_ar <- forecast_fun(ts, coef.AR = ar) mase_ar <- mase.fun(x = ts$x, xx = ts$xx, pred = forec_ar$mean)
mean(mase_ar)
all.equal(forec_sarima$mean, forec_ar$mean)


##### Xiaoqian Wang
###### Research Fellow

My research interests include time series analysis and statistical computing.