*

# Sequential Bayesian Predictive Synthesis

Riku Masuda111Graduate School of Economics, The University of Tokyo.
E-Mail<EMAIL_ADDRESS>and Kaoru Irie222Corresponding author. Faculty
of Economics, The University of Tokyo.
E-Mail<EMAIL_ADDRESS>
###### Abstract

Dynamic Bayesian predictive synthesis is a formal approach to coherently
synthesizing multiple predictive distributions into a single distribution. In
sequential analysis, the computation of the synthesized predictive
distribution has heavily relied on the repeated use of the Markov chain Monte
Carlo method. The sequential Monte Carlo method in this problem has also been
studied but is limited to a subclass of linear synthesis with weight
constraint but no intercept. In this study, we provide a custom, Rao-
Blackwellized particle filter for the linear and Gaussian synthesis,
supplemented by timely interventions by the MCMC method to avoid the problem
of particle degeneracy. In an example of predicting US inflation rate, where a
sudden burst is observed in 2020-2022, we confirm the slow adaptation of the
predictive distribution. To overcome this problem, we propose the
estimation/averaging of parameters called discount factors based on the power-
discounted likelihoods, which becomes feasible due to the fast computation by
the proposed method.

Key words and phrases: Dynamic Bayesian predictive synthesis, forecast density
combination, sequential Monte Carlo, Rao-Blackwellization, log discounting
framework.

##  1 Introduction

In modern forecasting problems, it is common that multiple predictive
distributions are available for an statistician. The source of these
predictive distributions can include human experts and deterministic computer
algorithms, but in data analysis, it is most likely statistical models
analyzed by the statistician. Predictive synthesis is a problem where the
statistician uses those predictive distributions to create a new, single
predictive distribution for an improved prediction. A formal approach to this
problem is Bayesian predictive synthesis (BPS); it is based on the series of
research on the coherent rule of updating the statistician’s predictive
distribution as “posterior” informed by the predictive distributions as “data”
(Genest and Schervish, 1985; West and Crosse, 1992; West, 1992) and includes
the traditional Bayesian model averaging as a special case. Its extension to
dynamic models has been proposed by McAlinn and West, (2019) (dynamic Bayesian
predictive synthesis, DBPS). Predictive analysis by (D)BPS has seen various
applications, including the prediction of macroeconomic indices (McAlinn et
al.,, 2020; McAlinn,, 2021), oil prices (Aastveit et al.,, 2023), real estate
prices (Cabel et al.,, 2022) and decision analysis (Tallman and West,, 2022).

Although the synthesized prediction by DBPS outperforms other predictive
methods in many scenarios (e.g., see McAlinn and West, 2019), its computation
relies heavily on the repeated use of the Markov chain Monte Carlo (MCMC)
method. In fact, in all the papers on DBPS cited above, the synthesized
predictive distributions are computed by re-doing the MCMC method. This is
because the synthesized predictive distribution is not analytically tractable,
even if we assume the linear and Gaussian synthesis structure. The computation
cost becomes more expensive if the statistician has to compute the synthesis
of the predictive distributions multiple times. Examples include real-time
monitoring, sequential forecasting, and calibrating the predictive synthesis
(or trying different ways of synthesizing the predictive distributions).
Especially in sequentially processing the streaming data, such repeated use of
the MCMC method is prohibitive, which has hindered the use of DBPS in
practice. In this context, seeking online computation by the sequential Monte
Carlo (SMC) method is natural.

In this study, we propose an SMC-based approach customized for DBPS and
evaluate its computational performance. The use of the particle filter (Gordon
et al.,, 1993) has been considered in the context of forecast density
pooling/combination (Billio et al.,, 2013; Casarin et al.,, 2015), and has
been limited to a particular class of synthesis functions. We, too, consider
the particle filter, but extend its application to the general class of
synthesis functions in the DBPS framework. Furthermore, we focus on the linear
and Gaussian synthesis and provide the custom, Rao-Blackwellized Bootstrap
particle filter for online predictive analysis (Chen and Liu, 2000; Doucet et
al., 2000). To handle the particle degeneracy problem, we monitor the
effective sample size (ESS) and, if the ESS is small and below a threshold, we
intervene in the sequential computation using the off-line MCMC method
instead. An appropriate and timely intervention by the MCMC method is crucial
for the accurate computation of the synthesized predictive distribution and
does not critically slow down the sequential computation.

The proposed method— the pair of the DBPS model and filtering algorithm— is
successful in providing almost the same predictive distribution as the MCMC
method does in most cases. However, the DBPS is critically insensitive to a
sudden burst of observations (or potential structural change) and could
perform poorly in predictions, regardless of the computational methods. This
sudden burst is typical in the analysis of inflation rates studied in the
literature, especially in the worldwide pandemic and economic crisis in
2020-2022. To tackle this problem, we introduce the loss discounting framework
(LDF, Bernaciak and Griffin, 2022) into the DBPS. This framework requires the
computation of the DBPS multiple times at every time point, which is a typical
example where the fast posterior/predictive computation by the proposed SMC
method is crucial.

The rest of this paper is organized as follows. Section 2 reviews the DBPS and
the linear and Gaussian synthesis functions, emphasizing properties related to
the SMC method. In Section 3, we develop the particle filter for the DBPS,
including its Rao-Blackwellization. We illustrate the proposed SMC method
using the dataset of the US inflation rates in Section 4, where the efficiency
and predictive performance of the proposed method are evaluated. In Section 5,
we discuss the introduction of the LDF into the DBPS and observe the improved
predictive performance in the inflation example. We comment on future research
in Section 6

Notational remarks: We denote the normal distribution with mean vector $\mu$
and variance matrix $\Sigma$ by $N(\mu,\Sigma)$. For integers $s<t$, we write
$s:t=(s,s+1,\dots,t-1,t)$.

##  2 Dynamic Bayesian predictive synthesis

###  2.1 Overview and terminology

Suppose that $K$ predictive distributions,
$\mathcal{H}_{t}=\\{h_{1t},\dots,h_{Kt}\\}$, are provided to forecast quantity
$y_{t}$. Bayesian predictive synthesis (or density combination if density
functions exist) discusses how we should combine the $K$ predictive
distributions into a single distribution as the final output for prediction.
The source of the predictive distributions can be anything, including reports
from human experts (e.g., Chernis, 2022). However, in most applications of the
BPS, those predictive distributions are obtained from statistical models. We
call those $K$ models and predictive distributions agent models and agent
predictive distributions, respectively. The resulting predictive distribution
of the BPS is named synthesized predictive distribution. In this study, we
focus on the one-step ahead prediction for simplicity. Still, the BPS is open
to multi-step ahead predictions and path forecasting, to which our study can
be easily generalized.

The dynamic BPS involves two steps of computation: model learning and
synthesis calibration. The former is the process of creating agent predictions
$\mathcal{H}_{t}$, and the latter is to revise (or calibrate) the way of
synthesizing $\mathcal{H}_{t}$. Since the calibration step depends on
$\mathcal{H}_{t}$ created by learning about the agent models, these steps
cannot be processed in parallel.

###  2.2 Agent predictive models

Let $\\{y_{t}\\}_{t=1:T}$ be univariate time series to be forecast. We have
$K$ agent models about $y_{t}$, denoted by $\\{\mathcal{M}_{k}\\}_{k=1:K}$. At
each time point $t-1$, we learn about each model $\mathcal{M}_{k}$ by using
observed $y_{1:(t-1)}$ and compute the one-step ahead agent predictive
distribution, namely, $p(y_{t}|y_{1:(t-1)},\mathcal{M}_{k})$. We denote this
agent predictive density by $h_{kt}(\cdot)$ as

$h_{kt}(x_{kt})=p(x_{kt}|y_{1:(t-1)},\mathcal{M}_{k}),$

where we use $x_{kt}$ as the argument of the density, not $y_{t}$, to
emphasize that this is the prediction made by agent model $\mathcal{M}_{k}$ at
time $t-1$. The set of agent predictive densities, $\mathcal{H}_{t}$, becomes
available at time $t-1$.

In this study, we use the conjugate dynamic models as the agent models whose
online posterior and predictive distributions are obtained analytically. In
general, the agent predictive distributions are not analytically tractable,
and are often approximated using simulation-based methods (SMC, MCMC, and
others). The computational methods we discuss below are also applicable to
such approximated predictive distributions, where the simulation from
(appoximated) $h_{kt}$ is feasible.

###  2.3 BPS, synthesis functions and priors

A formal approach to synthesizing the $K$ agent predictive distributions is to
use the conditional distribution $p(y_{t}|\mathcal{H}_{t})$, which is
difficult to compute directly using the Bayes rule. Instead, in dynamic
Bayesian predictive synthesis (DBPS, McAlinn and West, 2019), we synthesize
the $K$ agent predictive densities,
$\mathcal{H}_{t}=\\{h_{1t},\dots,h_{Kt}\\}$, by time-varying synthesis
function $\alpha_{t}(y_{t}|x_{t})$ at each time as follows:

$p(y_{t}|\mathcal{H}_{t},\alpha_{t})=\int\alpha_{t}(y_{t}|x_{t})h_{t}(x_{t})dx_{t},$

where $x_{t}=(x_{1t},\cdots,x_{Kt})^{\prime}$ and
$h_{t}(x_{t})=\prod_{k=1}^{K}h_{kt}(x_{kt})$. In fact, the distribution of
$y_{t}$ conditional on $\mathcal{H}_{t}$ must be of this form for the joint,
marginal, and conditional distributions to be consistently defined (Genest and
Schervish, 1985, West and Crosse, 1992, West, 1992). In practice, the
synthesis function, $\alpha_{t}(\cdot)$, must be specified by users. We
parametrize this synthesis function by using dynamic calibration parameters,
$\\{\Phi_{t}\\}_{t=1:T}$, and rewrite the synthesis equation as

$p(y_{t}|\mathcal{H}_{t},\Phi_{t})=\int\alpha(y_{t}|x_{t},\Phi_{t})h_{t}(x_{t})dx_{t},$

where the functional form of $\alpha(\cdot)$ is pre-determined and unchanged
over time. Note that, in predicting $y_{t}$, we condition not only set of
agent predictions $\mathcal{H}_{1:t}$ but also past observations
$y_{1:(t-1)}$. Thus, the synthesized predictive distribution, or the final
output for prediction, is obtained by marginalizing $\Phi_{t}$ out as

