Context-tree weighting for real-valued time series:
Bayesian inference with hierarchical mixture models

Consider a real-valued time series $(x_{1},x_{2},\ldots)$. The first key element of the present development is the use of an observable state for each $x_{n}$, based on discretised versions of some of the samples $(\ldots,x_{n-2},x_{n-1})$ preceding it. We refer to the string consisting of these discretised previous samples as the discrete context; it plays the role of a discrete-valued feature vector that can be used to identify useful nonlinear structure in the data. These contexts are extracted via simple quantisers $Q:\mathbb{R}\to A:=\{0,1,\dots,m-1\}$, of the form,

where in this section the thresholds $\{c_{1},\ldots,c_{m-1}\}$ and the resulting quantiser $Q$ are considered fixed. A systematic way to infer the thresholds from data is described in Section III-B.

This general framework can be used in conjunction with an arbitrary way of extracting discrete features, based on an arbitrary mapping to a discrete alphabet, not necessarily of the form in (1). However, the quantisation needs to be meaningful in order to lead to useful results. Quantisers as in (1) offer a generally reasonable choice although, depending on the application at hand, there are other useful approaches, e.g.,  quantising the percentage differences between samples.

Given a quantiser $Q$ as in (1), a maximum context length $D\geq 0$, and a proper $m$-ary context tree $T$, the context (or ‘state’) of each sample $x_{n}$ is obtained as follows. Let $t=(Q(x_{n-1}),\ldots,Q(x_{n-D}))$ be the discretised string of length $D$ preceding $x_{n}$; the context $s$ of $x_{n}$ is the unique leaf of $T$ that is a suffix of $t$. For example, for the context tree of Figure 1, if $Q(x_{n-1})=0$ and $Q(x_{n-2})=1$ then $s=01$, whereas if $Q(x_{n-1})=Q(x_{n-2})=1$ then $s=1$.

The leaves of the tree define the set of discrete states in the hierarchical model. So, for the example BCT-SSM of Figure 1, the set of states is $\mathcal{S}=\{1,01,00\}$. Equivalently, this process can be viewed as defining a partition of ${\mathbb{R}}^{2}$ into three regions indexed by the contexts ${\cal S}$ in $T.$

To complete the specification of the BCT-SSM, a different time series model $\mathcal{M}_{s}$ is associated with each leaf $s$ of the context tree $T$, giving a different conditional density for $x_{n}$. At time $n$, given the context $s$ determined by the past $D$ samples $(x_{n-1},\ldots,x_{n-D})$, the distribution of $x_{n}$ is determined by the model $\mathcal{M}_{s}$ assigned to $s$. Parametric models with parameters $\theta_{s}$ at each leaf $s$ are considered. Altogether, the BCT-SSM consists of an $m$-ary quantiser $Q$, an $m$-ary tree $T$ that defines the set of discrete states, and a collection of parameter vectors $\theta_{s}$ for the parametric models at its leaves.

Identifying $T$ with the collection of its leaves ${\cal S}$, and writing $x_{i}^{j}$ for the segment $(x_{i},x_{i+i},\ldots,x_{j})$, the likelihood is,

where $B_{s}$ is the set of indices $i\in\{1,2,\ldots,n\}$ such that the context of $x_{i}$ is $s$, and $\theta=\{\theta_{s}\;;\;s\in T\}$.

For the top level we consider collections of states represented by context trees $T$ in the class $\mathcal{T}(D)$ of all proper $m$-ary trees with depth no greater than $D$, where $T$ is proper if any node in $T$ that is not a leaf has exactly $m$ children.

Prior structure. For the trees $T\in{\mathcal{T}}(D)$ with maximum depth $D\geq 0$ at the top level of the hierarchical model, we use the Bayesian Context Trees (BCT) prior of [17],

where $\beta\in(0,1)$ is a hyperparameter, $\alpha=(1-\beta)^{1/(m-1)}$, $|T|$ is the number of leaves of $T$, and $L_{D}(T)$ is the number of leaves of $T$ at depth $D$. This prior penalises larger trees by an exponential amount, a desirable property that discourages overfitting. The default value $\beta=1-2^{-m+1}$ [17] is used. Given a tree model $T\in\mathcal{T}(D)$, an independent prior is placed on each $\theta_{s}$, so that $\pi(\theta|T)=\prod_{s\in T}\pi(\theta_{s})$.

Typically, the main obstacle in performing Bayesian inference with a time series $x$ is the computation of the normalising constant (or evidence) $p(x)$ of the posterior distriburion,

