Unusual use of time series: getting xreg from a given forecast value

Question

I recently fit an ARIMA model for some daily sales data. To account for seasonality I used various dummies in xreg for different days in the month, days in the week, holidays, etc. Since the sales would be related with the number of employees that day, one of the variables that went into model (in xreg) is number of employees that day.

Recall in simple linear regression, given an equation $y =$ $b_0 + b_1x $, we can calculate estimate of $y$ given $x$, and also estimate a value for $x$ that would give a prespecified $y$.

My question then would be, can we do this in time series in R? Basically what I want to do is the following:

1. Given the number of employees tomorrow, what would be the forecast for the sales? (usual scenario easily done in the forecast package)

2. Given I want tomorrow's sales to be $y$, what would be an estimate of the number of employees necessary?

To tackle this, I read the description of the package forecast and know that when xreg in not NULL, it is fitting an regression model with ARMA error. Does that mean I can use the regression coefficients, ARMA coefficients, and the ideal sales to get 2 above just like in a simple linear regression?

Answer 1

For point 1, you simply feed tomorrow's number of employees (along with all other covariates) into the xreg parameter of forecast.Arima (note the capital A). Here are a couple of SO questions on xreg.

Point 2 is a little more interesting. Recall that auto.arima() fits a regression with ARIMA errors, not an ARIMAX model (matters would be far more complicated for an ARIMAX model):

$$ y_t =(X\beta)_t + \epsilon_t, $$

where $\epsilon_t$ follows some ARIMA process.

Let's work with the AirPassengers dataset and use two random covariates:

library(forecast)
set.seed(1)
predictors <- data.frame(foo=rnorm(length(AirPassengers)),
  bar=rnorm(length(AirPassengers)))
model <- auto.arima(AirPassengers,xreg=predictors)

Now, say you are interested in what value of foo yields a forecast of 400.

target <- 400

You will of course need to specify what values the other covariates have. Let's say that $x_{\text{bar}}=1$. So you need to solve the following equation for $x_\text{foo}$:

$$ 400 = \hat{\beta}_0 + \hat{\beta}_\text{foo}\times x_\text{foo} + \hat{\beta}_\text{bar} \times 1 + \hat{\epsilon}_t. $$

You can predict $\hat{\beta}_0 + \hat{\beta}_\text{bar} \times 1 + \hat{\epsilon}_t$ with a call to forecast.Arima() with the appropriate xreg parameter with foo=0 and bar=1:

baz <- forecast(model,h=1,xreg=data.frame(foo=0,bar=1))
baz
#          Point Forecast    Lo 80    Hi 80    Lo 95   Hi 95
# Jan 1961       443.7337 429.0981 458.3694 421.3504 466.117

Solve the equation:

x.foo <- (target-baz$mean[1])/model$coef["foo"]
x.foo
#      foo 
# 52.50988 

Finally, test that we didn't make any errors:

forecast(model,h=1,xreg=data.frame(foo=x.foo,bar=1))
#          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
# Jan 1961            400 385.3643 414.6357 377.6167 422.3833

Hurray!

To conclude, I'll emphasize again that this relies crucially on the fact that we use a regression with ARIMA errors, not an ARIMAX model, where we would also need to specify when we want this particular target value.