2. M-step-ahead predictions

Assume we have a time series of observations $y=(y_1,y_2,\dots ,y_N)$ and let L be the minimum number of observations from the series that we will require before making predictions for future data. Depending on the application and how informative the data are, it may not be possible to make reasonable predictions for $y_{i+1}$ based on $(y_1,\dots ,y_i)$ until i is large enough so that we can learn enough about the time series to predict future observations. Setting L = 10, for example, means that we will only assess predictive performance starting with observation $y_{11}$, so that we always have at least 10 previous observations to condition on.

In order to assess M-SAP performance, we would like to compute the predictive densities (1) $p(y_{i+1:M}|y_{1:i})=p(y_{i+1},\dots ,y_{i+M}|y_1,\dots ,y_i)$(1) for each $i\in \{L,\dots ,N$-$M\}$, where we use $y_{1:i}=(y_1,\dots ,y_i)$ and $y_{i+1:M}=y_{(i+1):(i+M)}=(y_{i+1},\dots ,y_{i+M})$ to shorten the notation . ² As a global measure of predictive accuracy, we can use the expected log predictive density [ELPD; 8], which, for M-SAP, can be defined as (2) $\mathrm{ELPD}=\sum _{i=L}^{N-M}\int p_t({\tilde{y}}_{i+1:M})\log p({\tilde{y}}_{i+1:M}|y_{1:i})\mathrm{d}{\tilde{y}}_{i+1:M}.$(2) The distribution $p_t({\tilde{y}}_{i+1:M})$ describes the true data generating process for new data ${\tilde{y}}_{i+1:M}$. As these true data generating processes are unknown, we approximate the ELPD using LFO-CV of the observed responses $y_{i+1:M}$, which constitute a particular realization of ${\tilde{y}}_{i+1:M}$. This approach of approximating the true data generating process with observed data is fundamental to all cross-validation procedures. As we have no direct access to the underlying truth, the error implied by this approximation is hard to quantify but also unavoidable [c.f., 18].

Plugging in the realization $y_{i+1:M}$ for ${\tilde{y}}_{i+1:M}$ leads to [c.f., 7,18]: (3) ${\mathrm{ELPD}}_{\mathrm{LFO}}=\sum _{i=L}^{N-M}\log p(y_{i+1:M}|y_{1:i}).$(3) The quantities $p(y_{i+1:M}|y_{1:i})$ can be computed with the help of the posterior distribution $p(\theta |y_{1:i})$ of the parameters θ conditional on only the first i observations of the time series: (4) $p(y_{i+1:M}|y_{1:i})=\int p(y_{i+1:M}|y_{1:i},\theta )p(\theta |y_{1:i})\mathrm{d}\theta .$(4) In order to evaluate predictive performance of future data, it is important to predict $y_{i+1:M}$ only conditioning on the past data $y_{1:i}$ and not on the present data $y_{i+1:M}$. Including the present data in the posterior estimation, that is, using the posterior $p(\theta |y_{1:(i+M)})$ in (4), would result in evaluating in-sample fit. This corresponds to what Vehtari et al. [8] call log-predictive density (LPD), which overestimates predictive performance for future data [8].

Most time series models do not have conditionally independent observations, that is, $y_{i+1:M}$ depend on $y_{1:i}$ even after conditioning on θ. As such, we cannot simplify the integrand in (4) and always need to take $y_{1:i}$ into account when computing the predictive density of $y_{i+1:M}$. The concept of conditional independence of observations is closely related to the concept of factorizability of likelihoods. For the purpose of LFO-CV, we can safely use the time-ordering naturally present in time-series data to obtain a factorized likelihood even if observations are not conditionally independent. See Bürkner et al. [19] for discussion on computing predictive densities of non-factorized models and factorizability in general.

