I'm working through time-series forecasting, using models such as ETS, ARIMA, and vector autoregression as described in several texts (for example, Hyndman, R.J., & Athanasopoulos, G. (2021) Forecasting: principles and practice, 3rd edition, OTexts: Melbourne, Australia. OTexts.com/fpp3).

I've created a hypothetical where I assume I have only the first [12] months of time-series data and I forecast for months [13-24] based on the actuals for months [1-12]. I generate simulation paths for months [13-24], and distributions thereof. I then compare those forecasted simulation paths/distributions for months [13-24] with the actual data for months 13-24 in order to assess forecast reasonableness. Results using ETS and ARIMA have been fine, with some minor adjustment such as using logs.

However, these traditional time-series forecasting methods analyze/forecast essentially a single-line, depicted as the heavier trend line in the below image using my example data and labeled mean. In my data, that heavier trend line is simply an average of many underlying elements with disparate trends. The below is a simplified example of my actual data for the sake of post replicability and all of my actual curves take the form of nice smooth logarithmic functions. In the below example, there are elements v, w, x, y, and z, and their mean is mean in the example data frame. But the trends of the underlying elements in my actual data do look like this example data in terms of dispersion around the mean. Values never fall below zero.

For time-series forecasting such as for this form of example data, are there any other methods I should be considering, that take into account the additional information I have at hand for the many underlying elements? (In my actual data I have 48 months and 60,000 + elements trending over those 48 months).

enter image description here

Code to generate the above:


DF <- data.frame(
  mo = 1:24,
  v = c(rep(0,24)),
  w = c(0,0.1,rep(0.2,12),seq(0.2,0.5,length.out=10)),
  x = c(0,0,seq(0,0.5,length.out = 10),0.5,0.5,seq(0.5,0.98,length.out = 10)),
  y = seq(0, 1.5, length.out = 24),
  z = seq(0, 2.5, length.out = 24)

DF$mean <- rowMeans(DF[,2:6])

DF_reshape <- data.frame(
  x = DF$mo,                           
  y = c(DF$v, DF$w, DF$x,DF$y,DF$z,DF$mean),
  group = c(rep("v", nrow(DF)),
            rep("w", nrow(DF)),
            rep("x", nrow(DF)),
            rep("y", nrow(DF)),
            rep("z", nrow(DF)),
            rep("mean", nrow(DF))

ggplot(DF_reshape, aes(x, y, col = group)) +  
  geom_line() +
  geom_line(data = filter(DF_reshape,group == "mean"), linewidth = 2) +
  labs(x = "x axis = number of months elapsed") 
  • 1
    $\begingroup$ The first thing that comes to my mind would be ensembling: forecast each series separately, then take the average as an expectation forecast. (Incidentally, you write about using logs. I assume you use a bias correction when back-transforming?) $\endgroup$ Commented Mar 28, 2023 at 16:55

1 Answer 1


This issue is addressed in chapter 11 of the book you referenced (https://otexts.com/fpp3/hierarchical.html).

Basically, if you care about forecasts for the sub-series, you can (1) forecast the total and split it amongst the sub-series, or (2) forecast each sub-series individually and aggregate them to the total. Or, you can do both at once in the unified method discussed in that chapter.

If you go with approach (2), a simple starting place might be pooling all of the series together and estimating one common ARIMA model. Of course, this may not be a good idea if your series have very different ACFs. The opposite (separate ARIMA models for each series) would be the other extreme. Most likely, there is some way to group series and estimate their AR coefficient jointly.

One alternative to explicit grouping/pooling is to use a factor model (or principal components), extracting the first $k<N$ factors (where $N$ is the total number of series), and using those to form forecasts for the individual series. For two series ($y_1$ and $y_2$), with one factor ($f_t$) we would have: $$ \begin{bmatrix}y_{1,t} \\ y_{2,t} \end{bmatrix} = \begin{bmatrix}\alpha_1 \\ \alpha_2 \end{bmatrix} + \begin{bmatrix}\beta_1 \\ \beta_2 \end{bmatrix} f_t + \varepsilon_t \\ f_t = \phi f_{t-1} + v_t $$

Basically, the two series are driven by the same unobserved common factor. Once you estimate the system (either by MLE or Bayesian methods) you can forecast the factor and plug those forecasts into the top equation. A rough estimate of this model (I can't speak to it's consistency properties) would be to extract the first PC of $Y = [y_1 \ y_2]$ and plug that into the top equation. Then estimate an AR(1) on the first PC. Forecast the first PC according to the AR(1), and plug in the forecasts to the top equation.

If you were using the first $k$ factors in a larger context ($N>2$), you would simply extract the first $k$ factors, use them to get estimates of the parameters in the top equation, fit AR(1) models to each of the $k$ factors, forecast each of the $k$ factors, and then plug those forecasts into the estimated top equation to get forecasts for each of your series.


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.