I want to fit ARIMAX model in R. For simplicity, let's consider model: $Y_t = \theta Y_{t-1} + \beta X + \epsilon$.

I know function auto.arima(), but it fits ARIMA model or a regression with ARIMA errors (when we use xreg argument). This is not the same as ARIMAX model.

My question is what function should I use to fit that model? I could use linear regression to find $\theta$ and $\beta$ by using $lm(Y_t\sim Y_{t-1}+X$), but it fits model: $Y_t = \beta _1 Y_{t-1} + \beta _2 X + \epsilon$. These parametres are from OLS, while the $\theta$ should be from the Yule-Walker equation. Am I right? It's a big problem?

Other way is to fit separately the regression part and autoregressive, by using $lm(Y\sim X)$ and $ar/arima/auto.arima$ for $Y$, but how could I combine these results?


closed as off-topic by Firebug, Michael Chernick, Peter Flom Aug 25 '18 at 13:59

This question appears to be off-topic. The users who voted to close gave this specific reason:

  • "This question appears to be off-topic because EITHER it is not about statistics, machine learning, data analysis, data mining, or data visualization, OR it focuses on programming, debugging, or performing routine operations within a statistical computing platform. If the latter, you could try the support links we maintain." – Firebug, Michael Chernick, Peter Flom
If this question can be reworded to fit the rules in the help center, please edit the question.

  • $\begingroup$ Are you asking how to use lm() to write your own function to fit ARMA models or are you asking whether there is an ARIMAX package for R? $\endgroup$ – Digio Jul 21 '17 at 9:22
  • $\begingroup$ I don't want to fit ARMA/ARIMA model either a regression with ARIMA errors. In this link the best answer says:"4. Fit the model with the arima function in base R. This function can handle ARMAX models through the use of the xreg argument.", but it's not true. Here robjhyndman.com/hyndsight/arimax we can read that "The arima() function in R (and Arima() and auto.arima() from the forecast package) fits a regression with ARIMA errors. " so it's not the same. $\endgroup$ – Neokokaina Jul 21 '17 at 9:41
  • $\begingroup$ As I know, there isn't any auto.arimax function in R and the arimax() funtion dosn't fit the ARIMAX model (only transfer function), so I'm looking for any solutions. As I said, I have 2 options how to do this, but first is in my opinion not statistical correct. Additionally, I don't know how the coefficients in ARIMAX model (Yt=θYt−1+βX+ϵ) are estimated, by OLS? $\endgroup$ – Neokokaina Jul 21 '17 at 9:41
  • $\begingroup$ I'm familiar with the link you posted and I know exactly what you mean when you make a distinction between ARIMAX and regression with ARMA errors. However, since you can't find an ARMAX function, I'm guessing that you want to make your own. What you describe as a second option (fitting two separate models) does not sound as an ARMAX model at all since you're not adding the covariance on the right hand side. And if I understood well, the first option won't work either since fitting a MA() term is not as simple as a linear regression model. IMO, the best to do would be to find an ARMAX function. $\endgroup$ – Digio Jul 21 '17 at 10:18
  • $\begingroup$ I think this post describes a true ARMAX process in R. But why are you rejecting regression with ARMA errors in the first place? $\endgroup$ – Digio Jul 21 '17 at 10:24

auto.arima() is by far a wrap of arima(), hence will act as such in several directions.

In presence of covariates entered through xreg, the arima() Rd states "...If am xreg term is included, a linear regression (with a constant term if include.mean is true and there is no differencing) is fitted with an ARMA model for the error term...". Seems like it fits an LM on the x's, imposing an ARMA structure on the errors. This approach differs from that conveyed by the ARMAX framework, hence auto.arima() won't work.

