Optimal model-free prediction from multivariate time series

Forecasting a time series from multivariate predictors constitutes a challenging problem, especially using model-free approaches. Most techniques, such as nearest-neighbor prediction, quickly suffer from the curse of dimensionality and overfitting for more than a few predictors which has limited their application mostly to the univariate case. Therefore, selection strategies are needed that harness the available information as efficiently as possible. Since often the right combination of predictors matters, ideally all subsets of possible predictors should be tested for their predictive power, but the exponentially growing number of combinations makes such an approach computationally prohibitive. Here a prediction scheme that overcomes this strong limitation is introduced utilizing a causal preselection step which drastically reduces the number of possible predictors to the most predictive set of causal drivers making a globally optimal search scheme tractable. The information-theoretic optimality is derived and practical selection criteria are discussed. As demonstrated for multivariate nonlinear stochastic delay processes, the optimal scheme can even be less computationally expensive than commonly used suboptimal schemes like forward selection. The method suggests a general framework to apply the optimal model-free approach to select variables and subsequently fit a model to further improve a prediction or learn statistical dependencies. The performance of this framework is illustrated on a climatological index of El Niño Southern Oscillation.

Figure 1

Optimal prediction scheme. (i) Causal preselection in an example time-series graph (see text) for predicting $Y_{t+2}$. Even though $Z$ could have a high mutual information with $Y$, its influence is only indirect through the parents $X_t$ and $Y_{t+1}$ of $Y_{t+2}$. Among these the latter already lies in the unobserved future, but part of its information can still be recovered by measuring $Y_t$ and $X_{t-1}$ which share information along the paths marked with black arrows. These variables form the causal predictors ${\mathcal{P}}_{Y_{t+h}}$ (blue boxes) for $h=2$, which can be found by determining the Markov set. A suitable algorithm for this task will be discussed in Sec. 4a. All paths from nodes further in the past to $Y_{t+2}$ have to pass through this set. (ii) Selection of optimal subset from causal predictors. For all numbers of predictors $p$, the multivariate mutual information $\widehat{I}({\mathcal{P}}_{Y_{t+h}}^{\left(p\right)};Y_{t+h})$ of all possible subsets is estimated. (iii) The optimal predictors are those where the estimated MMI takes its maximum or can be obtained from cross-validation. For example, if the optimal predictors were ($Y_t,X_{t-1}$), only time series samples (blue shaded) of these predictors are used with a nearest-neighbor (Sec. IV C) or model-based (Sec. VIII) prediction scheme to forecast the future value $t+h$ of the time series of $Y$.

Figure 2

Multivariate mutual information (MMI, bottom red bars) and standardized root mean-squared prediction error (value proportional to lower end of gray bars at top) for all subsets of causal predictors for model (10), for details see Sec. 6. For each number of predictors $p$, the number of possible combinations varies according to the binomial coefficient $\left(\genfrac{}{}{0pt}{}{7}{p}\right)$. The predictor combinations are sorted by their prediction error. The maximum of MMI and also smaller values match the minimum prediction error very well. Note that MMI is estimated on the learning set while the prediction error is evaluated out-of-sample on the test set indicating that the bias of higher-dimensional MMIs here serves as a good proxy for overfitting as discussed in the text.

Figure 3

Scaling of computational complexity with the cardinality of the estimated causal predictors $|{\mathcal{P}}_{Y_{t+h}}|$ for ${\tau }_{max}=2$. Solid lines will be used for $N=10$, whereas dashed lines denote $N=30$. The curve for the causal CMI-selection scheme (black) differs only slightly from the constant curve for the MI-selection scheme (blue). The causal schemes have the advantage that their complexity rather depends on how many causal drivers are found, while the complexity of the MI- and CMI-selection schemes depends on the number of processes and the maximum lag. For the noncausal CMI-selection scheme (gray lines) we fixed the maximum number of predictors to $p_{max}=10$. To include the complexity due to the preselection step in the causal schemes we added 420 (1260) for $N=10$ ($N=30$) for these schemes according to Eq. (9) for $n_0=2$ and $n_i=3$ as in our numerical experiments (see Sec. 7). The plot shows that if the number of causal predictors can be reduced below 8 (10) here, the optimal scheme even takes less computation time than the noncausal CMI-selection scheme.

Figure 4