$\begin{split}p(y_{t}|\mathcal{H}_{1:t},y_{1:(t-1)})&=\int
p(y_{t}|\mathcal{H}_{t},\Phi_{t})p(\Phi_{t}|y_{1:(t-1)},\mathcal{H}_{1:(t-1)})d\Phi_{t}\\\
&=\iint\alpha(y_{t}|x_{t},\Phi_{t})h_{t}(x_{t})p(\Phi_{t}|\mathcal{H}_{1:(t-1)},y_{1:(t-1)})dx_{t}d\Phi_{t}.\end{split}$

As seen in this expression, we calibrate how to synthesize the multiple agent
predictive densities by learning about the calibration parameter, $\Phi_{t}$.
The set of calibration parameters includes the location and scale of each
forecast and adjusts the bias of the point prediction and predictive
uncertainty, as specified in Section 2.4 in detail.

To compute $p(y_{t}|\mathcal{H}_{1:t},y_{1:(t-1)})$ using the integral
expression above, we need to calculate
$p(\Phi_{t}|\mathcal{H}_{1:(t-1)},y_{1:(t-1)})$, or the online prior
distribution of calibration parameters. This distribution is further augmented
as

$\begin{split}&p(\Phi_{t}|\mathcal{H}_{1:(t-1)},y_{1:(t-1)})\\\ &=\int
p(\Phi_{t}|\Phi_{1:(t-1)},x_{1:(t-1)},y_{1:(t-1)})p(\Phi_{1:(t-1)},x_{1:(t-1)}|\mathcal{H}_{1:(t-1)},y_{1:(t-1)})d\Phi_{1:(t-1)}dx_{1:(t-1)}.\end{split}$

Here, the evolution of the calibration parameters,
$p(\Phi_{t}|\Phi_{1:(t-1)},x_{1:(t-1)},y_{1:(t-1)})$, depends on the
specification of synthesis function $\alpha$ and discussed in the next
subsection. The SMC method is applied to the computation of the online
posterior distribution of $(\Phi_{1:(t-1)},x_{1:(t-1)})$ and studied in
Section 3.

###  2.4 Dynamic linear model synthesis

Dynamic linear models (DLMs), or linear and Gaussian state space models, are
an important class of dynamic models (e.g., see Prado et al., 2021) and can
also be used in modeling synthesis function $\alpha(y_{t}|x_{t},\Phi_{t})$ and
the evolution of $\Phi_{t}$. Following McAlinn and West, (2019), we define the
DLM synthesis with discount factors $\beta,\delta\in(0,1]$ as

$\begin{split}&\alpha(y_{t}|x_{t},\Phi_{t})=N(y_{t}|F_{t}^{\prime}\theta_{t},\nu_{t})\\\
&\theta_{t}=\theta_{t-1}+\omega_{t}\ ,\ \omega_{t}\sim N(0,\nu_{t}W_{t})\\\
&\nu_{t}=\frac{\beta}{\gamma_{t}}\nu_{t-1}\ ,\ \gamma_{t}\sim
Be\left(\beta\frac{n_{t-1}}{2},(1-\beta)\frac{n_{t-1}}{2}\right)\end{split}$

where $F_{t}=(1,x_{1t},\cdots,x_{Kt})^{\prime}$,
$\theta_{t}=(\theta_{0t},\theta_{1t}\cdots,\theta_{Kt})^{\prime}$, and
$Be(a,b)$ denotes beta distribution. State variance $W_{t}$ is defined via a
standard, single discount factor specification with state evolution discount
factor $\delta\in(0,1]$. The residual variance $\nu_{t}$ follows the standard
beta-gamma random walk volatility model with discount factor $\beta\in(0,1]$
where the shape parameter is updated by $n_{t}=\beta n_{t-1}+1$. In this
setting, online posterior $p(\theta_{t},\nu_{t}|x_{1:t},y_{1:t})$ and prior
$p(\theta_{t+1},\nu_{t+1}|x_{1:t},y_{1:t})$ are normal-inverse-gamma
distributions whose sufficient statistics are deterministically updated.
Furthermore, predictive distribution $p(y_{t}|y_{1:(t-1)},x_{1:t})$ is
Student’s $t$-distribution, hence easily evaluated and simulated from. For
details about conjugacy and hyperparameter settings, see the Supplementary
Materials (Section S1).

In the DLM synthesis, $\Phi_{t}=\theta_{t}$ is the calibration parameter
($\nu_{t}$ is marginalized out). The conjugacy of the DLMs enables the
analytical computation of, and the direct simulation from, the online
posterior and prior distributions, $p(\Phi_{t}|x_{1:t},y_{1:t})$ and
$p(\Phi_{t+1}|x_{1:t},y_{1:t})$, and the predictive distribution,
$p(y_{t+1}|x_{1:t+1},y_{1:t})$. Utilizing this analytical property, the SMC
method for posterior computation can be customized accordingly, as discussed
in Secion 3.2.

In what follows, the discount factors, $\beta$ and $\delta$, are assumed to be
fixed. To calibrate the discount factors, one needs to evaluate the posterior
distribution of the calibration parameters multiple times for various values
of the discount factors, where the fast computation by the SMC method becomes
integral. We will come back to this problem in Section 5

#### Other synthesis functions

Models in the literature of forecast density combination can be viewed as the
synthesis functions in the DBPS framework. For example, Billio et al., (2013)
proposed the SMC method for the time-varying weight (TVW) models with the
corresponding synthesis function of the form,

$\alpha(y_{t}|x_{t},w_{t})=N(y_{t}|w_{t}^{\prime}x_{t},\sigma^{2}),$

where $\sigma^{2}>0$ is constant to be fixed or estimated (and could be
replaced with the stochastic volatility model). The weight parameters,
$w_{t}=(w_{1t},\dots,w_{Kt})^{\prime}$, are restricted to the simplex
($w_{kt}\geq 0$ for $k\in 1:K$ and $\sum_{k}w_{kt}=1$). While this constraint
on $w_{t}$ helps interpretation in the context of model averaging, it is not a
requirement from the original BPS theory and could be a hindrance to adaptive
forecasting and computation. In terms of model flexibility, the non-negativity
and sum-to-unity constraints, as well as the lack of the intercept term
($\theta_{0t}$ in the DLM synthesis), could limit the ability of bias
adjustment in prediction. In computation, unlike the DLM synthesis, the online
posterior and prior under the TVW synthesis are not tractable due to the
nonlinearity in the state variables, which limits the computational method
applicable to this synthesis model. In fact, Billio et al., (2013) used the
vanilla SMC for the online analysis of the TVW synthesis (or the bootstrap
filter in Section 3.1). The comparative analysis of the DLM and TVW syntheses
has been made in McAlinn and West, (2019), which shows the flexibility of the
DLM synthesis improves predictive performance in empirical studies.

The approach using the TVW synthesis has advanced by sophisticating the
synthesis function. Casarin et al., (2023) uses the additive noises with
stochastic volatilities for each agent forecast before taking their weighted
average. The sophisticated synthesis function requires more effort for
computation; the bootstrap filter for this model has been improved by
M-filtering (Baştürk et al.,, 2019). Our study differs from the series of
research on the TVW synthesis in using the DLM synthesis— a simple, linear and
Gaussian model without any hard constraint— and in fully utilizing its
analytical tractability in devising an efficient SMC algorithm.

##  3 Sequential Monte Carlo methods

To obtain the final output, $p(y_{t+1}|\mathcal{H}_{1:(t+1)},y_{1:t})$, we
have to compute the online posterior distribution
$p(\Phi_{1:t},x_{1:t}|\mathcal{H}_{1:t},y_{1:t})$. We first describe the
vanilla SMC method generally applicable to DBPS in the style of sequential
importance resampling (Rubin,, 1988). Then, we provide a custom SMC method for
the DLM synthesis. We also comment on the MCMC method for DBPS and its use as
a remedy for the problem of particle degeneracy.

###  3.1 SMC for DBPS

Denote the pair of calibration parameters and agent forecasts by
$Z_{t}=(\Phi_{1:t},x_{1:t})$. The goal of the SMC method is to update the
online posteriors from $p(Z_{t-1}|\mathcal{H}_{1:(t-1)},y_{1:(t-1)})$ to
$p(Z_{t}|\mathcal{H}_{1:t},y_{1:t})$. Suppose that, at time $t{-}1$, we have
particles $\\{Z_{t-1}^{i},W_{t-1}^{i}\\}_{i=1:M}$ to approximate the online
posterior $p(Z_{t-1}|\mathcal{H}_{1:(t-1)},y_{1:(t-1)})$ by

$\hat{p}^{M}_{t-1}(Z_{t-1})=\sum_{i=1}^{M}W_{t-1}^{i}\delta_{Z_{t-1}^{i}}(Z_{t-1}),$

where $\\{W_{t-1}^{i}\\}$ are weights and satisfy $W_{t-1}^{i}\geq 0$ for
$i=1:M$ and $\sum_{i}W_{t-1}^{i}=1$. Then, our objective is to obtain a
particle approximation of the online posterior distribution at the next time
point, or $p(Z_{t}|\mathcal{H}_{1:t},y_{1:t})$.

In DBPS, $x_{t}$ is independent of $(Z_{t-1},\Phi_{t})$ conditional on
$(y_{1:(t-1)},\mathcal{H}_{1:t})$. Thus, we have

$\begin{split}p(Z_{t}|\mathcal{H}_{1:t},y_{1:t})&\propto
p(y_{t}|Z_{t},\mathcal{H}_{1:t},y_{1:(t-1)})p(Z_{t}|\mathcal{H}_{1:t},y_{1:(t-1)})\\\
&=\alpha(y_{t}|x_{t},\Phi_{t})p(x_{t},\Phi_{t}|Z_{t-1},\mathcal{H}_{1:t},y_{1:(t-1)})p(Z_{t-1}|\mathcal{H}_{1:t},y_{1:(t-1)})\\\
&=\alpha(y_{t}|x_{t},\Phi_{t})p(\Phi_{t}|\Phi_{1:(t-1)},x_{1:(t-1)},y_{1:(t-1)})h_{t}(x_{t})p(Z_{t-1}|\mathcal{H}_{1:(t-1)},y_{1:(t-1)}).\end{split}$

Here, $p(\Phi_{t}|\Phi_{1:(t-1)},x_{1:(t-1)},y_{1:(t-1)})$ is dependent on the
choice of synthesis function $\alpha$, and online posterior
$p(Z_{t-1}|\mathcal{H}_{1:(t-1)},y_{1:(t-1)})$ is approximated by the
aforementioned particles. In approximating the online posterior at time $t$,
we use proposal distribution $q_{t}(x_{t},\Phi_{t})$ as

