Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
78,149 characters · 23 sections · 28 citation commands
Equivariant online predictions of non-stationary time series
\global\long\def\spacingset#1{ \global\long \global\long} \spacingset{1}
\if11 {
} \fi
\if01 {
} \fi
Keywords: Bayesian analysis; Exact minimax; Model misspecification; Time series; Ensemble methods.
\spacingset{1}
Prediction and decision making using time series data are central tasks for statistical methods in many fields, including economics, climatology, and epidemiology. For these problems, the goal of the statistical model is to use past data to characterize future data and its uncertainty in an online, sequential manner. Predictions are then used to inform and improve decision making. As these tasks are defined by their temporality-- past events and actions affect future events-- the success of a statistical model depends largely on its ability to capture temporal characteristics. To respond, many models have been developed and proposed, often categorized as dynamic models WestHarrison1997book2,Prado2010. Although the success of these models has been shown in practice, theoretical investigations have been limited, particularly for non-stationary data, despite its relevance in the domain.
There are several reasons why theoretical development in this area has been difficult, although the biggest reason is inherent in the data itself. As the data in time series tasks are, prima facie, non-stationary, defined by the lack of stationary assumptions, most of the theoretical apparatus used in statistics, which rely on assumptions, such as i.i.d., cannot be used. For example, in most time series data, we observe gradual and sudden shifts in trend and volatility as time progresses, and a model or covariate that works at one period often do not in another. This problem is compounded when we consider decisions, as these decisions themselves can affect future data (e.g. implementing strict mobility restrictions under a pandemic will likely affect future infection numbers). To partially get around this issue, many papers transform the data to appear stationary, often by taking the log difference several times. This, however, is often detrimental to the task itself, as these transformations do not guarantee that the data are actually stationary and, more importantly, they remove the signals and patterns that are critical to understand the data and make better predictions and decisions.
In developing theory for time series models, three points must be considered for it to be relevant to real world problems. First, given the non-stationary nature of the data, it cannot rely on asymptotics, nor the conditions required for it. This is not only required by the temporal characteristics of the data, but also reflects the need of the decision maker, since forecasting and decision making are always done with finite sample data and for finite horizons. For example, many economic data are measured at the monthly frequency, with only several decades in the past being relevant, at most (and much less so under economic shocks). This would not be enough data for the asymptotic results to be meaningful, even if the data can be assumed to be i.i.d., which they cannot. Second, since decision making is done under uncertainty, the theoretical results must be with regard to the whole predictive distribution, and not just the mean. Although the predictive literature has often focused on point metrics, such as mean squared forecast error, this is often insufficient for decision making. For example, in finance, investment portfolios are constructed by taking into account both the mean and covariance of the assets, and it is widely recognized that the latter quantification is critical for successful asset management. Third, all models must be assumed to be misspecified, a setting often referred to as an $\mathcal{M}$-open bernardo2009bayesian. This point reflects the reality that the “best" model at some point can change given time or even under some shock. It is simply unrealistic to assume that there exists a true model in time series contexts. Given all three points, thus, theoretical evaluations of statistical models for non-stationary time series must be done in finite sample, with regard to its distributional predictions, and where a true model cannot be assumed to exist, for it to be meaningful and provide insight for practice.
We contribute to this field by developing a theoretical framework that satisfies all three criteria. Specifically, we have three contributions. First, we define the Kullback-Leibler risk for predictions in non-stationary time series data (Section (ref)). This allows us to conduct theoretical analyses of statistical models for non-stationary time series data within a decision theoretic framework. Critical to this framework, we neither assume that the true model is nested nor any asymptotics. Second, using this decision theoretic framework, we show that a specific class of dynamic models produces exact minimax predictive distributions, under Gaussian assumptions (Section (ref)). This result provides finite sample predictive guarantees, making it a benchmark to compare other models against. We then generalize this result by relaxing the Gaussian assumption using semi-martingale processes, providing theoretical guarantees under more general settings (Section (ref)). Third, we extend the above results to the problem of combining several predictive distributions (i.e. ensemble methods), and show that dynamic Bayesian predictive synthesis mcalinn2019dynamic is exact minimax (Section (ref)). We highlight the theoretical results using three topical datasets: weekly average COVID-19 cases in Tokyo, monthly global mean sea level, and a high-dimensional monthly economic dataset (Section (ref)). Through all three applications (and simulation results in Appendix (ref)), we show that a method that achieves exact minimaxity is superior to methods that do not.
For the theoretical analysis, we first assume that the data generating process of both the target and covariates follow Gaussian processes. While this assumption is somewhat restrictive-- even if it is reasonable enough for most applications--, this is primarily done for ease of exposition, and to clarify the assumptions we make. This assumption of Gaussianity is relaxed and generalized using semi-martingale processes in Section (ref).
Let the $\mathbb{R}$-valued process to be predicted, $\left\{ y_{t}\right\} $, and the covariate $\mathbb{R}$-valued processes (which can include the lag of $y$), $\left\{ x_{j,t}\right\} $, be Gaussian processes. Thus, any finite linear combination of elements from a set of random variables, $\left\{ y_{s},x_{j,s};s=1,\cdots t,\ j\in J\right\} $, is a Gaussian random variable. Let the mean and variance of $\left\{ y_{t}\right\} $ be $\mu_{t}^{y}$ and $\sigma_{t}^{y}$, respectively, and (the elements of) the mean vector and covariance matrix of $\left\{ x_{j,t}\right\} $ be \[ \mu_{t}^{j}=\mathbb{E}\left[x_{j,t}\right],\ \sigma_{t}^{ij}=\mathbb{E}\left[\left(x_{j,t}-\mu_{t}^{j}\right)\left(x_{i,t}-\mu_{t}^{i}\right)\right]. \]
Let the $J+1$-dimensional multivariate normal distribution of $\left(y_{t},\boldsymbol{x}_{t}\right)$ be \[ \left[
\right]\sim N\left(\left[
\right],\left[
\right]\right), \] where $\boldsymbol{x}_{t}=\left[x_{1,t},\cdots,x_{J,t}\right]^{\top}$ $\boldsymbol{\mu}_{t}=\left[\mu_{t}^{1},\cdots,\mu_{t}^{J}\right]^{\top}$, $\boldsymbol{\sigma}_{t}^{yJ}=\left[\sigma_{t}^{y1},\cdots,\sigma_{t}^{yJ}\right]^{\top}$, $\sigma_{t}^{yj}=\mathbb{E}\left[\left(y_{t}-\mu_{t}^{y}\right)\left(x_{j,t}-\mu_{t}^{j}\right)\right]$, and $\varSigma_{t}^{J}$, which is a $J\times J$ matrix with $\sigma_{t}^{ij}$ as the $i,j$-th element. Then, the conditional distribution, $p_{t}\left(y_{t}\left|\boldsymbol{x}_{t}\right.\right)$, is, \[ p_{t}\left(y_{t}\left|\boldsymbol{x}_{t}\right.\right)=N\left(\mu_{t}+\left(\boldsymbol{\sigma}_{t}^{yJ}\right)^{\top}\left(\varSigma_{t}^{J}\right)^{-1}\left(x_{t}-\mu_{t}\right),\sigma_{t}^{y}-\left(\boldsymbol{\sigma}_{t}^{yJ}\right)^{\top}\left(\varSigma_{t}^{J}\right)^{-1}\boldsymbol{\sigma}_{t}^{yJ}\right), \] from the normal correlation theorem. If we write, \[ \mu_{t}+\left(\boldsymbol{\sigma}_{t}^{yJ}\right)^{\top}\left(\varSigma_{t}^{J}\right)^{-1}\left(\boldsymbol{x}_{t}-\mu_{t}\right)=\theta_{0,t+1}^{*}+\left\langle \boldsymbol{\theta}_{t+1}^{*},\boldsymbol{x}_{t+1}\right\rangle, \] where $\boldsymbol{\theta}_{t+1}^{*}=\left[\theta_{1,t+1}^{*},\cdots,\theta_{J,t+1}^{*}\right]^{\top}$, and $\sigma_{t+1}=\sigma_{t}^{y}-\left(\boldsymbol{\sigma}_{t}^{yJ}\right)^{\top}\left(\varSigma_{t}^{J}\right)^{-1}\boldsymbol{\sigma}_{t}^{yJ}$, we have the likelihood function of the best model,
Given the above formulation, the unknown parameters are $\left(\theta_{0,t+1}^{*},\boldsymbol{\theta}_{t+1}^{*},\sigma_{t+1}\right)$. By defining the best possible prediction, and identifying the unknown parameters, we are able to measure and analyze the predictive risk with regard to its Kullback-Leibler (KL) divergence. The specific definition of the KL loss and risk used in the statistical decision theory is given in Section (ref).
Note that our interest is in the predictive ability of time series models, and not in the geometric form of them. Because our interest is in the 1-step ahead forecasts, the model space that defines the KL risk is solely determined by the $J+2$ parameters, $\left(\theta_{0,t+1}^{*},\boldsymbol{\theta}_{t+1}^{*},\sigma_{t+1}\right)$, and they need not be seen as a function of $t$.
We consider an equivariant model for the decision problem. An equivariant model is defined here as a model that produces predictive distributions that are equivariant for the decision problem concerning the one-step ahead predictive KL loss function. Specifically, the equivariant model regarding the DGP (eq. (ref)) is given as follows (the proof of its equivariance is given in (ref)):
which is a random walk DLM, where $\boldsymbol{\theta}_{t+1}^{\prime}=\left[\theta_{0,t+1},\boldsymbol{\theta}_{t+1}\right]^{\top}$ evolves in time according to a linear, normal random walk with innovations variance matrix, $v_{t+1}\boldsymbol{W}_{t+1}$, at time $t+1$, and $v_{t+1}$ is the residual variance in predicting $y_{t+1}$ based on past information and the set of covariates. The residuals, $\nu_{t+1}$, and evolution innovations, $\boldsymbol{\omega}_{s+1}$, are independent over time and mutually independent for all $t,s.$
Specifically, the likelihood function is
where $\boldsymbol{x}_{t+1}^{\prime}=\left[1,\boldsymbol{x}_{t+1}\right]^{\top}$. The predictive distribution conditioned on the covariates, $\boldsymbol{x}_{t+1}^{\prime}$, are
Here, \[ \pi\left(\boldsymbol{\theta}_{t+1}^{\prime},v_{t+1}\left|\left\{ y_{s},\boldsymbol{x}_{s}\right\} _{s=1}^{t}\right.\right)=\int_{\left(0,\infty\right)}\int_{\mathbb{R}^{J+1}}\pi\left(\boldsymbol{\theta}_{t+1}^{\prime},v_{t+1}\left|\boldsymbol{\theta}_{t}^{\prime},v_{t}\right.\right)\pi\left(\boldsymbol{\theta}_{t}^{\prime},v_{t}\left|\left\{ y_{s},\boldsymbol{x}_{s}\right\} _{s=1}^{t}\right.\right)d\boldsymbol{\theta}_{t}^{\prime}dv_{t}, \] where $\pi\left(\boldsymbol{\theta}_{t}^{\prime},v_{t}\left|\left\{ y_{s},\boldsymbol{x}_{s}\right\} _{s=1}^{t}\right.\right)$ is the posterior distribution of $\left(\boldsymbol{\theta}_{t}^{\prime},v_{t}\right)$ at time $t$, after observing $y_{t}$ and $\boldsymbol{x}_{t+1}$. At $t=0$, we denote the (initial) prior distribution of $\left(\boldsymbol{\theta}_{0}^{\prime},v_{0}\right)$ as $\rho\left(\boldsymbol{\theta}_{0}^{\prime},v_{0}\right)$. Under this (initial) prior, the predictive distribution is denoted as $\hat{p}_{t}^{\rho}\left(y_{t+1}\left|\left\{ y_{s},\boldsymbol{x}_{s}\right\} _{s=1}^{t},\boldsymbol{x}_{t+1}\right.\right)$, and the posterior distribution as $\pi^{\rho}\left(\boldsymbol{\theta}_{t}^{\prime},v_{t}\left|\left\{ y_{s},\boldsymbol{x}_{s}\right\} _{s=1}^{t}\right.\right)$, both at time $t$.
Note that we can correspond eq. (ref) to eq. (ref), as follows:
In the next section, we show that the equivariant model also gives a minimax predictive density.
The goal is to show that the random walk DLM (eq. (ref)), denoted as $\hat{p}_{t}\left(y_{t+1}\left|\left\{ y_{s},\boldsymbol{x}_{s}\right\} _{s=1}^{t},\boldsymbol{x}_{t+1}\right.\right)$, is equivariant and minimax in terms of the KL risk with regard to the transition probability, $p_{t}\left(y_{t+1}\left|\left\{ y_{s},\boldsymbol{x}_{s}\right\} _{s=1}^{t},\theta_{0,t+1}^{*},\boldsymbol{\theta}_{t+1}^{*},\sigma_{t+1},\boldsymbol{x}_{t+1}\right.\right)$. Thus, for each time the covariate processes, $\left\{ \boldsymbol{x}_{t}\right\} $, are given, we can construct the predictive distribution, $\hat{p}_{t}\left(y_{t+1}\left|\left\{ y_{s},\boldsymbol{x}_{s}\right\} _{s=1}^{t},\boldsymbol{x}_{t+1}\right.\right)$, in a path-wise manner. We will show that the random walk DLM is a minimax estimate regarding the parameters, $\left(\theta_{0,t+1}^{*},\boldsymbol{\theta}_{t+1}^{*},\sigma_{t+1}\right)$, of the best model, $p_{t}\left(y_{t+1}\left|\left\{ y_{s},\boldsymbol{x}_{s}\right\} _{s=1}^{t},\theta_{0,t+1}^{*},\boldsymbol{\theta}_{t+1}^{*},\sigma_{t+1},\boldsymbol{x}_{t+1}\right.\right)$.
We first define the invariant decision problem. As the loss function, we employ KL divergence, and define it as follows:
Then, the KL risk is
Here, the probability measure, $\mu_{Y}\left(y_{t},\cdots,y_{0}\right)$, is a cylindrical measure of the process, $\left\{ y_{t}\right\} $. The sample space is the past observations, $\left\{ y_{s}\right\} _{s=1}^{t}$. As a transformation to the sample space, $\left\{ y_{s}\right\} _{s=1}^{t}$, we add the shift transform value, $-\left(\frac{1}{\sigma_{t+1}}\theta_{0,t+1}^{*}+\left\langle \frac{1}{\sigma_{t+1}}\boldsymbol{\theta}_{t+1}^{*},\boldsymbol{x}_{t+1}\right\rangle \right)$, and multiply each scale by $\frac{1}{\sigma_{t+1}}$, $\left(\frac{1}{\sigma_{t+1}}>0\right)$, to obtain,
This scale-shift transformation can transform the sample to an arbitrary value on $\mathbb{R}$. The parameter space is $\left(\theta_{0,t+1},\boldsymbol{\theta}_{t+1},v_{t+1}\right)\textrm{ or }\left(\theta_{0,t+1}^{*},\boldsymbol{\theta}_{t+1}^{*},\sigma_{t+1}\right)$. As with the sample, to transform the parameter, we scale-shift transform the parameters in the DLM model as,
or the parameters in the DGP as,
The decision space is the probability density function, $p\left(y_{t+1}\right)$, of $y_{t+1}$. For each time $t+1$, we consider the scale-shift transform, \[ y_{t+1}\rightarrow\tilde{y}_{t+1}=\frac{1}{\sigma_{t+1}}y_{t+1}-\left(\frac{1}{\sigma_{t+1}}\theta_{0,t+1}^{*}+\left\langle \frac{1}{\sigma_{t+1}}\boldsymbol{\theta}_{t+1}^{*},\boldsymbol{x}_{t+1}\right\rangle \right), \] for the target variable, $y_{t+1}$. The transformation of the probability density is denoted as \[ p\left(y_{t+1}\right)\rightarrow p\left(\tilde{y}_{t+1}\right). \]
Further, since \[ \widetilde{y}_{t+1}=\widetilde{\theta}_{0,t+1}^{*}+\left\langle \widetilde{\boldsymbol{\theta}}_{t+1}^{*},\boldsymbol{x}_{t+1}\right\rangle +\frac{1}{\sigma_{t+1}}\varepsilon_{t+1},\ \frac{1}{\sigma_{t+1}}\varepsilon_{t+1}\sim N\left(0,\frac{v_{t+1}}{\sigma_{t+1}^{2}}\right), \] we have, \[ p_{t}\left(y_{t+1}\left|\left\{ y_{s},\boldsymbol{x}_{s}\right\} _{s=1}^{t},\tilde{\theta}_{0,t+1}^{*},\tilde{\boldsymbol{\theta}}_{t+1}^{*},\tilde{\sigma}_{t+1},\boldsymbol{x}_{t+1}\right.\right)=p_{t}\left(\tilde{y}_{t+1}\left|\left\{ y_{s},\boldsymbol{x}_{s}\right\} _{s=1}^{t},\theta_{0,t+1}^{*},\boldsymbol{\theta}_{t+1}^{*},\sigma_{t+1},\boldsymbol{x}_{t+1}\right.\right). \] Then, the KL, as a loss function, satisfies invariance:
The data transformation, parameter transformation, and loss invariance define the invariant decision problem.
The predictive distribution, $\hat{p}_{t}$, that achieves \[ \hat{p}_{t}\left(\tilde{y}_{t+1}\left|\left\{ \tilde{y}_{s},\boldsymbol{x}_{s}\right\} _{s=1}^{t},\boldsymbol{x}_{t+1}\right.\right)=\hat{p}_{t}\left(y_{t+1}\left|\left\{ y_{s},\boldsymbol{x}_{s}\right\} _{s=1}^{t},\boldsymbol{x}_{t+1}\right.\right) \] is called an equivariant predictive distribution, and is minimax in terms of the KL risk, which minimizes the max risk: \[ \max_{\theta_{0,t+1}^{*},\theta_{t+1}^{*},\sigma_{t+1}}R_{\mathsf{KL}}\left(\left(\theta_{0,t+1}^{*},\boldsymbol{\theta}_{t+1}^{*},\sigma_{t+1}\right),\hat{p}_{t}\right),\ \left(\textrm{for each }t\right). \] This means that the sequentially updated Bayes predictive distribution, $\hat{p}_{t}$, is sequentially minimax. We will now show that the predictive distribution of eq. (ref) is equivariant and minimax.
The choice of criteria-- finite sample predictive minimax risk-- is done to reflect the need and reality of predictive tasks in non-stationary, time series data. For example, considering the risk of the estimated parameters, and not the predictive risk, is infeasible because the parameters dynamically change as the data generating process itself changes. Assuming strong conditions on the data generating process, such that it is equivalent to i.i.d. data, defeats the purpose and uninformative for the applications of interest. Further, any asymptotic analysis will require similarly strong assumptions, even if the interest is in predictive risk, making the analysis irrelevant. Evaluating the predictive minimax risk against the best predictive density allows us to assess the predictive performance in finite sample, making it possible to use it for non-stationary time series data and be relevant to the task.
As we will see later in Lemma (ref), if $\left\{ \theta_{0,t},\boldsymbol{\theta}_{t},v_{t}\right\} $ is updated via a random walk, the transition distribution of the scale-shift transformed $\left\{ \theta_{0,t},\boldsymbol{\theta}_{t},v_{t}\right\} $ is invariant. Thus, the statistical decision theoretic problem is invariant with regard to the scale-shift transform. The predictive distribution is given as
Here, $\left\{ \widetilde{\boldsymbol{\theta}}_{t-1},\cdots,\widetilde{\boldsymbol{\theta}}_{1}\right\} $ is $\left\{ \boldsymbol{\theta}_{t-1},\cdots,\boldsymbol{\theta}_{1}\right\} $ transformed by adding $-\frac{1}{\sigma_{t+1}}\boldsymbol{\theta}_{t}^{*}$ and multiplying $\frac{1}{\sigma_{t+1}}$.
Here, the prior distribution of $\boldsymbol{\theta}_{t}^{\prime}$, at $t=0$, is a Lebesgue measure, $\mathsf{m}\left(\cdot\right)$, the prior distribution of $v_{t}$ at $t=0$ is a scale-invariant prior, $\frac{1}{v_{0}}$, and the state space evolves following a random walk (eqs. (ref) and (ref)). Then, the following lemma holds.
Since $\hat{p}_{t}\left(\widetilde{y}_{t+1}\left|\widetilde{\boldsymbol{\theta}^{\prime}}_{t+1},\widetilde{v}_{t+1},\boldsymbol{x}_{t+1}\right.\right)$ is given from eq. (ref), we have \[ \hat{p}_{t}\left(\widetilde{y}_{t+1}\left|\widetilde{\boldsymbol{\theta}^{\prime}}_{t+1},\widetilde{v}_{t+1},\boldsymbol{x}_{t+1}\right.\right)=\hat{p}_{t}\left(y_{t+1}\left|\boldsymbol{\theta}_{t+1}^{\prime},v_{t+1},\boldsymbol{x}_{t+1}\right.\right). \] Therefore, with Lemma (ref), we have \[ \pi\left(\widetilde{\boldsymbol{\theta}^{\prime}}_{t+1},\widetilde{v}_{t+1}\left|\left\{ \widetilde{y}_{s},\boldsymbol{x}_{s}\right\} _{s=1}^{t}\right.\right)=\pi\left(\boldsymbol{\theta}_{t+1}^{\prime},v_{t+1}\left|\left\{ y_{s},\boldsymbol{x}_{s}\right\} _{s=1}^{t}\right.\right) \] and, as a result, can see that this is equivariant: \[ \hat{p}_{t}\left(\widetilde{y}_{t+1}\left|\left\{ \widetilde{y}_{s},\boldsymbol{x}_{s}\right\} _{s=1}^{t},\boldsymbol{x}_{t+1}\right.\right)=\hat{p}_{t}\left(y_{t+1}\left|\left\{ y_{s},\boldsymbol{x}_{s}\right\} _{s=1}^{t},\boldsymbol{x}_{t+1}\right.\right). \]
Next, we will show that the equivariance predictive distribution, $\hat{p}_{t}\left(y_{t+1}\left|\left\{ y_{s},\boldsymbol{x}_{s}\right\} _{s=1}^{t},\boldsymbol{x}_{t+1}\right.\right)$, is minimax in finite sample.
Given the KL risk, we have the following theorem:
The full proof is given in Appendix (ref). The general strategy of the proof for Theorem. (ref) is done by proving the following two conditions, according to the equalizer rule berger1985statistical, following Corollary 1. of george2006improved:
The Bayes risk is defined as \[ B\left(\rho_{k},\hat{p}_{t}^{\rho_{k}}\right)=\int R_{\mathsf{KL}}\left(\left(\theta_{0,t+1}^{*},\boldsymbol{\theta}_{t+1}^{*},\sigma_{t+1}\right),\hat{p}_{t}^{\rho_{k}}\right)\rho_{k}\left(\theta_{0,t+1}^{*},\boldsymbol{\theta}_{t+1}^{*},\sigma_{t+1}\right), \] which is written with regard to the posterior parameters, $\left(\theta_{0,t+1}^{*},\boldsymbol{\theta}_{t+1}^{*},\sigma_{t+1}\right)$, and not the prior on $\boldsymbol{\theta}_{0}^{\prime},\sigma_{0}$. However, this Bayes risk is effectively equivalent to the prior specification of $\boldsymbol{\theta}_{0}^{\prime},\sigma_{0}$, due to recursive updating. More specifically, the posterior distribution, $\pi\left(\boldsymbol{\theta}_{t}^{\prime},\sigma_{t}\left|\left\{ y_{s},\boldsymbol{x}_{s}\right\} _{s=1}^{t}\right.\right)$, is given by the generalized Bayes formula for filtering (see Appendix C for notation), \[ \pi\left(\boldsymbol{\theta}_{t}^{\prime},\sigma_{t}\left|\left\{ y_{s},\boldsymbol{x}_{s}\right\} _{s=1}^{t}\right.\right)=\frac{\int_{\left(\mathbb{R}^{J}\right)^{t-1}}\beta_{t}\left(\boldsymbol{\theta}^{\prime},\boldsymbol{\sigma},y\right)d\mu_{\theta,\sigma}^{\rho}\left(\left(\boldsymbol{\theta}_{0}^{\prime},\sigma_{0}\right),\cdots,\left(\boldsymbol{\theta}_{t-1}^{\prime},\sigma_{t-1}\right)\right)}{\int_{\left(\mathbb{R}^{J}\right)^{t}}\beta_{t}\left(\boldsymbol{\theta}^{\prime},\boldsymbol{\sigma},y\right)d\mu_{\theta,\sigma}^{\rho}\left(\left(\boldsymbol{\theta}_{0}^{\prime},\sigma_{0}\right),\cdots,\left(\boldsymbol{\theta}_{t}^{\prime},\sigma_{t}\right)\right)}, \] though because we assume a random walk (eq. (ref)) for the process, $\left\{ \boldsymbol{\theta}_{t}^{\prime},\sigma_{t}\right\} $, the cylindrical measure, $\mu_{\theta,\sigma}^{\rho}\left(\left(\boldsymbol{\theta}_{0}^{\prime},\sigma_{0}\right),\cdots,\left(\boldsymbol{\theta}_{t}^{\prime},\sigma_{t}\right)\right)$, is determined by the prior distribution of $\boldsymbol{\theta}_{0}^{\prime},\sigma_{0}$: $\rho\left(\boldsymbol{\theta}_{0}^{\prime}\right),\rho\left(\sigma_{0}\right)$.
Since the predictive distribution of the random walk DLM, $\hat{p}_{t}^{\frac{1}{v_{0}}\mathsf{m}}$, is risk constant (Appendix (ref)) and extended Bayes (Appendix (ref)), from Liang-Barron_04 Theorem 1., it is exact minimax with regard to the KL risk. This result can be interpreted as follows: under the least favorable prior distribution, the Bayes solution is minimax. However, if the state stochastic process, $\left\{ \boldsymbol{\theta}_{t}^{\prime}\right\} $, is stationary, even if the initial prior distribution is a Lebesgue measure, the updated posterior distribution converges in law to a proper probability measure, due to the individual ergodic theorem. Under this setting, the predictive distribution cannot be minimax, since it is not shift invariant. Therefore, for minimaxity, the state stochastic process, $\left\{ \boldsymbol{\theta}_{t}^{\prime}\right\} $, must not be stationary.
Minimaxity is useful as a theoretical benchmark because it examines the predictive under no prior information, a situation that both Bayesians and frequentists can agree to be relevant. In a situation where someone has prior information that informs their model or prior, that model or prior should, at the very least, be equal to or improve over the minimax predictive distribution. This is why, for non-stationary time series data, the minimax model is the theoretical benchmark for predictive performance. In other words, a statistical model should at least achieve minimaxity for it to be considered worthwhile, particularly when a minimax model is known to exist.
While the theorem in Section (ref) assumes that the data generating process and covariates follow Gaussian processes, this assumption can be relaxed using semi-martingale processes.
We assume that the predictive process, $y_{t}$, and the covariate process, $x_{j,t}$, $j\in J$, are both square integralable processes and filtration $\mathcal{F}_{t}$-adapted ($\mathcal{F}_{t}=\sigma\left(\left\{ y_{s},\boldsymbol{x}_{s}\right\} _{s=1}^{t}\right)$). Then, both $y_{t}$, $x_{j,t}$, $j\in J$ are uniquely decomposed, from the Doob decomposition, as the following:
where $y_{0},$ $\left\{ x_{j,0}\right\} _{j\in J}$ is a $\mathcal{F}_{0}$-measurable stochastic variable, $M,$ $\left\{ M_{j}\right\} _{j\in J}$ is a square-integrable martingale, and $V,$ $\left\{ V_{j}\right\} _{j\in J}$ is square-integrable and predictable, thus $\mathcal{F}_{t-1}$-adapted. Further, the optimal approximate process of $y_{t+1}$, utilizing the covariate process, $x_{j,t+1},j\in J$, can be represented as schweizer1995minimal:
where $\left(\theta_{y_{t+1}}^{*},\left\{ \boldsymbol{\theta}_{s}^{*}\right\} _{1\leqq s\leqq t+1}\right)$ are the true parameters of the FS decomposition (due to $\boldsymbol{\theta}_{t+1}^{*}$ being determined at $t$) and $\varDelta\boldsymbol{x}_{s}=\boldsymbol{x}_{s+1}-\boldsymbol{x}_{s}$. When $\left\{ \theta_{y_{t+1}}^{*}\right\} $ is constant, but $\left\{ \boldsymbol{\theta}_{s}^{*}\right\} _{1\leqq s\leqq t+1}$ vary, then $z_{t+1}$ is a martingale process, which is orthogonal to the martingale term of $\boldsymbol{x}_{t}$. While it might seem like an issue that the difference, $\varDelta x$, appears as a covariate in the linear model, it is done to make the cross term orthogonal. Using this expression, we can rewrite the best predictive model of $y_{t+1}$ at $t$ as \[ y_{t+1}=y_{t}+\theta_{y_{t+1}}^{*}-\theta_{y_{t}}^{*}+\left\langle \boldsymbol{\theta}_{t+1}^{*},\varDelta\boldsymbol{x}_{t}\right\rangle +z_{t+1}-z_{t}. \] This can be seen by subtracting the FS decomposition at $t$, $y_{t}=\theta_{y_{t}}^{*}+\sum_{s=1}^{t-1}\left\langle \boldsymbol{\theta}_{s+1}^{*},\varDelta\boldsymbol{x}_{s}\right\rangle +z_{t}$, from eq. (ref).
Since the process, $\left\{ z_{t}\right\} $, is a square-integrable martingale, if we consider $\sigma_{t+1}^{2}$ to be the quadratic variation of $z_{t+1}$, then $\sigma_{t+1}^{2}$ is the conditional variance of $\varDelta z_{t}$: $\sigma_{t+1}^{2}=\mathbb{E}\left[\left.\left(\varDelta z_{t}\right)^{2}\right|\left\{ y_{s},\boldsymbol{x}_{s}\right\} _{s=1}^{t}\right]$.
Here, we add a new assumption on the distribution of $\varDelta z_{t}$: the probability distribution, $p_{t,\varDelta z}$, for $\varDelta z_{t}$ only has the scale parameter, $\sigma_{t+1}$. In other words, if we scale transform $\varDelta z_{t}\rightarrow c\varDelta z_{t}$ and $\sigma_{t+1}^{2}\rightarrow c^{2}\sigma_{t+1}^{2}$, $p_{t,\varDelta z}$ is scale-invariant:
This assumption justifies the approximate normality of $p_{t,\varDelta z}$, when the time interval, $\varDelta t$, is not large (i.e. we assume that the 1-step ahead is not in the distant future).
Let us rewrite the extra term as $\theta_{0,t+1}^{*}=\theta_{y_{t+1}}^{*}-\theta_{y_{t}}^{*}$, and represent the unknown parameter at $t$ with $\left(\theta_{0,t+1}^{*},\boldsymbol{\theta}_{t+1}^{*}\right)$. Then, noting that $y_{t+1}=y_{t}+\varDelta y_{t}$, the probability density function of the increment of $y_{t}$, can be written using the probability density function of its martingale process, $\varDelta z_{t}=z_{t+1}-z_{t}$, as \[ p_{t}\left(y_{t+1}\left|\left\{ y_{s},\boldsymbol{x}_{s}\right\} _{s=1}^{t},\theta_{0,t+1}^{*},\boldsymbol{\theta}_{t+1}^{*},\sigma_{t}^{2},\varDelta\boldsymbol{x}_{t}\right.\right)=p_{t,\varDelta z}\left(\frac{y_{t+1}-\left(y_{t}+\theta_{0,t+1}^{*}+\left\langle \boldsymbol{\theta}_{t+1}^{*},\varDelta\boldsymbol{x}_{t}\right\rangle \right)}{\sigma_{t}^{2}}\right). \]
We consider some assumptions on $y_{t}$, in preparation for minimax analysis. For the decomposition of $y_{t}$ (eq. (ref)), we assume the martingale component of $\boldsymbol{x}_{t}$ and the increment of its orthogonal martingale process, $z_{t}$, are square-integrable for all $t$. The probability distribution function of $z_{t}$ is denoted as $p_{z_{t}}\left(\cdot\right)$, where $p_{z_{t}}\left(\cdot\right)$ can depend on $t$, given the subscript $t$, and further assume eq. (ref). Note that $p_{z_{t}}\left(\cdot\right)$ can be considered as a probability density conditional on $\boldsymbol{x}_{t}$, from williams1991probability, Section 9.5, Section 14.14.
For $\min_{\hat{p}_{t}}$ in the minimax theorem, we construct the Bayes predictive distribution of the DLM model following Section (ref) by corresponding the DLM parameters to the unknown parameters, $\{\theta_{0,t+1}^{*},\boldsymbol{\theta}_{t+1}^{*},\sigma_{t+1}\}$. This allows us to construct the minimum Bayes risk distribution.
Given the above specification, the proof for this theorem follows in a similar manner.
While the corresponding of DLM parameters to the unknown parameters of the DGP was straightforward for the case where we assumed Gaussian processes, the semi-martingale case is less straightforward.
Consider the random walk DLM,
where $\boldsymbol{\theta}_{t}^{\prime}$ is the vector $\left[\theta_{0,t},\boldsymbol{\theta}_{t}\right]^{\top}$, with intercept, $\theta_{0,t}$, $\boldsymbol{x}_{t}^{\prime}$ is $\left[1,\boldsymbol{x}_{t}\right]^{\top}$, and $\boldsymbol{\omega}_{t+1}$ is a $J+1$-dimensional vector of normal random numbers. We will now correspond the observational equation (eq. (ref)) to the FS decomposition in eq. (ref). First, noting that,
we have \[ y_{t+1}=\theta_{0,t+1}+\left\langle \boldsymbol{\theta}_{t},\boldsymbol{x}_{t}\right\rangle +\left\langle \Delta\boldsymbol{\theta}_{t},\boldsymbol{x}_{t}\right\rangle +\left\langle \boldsymbol{\theta}_{t+1},\Delta\boldsymbol{x}_{t}\right\rangle +\varepsilon_{t+1}, \] by setting
where $\theta_{0,0}$ is an initial value of $\left\{ \theta_{0,t}\right\} $. From this, we have
where $\theta_{0,t}$ follows a random walk. Therefore, eq. ((ref)) is transformed to
where the parameters considered in eq. ((ref)), $\left(\theta_{0,t+1}^{*},\boldsymbol{\theta}_{t+1}^{*}\right)$, correspond to $\left(\boldsymbol{\theta}_{t+1}^{\prime}\right)$ in eq. ((ref)) through
Further, we assume the error, $\varepsilon_{t+1}$, in eq. ((ref)), follows a normal distribution. In general, the distribution of $\varepsilon_{t+1}$ and $z_{t+1}$, obtained from the data generating process, need not be the same. Since $z_{t+1}$ is a process that includes omitted variables, it is unrealistic to exactly identify it.
As an extension and application of the above theoretical results, we consider the problem of ensembling multiple predictive distributions. In time series contexts, model misspecification and uncertainty pose a significant problem, particularly for predictive tasks. Ensemble methods, which include model averaging and forecast combination, have been used extensively to mitigate this uncertainty. In the Bayesian context, some popular methods include Bayesian model averaging hoeting1999bayesian and Bayesian stacking yao2018using. Recent interest in combining distributional information to account for better uncertainty quantification has spurred development in ensembling predictive densities HallMitchell2007,Geweke2011,Aastveit2014,diebold2019machine.
While ensembling strategies have shown to be effective in practice, improving over individual forecasts, the theoretical properties of the ensembled predictions have not been investigated in a non-stationary time series context. We are specifically interested in the Bayesian predictive synthesis framework proposed in mcalinn2019dynamic, which extended the work on expert opinion analysis GenestSchervish1985,West1992c,West1992d and generalizes other ensembling strategies as a special case. In particular, mcalinn2019dynamic developed a specific form of dynamic BPS for time series data, utilizing a random walk DLM as its synthesis function. Given its empirical success bianchi2018large,mcalinn2017multivariate,capek2020macroeconomic,mcalinn2017dynamic,aastveit2022quantifying,bernaciak2022loss, and relation to the random walk DLM, a natural question is whether the results in Sections (ref) and (ref) also hold in this context.
In the ensemble prediction set-up, a decision maker, $\mathcal{D}$, is interested in forecasting a quantity $y\in\mathbb{R}$, and solicits $J$ agents, where agents encompass models, forecasters, institutions, etc. Agents are denoted as $\mathcal{A}_{j}$, $j\in J.$ Each agent, $\mathcal{A}_{j}$, produces a predictive distribution, $h_{j}(\hat{y}_{j})$, which comprises the set, $\mathcal{H}=\left\{ h_{1}(\cdot),\cdots,h_{J}(\cdot)\right\} $.
In the BPS framework, the set, $\mathcal{H}$, is synthesized via Bayesian updating with the posterior of the form,
where $\hat{\boldsymbol{y}}=\hat{y}_{\seq1J}=(\hat{y}_{1},\ldots,\hat{y}_{J})^{\top}$ is a $J-$dimensional latent vector and $\alpha(y|\hat{\boldsymbol{y}})$ is a conditional probability density function for $y$ given $\hat{\boldsymbol{y}}$, called the synthesis function. Eq. (ref) is only a coherent Bayesian posterior if it satisfies the consistency condition\footnote{Specifically, the consistency condition states that given $\mathcal{D}$'s prior, $p(y)$, and her prior expectation of what the agents will produce before observing the forecasts, $m(\hat{\boldsymbol{y}})$, her priors have to be consistent: $p(y)=\int_{\hat{\boldsymbol{y}}}\alpha(y|\hat{\boldsymbol{y}})m(\hat{\boldsymbol{y}})d\hat{\boldsymbol{y}}$, for eq. (ref) to be a coherent Bayesian posterior.} GenestSchervish1985,West1992c,West1992d,mcalinn2019dynamic.
For time series data, mcalinn2019dynamic developed a dynamic specification of BPS. In the dynamic setting, $\mathcal{D}$ is sequentially predicting a time series $y_{t},t=1,2,\ldots,$ and receives, for all periods, forecast densities from $\mathcal{A}_{j}$, $j\in J.$ At each time $t,$ $\mathcal{D}$ aims to forecast $y_{t+1}$ and receives the set $\mathcal{H}_{t+1}=\left\{ h_{1,{t+1}}(\hat{y}_{1,{t+1}}),\cdots,h_{J,{t+1}}(\hat{y}_{J,{t+1}})\right\} $, a collection of agent predictive densities produced at $t$, from the agents. The full information set used by $\mathcal{D}$ is thus $\{\mbox{\boldmath$y$}_{\seq1{t}},\ \mathcal{H}_{\seq1{t+1}}\}.$ As time passes, $\mathcal{D}$ learns about the characteristics of the agents (bias, dependencies, etc.). Thus, the Bayesian model will involve parameters for which $\mathcal{D}$ updates information over time. From eq. (ref), $\mathcal{D}$ has a time $t$ distribution for $y_{t+1}$ of the form
where $\mbox{\boldmath$\Phi$}_{t+1}=(\theta_{0,{t+1}},\mbox{\boldmath$\theta$}_{t+1},v_{t+1},\boldsymbol{W}_{t+1})$, and the synthesis function, $\alpha_{t}(y_{t+1}|\hat{\boldsymbol{y}}_{t+1},\mbox{\boldmath$\Phi$}_{t+1})$, in mcalinn2019dynamic is specified as a standard dynamic linear model eqs. (ref)-(ref).
Finally, the transition probability can be obtained by integrating out the agent processes, $\hat{\boldsymbol{y}}_{t+1}$:
where $\hat{\boldsymbol{y}}_{t}^{\prime}$ is $\left[1,\hat{\boldsymbol{y}}_{t}\right]^{\top}$.
Given that the problem setting now involves densities to be synthesized, we show that Theorem. (ref) holds for the convolution case:
The Gaussianity assumptions in the above theorem can also be relaxed using semi-martingale processes, as in Section (ref). This result proves that dynamic BPS, as proposed in mcalinn2019dynamic, produces exact minimax predictive distributions.
Additionally, the above results show that linear combinations of forecasts (e.g. equal weight averaging, Bayesian model averaging, etc.) cannot produce predictive distributions that are exact minimax. This is because $\max_{\theta_{0,t+1}^{*}}\mathsf{KL}\left(\left(\theta_{0,t+1}^{*},\boldsymbol{\theta}_{t+1}^{*}\right),\hat{p}_{t}\right)$ can be infinite due to the BPS intercept, $\theta_{0,t+1}^{*}$, being under-parameterized for linear combination methods. This not only shows that BPS is theoretically better (in terms of being minimax), but that no linear combination of forecasts can be better under this criterion.
We consider three datasets from epidemiology, climatology, and economics (for a simulation study, see Appendix (ref)). All three datasets have common characteristics that highlight our theoretical results: they are all non-stationary time series data, a “true” model is not something that can be expected to obtain, and future predictions are tied to some decision making process, where uncertainty quantification is critical.
The first dataset concerns predicting the weekly average new cases of COVID-19 in Tokyo, Japan. The prediction of positive cases of COVID-19 has become a critical task in tackling the pandemic, drawing considerable interest awan2020prediction,gecili2021forecasting. This is because accurate predictions of new cases, as well as their variation, have been shown to be imperative for resource management in hospitals and local governments. For this application, we use daily data on new cases to predict the next week's average new cases, as the decision maker's interests are not necessarily in new cases “tomorrow,” but in the coming week, in order to, e.g., prepare medical resources. The dataset, taken from the Tokyo Metropolitan Government, begins on 2020/1/16, when the collection of data began, and ends at the end of 2021. Our interest is to predict over the entirety of 2021, where Tokyo experienced three major waves of new cases, including the Delta variant, and began to experience the Omicron variant.
The second dataset aims to predict the monthly global mean sea level from TOPEX and Jason Altimetry nerem2010estimating,masters2012comparison. The global mean sea level is a key indicator for climate change, reflecting the ocean's thermal expansion, meltwater from mountain glaciers, and discharge from the Greenland and Antarctic ice sheets. As such, the prediction of the global mean sea level is essential in understanding future trends in climate change. The satellite data is measured every ten days, though the focus is on the monthly (30-day cycle) trends of sea level, and its prediction. This is done, as with the COVID dataset, because shorter frequency (i.e. 10-day cycle) is prone to local fluctuations that may not be relevant for climate change, and might not be of interest to predict. We consider the monthly frequency to be of interest in this application in order to gauge the accuracy of ensemble strategies.
The third is a large-scale macroeconomic dataset to forecast monthly U.S. inflation, a context of interest in the field Cogley2005,Nakajima2010. Specifically, forecasting inflation is one of the central bank's key mandates, where that information is used to set policy. Typical horizon of interest is 1-, 3-, 12-, and 24-months ahead. Although long term forecasts can be, and are, cast as one step ahead forecasts within the BPS framework, as BPS can directly calibrate to the horizon of interest, we focus on the shortest horizon of interest, 1-month, as this exemplifies the theoretical results the most. For this study, we utilize a high-dimensional panel of $N=128$ monthly macroeconomic and financial covariates from mccracken2016fred (details of each variable can be found therein). This context of using a large dataset for economic forecasting is of topical interest, as leveraging large datasets for decision making is crucial for central banks setting their policy.
All three applications are conducted in a similar manner, with some variation to suit the context. We first divide the dataset equally into the training and evaluation set. The first half of the training set is used to build the individual models, in parallel. This is 2020/1/16--2020/6/27 for the COVID dataset, 1993:01--2002:05 for the sea level dataset, and 1986:01--1993:06 for the economic dataset. Specifically, we use a random walk DLM, for each agent model, $j=1{:}J$,
where the coefficients follow a random walk and the observation variance evolves with discount stochastic volatility. Priors for each predictive regression are set as $\mbox{\boldmath$\beta$}_{0j}|v_{0j}\sim N(\mbox{\boldmath$m$}_{0j},(v_{0j}/s_{0j})\mbox{\boldmath$I$})$ with $\mbox{\boldmath$m$}_{0j}={\bf 0}'$ and $1/v_{0j}\sim G(n_{0j}/2,n_{0j}s_{0j}/2)$ with $n_{0j}=10,s_{0}=0.01$.
In the COVID dataset, each model is built with different autoregressive lags, covariates, and day of the week indicators to capture the variation in testing capabilities. Since the data are daily, and the interest is the coming weekly average, we consider four models: AR(1)+Policy (Policy: two indicators for the two levels of state of emergency), AR(1)+Vax (Vax: six covariates for the vaccination rates for the first, second, and third shots for the entire population and 65 and above), AR(1)+Mob+Temp (Mob: six covariates for weekly average Google mobility for retail and recreation, grocery and pharmacy, parks, stations, workplaces, and residential; Temp: three covariates of weekly average rainfall, temperature, and humidity), AR(3) and AR(7).
For the sea level dataset, each model is built with different autoregressive lags: AR(1), AR(3), AR(9), and AR(18), as the data are monthly (30-day cycle).
For the economic dataset, the covariates are the $N=128$ monthly macroeconomic and financial covariates, partitioned into eight groups based on the existing qualitative classification. The eight main categories are: Output and Income, Labor Market, Consumption and Orders, Orders and Inventories, Money and Credit, Interest Rate and Exchange Rates, Prices, and Stock Market.
For the next half of the training data, while the individual DLMs are sequentially updated in parallel, we also begin calibrating and training the ensemble methods, including dynamic BPS. This is 2020/6/28--2020/12/31 for the COVID dataset, 2002:06--2010:12 for the sea level dataset, and 1986:01-2000:12 for the economic dataset.
Finally, the latter half of the entire dataset is used as the evaluation period, where we compare predictive performances of the ensemble methods in an online manner, mirroring real world situations. This is 2021/1/1--2021/12/31 for the COVID dataset, 2011:01--2020:12 for the sea level dataset, and 2001:01-2015:12 for the economic dataset. Updating of the individual models, as well as competing strategies, is done sequentially for each $t$ during this period.
We compare the mean squared forecast error (MSFE) and the log predictive density ratio (LPDR). If we assume a Gaussian process, the performance relation in KL risk is equivalent to the performance relation in the MSE risk. The log predictive density ratios (LPDR) for each $t$ is
where $p_{*}(y_{t+1}|y_{1{:}t})$ is the predictive density of the model being compared with. This metric directly corresponds to the KL divergence. Comparing both the MSFE and LPDR provides a more holistic assessment of the predictive performance.
We compare our framework against each agent DLM and benchmark ensemble strategies. Our goal is to show that dynamic BPS is, not only superior to other ensemble methods, but superior to the forecasts that it is ensembling; what would be considered a necessary requirement for ensemble methods to achieve.
Specifically, we compare dynamic BPS against equal weight averaging, BMA, Mallows $C_{p}$ averaging of hansen2007least, and the exponential weights algorithm of littlestone1994weighted. Equal weight averaging is standard in much of the forecast combination literature and ensemble learning literature, being almost ubiquitous in the latter. While the approach of taking the arithmetic mean of forecasts is very simple, it is, nonetheless, considered a benchmark that is hard to beat in real data applications Genre2013. On the other hand, BMA has the benefit of asymptotically converging to the true model, when the true model is nested ($\mathcal{M}$-closed), though it tends to converge to the “wrong" model in an $\mathcal{M}$-open setting. This is because, despite its name, BMA is more of a selection strategy rather than an ensemble strategy. In almost all real applications, including this one, we cannot assume that the true model is nested, which guarantees that BMA will converge to the “wrong" model. Nonetheless, it is standard in much of the Bayesian literature. Mallows $C_{p}$ averaging has the property of being asymptotically optimal, achieving the lowest possible squared error in a class of discrete model average estimators. The exponential weights algorithm, otherwise called the weighted majority algorithm, bounds the regret of not choosing the best agent.
The prior specification for dynamic BPS is $\theta_{0,0}|v_{0}\sim N(0,v_{0}/s_{0})$, $\mbox{\boldmath$\theta$}_{0}|v_{0}\sim N(\mbox{\boldmath$m$}_{0},(v_{0}/s_{0})\sigma^{2})$ with $\mbox{\boldmath$m$}_{0}=(0,\mbox{\boldmath$1$}'/J)'$ and $1/v_{0}\sim G(n_{0}/2,n_{0}s_{0}/2)$ with $n_{0}=10,s_{0}=0.002$ (the variance is scaled by 1/100 for the sea level dataset to reflect the variation in the data). The discount factors are set to $(\beta,\delta)=(0.95,0.99)$.
Panels A, B, and C of Table (ref) show the predictive results for the COVID dataset, sea level dataset, and macroeconomic dataset, respectively. Over the three datasets, dynamic BPS improves the out-of-sample forecasting accuracy compared to everything we consider, including all individual models and the four ensemble methods. This result is consistent for all datasets, though the performance gains for BPS differ amongst datasets. Most notably, the sea level datasets display lower gains compared to the other two, although the gains are still significant. This is likely due to the consistent trend exhibited in the data, where large fluctuations are not as clear as the other two. Nonetheless, dynamic BPS shows improved performance for both point and density evaluations.
In particular, for the COVID dataset and the macroeconomic dataset, dynamic BPS consistently improves over the other ensemble methods by large margins for both point and density forecasts. For the COVID data, specifically, dynamic BPS has an LPDR of approximately -40 to -100 across all ensemble methods compared against. The difficulty of standard ensemble methods in dealing with non-stationary data is also clear, with some even underperforming the agent models. For RMSE, Mallows $C_{p}$ edges out only in the sea level dataset, and for LPDR, Mallows $C_{p}$ slightly outperforms the best agent model in the COVID dataset and BMA is practically equivalent to the best model in the sea level dataset. Notably, the best performing agent model in terms of RMSE is not the same as the LPDR, except for the macroeconomic dataset. The standard ensemble methods seemingly struggle with this dissonance in performance, though dynamic BPS is able to successfully synthesize the information to outperform all in both metrics. This shows that, while the theory we presented is with regard to KL risk, dynamic BPS can outperform other ensemble methods for different decision making problems.
Finally, we look at the interpretability of dynamic BPS through its on-line coefficients. Figure. (ref) shows the on-line mean coefficients of BPS for the weekly COVID prediction dataset. Most notably, the intercept is quite volatile compared to the other coefficients. As the intercept can be viewed as the model set uncertainty (uncertainty not captured by the agent models), this shows how the different lags and covariates are not sufficient to predict weekly COVID cases. The intercept also plays a critical role in the theoretical results, and the fact that the intercept adapts over time at this magnitude, as seen in the figure, shows this in action.
Looking at the other coefficients, several points stand out. First, as an overall trend, the model with vaccination rates increases in importance over time. This suggests the importance of vaccination information to predict COVID cases, with its importance increasing as the vaccination rate increases. Second, the AR(7) model, which captures longer trends, spikes during increased COVID cases. This could be explained by the lagged nature of COVID diagnosis. Notably, this model is inversely correlated with the AR(1)+Mob+Temp model, which models mobility and weather data. Third, the model with policy indicators is effective during the first semi-emergency measure, though its importance diminishes throughout the latter emergencies. This is in line with the idea that these measures are not as effective in the latter iterations, due to people adapting or being complacent towards these measures.
As these results show, the on-line coefficients provide critical insight into the predictive process, as well as highlight the theoretical findings of this paper. In particular, the trend exhibited in the intercept supports the theoretical results that the intercept plays a crucial role. Added to the fact that dynamic BPS improves over other ensemble methods, the result empirically supports the theoretical findings.
Despite the importance of non-stationary time series forecasts in many domains, theoretical results are virtually non-existent. This is due to the dynamic nature of the data, where assumptions-- such as i.i.d.--, used in the existing literature to show asymptotic results, do not hold. This has stymied theoretical evaluations of predictive performances, especially when all models are misspecified. To evaluate the performance of statistical methods in this setting, we define the Kullback-Leibler risk for non-stationary data. Using this criterion, we prove that a random walk DLM is exact minimax, providing finite sample theoretical support for the method over other existing methods. We then show that, extending the problem to predictive density synthesis, dynamic BPS is also exact minimax. Through three relevant applications from epidemiology, climatology, and economics, we show that BPS; (i) outperforms individual models that it synthesizes, (ii) outperforms other ensemble methods, as the theory indicates, and (iii) provides useful information for decision making. All three real data applications reinforce the validity and applicability of the theoretical result.
We believe the results presented in this paper open up several avenues of research related to time series forecasting and decision making, as well as causal inference with the potential outcomes framework. If, for example, the interest is in predicting potential outcomes for time series data, this framework can provide some theoretical guarantee, even if the data cannot be assumed to be i.i.d. Different problem settings will require further development of the framework, though having the Kullback-Leibler risk defined is the first step towards solving each problem.
One notable direction for future development is in the direction of multivariate models, which is also relevant for policy decision making. However, extending the proposed framework to a multivariate setting is not trivial, due to the shift transform only holding for univariate or very specific multivariate data. For this purpose, a different approach to replace the shift transform is necessary, which will be future work. Further connections can be made to the generalized Bayes framework bissiri2016general, as seen in bernaciak2022loss,tallman2022bayesian, where the generalized Bayes framework is formulated as a special case of BPS (i.e. the synthesis function is a loss function). The minimax result, which is decision theoretic, could potentially produce fruitful results in this direction as well.