Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
119,117 characters · 27 sections · 99 citation commands
-0.2in The Bayesian Context Trees State Space Model for time series modelling and forecasting
A hierarchical Bayesian framework is introduced for developing tree-based mixture models for time series, partly motivated by applications in finance and forecasting. At the top level, meaningful discrete states are identified as appropriately quantised values of some of the most recent samples. At the bottom level, a different, arbitrary `base’ model is associated with each state. This defines a very general framework that can be used in conjunction with any existing model class to build flexible and interpretable mixture models. We call this the Bayesian Context Trees State Space Model, or the BCT-X framework. \textcolor{black}{Appropriate algorithmic tools are described, which allow for effective and efficient Bayesian inference and learning; these algorithms can be updated sequentially, facilitating online forecasting}. The utility of the general framework is illustrated in the particular instances when AR or ARCH models are used as base models. The latter results in a mixture model that offers a powerful way of modelling the well-known volatility asymmetries in financial data, revealing a novel, important feature of stock market index data, in the form of an enhanced leverage effect. In forecasting, the BCT-X methods are found to outperform several state-of-the-art techniques, both in terms of accuracy and computational requirements.
Keywords. Time series, Interpretable mixture models, Bayesian inference, Context-tree models, Forecasting, AR models, Conditional heteroscedastic models, Volatility asymmetries.
Time series modelling, inference and forecasting are well-studied tasks in statistics and machine learning (ML), with critical applications throughout finance, the sciences, and engineering. A wide range of approaches exist, from classical statistical methods box:book,durbin:book2, to modern ML techniques: notably neural networks benidis:20,alexandrov:20,zhang:98, Gaussian processes (GPs) rasmussen:book,roberts:13,frigola:15, and matrix factorisations yu:16,faloutsos:18. { Recent reviews that provide a thorough comparison between ML and classical methods specifically in the time series setting can be found in makridakis:18b,makridakis:18,ahmed:10. Motivated in part by important applications in financial time series and the inherent limitations of neural network models -- in that the they lack interpretability and have large training-data requirements makridakis:18b,benidis:20 -- in this work we propose a new, general class of flexible hierarchical Bayesian models, which are both naturally {\em interpretable} and suitable for applications with limited training data. For these models, we provide computationally efficient -- linear complexity -- algorithms for \textcolor{black}{effective} inference, forecasting, and posterior exploration. }
Roughly speaking, time series models can be broadly categorised in two classes, depending on the absence or presence of an underlying hidden state process. The first class includes the family of linear autoregressive (AR) and autoregressive integrated moving average (ARIMA) models box:book,tsay:05,tong:book, which directly model the stochastic mapping from previous to current observations, as well as numerous AR extensions and generalisations, including nonlinear AR models that employ GPs roberts:13,girard:02,turner:12,murray:01 and neural networks to model nonlinear dynamics. Feed-forward neural networks, including the Neural Network AR (NNAR) model, have been used since the 1990s in this setting zhang:98, and recurrent neural networks (RNN) have also been employed graves:13,salinas:20,oreshkin:19, with the DeepAR model of salinas:20 that uses Long Short-Term Memory (LSTM) cells hochreiter:97 being one of the most successful approaches.
The second class consists of State Space Models (SSM) durbin:book2,hyndman:book, which can describe more complex dynamics by combining an underlying hidden state process (`state-transition') with an observation model (`emission'). Classical approaches here include hidden Markov models (HMMs) cappe:book06, the linear Gaussian SSM kalman:60, exponential smoothing methods hyndman:book, \textcolor{black}{and dynamic factor models stock2011dynamic,forni2000generalized,stock2002macroeconomic}. Again, various extensions have been proposed by making the transition or emission equations nonlinear, using GPs frigola:13,frigola:14,eleftheriadis:17,turner:12 or neural networks krishnan:15,krishnan:17,zheng:17,karl:17, rangapuram:18. Inference in these settings is often a challenging task, requiring sophisticated and computationally intense approximate techniques, including Particle MCMC doucet:book,andrieu:10 and variational methods like stochastic gradient variational Bayes (SGVB) kingma:13,rezende:14.
{\bf The BCT-X framework.} In this work, we introduce a novel Bayesian modelling approach that combines important features of both the above classes. First, meaningful discrete states are identified, but these are {\em observable} rather than hidden: Each {\em state} is a short discrete sequence given by the appropriately quantised values of some of the most recent samples. This collection of states is described by a discrete {\em context-tree model}. Then a different time series model -- a {\em base model} -- is associated with each possible current state, generating the next sample.
In technical terms, we define a Bayesian hierarchical model. The top level consists of a set of relevant states, naturally represented as a discrete context-tree model, which admits a natural interpretation and captures important aspects of the structure present in the data. And the bottom level associates an arbitrary time series model to each state. \textcolor{black}{Equivalently, at the top level we can think of the context-tree model as defining a partition of the state space: Every discrete state actually corresponds to a different region of this partition, for which at the bottom level we associate a different base model (that is to be used in that region). } We call the resulting model class the {\em Bayesian Context Trees State Space Model}, or BCT-X. The `BCT' part refers to the discrete context-tree models, and the `X' indicates that any existing time series family of models can be employed as base models for the different states.
\textcolor{black}{Bayesian Context Trees (BCT) were originally introduced BCT-JRSSB:22 as a Bayesian framework for modelling variable-memory Markov chains. So far, the BCT framework has only been used in the restricted setting of {\em{discrete-valued}} time series -- which of course considerably limits its practical applicability. The main conceptual novelty of this work is that, although context trees are naturally suited for discrete data (and have only been used in this setting), we show that they can also be utilised in a very effective way for {\em{real-valued}} time series. Specifically, by introducing an appropriate quantiser, context-tree models can be used to define discrete states (or equivalently, partitions of the state space) as explained above, and hence provide a general way of building flexible and interpretable mixture models, by associating a different, arbitrary base-model to each state/state-space region. }
Although at times we refer to BCT-X as a `model', it is in fact a general framework for building Bayesian mixture models for time series, which can be used in conjunction with any existing model class. The resulting model family is rich, flexible, and much more general than the class one starts with. For example, using any of the standard linear families (like AR or ARIMA) leads to much more general mixture models that can capture highly nonlinear trends and are also easily interpretable. Moreover, this type of observable state process facilitates effective Bayesian inference. This is achieved by exploiting the structure of context-tree models, in a similar spirit to the methods developed for discrete-valued time series in BCT-JRSSB:22,ctw-isit:21.
{\bf Algorithmic tools for inference.} \textcolor{black}{ A family of important algorithms are introduced in conjunction with the BCT-X framework, by adapting the algorithms described in BCT-JRSSB:22 for the discrete-valued setting to the much more general setting considered here. First, the `Generalised Context-Tree Weighting' algorithm (GCTW) computes the evidence mackay:92 exactly, with all models and parameters integrated out. Second, the `Generalised Bayesian Context Trees' (GBCT) and $k$-GBCT algorithms identify the $k$ {\em a posteriori} most likely (MAP) context-tree models, along with their posterior probabilities. Lastly, a direct Monte Carlo procedure is developed for further exploring the posterior. }
Importantly, using this Bayesian approach, the set of relevant states -- namely, the context tree model -- is identified directly from the data, and does not need to be specified {\em a priori}. This avoids common overfitting problems that lead to lack of interpretability and poor out-of-sample performance. Moreover, the above algorithms have only linear complexity and allow for sequential updates. These properties are ideally suited for online forecasting and provide an important practical advantage compared to standard ML-based time series approaches.
The BCT-X modelling framework along with these algorithmic tools provide a powerful Bayesian framework for modelling and for \textcolor{black}{effective} -- and computationally efficient -- Bayesian inference. \textcolor{black}{The application of the general framework is illustrated by introducing two particular model classes that are found to be useful in practical applications, described next.}
{\bf BCT-AR models.} First, we examine the case where AR models are used as base models for BCT-X, with a different AR model associated to each state. We refer to the resulting model class as the {\em Bayesian Context Trees Autoregressive} (BCT-AR) model. This is shown to be a flexible, nonlinear AR mixture model that generalises popular AR mixtures, including the threshold AR (TAR) models tong:11,tong:80,tong:book and the mixture AR (MAR) models of wong:00. The BCT-AR model is found to outperform state-of-the-art methods in experiments on both simulated data and real-world data from standard applications of nonlinear time series in economics and finance, both in terms of forecasting accuracy and computational requirements.
{\bf BCT-ARCH models.} Second, we employ autoregressive conditional heteroscedastic (ARCH) models as base models. This results in yet another flexible mixture model class, BCT-ARCH, that provides a systematic and powerful way of modelling the well-known asymmetric response in volatility due to positive and negative shocks, which is an important feature of financial time series tsay:05. \textcolor{black}{In fact, modelling volatility asymmetries with the BCT-ARCH model serves as an important motivating example for illustrating the purpose and the practical utility of the general BCT-X methodology. Specifically, in contrast with traditional modelling of these asymmetries, the BCT-ARCH model identifies that it is not only the sign of the most recent value-change that is relevant in the volatility response, but that the exact pattern of recent ‘ups’ and ‘downs’ is also important; this is perfectly captured by a collection of discrete states in the BCT-ARCH model.} As a result, the BCT-ARCH model is found to outperform previous approaches that were developed to describe this effect, and is able to identify this interesting and newly observed structure in the data, in the form of an {\em enhanced leverage effect}.
Finally, we mention that a number of earlier approaches, e.g. alvarez:10,alvisi:07,berndt:94,fu:07,hu:14,liu:11, ouyang:10,hero:20, have employed discrete patterns in the analysis of real-valued time series. These works illustrate the fact that useful and meaningful information can indeed be extracted from discrete contexts. However, in most cases the methods developed have been either application- or task-specific, and often need to resort to {\em ad hoc} considerations for inference. In contrast, in this work discrete states are used in a natural manner, by defining a hierarchical Bayesian modelling structure upon which principled Bayesian inference is performed.
Before outlining the general construction of the BCT-X framework for real-valued time series $(x_n)$, we briefly illustrate its structure through a particular example of the BCT-AR model with context depth $D=2$ and AR order $p$. Given the past values $(\ldots,x_{n-2},x_{n-1})$ and the context-tree model $T$ shown in Figure (ref), the distribution of $x_n$ is determined as follows.
First, the {\em context} of $x_n$ is defined as the binary string given by $t:=({\mathbb I}_{\{x_{n-1}\geq 0\}},{\mathbb I}_{\{x_{n-2}\geq 0\}})\in\{0,1\}^2$ of length $D=2$. Then the corresponding {\em state} at time $n$ is the unique suffix $s$ of $t$ that appears as a leaf of the model $T$. For example, the context $t=(0,0)$ corresponds to the state $s=00$, while the context $t=(1,1)$ corresponds to the state $s=1$. Finally, the distribution of $x_n$ given $(x_{n-p},\ldots,x_{n-1})$ and $s$ is determined by the AR($p$) model with parameters $\theta_s$.
Although at first glance this may seem quite similar to the constructions of the mixture autoregressive (MAR) wong:00 or threshold autoregressive (TAR) tong:80 models, the BCT-X class is more general and critically different. First, the set of relevant states is not defined by the modeller but it is determined by data: A prior distribution is placed on {\em all possible models $T$} of maximum depth $D$, and Bayesian inference is performed by considering the resulting posterior. Second, any model class can be used for the {\em base models} in place of the AR models in this example.
Let $D\geq 0$ be a fixed, maximum context depth, and let $x$ denote the sequence of observations, consisting of the time series $x_1^n=(x_1,x_2,\ldots,x_n)$ together with the initial context $x_{-D+1}^0=(x_{-D+1},\ldots,x_0)$, with all $x_i\in{\mathbb R}$. Throughout, we write $x_i^j=(x_i,x_{i+1},\ldots,x_j)$, for $i<j$.
The first step in the development of the BCT-X framework is the construction of discrete, observable states. To that end, we consider piecewise constant quantisers from $\mathbb{R}$ to a finite alphabet $A = \{0,1,\dots , m-1 \}$, of the form,
where the alphabet size $m$, the thresholds $\{c_i\}$, and the resulting quantiser $Q$ are considered to be fixed for the rest of this section. \textcolor{black}{ A core part of our general methodology is the development of a systematic and computationally efficient way to infer the thresholds $\{c_i\}$ from data, as described in Section (ref). The quantiser function $Q$ plays a key role in the formulation of the BCT-X framework, as its role is crucial in describing the relevant states and state-space partitions.}
Given a maximum context depth $D\geq 0$, we consider the class ${\cal T}(D)$ of all proper $m$-ary trees of depth no more than $D$, where a tree $T$ is {\em proper} if any node in $T$ that is not a leaf has exactly $m$ children. The elements of $T\in{\cal T}(D)$ are the {\em context-tree models}: Given $T\in{\cal T}(D)$, the sequence $t=(Q(x_{i-1}),\ldots,Q(x_{i-D}))$ is the {\em context} at time $i\geq 1$. Since $T$ is proper, for any context $t$ there is a unique suffix $s$ of $t$ that is a leaf of $T$. This $s$ is the {\em state} at time $i$, and it plays the role of a discrete feature vector that can be used to identify useful structure in the data. The leaves of $T$ define the set ${\cal S}$ of discrete {\em states} in our hierarchical model. For example, for the tree of Figure (ref) and a binary quantiser with threshold $c=0$, we get $Q(x)={\mathbb I}_{\{x\geq 0 \}}$ and $\mathcal {S} = \{ 1, 01, 00\}$.
To complete the specification of the BCT-X framework, we associate a different time series model $\mathcal {M}_s$ to each state $s$, i.e., to each leaf $s$ of the context-tree model $T$, giving a different conditional distribution for the next sample. We refer to the class of all such models ${\cal M} _s$ as the {\em base model class}. At time $i$, given the current state $s$ determined by the context $t=(Q(x_{i-1}),\ldots,Q(x_{i-D}))$, the distribution of $x_i$ is given by the model $\mathcal {M}_s$ assigned to $s$. Although general non-parametric models could also be used, for the rest of this paper we consider parametric models with parameters $\theta_s$ for each state $s$. Altogether, a specific instance of the BCT-X framework consists of an $m$-ary quantiser $Q$, a context-tree model $T\in{\cal T}(D)$ that defines the set of discrete states ${\cal S}$, a collection $\theta=\{\theta_s\}$ of parameter vectors $\theta_s$ at the states $s\in{\cal S}$, and a corresponding collection $\{{\cal M}_s\}=\{{\cal M}(\theta_s)\}$ of base models.
Let $x$ denote a time series $x_1^n$ along with its initial context $x_{-D+1}^0$. Identifying $T$ with the collection of its leaves ${\cal S}$, the likelihood induced by the resulting BCT-X model is,
where $B_s$ is the set of indices $i\in\{1,2,\ldots ,n\}$ such that the state at time $i$ is $s$.
For the top level of the BCT-X framework, we consider context-tree models $T$ in the class $\mathcal T (D)$, and employ the Bayesian Context Trees (BCT) prior of BCT-JRSSB:22 on ${\cal T}(D)$,
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 models by an exponential amount to avoid overfitting; see BCT-JRSSB:22 for an extensive discussion of the properties of $\pi(T)$ and guidelines regarding the choice of $\beta$. Given a context-tree model $T \in \mathcal T (D)$, the parameter prior on $\theta=\{\theta_s\}$ is a product of independent priors for each $\theta_s$, $\pi (\theta | T) = \prod _ {s \in T} \pi (\theta_s)$. The exact form of $\pi(\theta_s)$ of course depends on the choice of the base models ${\cal M}_s$.
A typical obstacle to performing Bayesian inference is the difficulty in computing the normalising constant $p(x)$ of the posterior density (sometimes referred to as the {\em prior predictive likelihood} or simply as the {\em evidence}), which in this case is given by,
\textcolor{black}{ The power of the proposed Bayesian structure stems, in part, from the fact that, although $ \mathcal T (D)$ is enormously rich, consisting of doubly-exponentially many models in $D$, it is actually possible to perform \textcolor{black}{effective} Bayesian inference within ${\cal T}(D)$ very efficiently. By appropriately modifying and extending the corresponding algorithms for discrete time series from BCT-JRSSB:22, we introduce the Generalised Context-Tree Weighting (GCTW) algorithm and the Generalised Bayesian Context Tree (GBCT) algorithm, which can be used in the present general setting of real-valued time series. We show that GCTW computes $p(x)$ in ((ref)) (Theorem (ref)), and GBCT identifies the MAP context-tree model (Theorem (ref)). Theorems (ref) and (ref) are proved in Section (ref) of the Supplemental Material, where we also introduce the $k$-GBCT algorithm that obtains the top-$k$ a posteriori most likely context-tree models; the details are omitted here. }
Let $x$ be a time series $x_1^n$ with initial context $x_{-D+1}^0$. For a model $T$ with associated parameters $\theta=\{\theta_s\}$, the generalised estimated probabilities $P_e(s,x)$ defined at every node $s$ of $T$ play a central role in the algorithms' descriptions and in the proofs of the theorems. These are,
Write $q_i=Q(x_i)$ for the quantised samples. Given a maximum context depth $D\geq 0$ and the value of the hyperparameter $\beta$, the GCTW and GBCT algorithms operate as follows.
GCTW: The generalised context-tree weighting algorithm
GBCT: The generalised Bayesian context tree algorithm
{\color{black} By introducing an appropriate quantiser $Q$ and using the discretised versions $q_i = Q(x_i)$ of the samples $x_i$, those steps of the generalised versions of the algorithms that are associated with building the tree $T_{\rm MAX}$ and performing recursive operations on it remain the same as before. The important difference in the new algorithms comes from the generalised estimated probabilities $P_e(s,x)$ in ((ref)), computed in Step 2 for every node of $T_{\rm MAX}$. The necessary values that need to be stored at every node $s$ in order to compute $P_e(s,x)$ heavily depend on the choice of the base models; see for example Section (ref) and Section (ref) for the specific quantities that we need to store for the BCT-AR and BCT-ARCH models, respectively.
In fact, the discrete BCT framework of BCT-JRSSB:22 can be viewd as a simple special case of the general BCT-X framework, when simple categorical distributions are chosen as the base models associated to the leaves of the trees; other specific examples are the BCT-AR (Section (ref)) and BCT-ARCH (Section (ref)) models examined in this paper. What remains true in all cases, and is perhaps somewhat remarkable, is that in the general BCT-X setting, the new versions of the algorithms can be used to perform Bayesian inference in a dynamic programming fashion, for arbitrary base models associated at the leaves of the trees.}
{\bf Sampling from the posterior.} The GCTW, GBCT and $k$-GBCT algorithms compute the evidence $p(x)$ and identify the top-$k$ MAP context-tree models. Furthermore, it is in fact possible to obtain independent samples from the model posterior $\pi(T|x)$ as follows. For every node $s$, let $P_{b,s} = \beta P_e(s,x,) / P_{w,s}$. Start with the tree consisting of only the root node $\lambda$, and with probability $P_{b,\lambda}$ mark it as a leaf and stop, or, with probability $(1-P_{b,\lambda})$ add all $m$ of its children to $T$. Then, examine every new node $s$ and either mark it as a leaf with probability $P_{b,s}$, or add all $m$ of its children to $T$ with probability $(1-P_{b,s})$. Examining all non-leaf nodes at depths $<D$ recursively until no more eligible nodes remain, produces a random tree $T \in \mathcal {T} (D)$. Theorem (ref), proved in the Supplemental Material, states that $T$ is a sample from the posterior.
{\color{black}
The methodological tools developed in this section are based on the generalised form of the estimated probabilities $P_e(s,x)$ of ((ref)), whose exact form depends on the particular choice of base models. In the next section, the general principle is illustrated in an interesting case where $P_e(s,x)$ can be computed explicitly, and the resulting mixture model is a flexible nonlinear model of practical interest. Specifically, an AR model ${\cal M}_s$ is associated to each leaf $s$, and we hence refer to the resulting model class as the BCT-AR model; the BCT-AR model is only one particular instance of the general BCT-X framework.
However, even when the the integrals in ((ref)) are not available in closed form, the fact that the estimated probabilities are in the form of standard marginal likelihoods makes it possible to compute them approximately by using standard methods; see, e.g., kass:95,chib:95,chib:01,friel:08,wood:11. Then, the GCTW and GBCT algorithms can be used in exactly the same manner as before, with these approximations $\widehat {P_e} (s,x)$ in place of their exact values, and therefore still facilitate effective Bayesian inference. This is illustrated in detail via the BCT-ARCH version of BCT-X, where the estimated probabilities are not available in closed form; see Section (ref). }
Context-tree models -- introduced as “tree sources” in the information-theoretic literature by Rissanen rissanen:83b,rissanen:83,rissanen:86 in the 1980s -- have been employed widely in connection with various problems on {\em discrete} data, especially for data compression weinbergeretal:94,willems-shtarkov-tjalkens:95, willems:98,matsushima:94,matsushima:09. In the statistics literature they were first popularised as “variable-length Markov chains” buhlmann:00,buhlmann:04, and more recently in connection with the Bayesian Context Trees (BCT) framework BCT-JRSSB:22,ctw-isit:21, which has been found to be very effective in a range of statistical tasks, including model selection, prediction, change-point detection and entropy estimation papag-K-pre:23,branch-isit:22,lungu-arxiv:22, lungu-pap-K:22,papag-K-ITW:23. A central conceptual novelty of the present work is in showing that {\em discrete} context-tree models can in fact also be effectively utilised for modelling real-valued time series, by representing meaningful context-based discrete states that are used to build flexible and interpretable mixture models of practical interest.
{\color{black} In more technical terms, apart from this crucial observation, the major contributions of this work compared to existing BCT papers are summarised below:
}
In a different but related direction, starting with the algorithm of breiman:book2, the classification and regression trees (CART) procedures have been widely used for regression and classification tasks. Typically, these methods either rely on greedy algorithms to a grow a tree with some stopping criteria -- sometimes with additional penalties for pruning -- or they adopt a Bayesian CART approach and place a prior on trees chipman:98,bart; see loh:14 and linero:17 for recent reviews. The main difficulties in the use of these procedures are that greedily-constructed trees tend to overfit the data, while the corresponding Bayesian approaches require the use of MCMC sampling, which is both computationally expensive and not guaranteed to be effective for inference. These reasons perhaps partly explain why, in applications, approaches exploiting CARTs for time series data audrino:01,meek:02,dellaportas:07,taddy2011dynamic,deshpande2020vcbart have not been used as widely as much simpler and more restricted model classes like, e.g., threshold models tong:11,zakoian:94.
In contrast, in this work we take advantage of the special sequential nature of time-series data which, together with the desirable properties of the BCT-X framework, allow for exact Bayesian inference in selecting the tree model -- in a very efficient manner. The BCT-X family of algorithms can identify not only the MAP tree model (GBCT algorithm), but also the top-$k$ a posteriori most likely trees, for moderate $k$ ($k$-GBCT algorithm). Moreover, it is possible to further explore the posterior on model space by easily obtaining i.i.d.\ samples as described in Section (ref), without the need to resort to MCMC; which of course offers another important practical advantage compared to Bayesian CART methods.
Given a quantiser $Q:{\mathbb R}\to A=\{0,1,\ldots,m-1\}$, a fixed AR order $p$, and a maximum context depth $D\geq 0$, the BCT-AR model is the version of BCT-X with AR$(p)$ models as the {\em base model} class. { \color{black} With each possible state, that is, with each potential leaf $s$ in a context-tree model $T\in{\cal T}(D)$, we associate an AR$(p)$ model,
where $\boldsymbol \phi _ s = ( \phi _ {s,0}, \dots , \phi _ {s,p} ) ^ {\text{T}}$ and $ \mathbf{ \widetilde{ x} } _{n-1} = (1, x_{n-1},\dots, x_{n-p})^{\text{T}}$. } The parameters $\theta_s = (\boldsymbol \phi _ s , \sigma _ s ^2)$ of the model consist of the AR coefficients $\boldsymbol\phi_s$ and the noise variance $\sigma_s^2$.
We place an inverse-gamma prior on $\sigma_s^2$ and a Gaussian prior on $\boldsymbol\phi_s$, so that the joint prior on the parameters is $\pi(\theta _s) =\pi (\boldsymbol \phi _s | \sigma _s ^2) \pi (\sigma _ s ^2)$, with,
where $(\tau, \lambda, \mu _o , \Sigma _o )$ are the prior hyperparameters. This prior specification allows the exact computation of the estimated probabilities of ((ref)), and also gives closed-form posterior densities for the AR coefficients and the noise variance. These are given in Lemmas (ref) and (ref), whose their proofs are given in Section (ref) of the Supplemental Material. Importantly, since here the estimated probabilities $P_e(s,x)$ are available in closed form, the GCTW, GBCT and $k$-GBCT algorithms can be used directly for exact Bayesian inference.
Consider executing the GCTW algorithm for a time series $x_1^n$ of length $n$. For each observation $x_i$, $1 \leq i \leq n$, exactly $D+1$ nodes of $T_{\text{MAX}}$ need to be updated (or created if they do not already exist in $T_{\text{MAX}}$), corresponding to the discrete contexts of length $0, 1,\ldots , D$ preceding $x_i$. 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. After repeating this step for every $1 \leq i \leq n$, using Lemma (ref) the estimated probabilities $P_e (s,x)$ can be computed for all nodes of $T_{\text{MAX}}$, whose size is bounded as a function of the number of observations $n$. And since the recursive step of GCTW only performs operations on $T_{\text{MAX}}$, it follows that its overall complexity as a function of $n$ is only $\mathcal{O}(n)$, i.e., {\em linear} in the length of the time series. The same argument applies to GBCT as well, showing that both algorithms are computationally very efficient and scale well with large time series.
{\color{black} Since, as described above, the size of $T_{\text{MAX}}$ is bounded by $nD+1$, it is easy to see (taking into account the number of operations required to compute $P_e (s,x)$ for each node and the recursive steps of the algorithms) that the complexity of GCTW and GBCT is in fact $\mathcal{O}\left (nD (m+p^3) \right )$. Importantly, this is linear in each of $n$, $m$, and $D$, and it is actually only slightly higher than the $\mathcal{O}\left (np^2\right )$ complexity of fitting a single AR model using least squares. This means that exact Bayesian inference can be performed within this much richer model class at an only marginally higher computational cost.
The above argument also shows that both GCTW and GBCT can be updated {\em sequentially}, as for every additional observation $x_{n+1}$ only $D+1$ nodes of $T_{\text{MAX}}$ need to be updated. In particular, sequential prediction can be performed efficiently. Empirical running times for all forecasting experiments are reported in Section (ref) of the Supplemental Material, showing that the BCT-X methods are much more efficient than essentially all the alternatives examined. The difference is quite large, especially compared to ML models that require heavy training and cannot be efficiently updated sequentially, giving empirical running times that are typically larger by several orders of magnitude; a general review or related issues is given in makridakis:18b.
Finally, it is possible to show that the memory requirements of both algorithms are also linear in the length $n$ of the time series. Arguing as before, since we only need to store the tree $T_{\text{MAX}}$ and the quantities $\{|B_s|, s_1, \mathbf s_2 , S_3\}$ for each node of $T_{\text{MAX}}$, it is easy to see that the memory requirements of both algorithms are $\mathcal{O}\left (nDp^2 \right )$, again linear in both $n$ and $D$. }
The posterior distributions of $\boldsymbol \phi _s$ and $\sigma_ s ^2$ are not particularly sensitive to the prior hyperparameters. In all the experimental results below, the simple choice $\mu_o =0$ and $\Sigma_o =I$ in the AR coefficients' prior is made. In view of ((ref)), $\tau$ and $\lambda$ should be chosen to be relatively small in order to minimise their effect on the posterior while keeping the mode of the inverse-gamma prior, $\lambda/(\tau + 1)$, reasonable; setting $\lambda = \tau =1$ is a sensible choice when no additional prior information is available. For the context-tree model prior, the default value $\beta= 1-2^{-m+1}$ BCT-JRSSB:22 and a maximum context depth $D=10$ are adopted. These default values are used in all the experiments, without involving any hyperparameter tuning from data.
Finally, a principled Bayesian approach is taken for selecting the quantiser thresholds $\{c_i\}$ of ((ref)) and the AR order $p$. Viewing them as extra parameters on an additional layer above everything else, we place uniform priors on $\{c_i\}$ and $p$, and perform Bayesian model selection rasmussen:00,mackay:92 to obtain their MAP values. The resulting posterior $p(\{c_i\},p|x)$ is proportional to the evidence $p(x|\{c_i\},p)$, which can be computed exactly using the GCTW algorithm. Specifically, a suitable range of possible $\{c_i\}$ and $p$ is 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). In all forecasting experiments, the AR order and the quantiser thresholds are selected using the above procedure at the end of the training set.
{\em Threshold AR models.} Threshold autoregressive (TAR) models tong:80 have been used extensively in the analysis of nonlinear time series; see, e.g., tong:11,hansen:11 and the texts cryer:book,tong:book.
{\color{black} The most commonly used version of TAR models is the self-exciting TAR (SETAR) model, which considers partitions of the state space based on the quantised value of $x_{n-d}$:
where $ e_n \sim \mathcal N (0, 1)$, $p$ is the autoregressive order, $Q : \mathbb {R} \to A = \{0,\dots , m-1 \}$ is an $m$-ary quantiser of the form in ((ref)), and $d$ is called the delay parameter. In other words, the SETAR model class considers partitions of the state space based (only) on the value of $x_{n-d}$, with different parameters $(\boldsymbol \phi ^ {(j)}, \sigma ^ {(j)} )$ associated to each region.
It is clear from the above description that the BCT-AR model class is a strict generalisation of SETAR. In specific, the BCT-AR model class always contains the SETAR models as specific instances, corresponding to particular trees in ${\cal T} (D)$; for example, any threshold model with delay parameter $d=1$ can be represented as a BCT-AR model with respect to the full tree of depth $d=1$. However, the BCT-AR class also contains other, more complicated tree-based partitions of the state space that cannot be represented as simple threshold models. For example, asymmetric BCT-AR models with leaves at various depths define more complicated partitions that cannot be represented as linear combinations of threshold models. In fact, the BCT-X methodology can be viewed as a natural conceptual extension of the family of threshold models, which allows for a more systematic and powerful Bayesian way of breaking up the state space in possibly more -- and more complex -- regions, and then fitting a different time series model to each one of these regions. }
{\em Mixture AR models.} The mixture autoregressive (MAR) models of wong:00 consist of a simple linear mixture of $K$ Gaussian AR components. When the BCT-AR posterior distribution essentially concentrates on $K$ models, $T_1, \dots , T_K$, (which is both theoretically `allowed' and is commonly observed in practice), the posterior predictive distribution can be expressed as, \[ p(x_{n+1}|x) = \sum _{k=1} ^ K \pi (T_k|x) \ p(x_{n+1}|T_k, x) . \vspace*{-0.05cm} \] In this sense, BCT-AR can be viewed as a generalised MAR model, with components corresponding to the AR models at the leaves of each $T_k$, and Bayesian weights given by $\pi (T_k|x) $. Therefore, the BCT-AR model class is a strict generalisation of both MAR and SETAR.
{\em Markov switching AR models.} The Markov switching autoregressive model (MSA) hamilton:89 is a simple HMM, where the hidden state process is a discrete-valued first-order Markov chain (usually binary or ternary), and a different AR model is associated to each state. As the discrete states here are not observable tsay:05, the main difficulty in using the MSA is in estimating the model; the EM algorithm is usually employed for this task. This difficulty was also observed in practice in our experiments, where the MSA gave much larger running times compared to all other classical methods that do not involve neural networks (Supplemental Material, Section (ref)).
{\color{black} In this section, the performance of the BCT-AR model is evaluated on simulated and real-world datasets from standard applications of nonlinear time series in economics and finance. The complete descriptions of all datasets can be found in Section (ref) of the Supplemental Material. In all forecasting experiments, the training set consists of the first 50% of the observations in each dataset, and we consider out-of-sample 1-step ahead forecasts, allowing every model to be updated at every timestep. For the BCT-AR model, the current MAP tree model with its MAP estimated parameters is used at every timestep for forecasting.
Since, in our experiments, the posterior $\pi(T|x)$ on model space was typically found to concentrate on the MAP tree model with high posterior probability, this simple prediction rule already gives a reasonably good -- and computationally efficient -- approximation to the full Bayesian predictive distribution. In general, however, it is noted that prediction based on model-averaging over trees $T$ and over parameters $\theta$, may also be easily implemented.
The BCT-AR model is compared with the most successful previous methods for these applications, considering both classical and modern ML methods. Useful resources include the R package forecast hyndman:08 and the Python library `GluonTS' alexandrov:20, containing implementations of state-of-the-art classical and ML methods, respectively. We briefly discuss the methods used, and refer to the packages' documentation and Section (ref) of the Supplemental Material for more details on the methods and the training procedures carried out. Among classical statistical approaches, we compare with ARIMA and Exponential smoothing state space (ETS) models (implemented in forecast), with SETAR -- using the conditional least squares (CLS) method implemented in the R package TSA tsa -- and with MAR and MSA models, using the \texttt{R} packages \texttt{mixAR} mixar and \texttt{MSwM} mswm. Among ML-based techniques, we compare with the Neural Network AR (NNAR) model (implemented in \texttt{forecast}), and with the most-successful RNN-based approaches, deepAR salinas:20 and N-BEATS oreshkin:19 -- both implemented in GluonTS. }
Here we use simulated data to illustrate that the BCT-X methods are consistent and effective with data actually generated by BCT-X models. A time series $x$ is simulated from the BCT-AR model with the context-tree model $T$ in Figure (ref), a binary quantiser with threshold $c=0$, and the following AR(2) base models at the three states defined by $T$:
We first examine the posterior distribution $\pi(T|x)$ over $T\in{\cal T}(D)$. With $n=100$ observations, the MAP context-tree model identified by the GBCT 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 structure with multiple states. With $n=300$ observations, the MAP tree model is now the true underlying model, with posterior probability 57%. And with $n=500$ observations, the posterior of the true model is 99.9%. The complete BCT-AR model fitted from $n=1000$ observations with its MAP parameter estimates is shown in ((ref)). In Section (ref) of the Supplemental Material we also report values of the evidence $p(x|c,p)$, which is maximised at the true values of $c=0$ and $p=2$.
{\bf Interpretation.} Clearly, the three relevant states, the optimal quantiser, and the correct AR order were all identified without any prior training, based on only $n=1000$ observations. Also, the posterior distribution of the AR parameters appears to be well-concentrated around the true underlying values, verifying that all our inferential procedures are effective.
{\bf Forecasting.} The performance of the BCT-AR methods is evaluated in the task of out-of-sample 1-step ahead forecasts, and it is compared with state-of-the-art approaches in three simulated and three real datasets. The first simulated dataset (sim_1) consists of $n=600$ observations generated from the BCT-AR model used above, the second (sim_2) has $n=500$ observations generated by a BCT-AR model with a ternary context-tree model of depth 2, and the third (sim_3) consists of $n=200$ observations generated from a SETAR model of order $p=5$; see Section (ref) of the Supplemental Material.
The results in Table (ref) indicate that the BCT-AR model outperforms all the alternatives in all three simulated experiments, achieving a mean-squared error (MSE) that is lower by between 7% and 37% compared to the second-best method. As discussed in Section (ref) (see also Section (ref) of the Supplemental Material for more detailed results), the BCT-X methods also outperform the alternatives in terms of empirical running times, by anywhere between one and three orders of magnitude. Next, we examine the performance of the BCT-AR methods in real-world applications from economics and finance.
An important application of SETAR models is in modelling the US unemployment rate hansen:11,montgomery:98,tsay:05,rothman:98,koop:99. As described in montgomery:98,tsay:05, the unemployment rate moves countercyclically with business cycles, and rises quickly but decays slowly, indicating nonlinear behaviour. Here, we examine the quarterly US unemployment rate in the period from 1948 to 2019 (dataset unemp, 288 observations). Following montgomery:98, we consider the difference series ${\Delta x}_n = x_n-x_{n-1}$, and also include a constant term in the AR model. For the quantiser alphabet size, $m=2$ is a natural choice here, as will become apparent below. The threshold selected using the procedure of Section (ref) is $c=0.15$, and the resulting MAP context-tree model is the tree of Figure (ref), with states $\mathcal {S} = \{ 1, 01, 00\}$ and posterior probability 91.5%. The complete BCT-AR model with its MAP parameters is given below, where $e_n \sim \mathcal N (0, 1)$, \[ {\Delta x}_n \!=\! \left\{
\right. \]
Interpretation. The MAP BCT-AR model finds significant structure in the data, providing a natural interpretation. It identifies 3 meaningful states: First, jumps in the unemployment rate higher than 0.15 signify economic contractions ($s=1$). If there is not a jump at the most recent time-point, the model looks further back to determine the state. The state $s=00$ signifies a stable economy, as there are no jumps in the unemployment rate for two consecutive quarters. Finally, $s=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. It is higher in contractions ($\sigma = 0.42$), smaller in stable economy regions ($\sigma = 0.20$), and in-between for state 01 ($\sigma = 0.32$). This is an important finding, verifying the economists' understanding that there is much higher uncertainty during economic contractions.
{\bf Forecasting.} In addition to its appealing interpretation, the BCT-AR model outperforms all benchmarks in forecasting (together with MSA), achieving the lowest MSE (Table (ref)). In terms of empirical running times, the BCT-AR model vastly outperforms MSA; see Section (ref) of the Supplemental Material.
Another important example of nonlinear time series in economics is the US Gross National Product (GNP) potter:95,hansen:11. We examine the quarterly US GNP in the period from 1947 to 2019 (dataset gnp, 291 observations). Following potter:95, here we consider the difference in the logarithm of the series, $y_n = \log x_n - \log x_{n-1}$. 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 context-tree model is shown in Figure (ref). It has maximum depth 3, the states are $\mathcal {S} = \{0, 10, 110, 111 \}$, and its posterior probability is 42.6%. The complete set of MAP parameters at the leaves are shown in Section (ref) of the Supplemental Material.
{\bf Interpretation.} Compared with the previous example, the MAP BCT-AR model identifies even richer structure in this dataset, with 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 increases in expansions and decreases in contractions. And again, the model does not look further back whenever a contraction is detected -- it is a {\em renewal event}. In this example, the model shows that the effect of a contraction is still present even after {\em 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. An important feature captured by the model is again that the volatility is different in each case and, enhancing previous findings, it is found to decrease with the distance from the last contraction. It starts at $\sigma =1.23$ for $s=0$ and decreases to $\sigma = 0.75$ for $s=111$ (see Section {(ref)} of the Supplemental Material).
{\bf Forecasting.} As shown in Table (ref), the BCT-AR model is found to outperform all benchmarks in forecasting in this dataset, likely due to the additional structure identified, giving a significantly lower MSE than the second-best method (by 9%).
We revisit the daily IBM common stock closing price from May 17, 1961 to November 2, 1962 (dataset ibm, 369 observations). This is a well-studied dataset, with box:book fitting an ARIMA model, tong:book fitting a SETAR model, and wong:00 fitting a MAR model to the data. Following previous approaches, we consider the first-difference series, ${\Delta x}_n = x_n-x_{n-1}$. For the alphabet size of the quantiser we choose $m=3$, with the values $\{ 0,1,2\}$ naturally corresponding to the stock price changes \{down, steady, up\}. Using the procedure of Section (ref), the resulting quantiser thresholds are $\{-7,+7\}$ and the AR order $p=2$. The MAP context-tree model is shown in Figure (ref). It has maximum depth equal to 2, and it identifies five states, $\mathcal {S} = \{0,2,10,11,12 \}$. Its posterior is 99.3%, suggesting that there is very strong evidence in the data supporting this structure, even with only 369 observations. The complete BCT-AR model with its MAP parameters is given the Supplemental Material (Section (ref)).
{\bf Interpretation.} The BCT-AR model reveals important information about apparent structure in the data, which has not been identified before. Firstly, it admits a 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 the states 0, 2, 10, 12), or until reach the maximum depth of 2 is reached (state 11). Another important feature captured by this model is the commonly observed asymmetric response in volatility due to positive and negative shocks, sometimes called the {\em leverage effect} tsay:05,box:book. Indeed, the MAP model shows that negative shocks increase the volatility much more: the state $s=0$ has the highest volatility ($\sigma =12.3$), with $s=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 (between $\sigma = 5.17$ and $\sigma = 6.86$). Finally, we observe that when stabilising after a shock (states 10, 12), the latest value $x_{n-1}$ is not as important as $x_{n-2}$, whereas $x_{n-1}$ is dominant in all other cases.
{\bf Forecasting.} From Table (ref), it is observed that DeepAR outperforms all methods here, with BCT-AR giving a marginally higher MSE but much smaller empirical running times; see Section (ref) of the Supplemental Material.
A key aspect of financial time series analysis is modelling the dynamic evolution of volatility over time. To capture the well-known {\em volatility clustering} present in financial data, Engle's seminal work introduced the autoregressive conditional heteroscedasticity (ARCH) models engle:82, which, together with their numerous extensions, have been very widely used for modelling the volatility in financial time series. The perhaps more important limitation of ARCH models is their inability to describe the well-known asymmetric response in volatility due to positive and negative shocks -- the {\em leverage effect} tsay:05,box:book mentioned in the previous section. A number of approaches have attempted to incorporate this feature in ARCH models, notably including the works of nelson:91,glosten:93,zakoian:94,gray:96; see Section (ref) for more details.
In this section, as a second example of the general BCT-X class, we associate a different ARCH model to each state of the context-tree model, and refer to the resulting mixture model class as the {\em Bayesian Context Trees ARCH} (BCT-ARCH) model. The BCT-ARCH model is of high practical interest as it offers a systematic and powerful way of modelling volatility asymmetries. As shown in Section (ref), it is found to outperform previous approaches in real examples. Moreover, it reveals important structure present in the data that not been identified before, in the form of an {\em enhanced leverage effect}.
Given a quantiser $Q:{\mathbb R}\to A=\{0,1,\ldots,m-1\}$, a fixed ARCH order $p$, and a maximum context depth $D\geq 0$, the BCT-ARCH model is the version of BCT-X with ARCH$(p)$ models as base modes. With each potential state $s$ of a context-tree model $T\in{\cal T}(D)$ we associate an ARCH$(p)$ model,
where $\boldsymbol \alpha _ s = \left ( \alpha _ {s,0} , \alpha _{s,1} , \dots , \alpha _{s,p} \right )^{\text{T}}$ and $\mathbf{ { z } } _{n-1} = (1, x_{n-1} ^ 2,\dots, x_{n-p} ^ 2 )^{\text{T}}.$
Here, the parameters corresponding to each state $s$ are the ARCH coefficients, $\theta _ s = \boldsymbol \alpha _ s $. We follow vrontos:00 and use the following non-informative priors: $\pi ( \alpha _ {s,0}) = 1 / \alpha _ {s,0}$, and $\pi ( \alpha _ {s,j})\sim U(0,1)$ for $1\leq j\leq p$. Note that no hyperparameters need to be selected for this prior.
Our main tools for inference -- the GCTW and GBCT algorithms -- rely on the evaluation of the estimated probabilities $P_e(s,x)$ given in ((ref)). However, in this case these are not available in closed form. In fact, exact Bayesian inference is not possible even with a single ARCH model, and MCMC approaches are typically adopted for this task; see virbickaite:15 for a review. Therefore, we will employ effective approximations for the $P_e(s,x)$ and use them in GCTW and GBCT; the remaining methodology remains the same as before.
As noted earlier, the fact that $P_e(s,x)$ is in the form of marginal likelihoods allows us to use standard methods to approximate them. Specifically, at each node $s$ of $T_{\rm MAX}$ we use a variant of the Laplace method described in kass:95, which was found to be effective in the case of a multivariate GARCH model in vrontos:03. The resulting approximation is:
Here, $ L _ s ( {\theta} _ s )$ is the log-likelihood of data with context $s$,
$\widehat {\theta} _ s$ is its maximiser, which can be computed iteratively using the Fisher scoring algorithm,
and $\widehat {I} _ s$ the expected information matrix computed at $\widehat {\theta} _ s ^ { \ (k-1)}$. These quantities for the BCT-ARCH model are given in Lemma (ref), that is proven in Section (ref) of the Supplemental Material.
\textcolor{black}{ It is noted that another general alternative to the Laplace method for approximating the integrals of ((ref)) would be to employ MCMC as previously used in the literature to approximate marginal likelihoods; see, e.g., chib:95,chib:01. } Once the approximate values $\widehat {P_e} (s,x)$ of the estimated probabilities $P_e (s,x)$ have been computed as outlined, inference can be carried out exactly as in the case of the BCT-AR model. In particular, the GCTW, GBCT and $k$-GBCT algorithms can be used with these approximations for $P_e (s,x)$, and the selection of the ARCH order $p$ and the quantiser thresholds $\{c_i\}$ can be performed in the same way as before, using the method described in Section (ref) for the BCT-AR model.
In all the experiments presented in Section (ref), we take the number of iterations $M$ in ((ref)) to be equal to 10; see Section (ref) of the Supplemental Material for a discussion of this choice.
Consider, as before, executing the GCTW or GBCT algorithm for a time series $x_1 ^n$. Assuming for the moment that only one iteration of the Fisher scoring algorithm is performed for each node $s$ of $T_{\text{MAX}}$, the situation is the same with the BCT-AR case: For each observation $x_i$, $1 \leq i \leq n$, exactly $D+1$ nodes need to be updated, and for each one of these nodes, only the quantities in ((ref))--((ref)) need to be updated, which can be done efficiently by just adding one term to each sum. The only difference in the BCT-ARCH case is that this forward pass of the data for $1 \leq i \leq n$ now needs to be repeated for every iteration of the scoring algorithm, each time followed by the Fisher updates of ((ref)) for the nodes of $T_{\text{MAX}}$. Finally, the resulting estimates $\widehat {\theta} _ s$ can be used to compute the approximation $ \widehat {P_e} (s,x)$ of ((ref)) for all nodes $s$ of $T_{\text{MAX}}$, while the recursive step remains identical with before.
{ \color{black} Arguing as in Section (ref), and denoting the number of Fisher iterations as $M$, it is not hard to show that the overall complexity of the GCTW and GBCT algorithms for the BCT-ARCH model is $\mathcal {O} (nD (m+Mp^3))$. Importantly, this is still linear in all $n$, $m$, and $D$, and is actually only slightly higher than the $\mathcal {O} (nMp^2)$ complexity of fitting a single ARCH model, showing that, again, inference in this much richer model class can be performed at a negligible additional computational cost. Also, it is easy to see that the memory requirements of the algorithms are again $\mathcal {O} (nDp^2)$, i.e., linear in both $n$ and $D$. }
In Section (ref) we will compare the performance of the BCT-ARCH model with that of some of the most successful and commonly used alternative methods for modelling volatility. A brief summary of these methods is given here. Additional information can be found in the R packages rugarch rugarch, MSGARCH gray:96, and stochvol HK:21, and in Section (ref) of the Supplemental Material.
{\em GARCH and EGARCH models.} The most widely used extension of ARCH models are the Generalised ARCH (GARCH) models of bollerslev:86. The GARCH$(p,q)$ model is given by,
with the simple GARCH$(1,1)$ model being the most popular in practice. An important extension of the GARCH model is the Exponential GARCH (EGARCH) model of nelson:91, where the parametrisation is in terms of the logarithm of $\sigma _ n ^2$.
{\em GJR models.} A common alternative to GARCH, aimed at explicitly capturing volatility asymmetries, is the threshold model of Glosten, Jagannathan, and Runkle (GJR) glosten:93, given by,
where the indicator function forces negative returns to have higher volatility. A similar model parametrised in terms of the standard deviation instead of the variance is given in zakoian:94. \textcolor{black}{From the above description it is easy to see that the GJR model can be viewed as a particular case of a threshold GARCH model. So, the BCT-ARCH model class is more general than GJR, since it allows for more complicated partitions of the state space, and it is hence expected to perform better in practice; see also the corresponding and more detailed discussion and comparison with threshold models that is carried oout for the BCT-AR model in Section (ref).}
{\em MSGARCH models.} The structure of the Markov switching GARCH (MSGARCH) gray:96,haas:04 model is identical to that of MSA, except that GARCH models replace the AR models associated to each discrete hidden state. As with the MSA, performing inference in this setting is much more challenging compared to the previous approaches.
{\em SV models.} An alternative to the ARCH family -- which models the evolution of volatility deterministically -- is the family of {\em stochastic volatility} (SV) models, which model the volatility as a random process. This is usually done using an HMM where the hidden state is the logarithm of the variance, modelled as an AR process; see shephard:09,HK:21 for more details. As with MSGARCH, the randomness present in the hidden state process makes inference much more challenging and computationally more expensive compared to the other approaches.
We begin by examining the performance of the BCT-ARCH inference methods on data generated from a model within this class. Consider a time series $x$ consisting of observations simulated from a BCT-ARCH model where the context-tree model is the binary tree of depth 1 with states $\mathcal {S} = \{ 0,1 \}$, and with the associated ARCH(2) base models:
This is an intentionally difficult example, with two very similar ARCH models in the two regions. Let the maximum context depth $D=5$ and consider a binary quantiser with threshold $c$. With $n=1000$ observations the MAP context-tree model is the `empty' tree, corresponding to a single ARCH model, as there are not enough data to reveal a more complex structure. With $n=2500$ observations, the MAP context-tree model is now the true underlying model with a posterior probability of 42%, and with $n=5000$ observations the posterior probability of the true model becomes 90%. The parameter estimates from $n=5000$ observations are,
Here, the ARCH order $p=2$ and the quantiser threshold $c=0$ were selected as described in Section (ref). More extensive results on simulated data are shown in Section (ref) of the Supplemental Material, illustrating that the posterior probability of the true underlying model converges to 1, and that the estimates of the parameters converge to the true parameter values.
In this section, the BCT-ARCH framework is used to model the volatility on four real-world datasets corresponding to time series describing four major financial indices. In each example, the $n=7821$ daily values of a stock market index are examined, during a period of thirty years ending on 7 April 2023. Similarly to Section (ref), we examine the transformed time series, $y_n = 10 [\log x_ n -\log x_{n-1}]$. In order to explicitly capture the leverage effect and distinguish between positive and negative shocks, a binary quantiser with threshold $c=0$ is used. The two “memory length” parameters, namely, the maximum context depth $D$ and the ARCH order $p$, are both taken to be equal to 5, corresponding to a week of trading days.
{\bf Results: New structure.} The most important feature captured by the BCT-ARCH model used on stock market index data is an {\em enhanced leverage effect}. The fitted models suggest that, in understanding the asymmetries present in volatility, it is not only the sign of the most recent value-change that is relevant, but the exact pattern of recent `ups' and `downs' is also important. These relevant patterns are identified from the data as the {\em states} given by the leaves of the fitted MAP context-tree model. In all four cases examined here, that set of relevant states was found to be richer than $ \mathcal {S}=\{0,1\}$, which would correspond to just the sign of the previous change as in traditional leverage modelling. This strongly suggests the conclusion that modelling such enhanced leverage is required in practice, and that the BCT-ARCH model reveals this essential structure that has not been identified before.
{\bf FTSE 100.} First we examine the Financial Times Stock Exchange 100 Index, which is the most commonly used UK-based stock market indicator, including 100 companies listed on the London Stock Exchange. The fitted BCT-ARCH model exhibits the enhanced leverage effect, as it identifies the relevance of three meaningful states.
The MAP context-tree model given by the GBCT algorithm is shown in Figure (ref); it has depth 2 and three leaves, $\mathcal {S} = \{ 0, 10, 11\}$. Its posterior probability is 95.2%, signifying that there is very strong evidence in the data supporting this exact structure. State $s=0$ corresponds to a negative shock at the last timestep, $s=10$ to stabilising just after a negative shock whose effect is still present, and $s=11$ to a flourishing period. The complete BCT-ARCH model is given by,
Among the three states, the ARCH coefficients of state $s=11$ (where no negative shocks are present) are the smallest for all five lags. States $s=0$ and $s=10$ have very similar coefficients except for the one corresponding to the first lag, $y _{n-1} ^2$, which is essentially zero for $s=10$, i.e., when it corresponds to an increase in value. These observations are consistent with the common understanding that negative shocks lead to greater increases in volatility.
{\bf CAC 40 and DAX.} Next we examine CAC 40 and DAX, the main stock market indices in France and Germany, respectively, each one consisting of 40 major companies. The model fitted to the CAC and DAX time series finds very similar structure present in the two datasets: The same MAP context-tree model is identified by the GBCT algorithm in both cases (Figure (ref)). It has depth $3$ and four leaves, $\mathcal {S} = \{ 0, 10, 110, 111\} $. The estimated ARCH coefficients for the two models are given in Sections (ref) and (ref) of the Supplemental Material.
The interpretation of this result is similar to that in the case of FTSE, with negative shocks again playing the role of a `renewal' event: In order to determine the distribution of the current state, the model looks back into the past until the first time a negative shock is detected. The only difference is in the memory of the discrete state process, with negative shocks now being relevant even if they occur {\em three} timesteps in the past. This gives a total of four different states, meaning that some additional structure has been identified. In each case, the ARCH coefficients corresponding to lags with positive changes are small: For example, for $s=10$ the coefficient $\alpha _ 1 \approx 0$, while for $s=110$ both $\alpha _ 1 $ and $\alpha _ 2$ are small. Overall, the BCT-ARCH model again exhibits the enhanced leverage effect, and it gives strong evidence of a rich asymmetric response in volatility, with negative shocks having a stronger effect.
{\bf S&P 500.} The last index we examine is Standard and Poor's 500 (S&P 500), one of the most commonly followed indices worldwide, consisting of 500 of the largest companies in the US. Here, the MAP context-tree model (Figure (ref)) again displays the enhanced leverage effect, and it reveals even more structure compared to the previous cases. It has depth $3$, five leaves, $\mathcal {S} = \{ 00, 01, 10, 110, 111 \}$, and a posterior probability of $91.4 \%$. In fact, the MAP context-tree model is the same tree as for CAC and DAX, but with one additional branch added at $s=0$. Apart from identifying slightly more structure, the interpretation of the BCT-ARCH model is very similar with before. The values of the estimated ARCH coefficients (reported in Section (ref) of the Supplemental Material) are small when they correspond to lags with positive changes, again describing an asymmetric volatility response.
In this section, the performance of the BCT-ARCH model in forecasting is illustrated and contrasted with that of the alternative methods outlined in Section (ref). As the volatility is not directly observable, effectively comparing different models is known to be a challenging task tsay:05. Following standard approaches geweke2010comparing,vrontos:03b,dellaportas:07, in order to measure the relative predictive ability of the models, we consider the predictive distributions in 1-step ahead out-of-sample forecasting.
Following dellaportas:07, the last 130 observations are taken as the test set in each of the four datasets, corresponding to the trading days during a period of six months in each case; the first 7691 observations are used as the training set. \textcolor{black}{ At every timestep $i$ in the test set, all models are estimated using the entire past, and the resulting predictive density $\widehat {p} ( y _i | y_{1} ^ {i-1} )$ is evaluated at the next test datapoint, $y_i$. The complete details of the training process employed for each method are given in Section (ref) of the Supplemental Material. For the BCT-ARCH model, as in the case of the BCT-AR model earlier, we use the MAP tree model with its MAP estimated parameters for forecasting. Also, it is noted that, throughout this section, the threshold of the quantiser is still fixed at $c=0$ in order to explicitly differentiate between positive and negative shocks; but in terms of forecasting performance, the results of BCT-ARCH could possibly be further improved by also estimating the threshold from data in the usual way.}
{\color{black} As is standard practice geweke2010comparing,vrontos:03b, we examine the logarithm of the predictive density, i.e., $$ \vspace*{-0.2 cm} \mathcal {L} = \sum _ {i } \log \widehat {p} ( y _i | y_{1} ^ {i-1} ), \vspace*{-0.0 cm} $$ evaluated over all datapoints $y_i$ in the test set, in one-step ahead out-of-sample forecasts. }
The results for the four stock market indices are presented in Table (ref). The BCT-ARCH model is seen to outperform all the alternatives in all examples, the only exception being the S&P index data. Because of its low computational complexity and efficient sequential updates, the BCT-ARCH model also outperforms all the alternatives in terms of its computational requirements, giving empirical running times that anywhere between one and three orders of magnitude smaller than those of the alternatives; see Section (ref) of the Supplemental Material.
Finally, we further validate the findings of the previous section, namely that the BCT-ARCH model consistently outperforms the alternatives in forecasting the volatility of stock market indices -- mainly because of its ability to model asymmetries in a more flexible and systematic way. In this section we repeat the above experiment for a number of important stock market indices and test for statistical significance of the results.
Specifically, apart from the S&P 500 index we repeat the above experiment for 20 major European stock market indices (including the FTSE 100, CAC 40 and DAX as before); again, a complete set of the details for each dataset/index is given in Section (ref) of the Supplemental Material. In this section, we exclude from our comparisons the simple ARCH, GARCH and SV models, since they are more restricted model classes that do not account for volatility asymmetries, and as a result they were all found to perform consistently worse than the other methods in the previous section. {\color{black} The log predictive density results for these 21 major stock market indices are presented in Table (ref). }
The results of Table (ref) clearly show that the BCT-ARCH model outperforms the alternatives in volatility forecasting: It achieves the best performance in 15 out of 21 datasets, and it also achieves the best average ranking overall, which is 1.67 compared to 2.67 for the second-best EGARCH model. Further, we examine the statistical significance of these findings by implementing post-hoc Nemenyi tests hollander2013nonparametric,demvsar2006statistical using the R package tsutils tsutils; the results are presented in Figure (ref).
The Friedman test hollander2013nonparametric rejects the null hypothesis that all methods perform similarly, and the post-hoc Nemenyi tests shown in Figure (ref) further justify our findings. At the 90% confidence level, there is a statistically significant difference between the performance of the BCT-ARCH model and that of any other method. In particular, the Multiple Comparison with Best (MCB) plot koning2005m3 indicates that, in each case, the corresponding difference in average ranking is greater than the critical distance of the test. Moreover, the matrix plot kourentzes2019cross suggests that all other methods (EGARCH, GJR, MSGARCH) have comparable performance (among them), as there is not sufficient evidence of statistically significant differences between any pair of them.
{\color{black}
This work develops a general Bayesian framework for building flexible and interpretable mixture models for real-valued time series, that are based on context trees. The proposed framework can be combined with any existing model class as a base model, resulting in a much richer class of flexible mixture models, for which we provide algorithms that allow for Bayesian inference at a negligible additional computational cost compared to the original model. The utility of the proposed methodology has been illustrated by using AR and ARCH models as the base model, in both cases resulting in flexible mixture models of high practical interest.
The generality of the proposed framework leads to several possible directions for future work: For any given application, BCT-X can be combined with any existing state-of-the-art model to provide much greater modelling flexibility as well as potentially significant improvements in forecasting performance. As a few examples, candidate base model classes include ARIMA, EGARCH, general state space models, and MAR models, potentially leading to new and powerful ways to model feature-diverse time series datasets kang2020gratis. Similarly, BCT-X could be combined with modern ML models, including GPs, neural networks, and Deep Learning methods like DeepAR. In a different direction, forecasting performance might be improved by employing combination tools, ranging from simple averaging methods to modern ensemble learning techniques like bagging and boosting, which might lead to important practical improvements. Closing, we remark that the entire BCT-X framework can be extended to the multivariate time series setting, something that would greatly broaden the scope of potential applications. }
{
}