Comparison of prediction schemes for an ensemble of 500 realizations of model (10) for $a=0.4,b=2,c=0.4,\sigma =0.5$ with a time-series graph given in the inset [learning set length 500, test set length 125, nearest-neighbor prediction with $k=10$ neighbors, (C)MIs estimated from the learning set with parameter $k^{\mathrm{algo}}=50$ in the algorithm and $k^{\mathrm{CMI}/\mathrm{MMI}}=10$ for the optimization]. The box and whiskers plots give the ensemble median and the interquartile range of the standardized root-mean-squared prediction errors in the test set for each iteration step $p$ in the four schemes (from left to right: MI, CMI, causal-CMI, optimal scheme). The black line gives the median of the true minimal prediction error obtained by minimizing the out-of-sample prediction error for each number of predictors $p$ taken from the true causal drivers. For $p=3$, only the optimal scheme (red) selects the best (synergetic) predictors $Z_{t-1}^{\left(1\right)},Z_{t-1}^{\left(2\right)},Z_{t-1}^{\left(3\right)}$ and reaches the minimum possible error while the causal forward selection (black) first picks one of the less predictive $W_{t-1}^{(·)}$ and the pure forward selection (gray) and MI-based schemes first pick the two noncausal drivers $X_t^{\left(1\right)},X_t^{\left(2\right)}$. The nonoptimal schemes include the synergetic predictors only when the higher dimensionality is already worsening the prediction due to overfitting.

Figure 5

Results of numerical experiments for an ensemble of 500 trials of the synergetic model class (11) with time-series length $T=500$ (125 in the test set) for the four different prediction schemes (from left to right in the panels: MI in blue, CMI in gray, causal-CMI in black, optimal in red). (a) The four box plots show the ensemble interquartile range of (i) computational complexity (the green box on the right shows the added complexity due to the causal algorithm), (ii) the range of numbers of predictors $\widehat{p}$ selected by cross-validation (CV, here the green bar shows the number of causal predictors $|{\mathcal{P}}_{Y_{t+h}}|$ in the preselection step), (iii) the true positive rate (TPR), and (iv) false discovery rate (FDR). The latter are evaluated for each scheme at $p=8$, corresponding to the true number of causal drivers (or less if fewer causal drivers are estimated in the preselection step). (b) Box plots showing the median and interquartile range of the prediction error relative to the true minimal error obtained by minimizing the out-of-sample prediction error over all subsets of true causal drivers. On the left are the results if cross-validation is used to optimize $p$ for each scheme (whiskers show the 5% and 95% quantiles). The red box in the center shows the result for the heuristic optimality criterion. The range of boxes on the right shows the results for different thresholds $\lambda$ for the other schemes (only interquartile range). (c) Box and whiskers plots (showing the 5% and 95% quantiles) for the absolute prediction error of the optimal scheme at the cross-validated choice of $p$ for different numbers of nearest neighbors $k$.

Figure 6

Numerical experiments on a class of nonsynergetic but still nonlinearly coupled generalized additive models as analyzed in the supplement of Ref. [17]. Parameters as in Fig. 5 but with $N=4$ processes and $p_{max}=6$. The orange box plot in (b) shows the relative prediction error if the optimal predictors from the model-free selection scheme are used in conjunction with the linear autoregressive prediction model (12) (only the interquartile range shown).

Figure 7

Prediction of the El Niño Southern Oscillation (ENSO) index Nino3.4 in the period 2003–2014 (to December) using 1951–2002 as a learning set, causal algorithm run with significance threshold $I^{*}=0.03$ testing up to ${\tau }_{max}=12$ months. (a) Prediction error using nearest-neighbor prediction (solid red line) and linear prediction (dotted orange line) versus prediction step $h$. For both approaches the same optimal predictors obtained from the model-free scheme with cross-validated (fivefold within the learning set) choice of predictors are used. (b) Nino3.4 index with El Niño and La Niña events marked in red and blue, respectively. The black lines denote selected hindcasts and their $1\sigma$-prediction intervals (gray) using the linear prediction. The dots mark the starting time $t$ in May of each year, and the predicted values range from June ($h=1$) to October ($h=5$). The arrows mark correct (red), missed (gray), and false (black) hindcasts of the Nino3.4 index exceeding $0.5{}^{\circ }\mathrm{C}$ with more than 49% probability. The green line marks a real forecast starting in December 2014 giving a probability for the Nino3.4 index to stay above $0.5{}^{\circ }\mathrm{C}$ of 55–70% for the months until May 2015.

Figure 8

As in Figs. 6 and 6 but for a larger range of significance thresholds $I^{*}=0.001,0.002,0.004,0.005,0.006,0.008,0.009,0.01$ (row-wise from top left to bottom right).

Figure 9

As in Fig. 5 (left column) and Fig. 6 (right column) but for a larger range of time-series lengths $T=800,500,300$ (top to bottom, test set lengths $200,125,75$).