# Regularization for ARIMA models

I am aware of LASSO, ridge and elastic-net type of regularization in linear regression models.

Question:

1. Can this (or a similar) kind of penalized estimation be applied to ARIMA modelling (with a non-empty MA part)?

In building ARIMA models, it seems usual to consider a pre-selected maximum lag order ($p_{max}$,$q_{max}$) and then choose some optimal order $p \leqslant p_{max}$ and $q \leqslant q_{max}$ e.g. by minimizing AIC or AICc. But could regularization be used instead?

My further questions are:

1. Could we include all terms up to ($p_{max}$,$q_{max}$) but penalize the size of the coefficients (potentially all the way to zero)? Would that make sense?
2. If it would, has that been implemented in R or other software? If not, what was the trouble?

A somewhat related post can be found here.

• +1 for a very good question. Since P,Q are discrete values it may be more efficient to do a grid search to find optimal order of P,Q? – forecaster May 14 '15 at 13:22
• I am glad you liked it! Yes, a grid search is one of the options in the framework that I refer to as "the usual one". There one may search over a grid of possible combinations of $(p,q)$ from $(0,0)$ to $(p_{max},q_{max})$. However, this is still part of the "usual framework". As an alternative, I am interested in keeping all the lags but penalizing the size of the coefficients. – Richard Hardy May 14 '15 at 13:40
• columbia.edu/~sn2294/papers/forecast.pdf Supposedly LASSO works better as you can skip some lags instead of putting a maximum. Same can be done by AIC but then it gets computationally expensive. – Cagdas Ozgenc Jul 8 '15 at 14:52
• @CagdasOzgenc, I skimmed through the paper but it does not seem to be dealing with regularization applied on ARIMA models (although it mentions ARMA models in the context of information criteria). Could you please point out which part of the paper is relevant for my questions? – Richard Hardy Aug 2 '15 at 15:07
• 5.3 the table contains ARMAX models. The results apply to ARMA models. – Cagdas Ozgenc Aug 2 '15 at 15:16

# Answering Question 1.

Chen & Chan "Subset ARMA selection via the adaptive Lasso" (2011)* use a workaround to avoid the computationally demanding maximum likelihood estimation. Citing the paper, they

propose to find an optimal subset ARMA model by fitting an adaptive Lasso regression of the time series $y_t$ on its own lags and those of the residuals that are obtained from fitting a long autoregression to the $y_t$s. <...> [U]nder mild regularity conditions, the proposed method achieves the oracle properties, namely, it identifies the correct subset ARMA model with probability tending to one as the sample size increases to infinity, and <...> the estimators of the nonzero coefficients are asymptotically normal with the limiting distribution the same as that when the zero coefficients are known a priori.

Optionally, they suggest maximum likelihood estimation and model diagnostics for the selected subset ARMA model(s).

Wilms et al. "Sparse Identification and Estimation of High-Dimensional Vector AutoRegressive Moving Averages" (2017) do even more than I asked for. Instead of a univariate ARIMA model, they take a vector ARMA (VARMA) in high dimensions, and they use an $L_1$ penalty for estimation and lag order selection. They present the estimation algorithm and develop some asymptotic results.

In particular, they employ a two-stage procedure. Consider a VARMA model $$y_t = \sum_{l=1}^p \Phi_l y_{t-l} + \sum_{m=1}^q \Theta_m \varepsilon_{t-m} + \varepsilon_t$$ which needs to be estimated, but the lag orders $p$ and $q$ are uknown.

• In Stage 1, they approximate the VARMA model by a high-order VAR model and estimate it using a Hierarchical VAR estimator which places a lag-based hierarchical group-lasso penalty on the autoregressive parameters.
(The lag order is set to be $\lfloor 1.5\sqrt{T} \rfloor$. The model equations are estimated jointly and the Frobenius norm of the errors $||y-\hat y||_2^F$ is minimized with a hierarchical group-lasso penalty on the regression coefficients.)
They obtain residuals $\hat\varepsilon := y - \hat y$ to be used as proxies for the true errors in Stage 2.

• In Stage 2, they estimate a VARX model where X represents lagged residuals from Stage 1. That is, they minic a VARMA model but use estimated residuals in place of true errors, $$y_t = \sum_{l=1}^{\hat p} \Phi_l y_{t-l} + \sum_{m=1}^{\hat q} \Theta_m \hat\varepsilon_{t-m} + u_t,$$ which allows applying the same estimator (hierarchical group-lasso) again just like in Stage 1.
($\hat p$ and $\hat q$ are set to be $\lfloor 1.5\sqrt{T} \rfloor$.)

The approach of Wilms et al. is implemented in the R package "bigtime".

References

*Thanks to @hejseb for the link.

• This working paper is very fresh, posted on arXiv just yesterday. – Richard Hardy Aug 1 '17 at 13:55
• Is there any implementation in python or R? – David Masip Jul 11 '18 at 9:21
• @DavidMasip, see the updated post for an R implementation. – Richard Hardy Jul 19 '18 at 19:08