Consider a multivariate time series $(x_t,{\bf Y}_t)$ $1\le t \le n$ taking values in $\mathbb{R}^{d+1}$, and suppose that we wish to forecast $x_t$ using its own path as well as the "exogenous" series ${\bf Y}_t$, for instance using a linear time series model. For the sake of simplicity assume all series are (jointly) stationary. One approach that seems clear is to fit a VAR model, and then use the corresponding forecasts for $x_t$, but this seems sub-optimal, since the VAR model fit using ordinary least squares would not optimize the parameters to forecast $x_t$ alone.
Transfer function models also do not seem useful here, since, when forecasting $x_{n+h}$, they presume the corresponding "covariates" ${\bf Y}_s,\;\;s\le n+h$ are known when making the forecast, when we would only have observed ${\bf Y}_s,\;\;s\le n$. I guess one might forecast the ${\bf Y}_t$ series to input into the transfer function model, but again this doesn't seem optimal. Having done a fair bit of googling on the topic, I could not find any real guidance on what seems like a very standard problem.
Can anyone point me in the right direction? Or perhaps some of these simple ideas (VAR, transfer function model with forecasted covariates) is more optimal than I think?
In response to IrishStat's comment, I'll post a simple example using a transfer function model (at least a simple version of one using simple linear regression with ARIMA errors as implemented in auto.arima). Suppose that we wish to forecast the cardiovascular mortality in Los Angeles county, and we also have access to the daily temperature and particulate matter pollution concentration (data from the astsa package in R). I can do this in R as follows:
#Begin R code#
library(astsa)
library(forecast)
library(TSA)
\#taking data at monthly resolution (every four weeks, so that the approximate
\#seasonality/frequency is 13
cmort2=ts(lap[seq(1,508,by=4),3],frequency = 13)
temp2=ts(lap[seq(1,508,by=4),4],frequency = 13)
part2=ts(lap[seq(1,508,by=4),11],frequency = 13)
dat.mat=cbind(as.numeric(temp2),as.numeric(part2))
\#producing forecasts for temp2 and part2 using auto.arima to be
\#fed into the arimax model
temp.mod=auto.arima(temp2)
part.mod=auto.arima(part2)
temp.for=forecast(temp.mod,h=12)
plot(temp.for)
part.for=forecast(part.mod,h=12)
plot(part.for)
temp.for=ts(temp.for\\\$mean,frequency = 13)
part.for=ts(part.for\\\$mean,frequency = 13)
dat.mat.for=cbind(temp.for,part.for)
ar.regf=auto.arima(cmort2, xreg=dat.mat)
x=forecast(ar.regf,xreg=dat.mat.for,h=12)
autoplot(x)
\#for comparison, a simple SARIMA model excluding the covariates
ar.noregf=auto.arima(cmort2)
x.noreg=forecast(ar.noregf,h=12)
autoplot(x.noreg)
##end R code##
I guess my questions about this are:
- Is this really the best/a reasonable thing to do to forecast $x_t$=LA cardiac mortality? The sub-optimal part seems to be that we must input forecasts for the covariates, which rely on modeling the covariates separately from the response.
- I assume that the confidence bands produced in forecasting $x_t$ are not accurate since they do not account for the uncertainty in the forecast of the covariates. Does anybody know if this is the case? I could imagine fixing this myself by producing confidence bands via simulation, but one wonders how to automatically incorporate the forecast uncertainty into the confidence bands for the transfer function model.
