ARIMAX: only MA(1) model gives meaningful coefficients for exogenous variables

Question

Is it possible that the coefficients for the exogenous variables in (regression with ARIMA errors) are only meaningful with MA(1) models?

Here is my reasoning and a simulation. I would be grateful for any guidance. I have an engineering background with only undergrad-level stats.

Consider a dynamic process for a dependent variable y that responds to two independent variables x1, x2.

There is a steady state model that describes the value of y reached after a period of time. ys = 30 + 100 x1 + 200 x2.

There is a rate constant k that controls the approach to steady state. dy/dt = k(ys - y).

This describes any first-order physical process, for example,

y = the temperature of a bagel in a toaster
x1 = pressing down the lever to begin toasting
x2 = pressing the "bagel" button to further elevate the coil temperatures

A simple version of the heat equation is dT/dt = k(T0 - T) where T0 is the temperature of the coil and T the temperature of the bagel. Let's make a simple model of the coil temperature assuming it responds instantaneously to the inputs T0 = b0 + b1 x1 + b2 x2.

I simulated my data with room temperature as 30 degrees, pressing the lever heats to 130, and then pressing the button heats to 330 degrees. The bagel will gradually heat - this is shown by the solution for y in the below plot.

If x1 and x2 have values of zero or one for pressing the lever and the button, then b0 measures room temperature (30 degrees), b1 measures the additional degrees that pressing the lever heats the coil above room temperature (+100 degrees), and b2 the additional degrees that pressing the button gives (+200 degrees). The value of k has units of 1/sec and describes empirically the rate of heat transfer from the coil to the bagel.

I imagine this simple ODE system could model other dynamic processes approximately. I'm not a social scientist, but I imagine there are some dynamics after the application of an intervention and an eventual steady state.

Let's say I have observations of the bagel temperature and I want to measure the coefficients for x1 and x2 on dynamic data. It's natural to fit a regression model with ARIMA errors for y with exogenous regressors x1, x2.

In the equation dy/dt = k(ys - y), it seems to me that the driving force for the approach to steady state ys - y is related to the MA(1), x1, and x2 terms. It seems that there is no role for the AR terms or higher-order MA terms. And possibly, including these terms will destroy the relationship we would like to fit.

I simulated some data by integrating the above ODE and fit two regression models. I fit one model with auto.arima for automatic model selection. The coefficients for x1 and x2 are wrong and close to zero. b1 = 0.1, b2 = 0.2.

I fit a second model with only MA(1). The coefficients are much closer to the model that generated the data. b1 = 68, b2 = 181. I suppose there is some error created by the approximation of the ODE integration by the basic arithmetic of the MA(1) term.

edit: I wonder if the MA(1) model with exogenous regressors is equivalent to integrating an ODE with the Euler method, and we can expect some error in the estimate in a similar way that Euler method creates integration error that is addressed by Runge-Kutta.

library(tidyverse)
library(deSolve)
library(forecast)
#> Registered S3 method overwritten by 'quantmod':
#>   method            from
#>   as.zoo.data.frame zoo


# functions for generating dynamic data ----
  
  
fn <- function(x1, x2) {
  # predict steady state response from predictors
  30 + 100 * x1 + 200 * x2
}


calc_dy_dt <- function(t, y, p) {
  # return the derivative
  # first-order approach to steady state
  
  x1 <- approx(p$data$t, p$data$x1, t, rule = 2)$y
  x2 <- approx(p$data$t, p$data$x2, t, rule = 2)$y
  
  dy_dt <- p$k * (fn(x1, x2) - y)
  list(dy_dt)
}


integrate_y <- function(data, k) {
  # integrate the ode to produce dynamic response
  # and add as new column to the data frame
  
  soln <- ode(
    data$ys[1],
data$t,
    calc_dy_dt,
    lst(data, k)
  )
  
  data |> 
    add_column(y = soln[, 2])
}


plot_solution <- function(data) {
  # plot predictor, steady-state response, and dynamic response
  
  data |> 
    pivot_longer(-t) |> 
    mutate(
      label = case_match(name, "ys" ~ "steady state solution", "y" ~ "dynamic solution", .default = "input"),
      name = factor(name, c("x1", "x2", "ys", "y")),
      label = factor(label, c("input", "steady state solution", "dynamic solution"))
    ) |> 
    ggplot() + aes(t, value, color = label) + geom_line() +
    facet_wrap(~name, scales = "free_y")
}


# sample data ----

dat <- tibble(
  t = 1:100,
  x1 = if_else(t < 20, 0, 1),
  x2 = if_else(t < 50, 0, 1),
  ys = fn(x1, x2)
)

dat <- dat |> 
  integrate_y(0.1)


# plot ----

dat |> 
  plot_solution() |> 
  print()

# models ----

model <- auto.arima(dat$y, xreg = dat |> select(x1, x2) |> as.matrix())

model |> 
  print()
#> Series: dat$y 
#> Regression with ARIMA(5,1,0) errors 
#> 
#> Coefficients:
#>          ar1      ar2     ar3      ar4     ar5      x1      x2
#>       1.6899  -1.3043  0.9268  -0.5569  0.1946  0.1234  0.2463
#> s.e.  0.0979   0.1882  0.2084   0.1856  0.0959  0.4495  0.4497
#> 
#> sigma^2 = 1.47:  log likelihood = -157.72
#> AIC=331.44   AICc=333.04   BIC=352.2

model_ma <- Arima(dat$y, order = c(0, 0, 1), xreg = dat |> select(x1, x2) |> as.matrix())

model_ma |> 
  print()
#> Series: dat$y 
#> Regression with ARIMA(0,0,1) errors 
#> 
#> Coefficients:
#>          ma1  intercept       x1        x2
#>       0.7122    35.0462  68.6714  181.5979
#> s.e.  0.0546    10.0105  12.3623   10.1697
#> 
#> sigma^2 = 730.1:  log likelihood = -469.86
#> AIC=949.73   AICc=950.37   BIC=962.75

^{Created on 2024-06-03 with reprex v2.1.0}