Fitting an ARIMAX model with regularization or penalization (e.g. with the lasso, elastic net, or ridge regression) I use the auto.arima() function in the forecast package to fit ARMAX models with a variety of covariates. However, I often have a large number of variables to select from and usually end up with a final model that works with a subset of them.  I don't like ad-hoc techniques for variable selection because I am human and subject to bias, but cross-validating time series is hard, so I haven't found a good way to automatically try different subsets of my available variables, and am stuck tuning my models using my own best judgement.
When I fit glm models, I can use the elastic net or the lasso for regularization and variable selection, via the glmnet package. Is there a existing toolkit in R for using the elastic net on ARMAX models, or am I going to have to roll my own? Is this even a good idea?
edit: Would it make sense to manually calculate the AR and MA terms (say up to AR5 and MA5) and the use glmnet to fit the model?
edit 2: It seems that the FitAR package gets me part, but not all, of the way there.
 A: This is not a solution but some reflections on the possibilities and difficulties that I know of. 
Whenever it is possible to specify a time-series model as 
$$Y_{t+1} =  \mathbf{x}_t \beta + \epsilon_{t+1}$$
with $\mathbf{x}_t$ computable from covariates and time-lagged observations, it is also possible to compute the least-squares elastic net penalized estimator of $\beta$ using glmnet in R. It requires that you write code to compute $\mathbf{x}_t$ to form the model matrix that is to be specified in glmnet. This works for AR-models but not directly for ARMA-models, say. Moreover, the cross-validation procedures of glmnet are not sensible per se for time-series data. 
For more general models 
$$Y_{t+1} =  f(\mathbf{x}_t, \beta) + \epsilon_{t+1}$$
an implementation of an algorithm for computing the non-linear least-squares elastic net penalized estimator of $\beta$ is needed. To the best of my knowledge there is no such implementation in R. I am currently writing an implementation to solve the case where
$$Y_{t+1} =  \mathbf{x}_t g(\beta) + \epsilon_{t+1}$$ 
the point being that it is paramount for model selection that the lasso penalization is on $\beta$ and not $g(\beta)$. If I recall the ARIMA-parametrization correctly it also takes this form $-$ but I cannot offer any code at the moment. It is (will be) based on A coordinate gradient descent method for nonsmooth separable minimization.
Another issue is the selection of the amount of penalization (the tuning parameters). It will generally require a form of cross-validation for time-series, but I hope to be able to work out some less computationally demanding methods for specific models.  
A: I was challenged by a client to solve this problem in an automatic i.e. turnkey way. I implemented an approach that for each pair ( i.e. y and a candidate x ) , prewhiten , compute cross-correlations of the pre-whitened series, identify the PDL ( OR ADL AUTOREGRESSIVE DISTRIBUTED LAG MODEL including any DEAD TIME ) while incorporating Intervention Detection to yield robust estimates, develop a "measure" for this structure. After conducting this for ALL candidate regressors, rank them by the "measure" and then select the top K regressors based upon the "measure". This is sometimes referred to as Linear Filtering. We successfully incorporated this heuristic into our commercially available time series package. You should be able to "ROLL YOUR OWN" if you have sufficient time and statistical programming skills and/or available modules to implement some of these tasks. 