$p(Z_{t}|\mathcal{H}_{1:t},y_{1:t})\propto\alpha(y_{t}|x_{t},\Phi_{t})\frac{p(\Phi_{t}|\Phi_{1:(t-1)},x_{1:(t-1)},y_{1:(t-1)})h_{t}(x_{t})}{q_{t}(x_{t},\Phi_{t})}q_{t}(x_{t},\Phi_{t})p(Z_{t-1}|\mathcal{H}_{1:(t-1)},y_{1:(t-1)}).$

Based on this expression, we can construct the particle approximation as
follows. First, we sample $Z_{t-1}^{i}$’s from $\hat{p}^{M}_{t-1}(Z_{t-1})$,
by resampling $Z_{t-1}^{1:M}$ with probability $W_{t-1}^{1:M}$. Next, we
generate $(x_{t}^{i},\Phi_{t}^{i})$ from $q_{t}(x_{t},\Phi_{t})$ and calculate
the weight $W_{t}^{i}$ for particle
$Z_{t}^{i}=(Z_{t-1}^{i},x_{t}^{i},\Phi_{t}^{i})$ by

$W_{t}^{i}=\frac{w_{t}^{i}}{\sum_{j=1}^{M}w_{t}^{j}}\ \ \ \mathrm{and}\ \ \
w_{t}^{i}=\alpha(y_{t}|x_{t}^{i},\Phi_{t}^{i})\frac{p(\Phi_{t}^{i}|Z_{t-1}^{i},y_{1:(t-1)})h_{t}(x_{t}^{i})}{q_{t}(x_{t}^{i},\Phi_{t}^{i})}\
\ \text{for }i\in 1:M.$

Finally, the target distribution, $p(Z_{t}|\mathcal{H}_{1:t},y_{1:t})$, can be
approximated by

$\hat{p}_{t}^{M}(Z_{t})=\sum_{i=1}^{M}W_{t}^{i}\delta_{Z_{t}^{i}}(Z_{t}).$

In using this algorithm, the choice of proposal distribution
$q_{t}(x_{t},\Phi_{t})$ is the key. Bootstrap particle filter (BPF, Gordon et
al., 1993) uses
$q_{t}(x_{t},\Phi_{t})=h_{t}(x_{t})p(\Phi_{t}|\Phi_{1:(t-1)},x_{1:(t-1)},y_{1:(t-1)})$
and is applicable to a wide class of synthesis functions. Billio et al.,
(2013) uses BPF for the posterior computation under the TVW model, where
$p(\Phi_{t}|\Phi_{1:(t-1)},x_{1:(t-1)},y_{1:(t-1)})$ is simply a Gaussian
random walk and easy to simulate from. There is a series of research on a
better choice of $q_{t}(x_{t},\Phi_{t})$ including, for example, the auxiliary
particle filter (Pitt and Shephard,, 1999) and the resample-move strategy
(Gilks and Berzuini,, 2001). In the DLM synthesis, we can utilize the
analytical property of the model and derive the custom proposal distribution,
as we will see below.

###  3.2 Rao-Blackwellization for DLM synthesis

In the previous subsection, we approximate the joint online posterior
$p(Z_{t-1}|\mathcal{H}_{1:(t-1)},y_{1:(t-1)})$ to compute the online prior
$p(\Phi_{t}|\mathcal{H}_{1:(t-1)},y_{1:(t-1)})$. In some cases, this
distribution can be evaluated without using the joint online posterior as

$p(\Phi_{t}|\mathcal{H}_{1:(t-1)},y_{1:(t-1)})=\int
p(\Phi_{t}|x_{1:(t-1)},\mathcal{H}_{1:(t-1)},y_{1:(t-1)})p(x_{1:(t-1)}|\mathcal{H}_{1:(t-1)},y_{1:(t-1)})dx_{1:(t-1)}.$

Thus, if the marginal online posterior
$p(x_{1:(t-1)}|\mathcal{H}_{1:(t-1)},y_{1:(t-1)})$ is available, we do not
need the joint online posterior
$p(Z_{1:(t-1)}|\mathcal{H}_{1:(t-1)},y_{1:(t-1)})$. By marginalizing out
$\Phi_{1:t}$, we can improve not only the efficiency of the algorithm but also
the computational time. This approach is known as the Rao-Blackwellized
particle filter (Doucet et al.,, 2000).

Suppose that we have particles $x_{t-1}^{1:M}$ and weights $W_{t-1}^{1:M}$ to
approximate $p(x_{1:(t-1)}|\mathcal{H}_{1:(t-1)},y_{1:(t-1)})$ by
$\hat{p}^{M}_{t-1}(x_{1:(t-1)})$. With the proposal distribution
$q_{t}(x_{t})$, the online posterior distribution of $x_{1:t}$ can be written
as

$\begin{split}p(x_{1:t}|\mathcal{H}_{1:t},y_{1:t})&\propto
p(y_{t}|x_{1:t},y_{1:(t-1)})\frac{h_{t}(x_{t})}{q_{t}(x_{t})}q_{t}(x_{t})p(x_{1:(t-1)}|\mathcal{H}_{1:(t-1)},y_{1:(t-1)}).\end{split}$

Then, generate $x_{t}^{i}\sim q_{t}(x_{t})$ for each $i=1:M$, so that we can
approximate the marginal online posterior,
$p(x_{1:t}|\mathcal{H}_{1:t},y_{1:t})$, by

$\hat{p}^{M}(x_{1:t})=\sum_{i=1}^{M}W_{t}^{i}\delta_{x_{1:t}^{i}}(x_{1:t}),$

where the weights are computed as

$W_{t}^{i}=\frac{w_{t}^{i}}{\sum_{j=1}^{M}w_{t}^{j}}\ \ \ \mathrm{and}\ \ \
w_{t}^{i}=p(y_{t}|x_{1:t}^{i},y_{1:(t-1)})\frac{h_{t}(x_{t}^{i})}{q_{t}(x_{t}^{i})}\
\ \text{for }\ i\in 1:M.$

For the Rao-Blackwellized particle filter to be feasible, density
$p(y_{t}|x_{1:t}^{i},y_{1:(t-1)})$ used in weight $w_{t}^{i}$ must be
evaluated fast.

In the case of the DLM synthesis, the weights can be obtained in closed form
and easily evaluated. As pointed out in Section 2.4, the density of interest,
$p(y_{t}|x_{1:t}^{i},y_{1:(t-1)})$, is the Student’s $t$-distribution; for
details about the functional form, see Section S1. Then, we use
$q_{t}(x_{t})=h_{t}(x_{t})$ as the proposal distribution, so that the weight
is computed simply by $w_{t}^{i}=p(y_{t}|x_{1:t}^{i},y_{1:(t-1)})$. Therefore,
for the DLM synthesis, we only have to draw $x_{t}$ from $h_{t}(x_{t})$ and
set weights $w_{t}^{i}=p(y_{t}|x_{1:t}^{i},y_{1:(t-1)})$.

###  3.3 Interventions by MCMC

####  3.3.1 Gibbs sampler

The Gibbs sampler has been the standard method to compute the online
posteriors in the context of DBPS. The algorithm tailored for the DLM
synthesis is provided in McAlinn and West, (2019) and summarized as follows.
To estimate the posterior distribution at time $T$,
$p(x_{1:T},\Phi_{1:T}|y_{1:T},\mathcal{H}_{1:T})$, we iteratively sample
particles from the following full conditionals:

$\begin{split}\Phi_{1:T}^{n}|x_{1:T}^{n-1},y_{1:T}&\sim
p(\Phi_{1:T}|x_{1:T}^{n-1},y_{1:T})\\\
x_{1:T}^{n}|\Phi_{1:T}^{n},y_{1:T}&\sim\prod_{t=1}^{T}\alpha(y_{t}|x_{t},\Phi_{t}^{n})h_{t}(x_{t}),\end{split}$

for $n\in 1:N$, where $N$ is the chain size of MCMC. In this study, we use all
the historical data, $y_{1:T}$, in implementing the MCMC method. For large
datasets, one could instead use the latest $T_{0}$ observations,
$y_{T-T_{0}-1:T}$, for computational feasibility (Practical filtering; e.g.,
Prado et al., 2021).

In the DLM synthesis, sampling $\Phi_{1:T}^{n}$ from
$p(\Phi_{1:T}|x_{1:T}^{n-1},y_{1:T})$ can be implemented easily by the FFBS
algorithm (Frühwirth-Schnatter, 1994; Carter and Kohn, 1994). In sampling
$x_{1:T}^{n}$, we utilize the fact that all agent forecasts are Student’s
$t$-distributions and expressed as the scale mixture of normal distributions;
$h_{kt}(x_{kt})=\int
N(x_{kt}|\mu_{kt},\sigma_{kt}^{2}H_{kt})p(\sigma_{kt}^{2})d\sigma_{kt}^{2}$,
where $\sigma_{kt}^{2}$ follows some inverse-gamma distribution. Conditional
on latent scale $\sigma_{kt}^{2}$, the posterior distribution of $x_{t}$ is
simply a normal distribution. Thus, we simulate latent scale $\sigma_{kt}^{2}$
together as follows:

$\begin{split}\Phi_{1:T}^{n}|x_{1:T}^{n-1},y_{1:T}&\sim
p(\Phi_{1:T}|x_{1:T},y_{1:T})\ \ \text{by FFBS}\\\
x_{t}^{n}|\Phi_{1:T}^{n},(\sigma_{1:T}^{2})^{n-1},y_{1:T}&\sim\text{Normal}\ \
\text{for}\ t=1:T\\\
(\sigma_{kt}^{2})^{n}|\Phi_{1:T}^{n},x_{1:T}^{n}&\sim\text{Inverse-Gamma}\ \
\text{for}\ k=1:K,t=1:T.\end{split}$

Note that $x_{t}$’s and $\sigma_{kt}^{2}$’s can be sampled in parallel. For
details, see Section S1.

####  3.3.2 MCMC intervention

Particle degeneracy is an inevitable problem in using the SMC methods
repeatedly. As $t$ increases, it is likely to have fewer particle variations,
or increase the variance of the weights, resulting in an inaccurate
approximation of the online posteriors. To avoid this problem, it is advised
that one should monitor an efficiency measure and, once the measure gets below
some pre-specified threshold, then intervene in the algorithm to recover the
accuracy of the approximation (e.g., Gilks and Berzuini, 2001; Chopin, 2002).
We take the effective sample size (ESS) as the monitoring measure which is
computed by using the particles generated by the SMC method as

$\text{ESS}_{t}=\frac{1}{\sum_{m=1}^{M}(W_{t}^{m})^{2}},$