In practice, we will not be able to directly solve the integral in (4), but instead have to use Monte-Carlo methods to approximate it. Having obtained S random draws $({\theta }_{1:i}^{(1)},\dots ,{\theta }_{1:i}^{(S)})$ from the posterior distribution $p(\theta |y_{1:i})$, we can estimate $p(y_{i+1:M}|y_{1:i})$ as (5) $p(y_{i+1:M}|y_{1:i})\approx \frac{1}{S}\sum _{s=1}^{S}p(y_{i+1:M}|y_{1:i},{\theta }_{1:i}^{(s)}).$(5) In this paper, we use ELPD as a measure of predictive accuracy, but M-SAP (and the approximations we introduce below) may also be based on other global measures of accuracy such as the root mean squared error (RMSE) or the median absolute deviation (MAD). The reason for our focus on ELPD is that it evaluates a distribution rather than a point estimate (like the mean or median) to provide a measure of out-of-sample predictive performance, which we see as favourable from a Bayesian perspective [7]. The code we provide on GitHub (https://github.com/paul-buerkner/LFO-CV-paper) is modularized to support arbitrary measures of accuracy as long as they can be represented in a pointwise manner, that is, as increments per observation. In Appendix C, we also provide additional simulation results using RMSE instead of ELPD.

2.1. Approximate M-step-ahead predictions

The equations above make use of posterior distributions from many different fits of the model to different subsets of the data. To obtain the predictive density $p(y_{i+1:M}|y_{1:i})$, a model is fit to only the first i data points, and we need to do this for every value of i under consideration (all $i\in \{L,\dots ,N$-$M\}$). Below, we present a new algorithm to reduce the number of models that need to be fit for the purpose of obtaining each of the densities $p(y_{i+1:M}|y_{1:i})$. This algorithm relies in a central manner on Pareto smoothed importance sampling [8,9], which we will briefly review next.

2.1.1. Pareto smoothed importance sampling

Importance sampling is a technique for computing expectations with respect to some target distribution using an approximating proposal distribution that is easier to draw samples from than the actual target. If $f(\theta )$ is the target and $g(\theta )$ is the proposal distribution, we can write any expectation of some function $h(\theta )$ with respect to f as (6) ${\mathbb{E}}_f[h(\theta )]=\int h(\theta )f(\theta )\mathrm{d}\theta =\frac{\int [h(\theta )f(\theta )/g(\theta )]g(\theta )\mathrm{d}\theta }{\int [f(\theta )/g(\theta )]g(\theta )\mathrm{d}\theta }=\frac{\int h(\theta )r(\theta )g(\theta )\mathrm{d}\theta }{\int r(\theta )g(\theta )\mathrm{d}\theta }$(6) with importance ratios (7) $r(\theta )=\frac{f(\theta )}{g(\theta )}.$(7) Accordingly, if ${\theta }^{(s)}$ are S random draws from $g(\theta )$, we can approximate (8) ${\mathbb{E}}_f[h(\theta )]\approx \frac{\sum _{s=1}^{S}h({\theta }^{(s)})r({\theta }^{(s)})}{\sum _{s=1}^{S}r({\theta }^{(s)})},$(8) provided that we can compute the raw importance ratios $r({\theta }^{(s)})$ up to some multiplicative constant. The raw importance ratios serve as weights on the corresponding random draws in the approximation of the quantity of interest.

The main problem with this approach is that the raw importance ratios tend to have large or infinite variance and results can be highly unstable. In order to stabilize the computations, we can perform the additional step of regularizing the largest raw importance ratios using the corresponding quantiles of the generalized Pareto distribution fitted to the largest raw importance ratios. This procedure is called Pareto smoothed importance sampling [PSIS; 8,9,20] and has been demonstrated to have a lower error and faster convergence rate than other commonly used regularization techniques [9].

In addition, PSIS comes with a useful diagnostic to evaluate the quality of the importance sampling approximation. The shape parameter k of the generalized Pareto distribution fit to the largest importance ratios provides information about the number of existing moments of the distribution of the weights (smoothed ratios) and the actual importance sampling estimate. When k<0.5, the weight distribution has finite variance and the central limit theorem ensures fast convergence of the importance sampling estimate with increasing number of draws. This implies that approximate LOO-CV via PSIS is highly accurate for k<0.5 [9]. For $0.5\le k<1$, a generalized central limit theorem holds, but the convergence rate drops quickly as k increases [9]. Using both mathematical analysis and numerical experiments, PSIS has been shown to be relatively robust for $k<0.7$ [8,9]. As such, the default threshold is set to 0.7 when performing PSIS LOO-CV [8,20].

2.1.2. PSIS applied to M-step-ahead predictions

We now turn back to the task of estimating M-step-ahead performance for time series models. First, we refit the model using the first L observations of the time series and then perform a single exact M-step-ahead prediction step for $p(y_{L+1:M}|y_{1:L})$ as per (4). Recall that L is the minimum number of observations we have deemed acceptable for making predictions (setting L = 0 means the first data point will be predicted only based on the prior). We define $i^{\star }=L$ as the current point of refit. Next, starting with $i=i^{\star }+1$, we approximate each $p(y_{i+1:M}|y_{1:i})$ via (9) $p(y_{i+1:M}|y_{1:i})\approx \frac{\sum _{s=1}^S{w}_i^{(s)}p(y_{i+1:M}|y_{1:i},{\theta }^{(s)})}{\sum _{s=1}^S{w}_i^{(s)}},$(9) where ${\theta }^{(s)}={\theta }_{1:i^{\star }}^{(s)}$ are drawn from the posterior distribution based on the first $i^{\star }$ observations and $w_i^{(s)}$ are the PSIS weights obtained in two steps. First, we compute the raw importance ratios (10) $\begin{array}{rl}{r}_i^{(s)}& =r_i({\theta }^{(s)})=\frac{f_{1:i}({\theta }^{(s)})}{f_{1:i^{\star }}({\theta }^{(s)})}\propto \frac{p({\theta }^{(s)})\prod _{j\in 1:i}p(y_j|y_{1:(j-1)},{\theta }^{(s)})}{p({\theta }^{(s)})\prod _{j\in 1:i^{\star }}p(y_j|y_{1:(j-1)},{\theta }^{(s)})}\\ & =\prod _{j\in (i^{\star }+1):i}p(y_j|y_{1:(j-1)},{\theta }^{(s)}),\end{array}$(10) and then stabilize them using PSIS as described in Section 2.1.1. The function $f_{1:i}$ denotes the posterior distribution based on the first i observations, that is, $f_{1:i}=p(\theta |y_{1:i})$, with $f_{1:i^{\star }}$ defined analogously. The index set $(i^{\star }+1):i$ indicates all observations which are part of the data for the model $f_{1:i}$ whose predictive performance we are trying to approximate but not for the actually fitted model $f_{1:i^{\star }}$. The proportional statement arises from the fact that we ignore the normalizing constants $p(y_{1:i})$ and $p(y_{1:i^{\star }})$ of the compared posteriors, which leads to a self-normalized variant of PSIS [c.f. 8].

Continuing with the next observation, we gradually increase i by 1 (we move forward in time) and repeat the process. At some observation i, the variability of the importance ratios $r_i^{(s)}$ will become too large and importance sampling will fail. We will refer to this particular value of i as $i_1^{\star }$. To identify the value of $i_1^{\star }$, we check for which value of i does the estimated shape parameter k of the generalized Pareto distribution first cross a certain threshold τ [9]. Only then do we refit the model using the observations up to $i_1^{\star }$ and restart the process from there by setting ${\theta }^{(s)}={\theta }_{1:i_1^{\star }}^{(s)}$ and $i^{\star }=i_1^{\star }$ until the next refit. An illustration of this procedure is shown in Figure 1.

Approximate leave-future-out cross-validation for Bayesian time series models

All authors

Paul-Christian Bürkner , Jonah Gabry & Aki Vehtari

https://doi.org/10.1080/00949655.2020.1783262

Published online:

25 June 2020

Figure 1. Visualization of PSIS approximated one-step-ahead predictions. Predicted observations are indicated by X. In the shown example, the model was last refit at the $i^{\star }=4$th observation.

Display full size

Figure 1. Visualization of PSIS approximated one-step-ahead predictions. Predicted observations are indicated by X. In the shown example, the model was last refit at the $i^{\star }=4$th observation.

This bears some resemblance to Sequential Monte Carlo, also known as particle or Monte Carlo filtering [e.g. 15,16,21,22], in that applying SMC to state space models also entails moving forward in time and using importance sampling to approximate the next state from the information in the previous states [16,22]. However, in our case we are assuming we can sample from the posterior distribution and are instead concerned with estimating the model's predictive performance. Unlike SMC, PSIS-LFO-CV also entails a full recomputation of the model via Markov chain Monte Carlo (MCMC) once importance sampling fails.

In some cases, we may only need to refit once and in other cases we will find a value $i_2^{\star }$ that requires a second refitting, maybe an $i_3^{\star }$ that requires a third refitting, and so on. We refit as many times as is required (only when $k>\tau$) until we arrive at observation $i=N$-M. A detailed description of the algorithm in the form of pseudo code is provided in Appendix A. If the data are comprised of multiple independent time series, the algorithm can be applied to each of the time series separately and the resulting ELPD values can be summed up afterwards. If the data are comprised of multiple dependent time series, the algorithm should be applied to the joint likelihood across all time-series for each observation i in order to take the dependency into account.

Instead of moving forward in time, it is also possible to do PSIS-LFO-CV moving backwards in time. However, in that case the target posterior is approximated by a distribution based on more observations and, as such, the proposal distribution is narrower than the target. This can result in highly influential importance weights more often than when the proposal is wider than the target, which is the case for the forward PSIS-LFO-CV described above. In Appendix B, we show that moving backwards indeed requires more refits than moving forward, and without any increase in accuracy. When we refer to the PSIS-LFO-CV algorithm in the main text, we are referring to the forward version.

The threshold τ is crucial to the accuracy and speed of the PSIS-LFO-CV algorithm. If τ is too large, then we need fewer refits but accuracy will suffer. If τ is too small, then accuracy will be higher but many refits will be required and the computation time may be impractical. Limiting the number of refits without sacrificing too much accuracy is essential since almost all of the computation time for exact cross-validation of Bayesian models is spent fitting the models (not calculating the predictions). Therefore, any reduction we can achieve in the number of refits essentially implies a proportional reduction in the overall time required for cross-validation of Bayesian models. We will come back to the issue of setting appropriate thresholds in Section 3.

3. Simulations

To evaluate the quality of the PSIS-LFO-CV approximation, we performed a simulation study in which the following conditions were systematically varied:

• The number M of future observations to be predicted took on values of M = 1 and M = 4.

• The threshold τ of the Pareto k estimates was varied between k = 0.5 to k = 0.7 in steps of 0.1.

• Six different data generating models were evaluated, with linear and/or quadratic terms and/or autoregressive terms of order 2 (see Figure 2 for an illustration).

Approximate leave-future-out cross-validation for Bayesian time series models

All authors

Paul-Christian Bürkner , Jonah Gabry & Aki Vehtari

https://doi.org/10.1080/00949655.2020.1783262

Published online:

25 June 2020

Figure 2. Illustration of the models used in the simulations. Black points are observed data. The blue line represents posterior predictions of the model resembling the true data-generating process with 90% prediction intervals shown in grey. More details are provided in the text.

Display full size

In all cases, the time series consisted of N = 200 observations and the minimal number of observations required before make predictions was set to L = 25. We ran 100 simulation trials per condition.

Autoregressive (AR) models are some of the most commonly used time series models. An AR(p) model – an autoregressive model of order p – can be defined as (11) $y_i={\eta }_i+{\epsilon }_i\mathrm{with}{\epsilon }_i=\sum _{k=1}^p{\phi }_k{\epsilon }_{i-k}+e_i,$(11) where ${\eta }_i$ is the linear predictor for the ith observation, ${\phi }_k$ are the autoregressive parameters on the residuals ${\epsilon }_i$, and $e_i$ are pairwise independent errors, which are usually assumed to be normally distributed with equal variance ${\sigma }^2$. The model implies a recursive formula that allows for computing the right-hand side of the equation for observation i based on the values of the equations computed for previous observations. Observations from an AR process are therefore not conditionally independent by definition, but the likelihood still factorizes because we can write down a separate contribution for each observation [see 19 for more discussion on factorizability of statistical models].

In the quadratic model with AR(2) effects (the most complex model in our simulations), the true data generating process was defined as (12) $y_i=b_0+b_{1}t+b_2{t}^2+{\epsilon }_i\mathrm{with}{\epsilon }_i={\phi }_1{\epsilon }_{i-1}+{\phi }_2{\epsilon }_{i-2}+e_i,$(12) where t is the time variable scaled to the unit interval, that is, t = 0 for the smallest time point (1 in our simulations) and t = 1 for the largest time point (200 in our simulations). Specifically, we set the true regression coefficients to the values of $b_0=0$, $b_1=17$, $b_2=25$, and the true autocorrelation parameters to ${\phi }_1=0.5$, and ${\phi }_2=0.3$ (see Figure 2 for an illustration). The choices of the regression coefficients were done so that neither the linear nor quadratic term dominates the other within the specified time frame. The values of the autocorrelation parameters were set to represent typical positively autocorrelated data. In the simulation conditions without linear and/or quadratic and/or AR(2) terms, the corresponding true parameters were simply fixed to zero. We always fit the true data generating model to the data. This is neither required for the validity of LFO-CV in general nor for the validity of the comparison between exact and approximate versions but simply a choice of convenience. For example, a linear model without autocorrelation is used when all but $b_0$ and $b_1$ were set to zero in the simulations.

In addition to exact and approximate LFO-CV, we also computed approximate LOO-CV for comparison. This is not because we think LOO-CV is a generally appropriate approach for time series models, but because, in the absence of any approximate LFO-CV method, researchers may have used approximate LOO-CV for time series models in the past simply because it was the only available option. Demonstrating that LOO-CV is a biased estimate of LFO-CV underscores the importance of developing methods better suited for the task.

All simulations were done in R [23] using the brms package [24,25] together with the probabilistic programming language Stan [11,26] for model fitting, the loo R package [20] for the PSIS computations, and several tidyverse R packages [27] for data processing. The full code and all results are available on Github (https://github.com/paul-buerkner/LFO-CV-paper).

3.1. Results

In this section we focus on the ELPD as a measure out-of-sample predictive performance for reasons outlined in Section 2. In Appendix C, we provide additional simulation results for the RMSE.

Results of the 1-SAP simulations are visualized in Figure 3. Comparing the columns of Figure 3, it is clearly visible that the accuracy of the PSIS approximation is independent of the threshold τ when τ is within the interval $[0.5,0.7]$ motivated in 2.1.1 [9, this would not be the case if τ was allowed to be larger;]. For all conditions, the PSIS-LFO-CV approximation is highly accurate, that is, both unbiased and low in variance around the corresponding exact LFO-CV value (represented by the dashed line in Figure 3). The proportion of observations at which refitting the model was required did not exceed $3\mathrm{\%}$ under any of the conditions and only increased minimally when decreasing τ (see Table 1). At least for the models investigated in our simulations, using $\tau =0.7$ seems to be sufficient for achieving high accuracy and as such there is no need to lower the threshold below that value. As expected, LOO-CV (the lighter histograms in Figure 3) is a biased estimate of the 1-SAP performance for all non-constant models, in particular for models with a trend in the time series. More precisely, LOO-CV is positively biased, which implies that it systematically overestimates M-SAP performance of time series models.

Approximate leave-future-out cross-validation for Bayesian time series models

All authors

Paul-Christian Bürkner , Jonah Gabry & Aki Vehtari

https://doi.org/10.1080/00949655.2020.1783262

Published online:

25 June 2020

Figure 3. Simulation results of 1-step-ahead predictions. Histograms are based on 100 simulation trials of time series with N = 200 observations requiring at least L = 25 observations to make predictions. The black dashed lines indicate the exact LFO-CV result.

Display full size

Figure 3. Simulation results of 1-step-ahead predictions. Histograms are based on 100 simulation trials of time series with N = 200 observations requiring at least L = 25 observations to make predictions. The black dashed lines indicate the exact LFO-CV result.

Approximate leave-future-out cross-validation for Bayesian time series models

All authors

Paul-Christian Bürkner , Jonah Gabry & Aki Vehtari

https://doi.org/10.1080/00949655.2020.1783262

Published online:

25 June 2020

Table 1. Mean proportions of required refits for PSIS-LFO-CV.

CSVDisplay Table

Results of the 4-SAP simulations are visualized in Figure 4. Comparing the columns of Figure 4, it is again clearly visible that the accuracy of the PSIS approximation is independent of the threshold τ. The proportion of observations at which refitting the model was required did not exceed $3\mathrm{\%}$ under any condition and only increased minimally when decreasing τ (see Table 1). In light of the corresponding 1-SAP results presented above, this is not surprising because the procedure for determining the necessity of a refit is independent of M (see Section 2.1). PSIS-LOO-CV is not displayed in Figure 4 as the number of observations predicted at each step (4 vs. 1) makes 4-SAP LFO-CV and LOO-CV incomparable.

Approximate leave-future-out cross-validation for Bayesian time series models

All authors

Paul-Christian Bürkner , Jonah Gabry & Aki Vehtari

https://doi.org/10.1080/00949655.2020.1783262

Published online:

25 June 2020

Figure 4. Simulation results of 4-step-ahead predictions. Histograms are based on 100 simulation trials of time series with N = 200 observations requiring at least L = 25 observations to make predictions. The black dashed lines indicate the exact LFO-CV result.

Display full size

Figure 4. Simulation results of 4-step-ahead predictions. Histograms are based on 100 simulation trials of time series with N = 200 observations requiring at least L = 25 observations to make predictions. The black dashed lines indicate the exact LFO-CV result.

4. Case studies

4.1. Annual measurements of the level of lake Huron

To illustrate the application of PSIS-LFO-CV for estimating expected M-SAP performance, we will fit a model for 98 annual measurements of the water level (in feet) of https://en.wikipedia.org/wiki/Lake_HuronLake Huron from the years 1875–1972. This data set is found in the datasets R package, which is installed automatically with R [23]. The time series shows rather strong autocorrelation and some downward trend towards lower water levels for later points in time. Figure 5 shows the observed time series of water levels as well as predictions from a fitted AR(4) model.

Approximate leave-future-out cross-validation for Bayesian time series models

All authors

Paul-Christian Bürkner , Jonah Gabry & Aki Vehtari

https://doi.org/10.1080/00949655.2020.1783262

Published online:

25 June 2020

Figure 5. Water Level in Lake Huron (1875–1972). Black points are observed data. The blue line represents mean predictions of an AR(4) model with 90% prediction intervals shown in grey.

Display full size

Figure 5. Water Level in Lake Huron (1875–1972). Black points are observed data. The blue line represents mean predictions of an AR(4) model with 90% prediction intervals shown in grey.

Based on this data and model, we will illustrate the use of PSIS-LFO-CV to provide estimates of 1-SAP and 4-SAP when leaving out all future values. To allow for reasonable predictions, we will require at least $L=20$ historical observations (20 years) to make predictions. Further, we set a threshold of $\tau =$ 0.7 for the Pareto k estimates that indicate when refitting becomes necessary. Our fully reproducible analysis of this case study can be found on GitHub (https://github.com/paul-buerkner/LFO-CV-paper).

We start by computing exact and PSIS-approximated LFO-CV of 1-SAP. The computed ELPD values are ${\mathrm{ELPD}}_{\mathrm{exact}}=$ −93.48 and ${\mathrm{ELPD}}_{\mathrm{approx}}=$ −93.62, which are almost identical. Not only is the overall ELPD estimated accurately but so are all of the pointwise ELPD contributions (see the left panel of Figure 6). In comparison, PSIS-LOO-CV returns ${\mathrm{ELPD}}_{\mathrm{loo}}=$ −88.9, overestimating the predictive performance and as suggested by our simulation results for stationary autoregressive models (see fourth row of Figure 3). Plotting the Pareto k estimates reveals that the model had to be refit 3 times, out of a total of N-$L=$ 78 predicted observations (see Figure 7). On average, this means one refit every 26.0 observations, which implies a drastic speed increase compared to exact LFO-CV.

Approximate leave-future-out cross-validation for Bayesian time series models

All authors

Paul-Christian Bürkner , Jonah Gabry & Aki Vehtari

https://doi.org/10.1080/00949655.2020.1783262

Published online:

25 June 2020

Figure 6. Pointwise exact vs. PSIS-approximated ELPD contributions for 1-SAP (left) and 4-SAP (right) for the Lake Huron model. A threshold of $\tau =0.7$ was used for the Pareto k estimates. M is the number of predicted future observations.

Display full size

Figure 6. Pointwise exact vs. PSIS-approximated ELPD contributions for 1-SAP (left) and 4-SAP (right) for the Lake Huron model. A threshold of $\tau =0.7$ was used for the Pareto k estimates. M is the number of predicted future observations.

Approximate leave-future-out cross-validation for Bayesian time series models

All authors

Paul-Christian Bürkner , Jonah Gabry & Aki Vehtari

https://doi.org/10.1080/00949655.2020.1783262

Published online:

25 June 2020

Figure 7. Pareto k estimates for PSIS-LFO-CV of the Lake Huron model. The dotted red line indicates the threshold at which the refitting was necessary.

Display full size

Figure 7. Pareto k estimates for PSIS-LFO-CV of the Lake Huron model. The dotted red line indicates the threshold at which the refitting was necessary.

Performing LFO-CV for 4-SAP, we obtained ${\mathrm{ELPD}}_{\mathrm{exact}}=$ −411.41 and ${\mathrm{ELPD}}_{\mathrm{approx}}=$ −412.78, which are again very similar. In general, as M increases, the approximation will tend to become more variable around the true value in absolute ELPD units because the ELPD increment of each observation will be based on more and more observations (see also Section 3). For this example, we see some considerable differences in the pointwise ELPD contributions of specific observations which were hard to predict accurately by the model (see the right panel of Figure 6). This is to be expected because predicting M steps ahead using an AR model will yield highly uncertain predictions if most of the autocorrelation happens at lags smaller than M (see also the bottom rows in Figure 4). For such a model, it may be ill-advised to evaluate predictions too far into the future, at least when using the approximate methods presented in this paper. Since, for a constant threshold τ, the importance weights are the same independent of M, the Pareto k estimates are the same for 4-SAP and 1-SAP.

4.2. Annual date of the cherry blossoms in Japan

The cherry blossom in Japan is a famous natural phenomenon occurring once every year during spring. As the climate changes so does the annual date of the cherry blossom [28,29]. The most complete reconstruction available to date contains data between 801 AD and 2015 AD [28,29] and is available online (http://atmenv.envi.osakafu-u.ac.jp/aono/kyophenotemp4/).

In this case study, we will predict the annual date of the cherry blossom using an approximate Gaussian process model [30,31] to provide flexible non-linear smoothing of the time series. A visualization of both the data and the fitted model is provided in Figure 8. While the time series appears rather stable across earlier centuries, with substantial variation across consecutive years, there are some clearly visible trends in the data. Particularly in more recent years, the cherry blossom has tended to happen much earlier than before, which may be a consequence of changes in the climate [28,29].

Approximate leave-future-out cross-validation for Bayesian time series models

All authors

Paul-Christian Bürkner , Jonah Gabry & Aki Vehtari

https://doi.org/10.1080/00949655.2020.1783262

Published online:

25 June 2020

Figure 8. Day of the cherry blossom in Japan (812–2015). Black points are observed data. The blue line represents mean predictions of a thin-plate spline model with 90% regression intervals shown in grey.

Display full size

Figure 8. Day of the cherry blossom in Japan (812–2015). Black points are observed data. The blue line represents mean predictions of a thin-plate spline model with 90% regression intervals shown in grey.

Based on this data and model, we will illustrate the use of PSIS-LFO-CV to provide estimates of 1-SAP and 4-SAP leaving out all future values. To allow for reasonable predictions of future values, we will require at least L = 100 historical observations (100 years) to make predictions. Further, we set a threshold of $\tau =$ 0.7 for the Pareto k estimates to determine when refitting becomes necessary. Our fully reproducible analysis of this case study can be found on GitHub (https://github.com/paul-buerkner/LFO-CV-paper).

We start by computing exact and PSIS-approximated LFO-CV of 1-SAP. We compute ${\mathrm{ELPD}}_{\mathrm{exact}}=$ −2345.7 and ${\mathrm{ELPD}}_{\mathrm{approx}}=$ −2344.9, which are highly similar. As shown in the left panel of Figure 9, the pointwise ELPD contributions are highly accurate, with no outliers. The approximation has worked well for all observations. PSIS-LFO-CV performs much better than PSIS-LOO-CV (${\mathrm{ELPD}}_{\mathrm{approx}}=$ −2340.3), which overestimates the predictive performance. Plotting the Pareto k estimates reveals that the model had to be refit 6 times, out of a total of N-$L=$ 727 predicted observations (see Figure 10). On average, this means one refit every 121.2 observations, which implies a drastic speed increase as compared to exact LFO-CV.

Approximate leave-future-out cross-validation for Bayesian time series models

All authors

Paul-Christian Bürkner , Jonah Gabry & Aki Vehtari

https://doi.org/10.1080/00949655.2020.1783262

Published online:

25 June 2020

Figure 9. Pointwise exact vs. PSIS-approximated ELPD contributions of 1-SAP (left) and 4-SAP (right) for the cherry blossom model. A threshold of $\tau =0.7$ was used for the Pareto k estimates. M is the number of predicted future observations.

Display full size

Figure 9. Pointwise exact vs. PSIS-approximated ELPD contributions of 1-SAP (left) and 4-SAP (right) for the cherry blossom model. A threshold of $\tau =0.7$ was used for the Pareto k estimates. M is the number of predicted future observations.

Approximate leave-future-out cross-validation for Bayesian time series models

All authors

Paul-Christian Bürkner , Jonah Gabry & Aki Vehtari

https://doi.org/10.1080/00949655.2020.1783262

Published online:

25 June 2020

Figure 10. Pareto k estimates for PSIS-LFO-CV of the cherry blossom model. The dotted red line indicates the threshold at which the refitting was necessary.

Display full size

Figure 10. Pareto k estimates for PSIS-LFO-CV of the cherry blossom model. The dotted red line indicates the threshold at which the refitting was necessary.

Performing LFO-CV of 4-SAP, we compute ${\mathrm{ELPD}}_{\mathrm{exact}}=$ −9348.3 and ${\mathrm{ELPD}}_{\mathrm{approx}}=$ −9345.5, which are again similar but not as close as the corresponding 1-SAP results. This is to be expected as the uncertainty of PSIS-LFO-CV increases for increasing M (see Section 3). As displayed in the right panel of Figure 9, the pointwise ELPD contributions are highly accurate in most cases, with a few small outliers in both directions. For constant threshold τ, the importance weights are the same independent of M, so the Pareto k estimates are the same for 4-SAP and 1-SAP.

5. Conclusion

We proposed, evaluated and demonstrated PSIS-LFO-CV, a method for approximating cross-validation of Bayesian time series models. PSIS-LFO-CV is intended to be used when the prediction task is predicting future values based solely on past values, in which case leave-one-out cross-validation is inappropriate. Within the set of such prediction tasks, we can choose the number M of future observations to be predicted. For a set of common time series models, we established via simulations that PSIS-LFO-CV is an unbiased approximation of exact LFO-CV if we choose the threshold τ of the Pareto k estimates to not be larger than $\tau =0.7$. That is, PSIS-LFO-CV does not require a smaller (stricter) threshold than PSIS-LOO-CV to achieve satisfactory accuracy.

By nature of the approximated M-step-ahead predictions, the computation time of PSIS-LFO-CV still increases linearly with the number of observations N. However, in our numerical experiments, we were able to reduce to computation time by a factor of roughly 25–100 compared to exact LFO-CV, which is enough to make LFO-CV realistic for many applications.

A limitation of our current approach is that the uncertainty in the approximate LFO-CV estimates is hard to quantify. There are at least three types of uncertainty which could be considered here. First, there is uncertainty induced by approximating exact LFO-CV using (Pareto smoothed) importance sampling. Based on theoretical considerations of the approximation and numerical experiments both presented in Vehtari et al. [9], any PSIS approximation will be very close to the exact value as long as the Pareto k diagnostic does not exceed the threshold of 0.7, which we used as the refit criterion in our approximate LFO-CV approach. Second, there is uncertainty caused by finite amounts of data. For 1-step-ahead predictions, we can use an analogous approach to what is done in approximate LOO-CV by computing the standard error across the pointwise estimates for each observation [8]. More generally, for M-step-ahead predictions, we can compute the standard error by using every Mth value which are then independent. Third, there is uncertainty induced by the finite number of posterior draws but this uncertainty tends to be negligible with just a few thousand draws compared to the second source of uncertainty [8]. Investigating uncertainty measures for (approximate) LFO-CV in more detail is left for future research.

Lastly, we want to briefly note that LFO-CV can also be used to compute marginal likelihoods. Using basic rules of conditional probability, we can factor the log marginal likelihood as (13) $\log p(y)=\sum _{i=1}^N\log p(y_i|y_{1:(i-1)}).$(13) This is exactly the ELPD of 1-SAP if we set L = 0, that is if we choose to predict all observations using their respective past (the very first observation is only predicted from the prior). As such, marginal likelihoods may be approximated using PSIS-LFO-CV. Although this approach is unlikely to be more efficient than methods specialized for computing marginal likelihoods [32–34, e.g. bridge sampling;], it may be a noteworthy option if for some reason other methods fail.