Incorporating Prior Information Into Time Series Prediction

Question

Suppose I have data on my child C's height measured every week. Presumably there is a positive trend, due to growth, and some noise due to measurement errors, and maybe even seasonality (winter boots or haircuts add extra height). If want to predict C's height into the future (at say t+s), I might use an ARIMA model, which would represent C's height next week as a function of his height in previous weeks.

Now suppose that in the past, I had 2 other children, A and B, who are older than C. Thus I have prior information about the distribution my children's heights at time t+s. How can I incorporate this into my prediction?

For instance, consider the following data. Observe that while C6 (the height of C at t=6), is unknown, A6 and B6 are known.

set.seed(1)
A = seq(1,6,.5) + rnorm(11)/4 
B = seq(1,6,.5) + rnorm(11)/4
C = c((seq(1,3,.5) + rnorm(5)/4), rep(NA, 6))
df = data.frame(A, B, C)

> df
           A        B        C
1  0.8433865 1.097461 1.018641
2  1.5459108 1.344690 1.002662
3  1.7910928 1.446325 2.154956
4  2.8988202 2.781233 2.485968
5  3.0823769 2.988767 2.961051
6  3.2948829 3.495952       NA
7  4.1218573 4.235959       NA
8  4.6845812 4.705305       NA
9  5.1439453 5.148475       NA
10 5.4236529 5.729744       NA
11 6.3779453 6.195534       NA

I considered using 2NN (2 nearest neighbors) to predict C6. 2NN estimates C6 as the avg of 2.485968 and 2.961051, which is not great (although we get better results if we difference the data).
Alternatively we might say C6 is the expected value of the prior observations so (3.2948829, 3.495952)/2 in this case.

What are some other methods I could use to incorporate this prior information into my forecast? Could I combine my ARIMA forecast with this prior information somehow to form an ensemble forecast?

I've also started looking into Bayesian time series, dynamic linear models, and state space models, but don't have much experience in these areas and could use a pointer.

Answer 1

Great Question !

"Could I combine my ARIMA forecast with this prior information somehow to form an ensemble forecast?"

I have been involved with a commercial time series forecasting package called AUTOBOX and have incorporated delphi-type predictor series where the user provides probabilities of intervals and this is then used to monte-carlo a family of possible values for future values of input series where the user normally delivers 1 point estimates based upon perfect knowledge.

The "realizations" developed this way are then inter-joined with the arima simulations providing a family of ensemble forecast values for the dependent series that might also be effected by possible anomalies identified in the analysis stage via Intervention Detection schemes.

You should be able to program this with this advice as I have done. This problem/opportunuties arises quite naturally when the predictor series distribution can be "pre-guessed" such as alternative hypotheses for the price of oil for the next period. Armed with priors like this one can select from alternative offerings those with the greatest expected reward.

Answer 2

I think a good way to do this is via Bayesian Structural Time Series (BSTS). I found out about this approach via these 2 sites (1, 2). I would still be interested in other approaches.

Here is the example done with the bsts package in R. I use a time series component and a regression component. The regression component incorporates the prior information. A stack and slab prior is used on the regression component.

Below is a plot of the time-series features. I intentionally made C unlike the other features, and the response variable (D), in order to test the feature selection ability of BSTS.

The plot below shows that the BSTS model correctly observes that feature C isn't useful for predicting D.

Actual (Blue) vs Predicted (Red)

Code is below:

library(ggplot2, bsts, ggplotly, data.table)

# generate some data
set.seed(1)
n = 20
train_size = 10
A = seq(1,n) + rnorm(n)
B = seq(1,n) + rnorm(n)
# this variable is not like the others
C = rnorm(n) + 5*sin(seq(1,n))
D = seq(1,n) + rnorm(n)
X = data.table(A, B, C, D)
# transform the data for ggplot
long_data = melt(X)
m[, t := seq_len(.N), by = variable]
g1 = ggplot(data=m, aes(x=t, y=value, colour=variable)) + geom_line() + labs(title="Evolution of Parameters over Time")
ggplotly(g1)
#break the data into training/testing data
train_ind = seq(1,train_size)
train_X = X[train_ind,]
test_X = X[-train_ind,]


ss <- AddLocalLinearTrend(list(), train_X$D)
model4 <- bsts(D ~ .,
               state.specification = ss,
               niter = 1000,
               data = train_X,
               expected.model.size = 3) 
plot(model4, "components")
# observe that the model can tell that C isn't strongly related to D, but A and B are. 
plot(model4, "coef")
pred4 <- predict(model4, newdata = test_X, horizon = 24)
# plot predictions, vs actual (in red)
plot(pred4, ylim=c(0,50))
lines((max(train_ind)+1):nrow(X), test_X$D, col="red")