To fit an ARMAX(p, q) or any sub-class of it, you may want to try with vector generalized linear models (VGLMs) applied to time series, in R, part of my PhD. Particularly, my family function ARXff() estimates ARXs as that one above, $Y_t - \theta Y_{t - 1} = \beta X_t + \varepsilon_t$ (...[1]), by MLE using Fisher scoring. The following gives an example of such, assuming normal errors:

 nn <- 140
 x2 <- rnorm(nn)
 y  <- numeric(nn); y[1] <- 0
 theta <- 0.25
 beta  <- 1.5
 for (ii in 2:nn) 
   y[ii] <- theta * y[ii - 1] + beta * x2[ii] + rnorm(1)

 # Remove warming - up values.
 ts.data <- data.frame(y = y[-c(1:100)], x2 = x2[-c(1:100)])

# Modelling function: vglm(); Family function: ARXff()
 fit1 <- vglm(y ~ x2, ARXff(order = 1, zero = c("coeff", "Var"),
                            type.EIM = "exact"),
         data = ts.data, trace = TRUE)
  VGLM    linear loop  1 :  loglikelihood = -67.168821
  VGLM    linear loop  2 :  loglikelihood = -61.822737
  VGLM    linear loop  3 :  loglikelihood = -60.963229
  VGLM    linear loop  4 :  loglikelihood = -60.943982
  VGLM    linear loop  5 :  loglikelihood = -60.943972
  VGLM    linear loop  6 :  loglikelihood = -60.943972

 Checks on stationarity / invertibility successfully performed. 
 No roots lying inside the unit circle. 
 Further details within the 'summary' output.
 > coef(fit1, matrix = TRUE)
              ARdrift1 loge(noiseVar1) ARcoeff11
(Intercept) -0.086297         0.20932   0.26196
x2           1.339758         0.00000   0.00000

*** This is what arima() returns:

 with(ts.data, arima(y, order =c(1, 0, 0), xreg = x2))

 arima(x = y, order = c(1, 0, 0), xreg = x2)

         ar1  intercept     x2
       0.386     -0.174  1.292
 s.e.  0.170      0.290  0.187

 sigma^2 estimated as 1.3:  log likelihood = -62.02,  aic = 132.05

Given the normality assumption,fit1 is the same as fitting a normal distribution with mean conditional on $x_2$ and $Y_{t - 1}$, similar to [1] above. To see this, use the family uninormal(), from the VGAM package, as follows:

 ts.data <- transform(ts.data, ARcoeff = WN.lags(cbind(y), lags = 1))
 > fit2 <- vglm(y ~ x2 + ARcoeff, uninormal(var.arg = TRUE),
                data = ts.data, trace = TRUE)
 VGLM    linear loop  1 :  loglikelihood = -69.732118
 VGLM    linear loop  2 :  loglikelihood = -62.620483
 VGLM    linear loop  3 :  loglikelihood = -61.012317
 VGLM    linear loop  4 :  loglikelihood = -60.944089
 VGLM    linear loop  5 :  loglikelihood = -60.943972
 VGLM    linear loop  6 :  loglikelihood = -60.943972
> coef(fit2, matrix = TRUE)
                  mean loge(sd)
  (Intercept) -0.086297  0.10466
  x2           1.339758  0.00000
  ARcoeff      0.261963  0.00000

But ARXff() is even broader, e.g., you can model $\theta$ using VGLM-link functions, such as logit(), if needed. Also, I have implemented ARXff() to work with the exact expected information matrices (See type.EIM = "exact") for the ARX model.

However, the covariate effects on the ARXs as above may not be simple to interpret, e.g, as with ordinary LMs. This is the downside of ARX models. For the ease of interpretation, that regression model with ARMA errors arises probably as the most convenient choice. This is how arima() seems to work in presence of covariates. BTW. I have too implemented this approach through my family function ARIMAX.errors.ff().

ARMAXff() and ARIMAXff() have also been implemented accordingly, so you can also fit ARMAX models, similar to fit1. But for regression models with ARMAX errors, ARIMAX.errors.ff() is the choice. These are incorporated in my package, called VGAMextra, an extension of VGAM in a few directions, including time series analysis. At the moment, VGAMextra concetrates on modelling and estimation. I will incorporate, e.g., automatic forecasting over time.


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