How to meaningfully compute the accuracy of a multi-step forecast produced by a model

Question

I am trying to measure the accuracy of my model in producing a multi-step forecast and I have read a lot of different opinions on the matter and am now rather confused.

The goal of my model is to produce a forecast $H$ weeks into the future. Let's say we are at week $t$ and $H=4$ my model will produce a forecast for weeks $\{t+1, t+2, t+3, t+4\}$. The forecasts will be denoted $\{\hat{y}_{t+1},\cdots,\hat{y}_{t+H}\}$ and the true future observations $\{y_{t+1},\cdots,y_{t+H}\}$ of the time series $\{y_1,\cdots,y_t\}$ and the forecast horizon $h\in \{1,\cdots,H\}$.

From this I would like an accuracy measure for each week so that I can say at week $t+1$ my model has an error of 80% at week $t+2$ it has an accuracy of 70%. To compute the error I perform the following calculation $e_{t+h}=\hat{y}_{t+h}-y_{t+h}$ which is pretty standard.

What I also would like to compute is the average error of my model over the whole forecasting horizon so the $\bar{e_H} = \frac{\sum^{H}_{h=1}e_{t+h}}{H}$. The idea is so that I can first pick parameters for my model as I am dealing with many time series based on an average error of the whole forecast horizon. After these have been selected I can further analyse that this models performance will deteriorate by $X$ for each step of the forecast for example.

Collecting these measurements isn't a problem what I am struggling with is collecting them in a way that gives me independent results so that I can get a statistically significant measurement of my results and make a meaningful parameter selection.

I found a method detailed in several places called Time Series Cross Validation (see here, here and here). All of these sources however deal with a single step forecast so that they are able to keep independence between the sub-samples (or folds). For example for the time series $T=1,\cdots,6$ using the following validation method I would get the following train and test splits respectively (assuming a minimum train size of 3):

[1, 2, 3],       [4]
[1, 2, 3, 4],    [5]
[1, 2, 3, 4, 5], [6].

The method can also be applied for an $h$ step ahead forecast lets say $h=2$ and the following train and test sets would be produced

[1, 2, 3],    [5]
[1, 2, 3, 4], [6].

These methods do not detail however how I should handle when my test set should consist of more than one value as I am interested in both the average forecast error over the whole horizon and the individual forecast errors. Up until now I had been using a method which produced train and test sets like the following

[1, 2, 3],       [4, 5]
[1, 2, 3, 4],    [5, 6]
[1, 2, 3, 4, 5], [6, 7].

It was recently pointed out to me that that my comparisons of the error between these sets is not independent and therefore cannot be considered meaningful. So I would like to know how do I perform a meaningful measurement of the accuracy of my model using Time Series Cross Validation when all the individual errors of the forecast horizon are important and not just the last one?

I had thought that maybe the following would work for a train and test set split

[1, 2, 3],             [4, 5]
[1, 2, 3, 4, 5],       [6, 7]
[1, 2, 3, 4, 5, 6, 7], [8, 9].

The idea being to just expand the window by the full forecast horizon each time or to use a sliding window instead of expanding

[1, 2, 3], [4, 5]
[3, 4, 5], [6, 7]
[5, 6, 7], [8, 9].

Would either of these methods be considered a valid way to make a valid measurement of the average accuracy of my model?

EDIT

This edit is a follow up to an answer given below in order to provide more information.

If I then decide to use normal k-fold cross validation to validate the accuracy of my model, is it also valid to have a multi-step test set instead of a test set with just one result in it. So to use your example I take the dataset below

[1, 2, 3 | 4]
[2, 3, 4 | 5]
[3, 4, 5 | 6]
[4, 5, 6 | 7]

Which I can split into the following

x_train = [[1,2,3],
           [2,3,4],
           [3,4,5]]

y_train = [4, 5, 6]

an x_test = [4, 5, 6] and a y_test = [7]? To put the problem in ML terms.

For my specific problem I am dealing with a multi-step forecast is it then also valid to use kfold cross validation when I have more than one target?

For example consider the following with a window width of three and two target variables (one for each output)

X1  X2  X3  | y1  y2
100 110 120 | 130 140 
110 120 130 | 140 150
120 130 140 | 150 160
130 140 150 | 160 170.

This could then also be split using kfold cross validation to train on

X1  X2  X3  | y1  y2
100 110 120 | 130 140 
110 120 130 | 140 150
120 130 140 | 150 160

test on

X1  X2  X3  | y1  y2
130 140 150 | 160 170.

then train on

X1  X2  X3  | y1  y2
110 120 130 | 140 150
120 130 140 | 150 160
130 140 150 | 160 170.

and test on

X1  X2  X3  | y1  y2
100 110 120 | 130 140,

and so on, correct?