We set some threshold $C>0$ and, if $\text{ESS}_{t}<C$, we discard all the
particles obtained by the SMC method at time $t$ and generate a new set of
particles using the MCMC method. That is, we re-do the approximation of the
online posterior at time $t$ by the Gibbs sampler with chain size $N$. Then,
using the generated chain $\\{x_{1:t}^{1},\cdots,x_{1:t}^{N}\\}$ (and ignoring
the chains of $\Phi_{t}$’s), we approximate the online posterior by

$\hat{p}_{t}^{N}(x_{1:t})=\sum_{n=1}^{N}\frac{1}{N}\delta_{x_{1:t}^{n}}(x_{1:t}).$

This distribution obtained by the MCMC method converges weakly to
$p(x_{1:T}|y_{1:T})$ as in the SMC method, so we use it as the approximation
of the online posterior at time $t$. At the next time point of $t+1$, we
implement the SMC method by using $\hat{p}_{t}^{N}(x_{1:t})$. Thus, the
proposed algorithm is viewed as the combination of the Rao-Blackwellized
particle filter and the Gibbs sampler, as summarized in Algorithm 1.

Note that the chain size $N$ of the MCMC method can differ from the particle
size $M$ of the SMC method. In our application, we set $N=M$ for simplicity.
Since the MCMC method is efficient but time-consuming, one can set $N<M$ in
practice. The threshold $C$ should be chosen carefully. The larger $C$ is, the
more interventions are expected, requiring longer computational time. It is
essential to keep $C$ moderate so that the sequential computation is feasible
while the accuracy of the particle approximation remains satisfactory. We
conducted the sensitivity analysis on the choice of $C$ in the real data
application, reported in Section S2.

Algorithm 1 Rao-Blackwellized Bootstrap Particle Filter with MCMC intervention

for $m=1:M$ do

Generate $x_{1}^{m}\sim h_{1}(x_{1})$

Set $w_{1}^{m}=p(y_{1}|x_{1}^{m})$

end for

Set $W_{1}^{m}=w_{1}^{m}/\sum_{i}w_{1}^{i}$ for $m=1:M$

Set $\hat{p}(x_{1})=\sum_{i=1}^{M}W_{1}^{i}\delta_{x_{1}^{i}}(x_{1})$

for $t=2:T$ do

for $m=1:M$ do

Generate $x_{1:(t-1)}^{m}\sim\hat{p}(x_{1:(t-1)})$

Generate $x_{t}^{m}\sim h_{t}(x_{t})$

Set $w_{t}^{m}=p(y_{t}|x_{1:t}^{m},y_{1:(t-1)})$

end for

Set $W_{t}^{m}=w_{t}^{m}/\sum_{i}w_{t}^{i}$ for $m=1:M$

Set $\hat{p}(x_{1:t})=\sum_{i=1}^{M}W_{t}^{i}\delta_{x_{1:t}^{i}}(x_{1:t})$

if $\text{ESS}_{t}<C$ then

Generate a chain $\\{x_{1:t}^{1},\cdots,x_{1:t}^{N}\\}$ from MCMC with chain
size $N$.

Set $\hat{p}(x_{1:t})=\sum_{i=1}^{N}\delta_{x_{1:t}^{i}}(x_{1:t})/N$

end if

end for

##  4 Numerical evaluation of the proposed SMC methods

In this section, we apply the proposed SMC method to the predictions of the
quarterly US inflation rates. This numerical study aims to confirm that the
predictive performance of the SMC method is as competitive as the MCMC method
with a sufficiently large chain size. Other properties of the SMC method, such
as the necessity of occasional MCMC interventions, are also illustrated in
this example. The settings of this numerical example, including the dataset
and agent models, have been used in McAlinn and West, (2019) to evaluate the
predictive performance of the DBPS.

###  4.1 Settings

#### Data description

