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


  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
    $\begingroup$ +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? $\endgroup$
    – forecaster
    May 14, 2015 at 13:22
  • 2
    $\begingroup$ 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. $\endgroup$ May 14, 2015 at 13:40
  • 1
    $\begingroup$ 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. $\endgroup$ Jul 8, 2015 at 14:52
  • 1
    $\begingroup$ @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? $\endgroup$ Aug 2, 2015 at 15:07
  • 1
    $\begingroup$ 5.3 the table contains ARMAX models. The results apply to ARMA models. $\endgroup$ Aug 2, 2015 at 15:16

1 Answer 1


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".


*Thanks to @hejseb for the link.

  • 3
    $\begingroup$ This working paper is very fresh, posted on arXiv just yesterday. $\endgroup$ Aug 1, 2017 at 13:55
  • $\begingroup$ Is there any implementation in python or R? $\endgroup$ Jul 11, 2018 at 9:21
  • $\begingroup$ @DavidMasip, see the updated post for an R implementation. $\endgroup$ Jul 19, 2018 at 19:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.