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.

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.


2 Answers 2


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.

  • $\begingroup$ Thanks! I'll investigate this approach. By the way, I have very much enjoyed reading your insightful answers to other questions on this site. Do you have an opinion on the bsts approach I posted below? $\endgroup$ Commented Aug 3, 2019 at 20:12
  • $\begingroup$ I am not a USER OF bsts biur wouldn't mind learning more about what you like ... please contact me at my land line as I have many questions and little patience for keyed (i.e. non-voice) interaction. perhaps both of us jointly learn. $\endgroup$
    – IrishStat
    Commented Aug 3, 2019 at 20:22
  • $\begingroup$ If you are happy with my response please accept it to close the question $\endgroup$
    – IrishStat
    Commented Aug 15, 2019 at 8:19

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 Features

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

Feature Selection

Actual (Blue) vs Predicted (Red) Actual (Blue) vs Predicted (Red)

Code is below:

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

# generate some data
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")
#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")
  • $\begingroup$ What is 'm' in 'm[, t := seq_len(.N), by = variable]'? $\endgroup$
    – Masoud
    Commented Nov 3, 2021 at 15:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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