The power of the proposed Bayesian structure is that, although $\mathcal{T}(D)$ is enormously rich, consisting of doubly-exponentially many models in $D$, it is actually possible to perform exact Bayesian inference efficiently. To that end, we introduce the Continuous Context Tree Weighting (CCTW) algorithm, and the Continuous Bayesian Context Tree (CBCT) algorithm, generalising the corresponding algorithms for discrete time series in [17]. It is shown that CCTW computes the normalising constant $p(x)$ exactly (Theorem 1), and CBCT identifies the MAP tree model (Theorem 2). The main difference from the discrete case in [17], both in the algorithmic descriptions and in the proofs of the theorems (given in Appendix A), is that a new generalised form of estimated probabilities is introduced and used in place of their discrete versions. For a time series $x=x_{-D+1}^{n}$, these are,

Let $x=x_{-D+1}^{n}$ be a time series, and let $y_{i}=Q(x_{i})$ denote the corresponding quantised samples.

CCTW: Continuous context-tree weighting algorithm

1. 1.

   Build the tree $T_{\text{MAX}}$, whose leaves are all the discrete contexts $y_{i-D}^{i-1},\ i=1,2,\ldots,n$. Compute $P_{e}(s,x)$ as given in (4) for each node $s$ of $T_{\text{MAX}}$.