We predict the annual percentage change in a chain-weighted GDP price index
$(y_{t})$ as an inflation index. We also use the unemployment rate $(u_{t})$
and the yield on three-month Treasury bills as short-term nominal interest
rate $(r_{t})$ as covariates. We use the data observed from 1961/Q1 to 2022/Q4
($t=1:248$). The data is retrieved at the webpage of the Federal Reserve Bank
of St. Louis (https://research.stlouisfed.org/).

Our numerical study is largely based on McAlinn and West, (2019); we decided
to stick to the univariate data and simple agent models to focus on the
computational problem of our main interest. One notable change from the
experiment in McAlinn and West, (2019) is the use of recent data, especially
those observed in 2020-2022. In this period, when the COVID-19 pandemic hit
the US and world economies, a sudden burst of the inflation rate was observed,
as seen later in Figure 6 (For unemployment and interest rates, see Figure
S1). Sequential inference and prediction by the SMC method became challenging
in this period, as discussed below.

#### Agent models

We use four agent models ($K=4$) to create predictive distributions
$\mathcal{H}_{t}$. All agents are DLMs but with different predictors:
$\mathcal{M}_{1}$ \- $y_{t-1}$ ; $\mathcal{M}_{2}$ \-
$y_{t-3:t-1},u_{t-3:t-1},r_{t-3:t-1}$ ; $\mathcal{M}_{3}$ \- $y_{t-3:t-1}$ ;
$\mathcal{M}_{4}$ \- $y_{t-1},u_{t-1},r_{t-1}$. We use the same
hyperparameters for the four agent models:
$(m_{0},C_{0},n_{0},s_{0},\beta,\delta)=(\mathbf{0},\mathbf{I},2,0.01,0.99,0.95)$,
where $\mathbf{I}$ is the identity matrix. At each $t$, we learn about these
four models from the observed data to obtain the set of predictive densities,
$\mathcal{H}_{t}=\\{h_{1t},\cdots,h_{4t}\\}$, where
$h_{kt}(y_{t})=p(y_{t}|y_{1:t-1},u_{1:t-1},r_{1:t-1},\mathcal{M}_{k})$. Since
the agent models are all DLMs, the predictive densities, $h_{kt}(\cdot)$’s,
are simply the Student’s $t$-distributions and can be computed without any
simulation-based method.

#### DBPS

We use the DLM synthesis in Section 2.4 with hyperparameters
$(m_{0},C_{0},n_{0},s_{0},\beta,\delta)=((0,\mathbf{1}^{\prime}/4)^{\prime},\mathbf{I},$
$10,0.002,0.99,0.95)$ where $\mathbf{1}$ is the 4-dimensional vector with all
elements $1$.

#### Training and test datasets

We divide the whole dataset into the three sub-datasets: the first learning
period from 1961/Q1 to 1977/Q1 $(t=1:65)$, the second learning period from
1977/Q2 to 1989/Q4 $(t=66:116)$, and the evaluation period from 1990/Q1 to
2022/Q4 $(t=117:248)$. In the first learning period, we only learn about the
agent models and predictive distributions. In the second learning period, we
keep learning about the agent models but start to estimate the calibration
parameters in the DLM synthesis. However, we do not evaluate its predictive
performance yet. In the evaluation period, we still keep model learning and
calibration, and we assess the synthesized predictive distribution
$p(y_{t}|y_{66:t-1},\mathcal{H}_{66:t})$ at each time $t$, computing the point
and intervals estimations and other measures introduced below.

#### Computational details

For comparative analysis, we take three approaches to the sequential analysis
of the inflation data. First, we implement the Rao-Blackwellized particle
filter with MCMC interventions. We set the SMC particle size $M=10000$, the
intervention threshold $C=M/20=500$, and the MCMC chain size $N=M=10000$.
Second, we implement the same SMC method without MCMC interventions. In this
fully SMC method, we set the particle size $M=100000$. As we see later, the
problem of particle degeneracy is unavoidable even with this large particle
size. Third, we employ the Gibbs sampler repeatedly at every $t$ as practiced
in the literature. Depending on the purpose of the analysis, we use different
values for chain size $N$.

All the computations are implemented in R on a laptop computer with Intel Core
i7-7700HQ CPU @2.80GHz, 2.80GHz, RAM-8GB.

###  4.2 Computational time and efficiency

Figure 1 shows the raw computational time (log second) of the proposed SMC
method and the repeated use of Gibbs sampler with different chain sizes
($N=100$, $1000$, and $10000$) at each of $t=117:T$. Since the Gibbs sampler
processes all the historical data and sample $(t-66)$ variables, the
computational time increases linearly in $t$. By contrast, the SMC method only
processes the data and parameters at time $t$, whose computational time is
constant over different time points. The advantage of the SMC method in
computation can also be seen in the number of parameters to sample; the Gibbs
sampler generates $(x,\Phi,\sigma^{2})$, while the SMC method needs to sample
$x$ only. For details, see Section S2. The average time of $T-66$ computations
by SMC $(M=10000)$ is 0.56 seconds and the maximum time is 0.76 seconds
(standard deviation is 0.046). The computational time of the MCMC method
$(N=10000)$ at $t=248$ is 363.5 seconds, while it is 116.36 seconds at
$t=117$.

Figure 1: Logarithmic computational time (seconds) over the evaluation period
by the MCMC method with $N=100$, $1000$, $10000$ (dotted, from bottom to top)
and the SMC method with $M=10000$ (solid).

In Figure 2, the realized ESS are plotted over $t=66:248$. The figure also
indicates time points when $\text{ESS}_{t}<C=500$ and the MCMC intervention
with size $N=10000$ is triggered. With this threshold, we made the
intervention five times out of $183\ (t=66:248)$ time points. The total time
for this computation is 1266.79 seconds, mainly for the five MCMC
interventions. Sensitivity analysis about different values of threshold $C$
can be found in Section S2, where the number of interventions, computational
time, and predictive performance are reported.

Figure 2: Effective sample size (ESS) over $t=66:248$ (solid). The starting
point of the evaluation period ($t=117$) is indicated by the vertical solid
line. The five MCMC interventions with chain size $N=10000$ (time points when
$\text{ESS}_{t}<C$) are indicated by the vertical dotted lines. The threshold
is set to $C=500$ and indicated by the horizontal dashed line.

###  4.3 Approximation accuracy

We assess the accuracy of particle approximation made by the SMC methods for
the synthesized predictive distributions,
$p(y_{t}|y_{66:t-1},\mathcal{H}_{66:t})$, and the posterior distributions of
the calibration parameters, $p(\Phi_{t}|y_{66:t-1},\mathcal{H}_{66:t-1})$, for
$t=117:248$. First, we employ the MCMC method repeatedly, using chain size
$N=70000$ at each time, to obtain the accurate approximation of the posterior
and predictive distributions. Then, we implement the computation by the SMC
methods with and without the MCMC interventions and compare their results with
the repeated MCMC method.

In Figure 3, we plot the posterior distributions of the location calibration
parameter $\theta_{t,0}$, or
$p(\theta_{t,0}|y_{66,t-1},\mathcal{H}_{66,t-1})$, for $t=117:248$. Without
MCMC interventions, the posteriors computed by the generated particles are
biased and degenerated. The fully SMC-based approach is also sensitive to the
sudden burst of the inflation rate at $t=242$ (2021/Q2). At this time point,
both the agent models and the DLM synthesis significantly fail in prediction,
causing a severe particle degeneracy in the sequential analysis. With
interventions, the SMC method accurately approximates the target predictive
distribution; their median and 90% credible intervals are almost identical to
those computed by the MCMC method with the large chain size. We confirmed that
the SMC method with interventions can achieve the accurate particle
approximation in other posterior plots; see Section S2 for details.

Figure 3: The $5\%,50\%,$ and $95\%$ quantiles of posterior distributions
$p(\theta_{t,0}|y_{66:t-1},\mathcal{H}_{66:t-1})$ computed by MCMC (solid) and
SMC (dashed). The top panel shows the result of the fully SMC-based method,
while the bottom allows the MCMC interventions. Figure 4: Log predictive
density rations (LPDRs) of the SMC methods with and without MCMC interventions
(solid and dotted, respectively) and four agent models (dash-dotted) against
the MCMC method. The LPDR of the SMC method with interventions is close to
zero, showing its accuracy of particle approximation. The fully SMC-based
approach shows decreasing LPDRs due to the accumulated approximation errors.

To evaluate the SMC method from the viewpoint of predictive performance, we
compute the log predictive density ratios (LPDR) over $t=117:248$, which is
defined by

$\text{LPDR}_{t}=\sum_{s=117}^{t}\log(p(y_{s}|y_{66:s-1},\mathcal{H}_{66:s})/p_{\text{MCMC}}(y_{s}|y_{66:s-1},\mathcal{H}_{66:s}))$

where $p_{\text{MCMC}}(y_{s}|y_{66:s-1},\mathcal{H}_{66:s})$ is the predictive
density computed by using the MCMC method repeatedly, and
$p(y_{s}|y_{66:s-1},\mathcal{H}_{66:s})$ is the predictive density of
interest. If the LPDR is close to $0$, the predictive distribution of interest
is as performative as the synthesized predictive distribution computed by the
MCMC method. Here, we evaluate the LPDR for the particle approximation of the
synthesized predictive density computed by the SMC methods with or without
MCMC interventions and the four agent predictive densities. The LPDRs are
plotted in Figure 4. Again, without interventions, the approximation errors
accumulate gradually over time and critically at the sudden burst of the
inflation rate. The MCMC interventions improve approximation accuracy
regarding LPDRs at all time points.

##  5 Application to the sequential loss discounting

In the previous section, we confirmed that the proposed SMC method could
accurately approximate the posterior and predictive distributions. However,
this does not mean that the synthesized predictive distribution is always
useful. Rather, the synthesized predictive distribution performs poorly in
2020-2022, regardless of the computational methods. This is all due to the
sudden burst of the inflation rate, to which all four agent models are slow to
adapt. It is also true for the calibration parameters in the DLM synthesis
that adapting to this sudden change is difficult, due to the use of the
Gaussian random walk as the prior.

To make the predictive analysis more flexible, we should not fix the discount
factors, $\beta$ and $\delta$, in the DLM synthesis, but adaptively change
their values. In doing so, one can take a decision-theoretic approach and make
a manual intervention to the sequential analysis (West and Harrison, 1997,
Chapter 11). Although practiced in applied statistics (e.g., Chen et al.,
2018), this approach involves many subjective decisions, such as deciding the
criterion for interventions. Using priors for discount factors is a formal
Bayesian approach (e.g., Irie et al., 2022), but involves more intensive use
of the SMC methods (e.g., Liu and West, 2001). As an alternative easily
implemented with fewer manual interventions, one can calibrate discount
factors— trying different values of discount factors and comparing the
resulting synthesized predictive distributions. This approach requires the
computation of the synthesized predictive distribution multiple times, for
which fast computation by the proposed SMC method is crucial.

###  5.1 Loss Discounting Framework

In the presence of possible structural changes, power discounting has been
utilized in combining the predictive densities (e.g., Zhao et al., 2016). This
idea is formalized as the loss discounting framework (LDF) in the context of
forecast combination (Bernaciak and Griffin,, 2022). Suppose that we have $K$
agent models and their predictive densities,
$h_{kt}(y_{t})=p(y_{t}|y_{1:t-1},\mathcal{M}_{k})$, for $k=1:K$. In evaluating
the historical performance of the agent models, the LDF uses log-discounted
predictive likelihood (LDPL) defined by

$\text{LDPL}_{k,t}(\gamma_{1})=\sum_{s=1}^{t}\gamma_{1}^{(t-s)}\log(p(y_{s}|y_{1:s-1},\mathcal{M}_{k})),$

where $\gamma_{1}\in(0,1]$ is called the first layer discount factor. If
$\gamma_{1}=1$, then the LDPL reduces to the log-marginal likelihood of agent
model $\mathcal{M}_{k}$, which is used in the standard Bayesian analysis. The
LDPL rates the recent predictive performance of the agent model higher than
those in the past, making the model evaluation sensitive to sudden structural
changes. Using the LDPLs, the $K$ predictive densities are averaged as

$\begin{split}p_{\mathrm{LDF}}(y_{t}|y_{1:t-1},\mathcal{H}_{1:t},\gamma_{1})=\sum_{k=1}^{K}w_{k}(\text{LDPL}_{1:K,t-1}(\gamma_{1}))p(y_{t}|y_{1:t-1},\mathcal{M}_{k}),\end{split}$

where the weight function, $w_{k}(\cdot)$, depends on the LDPLs. If $w_{k}$ is
the expit function, the LDF reduces to the dynamic model averaging (Raftery et
al.,, 2010) and, in the case of $\gamma_{1}=1$, it is the Bayesian model
averaging. Bernaciak and Griffin, (2022) proposes two weight functions:
softmax and argmax weights defined by

$\begin{split}\text{softmax-
weights}(a_{1},\cdots,a_{K})&=\left(\frac{\exp(a_{1})}{\sum_{k}\exp(a_{k})},\cdots,\frac{\exp(a_{K})}{\sum_{k}\exp(a_{k})}\right),\\\
\text{argmax-
weights}(a_{1},\cdots,a_{K})&=\Big{(}0,\cdots,\underbrace{1}_{k_{\max}\text{-th}},\cdots,0\Big{)},\
\ \ \ \ a_{k_{\max}}\geq a_{k},\ \ \ \ \mathrm{for\ all\ }k.\end{split}$

The softmax weights utilize all the $K$ models as in the model averaging,
while the argmax weights choose the best model as the final output for
prediction. Denote the LDFs with the softmax and argmax weight functions by
$\text{LDF}_{s}$ and $\text{LDF}_{a}$, respectively.

The first-layer discount factor, $\gamma_{1}$, is an important tuning
parameter to control how much we discount the past predictive performance of
the agent models. To emphasize the dependence on $\gamma_{1}$, we write the
predictive distribution in the LDF as
$p_{\mathrm{LDF}_{\ast}}(y_{t}|y_{1:t-1},\mathcal{H}_{1:t},\gamma_{1})$, where
$\ast\in\\{s,a\\}$. To compare different values of $\gamma_{1}$, we have to
obtain multiple predictive distributions with different values of
$\gamma_{1}$. These distributions need to be combined in the LDF again. To be
precise, we introduce the second-layer discount factor $\gamma_{2}\in[0,1]$,
evaluate the LDPL of
$p_{\mathrm{LDF}_{\ast}}(y_{t}|y_{1:t-1},\mathcal{H}_{1:t},\gamma_{1})$ for
each value of $\gamma_{1}$, and combine the predictive densities using a
weight function that is not necessarily the same as one used in the first
layer combination. We denote this forecast combination by
$\mathrm{LDF}_{\ast_{1},\ast_{2}}$, where $\ast_{1},\ast_{2}\in\\{s,a\\}$. For
example, if we use the softmax weights to combine the predictions made by
$\mathrm{LDF}_{a}$, then it is $\mathrm{LDF}_{a,s}$.

Finally, note that we still have to choose the second-layer discount factor,
$\gamma_{2}$ carefully. In this study, we set $\gamma_{2}$ to a fixed value.
Having extra layers is possible in the general LDF, but is expected to have
less impact on the final predictive result.

###  5.2 LDF for BPS

The LDF can be implemented quickly and fast, but its flexibility is limited
since the resulting density is simply the weighted average of $K$ agent
predictive densities. Here, we use the idea of the LDF in calibrating the
discount factors of the DLM synthesis function of the DBPS, namely, $\beta$
and $\delta$.

Denote the set of $J$ possible values of discount factors by
$S=\\{(\beta_{1},\delta_{1}),\dots,(\beta_{J},\delta_{J})\\}$. At each time
$t-1$, for each pair $(\beta_{j},\delta_{j})\in S$, we compute the synthesized
predictive distribution,
$p(y_{t}|\mathcal{H}_{1:t},y_{1:(t-1)},\beta_{j},\delta_{j})$. Consequently,
we have $J$ synthesized predictive distributions, which are combined by using
the LDF discount factor, $\gamma$, and the softmax or argmax weight function
with argument $\text{LDPL}_{j,t}(\gamma)$ computed by
$p(y_{t}|\mathcal{H}_{1:t},y_{1:(t-1)},\beta_{j},\delta_{j})$. We denote this
approach by $\text{LDF}_{B,\ast}$, where $\ast\in\\{a,s\\}$. From the
viewpoint of the two-layer LDF, this approach replaces the first-layer density
combination with the DBPS.

In this approach, we must evaluate the $J$ synthesized predictive
distributions at each time point, which could be costly when the computational
resource is limited. If these computations are not fully parallelized, then
using the MCMC method in a butch process could increase the computational
cost, making this approach infeasible. The fast computation by the proposed
SMC method is beneficial in this context.

###  5.3 Real Data analysis

Here, we again analyze the US inflation rates using the LDF to calibrate the
discount factors. The same agent models and hyperparameters are used, but this
time we consider $J=35$ values for the discount factors, $(\beta,\delta)$, in
the DLM synthesis function. The 35 synthesized predictive densities are
combined in the LDF with discount factor $\gamma=0.98$. Likewise, we apply the
two-layer LDF approach to the same data/agent models with 15 values for
$\gamma_{1}$ in the first layer and $\gamma_{2}=0.98$ in the second layer. For
the complete list of the values used for the discount factors, see Section S3.

In computation, we set particle/chain sizes to $N=M=10000$ and the
intervention threshold to $C=700$. The synthesized predictive distributions
are computed in the same computational environment described in Section 4.1.
Note that we cannot fully parallelize the computation of $J=35$ synthesized
predictive distributions by using the 8-core CPUs. The raw computational time
is about 5 hours in total. For faster computation, one could use fewer values
of the discount factor (smaller $J$), especially in using the MCMC method as
the intervention.

Figure 5 shows the LPDRs of the $\text{LDF}_{B,a}$, $\text{LDF}_{B,s}$ and
$\text{LDF}_{s,a}$ against the benchmark BPS with fixed discount factors
$(\beta,\delta)=(0.99,0.95)$. Both the $\text{LDF}_{B,a}$ and
$\text{LDF}_{B,s}$ improve the predictive performance of the BPS with the
fixed discount factors, especially in the period of the sudden burst in
2020-2022. The $\text{LDF}_{B,a}$ is slightly better than the
$\text{LDF}_{B,s}$; in our example, using all of the 35 discount factors might
have overestimated the predictive uncertainty. Even with different choices of
$\gamma$, the $\text{LDF}_{B,a}$ and $\text{LDF}_{B,s}$ consistently
outperform the benchmark DBPS, but an extreme value of $\gamma$, such as
$\gamma=1$, is not recommended. For details, see Section S3.

Figure 5: LPDRs in the three loss discounting frameworks against the benchmark
DBPS with the fixed discount factors: $\text{LDF}_{B,a}$ (solid),
$\text{LDF}_{B,s}$ (dashed) and $\text{LDF}_{s,a}$ (dotted). Both
$\text{LDF}_{B,a}$ and $\text{LDF}_{B,s}$ improve the prediction by the
benchmark DBPS. $\text{LDF}_{s,a}$ is worse than the benchmark, although it
shows an improvement of the predictive performance in 2020-2022.

By contrast, the $\text{LDF}_{s,a}$ is far less performative than the
benchmark because the original LDF is simply the weighted average of the agent
models and less flexible than the DBPS. However, the advantage of the LDF in
discounting the agent predictions in the past can be seen clearly in the
increase of its LPDR in 2020-22. Other LDFs provide similar predictions, which
are reported in Section S3.

The resulting predictive distributions of the $\text{LDF}_{B,a}$ and benchmark
DBPS are shown in Figure 6. The $\text{LDF}_{B,a}$ benefits from the slightly
narrower predictive intervals when the dynamics of the inflation rates are
relatively stable (e.g., in 1990-2000) and from the wider intervals when the
inflation rates are volatile (e.g., after 2010).

Figure 6: Predictive medians and 90% predictive intervals of DBPS and
$\text{LDF}_{B,a}$ with $\gamma=0.98$ laid over the data in the period
$t=117:248$ (top) and $t=117:236$ (bottom). $\text{LDF}_{B,a}$ provides
diffuse predictive distributions after 2020 due to lower discount factors used
for the adaptation to a sudden structural change. The predictive intervals in
the $\text{LDF}_{B,a}$ are slightly narrower in 1990-2000 and wider in
2010-2020 than the benchmark DBPS.

##  6 Future research

The proposed SMC method for DBPS with appropriate interventions by the Gibbs
sampler enables the fast computation of the online posterior and synthesized
predictive distribution. The fast computation by the SMC method is crucial in
not only real-time monitoring but also situations where the online posteriors
must be computed multiple times. As an example of the latter cases, we
considered the loss discounting DBPS, where the synthesis function has the
discount factors to be tuned at each time point. Another important example not
studied in this paper is the synthesis of many predictive densities, or the
“large $K$” problem. In this situation, the synthesized predictive
distribution could be less performative since we have many calibration
parameters to estimate. In the DLM synthesis, selecting predictive densities
in the synthesis function is equivalent to the variable selection problem in
the high-dimensional regression. The study of this “model selection” problem
in the context of BPS has been hindered due to the computational cost of
evaluating the posterior of the calibration parameters, limiting the feasible
approaches to using shrinkage priors (Chernis,, 2022) and clustering
(Kobayashi et al.,, 2023). The variable selection approach is interesting
future research that we believe is enabled by the proposed SMC method.

The Rao-Blackwellized particle filter proposed in this study is a custom SMC
method for the DLM synthesis. Although many other SMC methods could be used
for sequential forecasting by the DBPS, the problem of particle degeneracy is
most likely unavoidable, especially in the period of sudden bursts/structural
changes. Developing custom computational algorithms for other synthesis
functions, such as mixture synthesis (Johnson and West,, 2023), is undoubtedly
important future research.

## Acknowledgments

We thank Naoki Awaya for his comments and suggestions. The second author’s
research was partly supported by JSPS KAKENHI Grant Number 22K20132 from Japan
Society for the Promotion of Science.

## References

  * Aastveit et al., (2023) Aastveit, K. A., Cross, J. L., and van Dijk, H. K. (2023).  Quantifying time-varying forecast uncertainty and risk for the real price of oil.  Journal of Business & Economic Statistics, 41(2):523–537.
  * Baştürk et al., (2019) Baştürk, N., Borowska, A., Grassi, S., Hoogerheide, L., and van Dijk, H. K. (2019).  Forecast density combinations of dynamic models and data driven portfolio strategies.  Journal of Econometrics, 210(1):170–186.
  * Bernaciak and Griffin, (2022) Bernaciak, D. and Griffin, J. E. (2022).  A loss discounting framework for model averaging and selection in time series models.  arXiv preprint arXiv:2201.12045.
  * Billio et al., (2013) Billio, M., Casarin, R., Ravazzolo, F., and Van Dijk, H. K. (2013).  Time-varying combinations of predictive densities using nonlinear filtering.  Journal of Econometrics, 177(2):213–232.
  * Cabel et al., (2022) Cabel, D., Kato, M., McAlinn, K., Sugasawa, S., and Takanashi, K. (2022).  Spatially-varying Bayesian predictive synthesis for flexible and interpretable spatial prediction.  arXiv preprint arXiv:2203.05197.
  * Carter and Kohn, (1994) Carter, C. K. and Kohn, R. (1994).  On Gibbs sampling for state space models.  Biometrika, 81(3):541–553.
  * Casarin et al., (2015) Casarin, R., Grassi, S., Ravazzolo, F., van Dijk, H., et al. (2015).  Parallel sequential Monte Carlo for efficient density combination: The DeCo MATLAB toolbox.  JOURNAL OF STATISTICAL SOFTWARE, 68.
  * Casarin et al., (2023) Casarin, R., Grassi, S., Ravazzolo, F., and van Dijk, H. K. (2023).  A flexible predictive density combination for large financial data sets in regular and crisis periods.  Journal of Econometrics.
  * Chen and Liu, (2000) Chen, R. and Liu, J. S. (2000).  Mixture kalman filters.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 62(3):493–508.
  * Chen et al., (2018) Chen, X., Irie, K., Banks, D., Haslinger, R., Thomas, J., and West, M. (2018).  Scalable Bayesian modeling, monitoring, and analysis of dynamic network flow data.  Journal of the American Statistical Association, 113(522):519–533.
  * Chernis, (2022) Chernis, T. (2022).  Combining large numbers of density predictions with Bayesian predictive synthesis.  Preprint.
  * Chopin, (2002) Chopin, N. (2002).  A sequential particle filter method for static models.  Biometrika, 89(3):539–552.
  * Doucet et al., (2000) Doucet, A., Godsill, S., and Andrieu, C. (2000).  On sequential Monte Carlo sampling methods for Bayesian filtering.  Statistics and computing, 10:197–208.
  * Frühwirth-Schnatter, (1994) Frühwirth-Schnatter, S. (1994).  Data augmentation and dynamic linear models.  Journal of time series analysis, 15(2):183–202.
  * Genest and Schervish, (1985) Genest, C. and Schervish, M. J. (1985).  Modeling expert judgments for Bayesian updating.  The Annals of Statistics, pages 1198–1212.
  * Gilks and Berzuini, (2001) Gilks, W. R. and Berzuini, C. (2001).  Following a moving target—Monte Carlo inference for dynamic Bayesian models.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 63(1):127–146.
  * Gordon et al., (1993) Gordon, N. J., Salmond, D. J., and Smith, A. F. (1993).  Novel approach to nonlinear/non-Gaussian Bayesian state estimation.  In IEE proceedings F (radar and signal processing), volume 140(2), pages 107–113. IET.
  * Irie et al., (2022) Irie, K., Glynn, C., and Aktekin, T. (2022).  Sequential modeling, monitoring, and forecasting of streaming web traffic data.  The Annals of Applied Statistics, 16(1):300–325.
  * Johnson and West, (2023) Johnson, M. C. and West, M. (2023).  Bayesian predictive synthesis with outcome-dependent pools.  arXiv preprint.  arXiv:1803.01984.
  * Kobayashi et al., (2023) Kobayashi, G., Sugasawa, S., Kawakubo, Y., Han, D., and Choi, T. (2023).  Clustering and predicting multiple time series count data via mixture of Bayesian predictive syntheses: Analysis of COVID-19 hospitalisation in Japan and Korea.  arXiv preprint.  arXiv:2308.06134.
  * Liu and West, (2001) Liu, J. and West, M. (2001).  Combined parameter and state estimation in simulation-based filtering.  In Sequential Monte Carlo methods in practice, pages 197–223. Springer.
  * McAlinn, (2021) McAlinn, K. (2021).  Mixed-frequency Bayesian predictive synthesis for economic nowcasting.  Journal of the Royal Statistical Society: Series C (Applied Statistics).
  * McAlinn et al., (2020) McAlinn, K., Aastveit, K. A., Nakajima, J., and West, M. (2020).  Multivariate Bayesian predictive synthesis in macroeconomic forecasting.  Journal of the American Statistical Association, 115(531):1092–1110.
  * McAlinn and West, (2019) McAlinn, K. and West, M. (2019).  Dynamic Bayesian predictive synthesis in time series forecasting.  Journal of Econometrics, 210(1):155–169.
  * Pitt and Shephard, (1999) Pitt, M. K. and Shephard, N. (1999).  Filtering via simulation: Auxiliary particle filters.  Journal of the American statistical association, 94(446):590–599.
  * Prado et al., (2021) Prado, R., Ferreira, M. A., and West, M. (2021).  Time series: modeling, computation, and inference.  CRC press.
  * Raftery et al., (2010) Raftery, A. E., Kárnỳ, M., and Ettler, P. (2010).  Online prediction under model uncertainty via dynamic model averaging: Application to a cold rolling mill.  Technometrics, 52(1):52–66.
  * Rubin, (1988) Rubin, D. (1988).  Using the sir algorithm to simulate posterior distributions.  Bayesian Statistics, 3:395–402.
  * Tallman and West, (2022) Tallman, E. and West, M. (2022).  Bayesian predictive decision synthesis.  Submitted for publication.  arXiv:2206.03815.
  * West, (1992) West, M. (1992).  Modelling agent forecast distributions.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 54:553–567.
  * West and Crosse, (1992) West, M. and Crosse, J. (1992).  Modelling probabilistic agent opinion.  Journal of the Royal Statistical Society: Series B (Methodological), 54(1):285–299.
  * West and Harrison, (1997) West, M. and Harrison, P. J. (1997).  Bayesian Forecasting and Dynamic Models.  Springer Verlag, 2nd edition.
  * Zhao et al., (2016) Zhao, Z. Y., Xie, M., and West, M. (2016).  Dynamic dependence networks: Financial time series forecasting and portfolio decisions.  Applied Stochastic Models in Business and Industry, 32(3):311–332.

Supplementary Materials for
“Sequential Bayesian Predictive Synthesis”

##  Appendix S1 DLM synthesis

Here, we describe the details of DLMs, including the forward filter to compute
the online posterior and predictive distributions. Note that the DLMs are used
as both the agent models and synthesis functions. For further discussions, see
Prado et al., (2021).

A DLM with discount factors $\beta,\delta\in(0,1]$ is defined as

$\begin{split}&y_{t}=F_{t}^{\prime}\theta_{t}+\epsilon_{t}\ ,\
\epsilon_{t}\sim N(0,\nu_{t})\\\ &\theta_{t}=\theta_{t-1}+\omega_{t}\ ,\
\omega_{t}\sim N(0,\nu_{t}W_{t})\\\ &\nu_{t}=\frac{\beta}{\gamma_{t}}\ ,\
\gamma_{t}\sim
Be\left(\beta\frac{n_{t-1}}{2},(1-\beta)\frac{n_{t-1}}{2}\right)\end{split}$

where
$F_{t}=(1,x_{1t},\cdots,x_{Kt})^{\prime},\theta_{t}=(\theta_{0t},\theta_{1t}\cdots,\theta_{Kt})^{\prime},$
and $Be(a,b)$ denotes beta distribution. $W_{t}$ is defined implicitly so that
the following relationship holds;

$\begin{split}&p(\theta_{t-1}|\nu_{t-1},x_{1:(t-1)},y_{1:(t-1)})=N\left(\theta_{t-1}\left|m_{t-1},\frac{\nu_{t-1}}{s_{t-1}}C_{t-1}\right)\right.\\\
&p(\nu_{t-1}|x_{1:(t-1)},y_{1:(t-1)})=IG\left(\nu_{t-1}\left|\frac{n_{t-1}}{2},\frac{n_{t-1}s_{t-1}}{2}\right)\right.\\\
&p(\theta_{t}|\nu_{t},x_{1:(t-1)},y_{1:(t-1)})=N\left(\theta_{t}\left|a_{t},\frac{\nu_{t}}{s_{t-1}}R_{t}\right)\right.\\\
&p(\nu_{t}|x_{1:(t-1)},y_{1:(t-1)})=IG\left(\nu_{t}\left|\frac{r_{t}}{2},\frac{r_{t}s_{t-1}}{2}\right)\right.\\\
\end{split}$

where $IG(a,b)$ denotes inverse-gamma distribution. Sufficient statistics in
the above equations are calculated recursively as follows;

$\begin{split}\text{posterior to prior}&\\\ a_{t}&=m_{t-1}\ ,\
R_{t}=\frac{1}{\delta}C_{t-1}\ ,\ r_{t}=\beta n_{t-1}\\\ \text{prior to
posterior}&\\\ m_{t}&=a_{t}+A_{t}\ ,\
C_{t}=(R_{t}-A_{t}A_{t}^{\prime}q_{t})z_{t}\ ,\ n_{t}=r_{t}+1\ ,\
s_{t}=s_{t-1}z_{t}\\\ \text{where}&\quad
q_{t}=s_{t-1}+F_{t}^{\prime}R_{t}F_{t}\ ,\
z_{t}=(r_{t}+(y_{t}-F_{t}^{\prime}a_{t})^{2}/q_{t})/(r_{t}+1)\ ,\
A_{t}=R_{t}F_{t}/q_{t}\end{split}$

The above procedure is called Forward Filtering (e.g. Prado et al., 2021). To
implement this model, it is needed to set the initial distribution of the
latent variables;

$p(\theta_{1},\nu_{1})=N\left(\theta_{1}\left|a_{1},\frac{\nu_{1}}{s_{0}}R_{1}\right)\right.IG\left(\nu_{1}\left|\frac{r_{1}}{2},\frac{r_{1}s_{0}}{2}\right)\right.,$

and it can be seen that we have to set $(m_{0},C_{0},n_{0},s_{0})$. Therefore,
the all hyperparameters in this DLM are
$(m_{0},C_{0},n_{0},s_{0},\beta,\delta)$. The predictive distribution is
Stundent’s t-distribution with the degree of freedom $r_{t}$;

$p(y_{t}|x_{1:t},y_{1:(t-1)})=t(y_{t}|r_{t},f_{t},q_{t})$

where $f_{t}=F_{t}^{\prime}a_{t}$ is the location parameter and $\sqrt{q_{t}}$
is the scale parameter. The density of this distribution is given by

$p(y_{t}|x_{1:t},y_{1:(t-1)})=(r_{t}s_{t-1})^{-\frac{1}{2}}\
B\left(\frac{1}{2},\frac{r_{t}}{2}\right)^{-1}\left(1+\frac{1}{s_{t-1}}F_{t}^{\prime}R_{t}F_{t}\right)^{-\frac{1}{2}}\left(1+\frac{1}{r_{t}}\frac{(y_{t}-a_{t}^{\prime}F_{t})^{2}}{s_{t-1}+F_{t}^{\prime}R_{t}F_{t}}\right)^{-\frac{1+r_{t}}{2}},$

where $B(a,b)$ is the beta function. Conditional on forecasters $x_{1:T}$, all
online posterior, online prior, and predictive distributions are computed
analytically by Forward Filtering.

The offline posterior distribution $p(\theta_{1:t},\nu_{1:t}|y_{1:t},x_{1:t})$
can be approximated by particles. Forward Filtering and Backward Sampling
(FFBS) algorithm (e.g., see Prado et al., 2021) generates independent
particles directly from this distribution. FFBS algorithm at time $T$ is
summarized as follows;

  * 1.

Forward Filtering: Calculate sufficient statistics $(m_{t},C_{t},n_{t},s_{t})$
for $t=1:T$ in sequence.

  * 2.

Backward sampling: Sample $\Phi_{T}^{1:N}$ from the online posterior
distribution;

$\begin{split}\nu_{T}^{n}&\sim
IG\left(\frac{n_{T}}{2},\frac{n_{T}s_{T}}{2}\right)\\\ \theta_{T}^{n}&\sim
N\left(m_{T},\frac{\nu_{T}^{n}}{s_{T}}C_{T}\right)\ \ \text{for
}n=1:N,\end{split}$

and after that, simulate $\Phi_{t}^{1:N}$ for $t=T-1:1$ as follows;

$\begin{split}(\nu_{t}^{n})^{-1}&=\beta(\nu_{t+1}^{n})^{-1}+\gamma_{t}^{n}\ ,\
\gamma_{t}^{n}\sim
Ga\left((1-\beta)\frac{n_{t}}{2},\frac{n_{t}s_{t}}{2}\right)\\\
\theta_{t}^{n}&\sim
N\left(m_{t}+\delta(\theta_{t+1}^{n}-m_{t}),(1-\delta)\frac{\nu_{t}^{n}}{s_{t}}C_{t}\right)\
\ \text{for }n=1:N,\end{split}$

where $Ga(a,b)$ denotes gamma distribution. Particles $\Phi_{1:T}^{1:N}$
generated by the FFBS algorithm approximate the offline posterior
distribution.

Since agent predictive distributions are all Student’s t-distributions, they
can be expressed as the scale mixture of normal distributions;

$\begin{split}h_{kt}(x_{kt})&=t(x_{kt}|e_{kt},\mu_{kt},H_{kt})\\\ &=\int
N(x_{kt}|\mu_{kt},\sigma_{kt}^{2}H_{kt})IG\left(\sigma_{kt}^{2}\left|\frac{e_{kt}}{2},\frac{e_{kt}}{2}\right)\right.d\sigma_{kt}^{2}\end{split}$

where $e_{kt}$ is the degree of freedom, $\mu_{kt}$ is the location parameter,
and $\sqrt{H_{kt}}$ is the scale parameter. In this setting, we can implement
Gibbs Sampling at time $T$ as follows;

$\begin{split}\Phi_{1:T}&\sim p(\Phi_{1:T}|x_{1:T},y_{1:T})\ \ \text{by
FFBS}\\\ x_{t}&\sim N(\mu_{t}+b_{t}c_{t},H_{t}-b_{t}b_{t}^{\prime}g_{t})\ \
\text{for}\ t=1:T\\\ \sigma_{kt}^{2}&\sim
IG\left(\frac{e_{kt}+1}{2},\frac{e_{kt}+d_{kt}}{2}\right)\ \ \text{for}\
k=1:K,t=1:T,\end{split}$

where
$H_{t}=\text{diag}(\sigma_{1t}^{2}H_{1t},\cdots,\sigma_{Kt}^{2}H_{Kt}),\mu_{t}=(\mu_{1t},\cdots,\mu_{Kt})^{\prime},c_{t}=y_{t}-\theta_{t0}-\mu_{t}^{\prime}\theta_{t,1:K},g_{t}=\nu_{t}+\theta_{t,1:K}^{\prime}H_{t}\theta_{t,1:K},b_{t}=H_{t}\theta_{t,1:K}/g_{t}$
and $d_{kt}=(x_{kt}-\mu_{kt})^{2}/H_{kt}$.

##  Appendix S2 Additional Results for Section 4

###  S2.1 Data description

The data on the US inflation, unemployment, and interest rates is retrieved at
the webpage of the Federal Reserve Bank of St. Louis
(https://research.stlouisfed.org/series/). The unemployment and interest rates
are recorded monthly, so we create the quarterly data by using the values of
the last month in each quarter. Figure S1 shows the time series of the three
variables.

Figure S1: US inflation rate (top), unemployment rate (middle) and interest
rate (bottom).

###  S2.2 Additional posterior plots

#### Comparison with the MCMC method

In the main text, we illustrated the approximation accuracy of the proposed
SMC method by confirming that the marginal posteriors computed by the SMC and
MCMC methods become similar. More marginal plots are shown in Figure S2,
showing the accuracy of the proposed SMC method with interventions. The
posteriors approximated without MCMC interventions are biased, especially in
estimating variance $\nu_{t}$, and tend to underestimate the posterior
uncertainty.

Figure S2: The $5\%,50\%,$ and $95\%$ quantiles of posterior distributions
$p(\Phi_{t}|y_{66:t-1},\mathcal{H}_{66:t-1})$ computed by MCMC (solid), SMC
without interventions (dotted) and SMC with interventions (dashed).

Here, we show in Figure S3, the two-dimensional scatter plots of generated
$(\theta_{2t},\theta_{4t})$ at $t=247$. These joint posterior distributions
are skewed in a similar way, which supports our claim on the accuracy of the
SMC method.

Figure S3: Scatter plots of the joint posterior distribution
$p(\theta_{2t},\theta_{4t}|y_{66:t-1},\mathcal{H}_{66:t-1})$ for $t=247$
computed by MCMC (top) and SMC with interventions (bottom). The horizontal
axis is $\theta_{2t}$ and the vertical axis is $\theta_{4t}$.

#### Posteriors of $\theta_{t}$

Figure S4 shows the posterior means of calibration coefficients $\theta_{t}$.
Note that this is not a weight vector, so does not sum to unity and can be
negative. The highest coefficient is $\theta_{2t}$, the coefficient of
$\mathcal{M}_{2}$ that is most flexible with 9 covariates. In 2020-2022, the
coefficient of $\theta_{4t}$ increases in response to the sudden burst of
inflation rates. Also, in the same period, the intercept, $\theta_{0t}$,
frequently changes to adjust the increased predictive bias.

Figure S4: Means of $p(\theta_{t}|y_{66:t-1},\mathcal{H}_{66:t-1})$ for
$t=117:248$ computed by MCMC $(N=70000)$. The posterior mean of $\theta_{3t}$
decreases gradually in 2005-2010, while $\theta_{2t}$ becomes more dominant in
the same period.

###  S2.3 Rough estimate of raw computational time

We can estimate the computational time of the SMC/MCMC methods by counting the
number of parameters to be sampled. Remember that we start the BPS at $t=66$
until $t=T=248$, and use $M$ particles in the SMC method and the chain of
length $N$ in the MCMC method. In implementing the Gibbs Sampler, we generate
three latent variables $\\{x_{t},\Phi_{t},\sigma_{t}^{2}\\}$ at each time
$t=66:T$ where $\sigma_{t}^{2}=(\sigma_{t,1}^{2},\cdots,\sigma_{t,K}^{2})$.
Thus, the total number of generated particles in the Gibbs sampler is
$3(T-65)N$. In the SMC method, at each time $t$, we need only $M$ particles to
obtain the posterior at time $t$. Then, we estimate the computational time in
MCMC by the ratio of the numbers of particles; the MCMC method takes
$3(T-65)N/M$ times longer than the SMC method.

To see this estimate in the real application, consider the actual raw
computational times at $T=248$ with particle size $N=M=10000$, which is also
reported in Figure 1. The SMC method takes $0.53$ seconds to complete, while
the MCMC method takes $365.5$ seconds, or $\log(363.5)\approx 5.89$ in the
scale of log-seconds. The computational time to implement the MCMC method is
estimated as $\log((248-65)\times 3)+\log(0.53)\approx 5.67$, which is close
to the actual computational time. This estimate could be useful to infer about
the computational time for the MCMC method before its implementation and to
decide the number of iterations. However, note that both methods can be speed-
up by utilizing parallel computation.

###  S2.4 Thresholds of ESS for MCMC intervention

The intervention by the MCMC method occurs when $\mathrm{ESS}_{t}<C$. In the
main text, we set $C=500$. Here, for different values of threshold $C$, we
show the number of interventions, raw computational time and LPDRs in Table
S1. As easily expected, the higher the threshold $C$ is, the more
interventions are made, and the more computational time is needed. It is
interesting to see that both $C=100$ and $C=250$ have two interventions, but
$C=100$ takes longer to compute the posteriors. This is because the
interventions occur at later time points (large $t$) when using $C=100$,
meaning that the MCMC method is applied to the longer time series and needs to
sample many latent variables. The patterns of LPDRs are difficult to interpret
as the function of $C$, for the timing of interventions could abruptly change
as $C$ increases/decreases. Hence, choosing an optimal $C$ is difficult and
requires manual, trial-and-error learning.

$C=$ |  | the number of
---
interventions
| computational
---
time (sec)
LPDR $(t=200)$ | LPDR $(t=248)$
$100$ | 2 | 831.53 | -1.045 | -2.367
$250$ | 2 | 535.21 | -0.622 | 1.176
$500$ | 5 | 1266.79 | 0.875 | 1.279
$1000$ | 9 | 2456.86 | 0.190 | -0.717
$2000$ | 25 | 6316.75 | -0.019 | 0.042
Table S1: Effective sample sizes (ESSs), raw computational time (in seconds)
and log predictive density ratios (LPDRs) against the MCMC method $(N=70000)$
at two points for different values of intervention threshold $C$. The raw
computational time with $C=100$ is longer than that with $C=250$ because, with
$C=100$, the intervention occurs at a later time, which makes computational
more costly.

##  Appendix S3 LDF for BPS ; Real Data analysis setting

###  S3.1 Discount factors and LDPLs

In the real data analysis, we use $J=35$ different values for the DBPS
discount factors, $(\beta,\delta)$. Those 35 values are listed as set $S$
below:

$\begin{split}S=&\left\\{(r,r)|r=0.99-i\cdot 0.01\ \text{for}\
i=1:20\right\\}\\\ &\cup\left\\{(r,r-0.02)|r=0.99-i\cdot 0.02\ \text{for}\
i=0:4\right\\}\\\ &\cup\left\\{(r,r-0.04)|r=0.99-i\cdot 0.02\ \text{for}\
i=0:4\right\\}\\\ &\cup\left\\{(r,r-0.06)|r=0.99-i\cdot 0.02\ \text{for}\
i=0:4\right\\}.\end{split}$

For each $(\beta_{j},\delta_{j})\in S$, we can compute the synthesized
predictive density $p(y_{s}|y_{1:s-1},\mathcal{H}_{1:s},\beta,\delta)$ for
$s=66:t$ and the LDPL at time $t$ by

$\text{LDPL}_{j,t}(\gamma)=\sum_{s=66}^{t}\gamma^{t-s}\log
p(y_{s}|\mathcal{H}_{1:s},y_{1:(s-1)},\beta_{j},\delta_{j}),$

where $\gamma$ is the LDF discount factor and set as $\gamma=0.98$ in the main
text. We combine the 35 synthesized predictions based on their LDPLs and
obtain our prediction for $y_{t+1}$.

For each $\text{LDF}_{\ast_{1},\ast_{2}}$, we compute the LDPLs in the period
of $s=66:t$. For the first layer discount factor, $\gamma_{1}$, we consider
the following 15 values:

$\gamma_{1}\in\\{0.01,0.3,0.5,0.6,0.7,0.75,0.8,0.85,0.9,0.92,0.95,0.97,0.98,0.99,1\\}.$

The second layer discount factor is set as $\gamma_{2}=0.98$.

###  S3.2 Results about other DBPSs in the LDF

The predictive distributions of $\text{LDF}_{B,a}$ and $\text{LDF}_{B,s}$ are
shown in Figure S5. Both DBPSs provide similar predictive distributions. A
slight difference can be seen, for example, in the narrower predictive
intervals of the $\text{LDF}_{B,a}$ before 2020. This difference can be
explained by the use of the argmax function; the predictive distribution of
$\text{LDF}_{B,a}$ uses a single discount factor.

Figure S5:  Predictive medians and 90% predictive intervals of
$\text{LDF}_{B,a}$ and $\text{LDF}_{B,s}$ with $\gamma=0.98$ laid over the
data in the period $t=117:248$ (top) and $t=117:236$ (bottom). After the
sudden structural change at $t=242$, the $\text{LDF}_{B,a}$ has wider
predictive distributions due to the low discount factor it uses.
$\text{LDF}_{B,s}$ has slightly wider predictive distributions before that
change.

We consider three values for the LDF discount factors:
$\gamma\in\\{0.95,0.98,1\\}$. Figure S6 shows the LPDRs of $\text{LDF}_{B,a}$
and $\text{LDF}_{B,s}$ with three discount factors. If $\gamma$ is high, then
the prediction is less adaptive to the sudden change after 2020. If $\gamma$
is low, then the predictions become myopic and less performative before 2020.
We conclude that higher values of $\gamma$, but not unity, are preferred.

Figure S6: Log predictive density rations (LPDRs) of $\text{LDF}_{B,a}$ and
$\text{LDF}_{B,s}$ with LDF discount factor $\gamma=0.95,0.98,1$ again the
benchmark DBPS with the fixed discount factors.

Figure S7 shows the values of BPS discount factors $\beta$ and $\delta$ that
have the highest LDPLs in the $\text{LDF}_{B,a}$. In response to the sudden
burst, both discount factors to adjust bias in predictive location and
uncertainty. Furthermore, in cases of $\gamma=0.98$, the value of the discount
factor increases again after the burst, reducing the predictive uncertainty.

Figure S7: The BPS discount factors, $\beta$ (solid) and $\delta$ (dashed),
selected in the $\text{LDF}_{B,a}$ by having the highest LDPLs with the LDF
discount factor $\gamma=1$ (top), $\gamma=0.98$ (middle) and $\gamma=0.95$
(bottom).

###  S3.3 Results about other loss discounting approaches

Figure S8 shows the LPDRs computed in the LDF for four combinations of the
softmax and argmax weights. The $\text{LDF}_{s,a}$ used in the main text
performs satisfactorily both before and after 2020. All four models are not as
competitive as the benchmark DBPS.

Figure S8: LPDRs of the four LDFs with $\gamma_{2}=0.98$.

The first layer discount factor, $\gamma_{1}$, that has the highest LDPL is
shown in Figure S9. The sudden decrease in the discount factor can also be
seen in the LDF. In some cases, however, this decrease occurs earlier than
2020, responding to fluctuations in inflation rates in 2013-2016.

Figure S9:  The first layer discount factors, $\gamma_{1}$, selected in the
$\text{LDF}_{a,a}$ (solid) and $\text{LDF}_{s,a}$ (dashed) by having the
highest LDPLs with the LDF discount factor $\gamma_{2}=1$ (top),
$\gamma_{2}=0.98$ (middle) and $\gamma_{2}=0.95$ (bottom).