Similarly when framing my problem to be used by ARIMA say, I would simply follow this procedure as well? Is it correct to say that I will then have to implement my own ARIMA and and ETS methods as statsmodels doesn't support feeding time series data formatted this way?

Answer 1

First of all, I wouldn't use:

$\bar{e_H} = \frac{\sum^{H}_{h=1}e_{t+h}}{H}$

Because negative and positive error terms will cancel each other out.

Instead I would use:

$MSE = \frac{\sum^{H}_{h=1}e_{t+h}^2}{H}$ or $RMSE=\sqrt{MSE}$

Also frequently used for time series problems is the MAPE:

$\bar{e_H} = \frac{\sum^{H}_{h=1}|{e_{t+h}}|}{H}$.

I am struggling with is collecting them in a way that gives me independent results so that I can get a statistically significant measurement of my results and make a meaningful parameter selection.

For this, you need to specify first was it the purpose of your forecast in the first place. There are many ways of evaluating error and depending on the domain you are applying your forecasting approach to, one method might be more suitable than the other (See here for some details, and here for some discussion).

For example, if you are forecasting weekly demand for a product being shipped from a supplier, and your lead time for that supplier is 3 weeks then all you care about is the forecast values for $H=3$, so your cross validation split should look like:

[1, 2, 3],       [6]
[1, 2, 3, 4],    [7]
[1, 2, 3, 4, 5], [8]

etc... And your evaluation metric should be:

$|e_{t+3}|$ or $e_{t+3}^2$.

However, If you need to forecast values over multiple steps in the forecast horizon, then the above mentioned MSE and MAPE should work.

Note also that if you are performing recursive forecasting using traditional methods like ETS or ARIMA, then you can forgo CV all together and use the AIC or the BIC as a model selection criteria.

---

Regarding independence, CV is a little bit tricky in that regard but your overall approach is still correct, because you are preserving the order of the time series when you do a CV split of the type:

[1, 2, 3],       [4, 5]
[1, 2, 3, 4],    [5, 6]
[1, 2, 3, 4, 5], [6, 7].

Independence would have been an issue for some time series models like ARIMA and ETS if you had used normal CV instead of time series CV, for example, leave one out CV approach such as:

[1, 2, 4, 5, 6], [3]
[1, 2, 3, 5, 6], [4]
[1, 2, 3, 4, 6], [5]

Wouldn't work with most time series methods because it ignores the ordered nature of time series data, i.e. prediction the value at [3] using values from [4] and [5] means we are using the future to predict the past, which doesn't make sense. I think this maybe why independence was pointed out to you as an issue for time series.

But things get more interesting. Recently, Bergmeir, Hyndman and Koo, have shown that for purely auto-regressive models, such as AR(p) models, Neural Networks, or Support Vector Regression, even normal CV can be used, be used, as long as you format your training data so that it looks like a supervised machine learning problem instead of a time series problem.

To understand how this works, first you need to note that pure autoregressive method use only a fixed number of previous periods, so that instead of:

[1, 2, 3],          [4]
[1, 2, 3, 4],       [5]
[1, 2, 3, 4, 5],    [6]
[1, 2, 3, 4, 5, 6], [7]

You would have for example (for an order 3 (i.e. 3 lags) autoregressive model):

[1, 2, 3],   [4]
[2, 3, 4],   [5]
[3, 4, 5],   [6]
[4, 5, 6],   [7]

Then you need to think of your data in the ML format instead of the time series format:

[1, 2, 3 | 4]
[2, 3, 4 | 5]
[3, 4, 5 | 6]
[4, 5, 6 | 7]

So that that your data point is not just a single value like [4] but instead a vector with 3 inputs and one target value: [1, 2, 3 | 4]

In this case, using normal cross validation (instead of time series cross validation) is valid, because dependence is no longer an issue, and you can train your model on:

[2, 3, 4 | 5]
[3, 4, 5 | 6]
[4, 5, 6 | 7]

And test with:

[1, 2, 3 | 4]

then train the model on:

[1, 2, 3 | 4]
[3, 4, 5 | 6]
[4, 5, 6 | 7]

And test with:

[2, 3, 4 | 5]

And so on....

See "A Note on the Validity of Cross-Validation for Evaluating Autoregressive Time Series Prediction", by Christoph Bergmeir, Rob J Hyndman, and Bonsoo Koob for details.

---

Addressing the additional final part of the question:

Yes, you can use a similar K-fold CV approach for multistep forecasts, as long as you are using purely auto-regressive methods.

Which means AR models (or ARMA/ARIMA models with q=0) and ML models will work, but not ARIMA models (q ≠ 0) or Exponential Smoothing models.

Then again ARIMA models and Exponential Smoothing models are inherently sequential, so I don't know that you can use them for direct forecasting at all, and the question is probably moot.