2. 2.

   Starting at the leaves and proceeding recursively towards the root compute:

   $$P_{w,s}\!=\!\left\{\begin{array}[]{ll}P_{e}(s,x),\;\;\;\mbox{if $s$ is a leaf,}\\ \beta P_{e}(s,x)+(1-\beta)\prod_{j=0}^{m-1}P_{w,sj},\;\;\;\mbox{o/w,}\end{array}\!\!\right.\!\!\vspace*{-0.05 cm}$$

   where $sj$ is the concatenation of context $s$ and symbol $j$.

CBCT: Continuous Bayesian context tree algorithm

1. 1.

   Build the tree $T_{\text{MAX}}$ and compute $P_{e}(s,x)$ for each node $s$ of $T_{\text{MAX}}$, as in CCTW.

2. 2.

   Starting at the leaves and proceeding recursively towards the root compute:

   $$P_{m,s}\!=\!\left\{\begin{array}[]{ll}P_{e}(s,x),\;\;\;\mbox{if $s$ is a leaf at depth $D$,}\\ \beta,\;\;\;\mbox{if $s$ is a leaf at depth $<D$,}\\ \max\big{\{}\beta P_{e}(s,x),(1-\beta)\prod_{j=0}^{m-1}P_{m,sj}\big{\}},\ \mbox{o/w.}\end{array}\right.$$

3. 3.

   Starting at the root and proceeding recursively with its descendants, for each node: If the maximum is achieved by the first term, prune all its descendants from $T_{\text{MAX}}$.

Even in cases where the integrals in (4) are not tractable, the fact that they are in the form of standard marginal likelihoods makes it possible to compute them approximately using standard methods, e.g., [8, 9, 46]. The above algorithms can then be used with these approximations as a way of performing approximate inference for the BCT-SSM. However, this is not investigated further in this work. Instead, the general principle is illustrated via an interesting example where the estimated probabilities can be computed explicitly and the resulting mixture model is a flexible nonlinear model of practical interest. This is described in the next section, where AR models ${\cal M}_{s}$ are associated to each leaf $s$. We refer to the resulting model as the Bayesian context tree autoregressive (BCT-AR) model, which is just a particular instance of the general BCT-SSM.

For each symbol $x_{i}$ in a time series $x_{1}^{n}$, exactly $D+1$ nodes of $T_{\text{MAX}}$ need to be updated, corresponding to its contexts of length $0,1,\ldots,D$. For each one of these nodes, only the quantities $\{|B_{s}|,s_{1},\mathbf{s}_{2},S_{3}\}$ need to be updated, which can be done efficiently by just adding an extra term to each sum. Using these and Lemma 1, the estimated probabilities $P_{e}(s,x)$ can be computed for all nodes of $T_{\text{MAX}}$.

Hence, the complexity of both algorithms as a function of $n$ and $D$ is only $\mathcal{O}(nD)$: linear in the length of the time series and the maximum depth. Therefore, the present methods are computationally very efficient and scale well with large numbers of observations. [Taking into account $m$ and $p$ as well, it is easy to see that the complexity is $\mathcal{O}\left(nD(p^{3}+m)\right)$.]

The above discussion also shows that, importantly, all algorithms can be updated sequentially. When observing a new sample $x_{n+1}$, only $D+1$ nodes need to be updated, which requires $\mathcal{O}(D)$ operations, i.e., $\mathcal{O}(1)$ as a function of $n$. In particular, this implies that sequential prediction can be performed very efficiently. Empirical running times in the forecasting experiments show that the present methods are much more efficient than essentially all the alternatives examined. In fact, the difference is quite large, especially when comparing with state-of-the-art ML models that require heavy training and do not allow for efficient sequential updates, usually giving empirical running times that are larger by several orders of magnitude; see also [24] for a relevant review comparing the computational requirements of ML versus statistical techniques.

Finally, a principled Bayesian approach is introduced for the selection of the quantiser thresholds $\{c_{i}\}$ of (1) and the AR order $p$. Viewed as extra parameters on an additional layer above everything else, uniform priors are placed on $\{c_{i}\}$ and $p$, and Bayesian model selection [22] is performed to obtain their MAP values: The resulting posterior distribution $p(\{c_{i}\},p|x)$ is proportional to the evidence $p(x|\{c_{i}\},p)$, which can be computed exactly using the CCTW algorithm (Theorem 1). So, in order to select appropriate values, suitable ranges of possible $\{c_{i}\}$ and $p$ are specified, and the values with the higher evidence are selected. For the AR order we take $1\leq p\leq p_{\text{max}}$ for an appropriate $p_{\text{max}}$ ($p_{\text{max}}=5$ in our experiments), and for the $\{c_{i}\}$ we perform a grid search in a reasonable range (e.g., between the 10th and 90th percentiles of the data).

First an experiment on simulated data is performed, illustrating that the present methods are consistent and effective with data generated by BCT-SSM models. The context tree used is the tree of Figure 1, the quantiser threshold is $c=0$, and the AR order $p=2$. The posterior distribution over trees, $\pi(T|x)$, is examined. On a time series consisting of only $n=100$ observations, the MAP tree identified by the CBCT algorithm is the empty tree corresponding to a single AR model, with posterior probability 99.9%. This means that the data do not provide sufficient evidence to support a more complex state space partition. With $n=300$ observations, the MAP tree is now the true underlying model, with posterior probability 57%. With $n=500$ observations, the posterior of the true model is 99.9%. So, the posterior indeed concentrates on the true model, indicating that the BCT-AR inferential framework can be very effective even with limited training data.

Forecasting. The performance of the BCT-AR model is evaluated on out-of-sample 1-step ahead forecasts, and compared with state-of-the-art approaches in four simulated and three real datasets. The first simulated dataset (sim_1) is generated by the BCT-AR model used above, and the second (sim_2) by a BCT-AR model with a ternary tree of depth 2. The third and fourth ones (sim_3 and sim_4) are generated from SETAR models of orders $p=1$ and $p=5$, respectively. In each case, the training set consists of the first 50% of the observations; also, all models are updated at every time-step in the test set. For BCT-AR, the MAP tree with its MAP parameters is used at every time-step, which can be updated efficiently (Section III-A). From Table I, it is observed that the BCT-AR model outperforms the alternatives, and achieves the lowest mean-squared error (MSE) even on the two datasets generated from SETAR models. As discussed in Section III-A, the BCT-AR model also outperforms the alternatives in terms of empirical running times.

An important application of SETAR models is in modelling the US unemployment rate [13, 26, 42]. As described in [42], the unemployment rate moves countercyclically with business cycles, and rises quickly but decays slowly, indicating nonlinear behaviour. Here, the quarterly US unemployment rate in the time period from 1948 to 2019 is examined (dataset unemp, 288 observations). Following [26], the difference series ${\Delta x}_{n}=x_{n}-x_{n-1}$ is considered, and a constant term is included in the AR model. For the quantiser alphabet size, $m=2$ is a natural choice, as will become apparent below. The threshold selected using the procedure of Section III-B is $c=0.15$, and the MAP tree is the tree of Figure 1, with depth $d=2$, leaves $\{1,01,00\}$, and posterior probability 91.5%. The fitted BCT-AR model with its MAP parameters is,

with $e_{n}\sim\mathcal{N}(0,1)$, and corresponding regions $s=1$ if ${\Delta x}_{n-1}>0.15$, $s=01$ if ${\Delta x}_{n-1}\leq 0.15,{\Delta x}_{n-2}>0.15$, and $s=00$ if ${\Delta x}_{n-1}\leq 0.15,\ {\Delta x}_{n-2}\leq 0.15$.

Interpretation. The MAP BCT-AR model finds significant structure in the data, providing a very natural interpretation. It identifies 3 meaningful states: First, jumps in the unemployment rate higher than 0.15 signify economic contractions (context 1). If there is no jump at the most recent time-point, the model looks further back to determine the state. Context 00 signifies a stable economy, as there are no jumps in the unemployment rate for two consecutive quarters. Finally, context 01 identifies an intermediate state: “stabilising just after a contraction”. An important feature identified by the BCT-AR model is that the volatility is different in each case: Higher in contractions ($\sigma=0.42$), smaller in stable economy regions ($\sigma=0.20$), and in-between for context 01 ($\sigma=0.32$).

Forecasting. In addition to its appealing interpretation, the BCT-AR model outperforms all benchmarks in forecasting, giving a 6% lower MSE than the second-best method (Table I).

Another important example of a nonlinear time series is the US Gross National Product (GNP) [33, 13]. The quarterly US GNP in the time period from 1947 to 2019 is examined (291 observations, dataset gnp). Following [33], the difference in the logarithm of the series, $y_{n}=\log x_{n}-\log x_{n-1}$, is considered. As above, $m=2$ is a natural choice for the quantiser size, helping to differentiate economic expansions from contractions – which govern the underlying dynamics. The MAP BCT-AR tree is shown in Figure 2: It has depth $d=3$, its leaves are $\{0,10,110,111\}$ and its posterior probability is 42.6%.

Interpretation. Compared with the previous example, here the MAP BCT-AR model finds even richer structure in the data and identifies four meaningful states. First, as before, there is a single state corresponding to an economic contraction (which now corresponds to $s=0$ instead of $s=1$, as the GNP obviously increases in expansions and decreases in contractions). And again, the model does not look further back whenever a contraction is detected. Here, the model also shows that the effect of a contraction is still present even after three quarters ($s=110$), and that the exact ‘distance’ from a contraction is also important, with the dynamics changing depending on how much time has elapsed. Finally, the state $s=111$ corresponds to a flourishing, expanding economy, without a contraction in the recent past. An important feature captured by the BCT-AR model is again that the volatility is different in each case. More specifically, it is found that the volatility strictly decreases with the distance from the last contraction, starting with the maximum $\sigma=1.23$ for $s=0$ and decreasing to $\sigma=0.75$ for $s=111$.

Forecasting. The BCT-AR model outperforms all benchmarks in forecasting, giving a 12% lower MSE than the second-best method, as presented in Table I.

Finally, the daily IBM common stock closing price from May 17, 1961 to November 2, 1962 is examined (dataset ibm, 369 observations), taken from [7]. This is a well-studied dataset, with [7] fitting an ARIMA model, [40] fitting a SETAR model, and [45] fitting a MAR model to the data. Following previous approaches, the first-difference series, ${\Delta x}_{n}=x_{n}-x_{n-1}$, is considered. The value $m=3$ is chosen for the alphabet size of the quantiser, with contexts $\{0,1,2\}$ naturally corresponding to the states {down, steady, up} for the stock price. Using the procedure of Section III-B to select the thresholds, the resulting quantiser regions are: $s=0$ if ${\Delta x}_{n-1}<-7$, $s=2$ if ${\Delta x}_{n-1}>7$, and $s=1$ otherwise. The MAP tree is shown in Figure 3: It has depth $d=2$, and its leaves are $\{0,2,10,11,12\}$, hence identifying five states. Its posterior probability is 99.3%, suggesting that there is very strong evidence in the data supporting this exact structure, even with only 369 observations.

Interpretation. The BCT-AR model reveals important information about apparent structure in the data. Firstly, it admits a very simple and natural interpretation: In order to determine the AR model generating the next value, one needs to look back until there is a significant enough price change (corresponding to contexts 0, 2, 10, 12), or until the maximum depth of 2 (context 11) is reached. Another important feature captured by this model is the commonly observed asymmetric response in volatility due to positive and negative shocks, sometimes called the leverage effect [42, 7]. Even though there is no suggestion of that in the prior, the MAP parameter estimates show that negative shocks increase the volatility much more: Context 0 has the highest volatility ($\sigma=12.3$), with 10 being a close second ($\sigma=10.8$), showing that the effect of a past shock is still present. In all other cases the volatility is much smaller ($5.17\leq\sigma\leq 6.86$).

Forecasting. As shown Table I, in this example deepAR performs marginally better in 1-step ahead forecasts, with BCT-AR, ETS and MAR also having comparable performance.