Bayesian Forecasting of Dynamic Extreme Quantiles

Abstract

In this paper, we provide a novel Bayesian solution to forecasting extreme quantile thresholds that are dynamic in nature. This is an important problem in many fields of study including climatology, structural engineering, and finance. We utilize results from extreme value theory to provide the backdrop for developing a state-space model for the unknown parameters of the observed time-series. To solve for the requisite probability densities, we derive a Rao-Blackwellized particle filter and, most importantly, a computationally efficient, recursive solution. Using the filter, the predictive distribution of future observations, conditioned on the past data, is forecast at each time-step and used to compute extreme quantile levels. We illustrate the improvement in forecasting ability, versus traditional methods, using simulations and also apply our technique to financial market data.

1. Introduction

Forecasting extreme quantiles, which is a part of extreme value (EV) analysis, is an important problem in many applications. For example, in analyzing ocean storm severity it is often critical to forecast an extreme quantile, such as the 99% quantile, for the maximum wave height over a certain time period to facilitate adequate construction design [1]. In financial risk-management, forecasting extreme quantiles, or so-called VaR analysis, is critical for establishing adequate bank capital levels [2]. In fact, extreme value analysis, and quantile forecasting, has found wide application in fields as varied as target detection [3], communication systems [4], image analysis [5], power systems [6], and population studies [7].

In this paper, which is an extension of the work presented in [8], we derive a recursive method for Bayesian forecasting of extreme quantiles where the underlying process is non-stationary. This differs from our previous work given our focus here is on the predictive distribution of future observations, rather than simply the model parameters, and we also properly invert the predictive distribution to obtain efficient quantile forecasts. In our previous work, we obtained improved parameter estimates, versus traditional methods, but our prediction ability, as measured by predictive p-values, was hampered by a number of issues. We made a number of improvements to our previous work including modifying the support assumptions on the observed data as well as implementing an adaptive particle filtering algorithm both of which resulted in demonstrable improvement in quantile forecasting.

We model the maximum of a block of data, which has an asymptotic general extreme value (GEV) distribution [9], via a non-linear state-space model where the parameters are driven by a hidden stochastic process with unknown static parameters. To recursively estimate the model parameters, we utilize a particle filter (PF) [10] and, in particular, we derive a Rao-Blackwellized particle filter (RBPF) [11] to marginalize the unknown, static parameters. Importantly, we derive a recursive solution for the predictive state density, which is required for the particle filter, and we design an algorithm that bears similarity to the well-known Kalman filter [12] and is therefore readily implemented.

The work presented here, and in [8], is an extension of some previous studies. In [13], a deterministic trend is applied to a subset of the GEV parameters, while in [14,15,16], a GEV parameter subset is dynamic with known state equation. More recently, Ref. [17] developed a dynamic model for the tail-index and a Gaussian-mixture approximation to linearize the estimation problem. Our work has extended these previous studies in a number of ways. First, under a reasonable assumption, we reduce the number of GEV parameters and model the remaining set as a vector Markov process, with unknown system and covariance matrix. Second, we recursively compute the marginalized state density, eliminating the need to estimate unwanted nuisance parameters and, in the process, we derive recursive expressions for the necessary sufficient statistics. This allows for a fast, real-time implementation without the need to batch-process observations. Lastly, we derive the predictive distribution, for the block-maximum, that we use to forecast extreme quantiles.

The paper is organized as follows. In the next section, we formulate the problem for the time-varying block-maximum and propose a Bayesian approach by deriving the recursive solution for marginalizing the nuisance parameters-Rao-Blackwellization. In Section 3, we implement our solution and derive our simplified likelihood function. We describe our algorithm in detail and show the method’s ability to forecast extreme quantiles, outperforming traditional methods, using simulated data. In the last section, we discuss our results and use our approach to analyze S&P 500 stock market returns from 1928 to 2020. We also conclude the paper with ideas for further research.

2. Materials and Methods

The starting point for our analysis is a sequence of block-maximums denoted by $y_k\in \mathbb{R},k=1,\dots N$. The k’th block-maximum is the maximum of a strictly stationary time series, $s_t,$ so that

$$y_k=\underset{N_{k-1}<t\le N_k}{max}{s}_t,$$

(1)

where the index t typically represents time and the k’th block consists of data with time indices $(N_{k-1},N_k]$ and $N_0=0$. An example would be in financial markets where the underlying data is the daily return of a stock market index (e.g., S&P 500) and $y_k$ is the largest daily return loss experienced in year k (see Figure 1). In this case, $s_t=-r_t$, where $r_t=(p_t-p_{t-1})/p_{t-1}$ is the daily return defined as the percentage change in the price, $p_t$.

Similar to the central-limit theorem, there is a limiting distribution for the normalized block-maximum of i.i.d. random variables (RVs). The Fisher-Tippet-Gnedenko (FTG) theorem states that the only non-degenerate limiting distribution for the block-maximum, is in the generalized extreme value (GEV) family [18] with shape parameter, or tail-index, $\xi$. That is, the cumulative distribution function (cdf) of the block-maximum, $y_k$, properly normalized, converges to a distribution $H_{\xi }$ as $(N_k-N_{k-1})\to \infty$ which has the Jenkinson-von Mises representation [19],

$$H_{\xi }\left(y\right)=exp(-{(1+\xi y)}^{-1/\xi }).$$

(2)

When the tail-index is positive ($\xi >0$), the resulting distribution is the Fréchet distribution with $y\ge -1/\xi$. The Fréchet is the limiting distribution for many heavy-tailed underlying random variables [20]. Upon normalization, a three parameter family is obtained and, therefore, we asymptotically model each of the block-maximum $y_k~H_{\xi ,\beta ,\mu }\left(y_k\right)=H_{\xi }\left(\frac{y_k-\mu }{\beta }\right)$, with $\beta >0$ and support $y\in [\mu -\beta /\xi ,\infty )$. We note that the FTG theorem’s i.i.d. assumption can be relaxed to include most strictly stationary processes [21].

To capture the non-stationary effects and clustering often witnessed in real-world data (see Figure 1), we propose to model the parameters of the GEV distribution as time-varying and, to insure positivity for the shape and scale parameters, we define the unknown state vector as

$${\mathbf{x}}_k={[log\left({\xi }_k\right),log\left({\beta }_k\right),{\mu }_k]}^{\top }$$

(3)

and the state transition equation as

$${\mathbf{x}}_k=\mathbf{\eta }+\Theta {\mathbf{x}}_{k-1}+{\mathbf{C}}^{1/2}{\mathbf{u}}_k,$$

(4)

where ${\mathbf{x}}_k\in {\mathbb{R}}^3$, $k\in {\mathbb{N}}_0$, ${\mathbf{u}}_k~\mathcal{N}(0,{\mathbf{I}}_3)$, $\mathbf{\eta }\in {\mathbb{R}}^3$, and $\Theta ,\mathbf{C}\in {\mathbb{R}}^{3\times 3}$. We assume the matrices $\mathbf{\eta },\Theta ,$ and $\mathbf{C}$ are unknown and the Gaussian process noise, ${\mathbf{u}}_k$, was chosen for its analytical tractability and maximal entropy property. Modeling the state as a first-order auto-regressive process, as in Equations (3) and (4), is similar to stochastic volatility models [22].

We specify the observation, or measurement, equation as

$$y_k=x_{3_k}+e^{x_{2_k}}{w}_k,$$

(5)

where $e^{x_{2_k}}={\beta }_k,x_{3_k}={\mu }_k$, and $w_k$ is distributed according to (2) and, hence, a function of $e^{x_{1_k}}={\xi }_k$. Specifically, $y_k$ is a standard GEV RV that is scaled and translated and, combined, Equations (4) and (5) form a non-linear state-space model from which we wish to make Bayesian inferences.

The Bayesian inferences we wish to make come in the form of three probability density functions (pdf). The first is the state filtering density $p\left({\mathbf{x}}_k\right|y_{1:k})$ where $y_{1:k}\equiv \{y_1,y_2,\dots ,y_k\}$. The second is the state predictive density $p\left({\mathbf{x}}_{k+1}\right|y_{1:k})$. Lastly, and most germane to our goal, is the predictive distribution of the observation

$$p\left(y_{k+1}\right|y_{1:k})=\int p\left(y_{k+1}\right|{\mathbf{x}}_{k+1})p\left({\mathbf{x}}_{k+1}\right|y_{1:k})d{\mathbf{x}}_{k+1}.$$

(6)

Since our goal is to estimate extreme quantiles, we require the cdf of this predictive density,

$$F\left(y_{k+1}\right|y_{1:k})={\int }_{-\infty }^{y_{k+1}}p\left(\zeta \right|y_{1:k})d\zeta ,$$

(7)

which we can obtain under suitable conditions as

$$F\left(y_{k+1}\right|y_{1:k})=E_{{\mathbf{x}}_{k+1}|y_{1:k}}\left[H_{{\mathbf{x}}_{k+1}}\left(y_{k+1}\right)\right].$$

(8)

In words, the cumulative predictive distribution is the expected GEV distribution conditioned on the state predictive density $p\left({\mathbf{x}}_{k+1}\right|y_{1:k})$ and our goal is to forecast ${\eta }_{\alpha }=F^{-1}\left(\alpha \right|y_{1:k})$ or the $\alpha$%-quantile.

To accomplish our goal, we need to derive a recursive solution for the filtering density, $p\left({\mathbf{x}}_k\right|y_{1:k})$, and the state predictive density, $p\left({\mathbf{x}}_{k+1}\right|y_{1:k})$. This can be done analytically for the case of known linear state-space models with Gaussian noise. However, with a non-linear state-space model, or non-Gaussian disturbances, a numerical approach is needed.

From standard Bayesian analysis, we can write

$$p\left({\mathbf{x}}_{0:k}\right|y_{1:k})\propto p\left(y_k\right|{\mathbf{x}}_k)p\left({\mathbf{x}}_{0:k}\right|y_{1:k-1}),$$

(9)

where ${\mathbf{x}}_{0:k}\equiv \{{\mathbf{x}}_0,{\mathbf{x}}_1,\dots ,{\mathbf{x}}_k\}$ is the state stream which are the state vectors from the initial state, at time 0, up to the current state, at time k. The likelihood function, $p\left(y_k\right|{\mathbf{x}}_k)$, highlights the assumption that the current observation depends only on the current value of the state vector. If we further assume that the state vector process is Markov, independent of the observations, we can write

$$p\left({\mathbf{x}}_{0:k}\right|y_{1:k})\propto \left\{p\left(y_k\right|{\mathbf{x}}_k)p\left({\mathbf{x}}_k\right|{\mathbf{x}}_{0:k-1})\right\}p\left({\mathbf{x}}_{0:k-1}\right|y_{1:k-1}),$$

(10)

which is the classic Bayesian recursion equation. If we can further invoke that the state process is first-order Markov, so that $p\left({\mathbf{x}}_k\right|{\mathbf{x}}_{0:k-1})=p\left({\mathbf{x}}_k\right|{\mathbf{x}}_{k-1})$, then we get the integral form of the Bayesian recursion [23],

$$p\left({\mathbf{x}}_k\right|y_{1:k})\propto p\left(y_k\right|{\mathbf{x}}_k){\int }_{{\mathbf{x}}_{k-1}}p\left({\mathbf{x}}_k\right|{\mathbf{x}}_{k-1})p\left({\mathbf{x}}_{k-1}\right|y_{1:k-1})d{\mathbf{x}}_{k-1},$$

(11)

which can be implemented with a standard particle filter [24].

In our case, where the static parameters in the state transition equation (4) are unknown, we can not directly implement Equations (10) or (11) since the current state depends on its complete history through the unknown parameters. The tact we take is to derive a Rao-Blackwellized particle filter via marginalizing $p({\mathbf{x}}_k,\mathbf{\eta },\Theta ,\mathbf{C}|{\mathbf{x}}_{0:k-1})$,

$${\int }_{\mathbf{\eta },\Theta ,\mathbf{C}}p\left({\mathbf{x}}_k\right|\mathbf{\eta },\Theta ,\mathbf{C},{\mathbf{x}}_{0:k-1})dP(\mathbf{\eta },\Theta ,\mathbf{C}|{\mathbf{x}}_{0:k-1}),$$

(12)

to obtain $p\left({\mathbf{x}}_k\right|{\mathcal{T}}_{k-1}$), where ${\mathcal{T}}_{k-1}$ are a set of sufficient statistics dependent only on the state stream up to time $k-1$. In doing so, we need an efficient, recursive solution to ${\mathcal{T}}_{k-1}$ so that the filter can be implemented without the need for repeated, large matrix operations or the need to retain the complete state history.

To start, we note that Equation (4) can be written compactly as a general linear model for the complete state stream up to time k as

$${\mathbf{X}}_k=\mathsf{\Phi }{\mathbf{H}}_k+{\mathbf{U}}_k,$$

(13)

where ${\mathbf{X}}_k=[{\mathbf{x}}_k,{\mathbf{x}}_{k-1},\dots ,{\mathbf{x}}_1]\in {\mathbb{R}}^{p\times k}$, with p being the number of state variables. The unknown matrix $\mathsf{\Phi }=[\mathbf{\eta }\Theta ]$$\in {\mathbb{R}}^{p\times p+1}$ and ${\mathbf{H}}_k=[{\mathbf{h}}_k,{\mathbf{h}}_{k-1},\dots ,{\mathbf{h}}_1]\in {\mathbb{R}}^{q\times k}$ has columns given by ${\mathbf{h}}_k^{\top }=\left[1{\mathbf{x}}_{k-1}^T\right]$. The noise term ${\mathbf{U}}_k\in {\mathbb{R}}^{p\times k}$ is a random Normal matrix whose columns are i.i.d. $\mathcal{N}(\mathbf{0},\mathbf{C})$ with $\mathbf{C}$ unknown. For the Rao-Blackwellized particle filter, the columns of ${\mathbf{X}}_k$ are formed by the particles from 1 to k and the columns of ${\mathbf{H}}_k$ include the past particles from 0 to $k-1$. Therefore, Ref. (13) is the state transition equation that describes the evolution of the state stream and, in the filter, each particle stream operates under its own version of this multivariate regression model.

The predictive density for the general linear model, with a non-informative, Jeffreys prior [25], results in a multivariate Student-t distribution [26], which for a p-dimensional random vector is denoted as $t_p(v,\mathbf{\mu },\mathbf{\Sigma })$$\propto {[1+{(\mathbf{x}-\mathbf{\mu })}^T{\mathbf{\Sigma }}^{-1}(\mathbf{x}-\mathbf{\mu })/v]}^{-(v+p)/2}$ with v degrees of freedom, mean $\mathbf{\mu }$, and scaling matrix $\mathbf{\Sigma }$. Using the non-informative prior, $p(\Theta ,\Theta ,{\mathbf{C}}^{-1})\propto {\left|\mathbf{C}\right|}^{(p+1)/2}$, we can write the marginal predictive density for ${\mathbf{x}}_{k+1}$ as a multivariate Student-t distribution with $v_k=k-2p$ degrees of freedom, i.e.,

$$p\left({\mathbf{x}}_{k+1}\right|{\mathbf{X}}_k,{\mathbf{H}}_k,{\mathbf{h}}_{k+1})~t_p(v_k,{\widehat{\mathbf{x}}}_{k+1},{\widehat{\mathbf{\Sigma }}}_{k+1}).$$

The mean of this distribution is derived as

$${\widehat{\mathbf{x}}}_{k+1}={\mathbf{S}}_{XH}^k{\left({\mathbf{S}}_{HH}^k\right)}^{-1}{\mathbf{h}}_{k+1}={\widehat{\mathsf{\Phi }}}_k{\mathbf{h}}_{k+1}$$

(14)

and the scale matrix is

$${\widehat{\mathbf{\Sigma }}}_{k+1}=\frac{1+{\mathbf{h}}_{k+1}^{\top }{\left({\mathbf{S}}_{HH}^k\right)}^{-1}{\mathbf{h}}_{k+1}}{v_k}({\mathbf{S}}_{XX}^k-{\widehat{\mathsf{\Phi }}}_k{{\mathbf{S}}_{XH}^k}^{\top }),$$

(15)

with the matrices defined as

$${\mathbf{S}}_{XX}^k={\mathbf{X}}_k{\mathbf{X}}_k^{\top },{\mathbf{S}}_{XH}^k={\mathbf{X}}_k{\mathbf{H}}_k^{\top },{\mathbf{S}}_{HH}^k={\mathbf{H}}_k{\mathbf{H}}_k^{\top }.$$

(16)

As is, Equations (14)–(16) can be used to compute the state predictive density and to extrapolate particles from time k to time $k+1$. That said, the resulting matrix operations are expensive computationally and a recursive solution is desired. Letting ${\mathbf{P}}_k={\left({\mathbf{S}}_{HH}^k\right)}^{-1}$, we derived recursive expressions for the mean and the scale matrix as follows:

$$\begin{array}{cc}\hfill {\mathbf{K}}_k& ={\mathbf{h}}_k^{\top }{\mathbf{P}}_{k-1}{[1+{\mathbf{h}}_k^{\top }{\mathbf{P}}_{k-1}{\mathbf{h}}_k]}^{-1},\hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill {\mathbf{P}}_k& ={\mathbf{P}}_{k-1}[\mathbf{I}-{\mathbf{h}}_k{\mathbf{K}}_k],\hfill \end{array}$$

(18)

$$\begin{array}{cc}\hfill {\widehat{\mathsf{\Phi }}}_k& ={\widehat{\mathsf{\Phi }}}_{k-1}+({\mathbf{x}}_k-{\widehat{\mathbf{x}}}_k){\mathbf{K}}_k,\hfill \\ \hfill {\widehat{\mathbf{\Sigma }}}_{k+1}& =\frac{v_k-1}{v_k}\left\{\frac{1+{\mathbf{h}}_{k+1}^{\top }{\mathbf{P}}_k{\mathbf{h}}_{k+1}}{1+{\mathbf{h}}_k^{\top }{\mathbf{P}}_{k-1}{\mathbf{h}}_k}\right\}\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill & \times \left\{{\widehat{\mathbf{\Sigma }}}_k+({\mathbf{x}}_k-{\widehat{\mathbf{x}}}_k){({\mathbf{x}}_k-{\widehat{\mathbf{x}}}_k)}^{\top }/(v_k-1)\right\},\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill {\widehat{\mathbf{x}}}_{k+1}& ={\widehat{\mathsf{\Phi }}}_k{\mathbf{h}}_{k+1}.\hfill \end{array}$$

(21)

These equations bear a resemblance to the recursive equations found in the Kalman filter [27] and are readily implemented with minimum computational and memory requirements. The specific algorithm we use is discussed in the next section and we note that the above procedure can likewise be used in the case of a linear observation equation with unknown observation matrix.

3. Results

To implement our solution, we first make a simplifying assumption that reduces our state vector to ${\mathbf{x}}_k={[log\left({\xi }_k\right),log\left({\beta }_k\right)]}^{\top }$. Recall that for the GEV distribution, with $\xi >0$, the support for the observation is $y\in [\mu -\beta /\xi ,\infty )$, where we have dropped the subscript k for simplicity. To eliminate the parameter $\mu$, we constrain $H_{\xi ,\beta ,\mu }\left(y_0\right)={\alpha }_0$ allowing us to solve for $\mu =-H_{\xi ,\beta ,y_0}^{-1}\left({\alpha }_0\right)$. This is a generalization of the constraint used in [8] where the support was $y\in [0,\infty )$ or ${\alpha }_0=y_0=0$ resulting in $\mu =\beta /\xi$. The main issue in the previous study was that small values for the parameter $\xi$ resulted in a substantial portion of the left tail with negligible probability since $H_{\xi ,\beta ,\beta /\xi }^{-1}\left(\alpha \right)\to \infty$ as $\xi \to 0^{+}\forall \alpha >0$ resulting in sub-par prediction performance. While the parameter estimates in [8] were quite reasonable, and outperformed the maximum likelihood (ML) method, quantile estimates were biased.

Under our support constraint, we can write

$$\mu =y_0-\frac{\beta }{\xi }\left[{(-ln\left({\alpha }_0\right))}^{-\xi }-1\right]$$

(22)

and the now two-parameter GEV distribution is

$$H_{\xi ,\beta }\left(y\right)=exp\left[-{\left(\frac{\xi }{\beta }(y-y_0)+{(-ln\left({\alpha }_0\right))}^{-\xi }\right)}^{-1/\xi }\right]$$

(23)

with support $y\in [y_0-\frac{\beta }{\xi }{(-ln\left({\alpha }_0\right))}^{-\xi },\infty )$. The likelihood function, $p\left(y\right|\mathbf{x})$, can then be written as

$$p\left(y\right|\mathbf{x})=\frac{d}{dy}{H}_{\xi ,\beta }\left(y\right)=\frac{1}{\beta }{H}_{\xi ,\beta }\left(y\right){\left(\frac{\xi }{\beta }(y-y_0)+{(-ln\left({\alpha }_0\right))}^{-\xi }\right)}^{-(1+{\xi }^{-1})},$$

(24)

where the substitutions $\xi =e^{x_1}$ and $\beta =e^{x_2}$ can readily be made.

Armed with the likelihood in Equation (24) and the sufficient statistics, computed recursively in Equations (17)–(21), we are ready to implement our version of the RBPF (see Algorithm 1).

Algorithm 1 RBPF recursive implementation.

Initialization ($k=0$, $m=1:M$ particles): Generate $\tilde{p}>2p+1$ random particles, ${\mathbf{x}}_{0:\tilde{p}}^{\left(m\right)}$ from (13) with random $\mathsf{\Phi },\mathbf{C}$. Compute ${\mathbf{S}}_{XX}^{\tilde{p}}$, ${\mathbf{S}}_{XH}^{\tilde{p}}$, and ${\mathbf{S}}_{HH}^{\tilde{p}}$ from (16). Initialize ${\widehat{\mathsf{\Phi }}}_0={\mathbf{S}}_{XH}^{\tilde{p}}{\left({\mathbf{S}}_{HH}^{\tilde{p}}\right)}^{-1}$, ${\widehat{\mathbf{x}}}_1$ from (14) and ${\widehat{\mathbf{\Sigma }}}_1$ from (15). Let ${\mathbf{P}}_0={\left({\mathbf{S}}_{HH}^{\tilde{p}}\right)}^{-1}$ and $v_k=\tilde{p}-2p$. Set particle weights $w^{\left(m\right)}=1/M$ and $M_1=M$. Recursive Loop ($k=1:N$): Sample ${\mathbf{x}}_k^{\left(m\right)-}~t_p(v_k,{\widehat{\mathbf{x}}}_k,{\widehat{\mathbf{\Sigma }}}_k)$ for $m=1,\dots M_k$. (particle extrapolation) Compute predictive statistics: ${\widehat{{\mathbf{x}}_k}}^{-}={\sum }_{m=1}^M{w}^{\left(m\right)}{\mathbf{x}}_k^{\left(m\right)-}$ Compute quantile estimates, p-values, and $M_k$ (see Algorithm 2) Systematically resample particles, ${\mathbf{x}}_k^{\left(m\right)-}\to {\mathbf{x}}_k^{\left(m\right)+}$, using $p\left(y_k\right|{\mathbf{x}}_k)$ in (24) Compute contemporaneous statistics: ${\widehat{{\mathbf{x}}_k}}^{+}={\sum }_{m=1}^M{w}^{\left(m\right)}{\mathbf{x}}_k^{\left(m\right)+}$ Compute ${\mathbf{K}}_k$ from (17). Update ${\mathbf{P}}_{k-1}\to {\mathbf{P}}_k$ via (18) and ${\widehat{\mathsf{\Phi }}}_{k-1}\to {\widehat{\mathsf{\Phi }}}_k$ via (19) Update to ${\widehat{\mathbf{\Sigma }}}_{k+1}$ from (20) and ${\widehat{\mathbf{x}}}_{k+1}$ from (21) $k=k+1$, $v_k=v_k+1$, and repeat loop

In Algorithm 1, estimates for ${\xi }_k$ and ${\beta }_k$ are produced by a weighted sum of the particles ${\mathbf{x}}_k^{\left(m\right)+}$ available after processing the observation, $y_k$, at time k. To illustrate the performance of the algorithm, we simulated the stochastic GEV parameters from an uncoupled, mean-reverting, AR(1) process as follows

$$\begin{array}{cc}\hfill log{\xi }_k& =(1-0.95)\times log\left(0.45\right)+0.95log{\xi }_{k-1}+0.1u_{1,k},\hfill \end{array}$$

(25)

$$\begin{array}{cc}\hfill log{\beta }_k& =(1-0.9)\times log\left(0.012\right)+0.9log{\beta }_{k-1}+0.1u_{2,k},\hfill \end{array}$$

(26)

where $u_{1,k}$ and $u_{2,k}$ are independent $\mathcal{N}(0,1)$. The process was chosen with $\overline{\xi }=0.45$ and $\overline{\beta }=1.25\%$ to be representative of time-series encountered in heavy-tailed phenomena such as financial markets and we refer the reader to [8] for representative simulation examples.

We produced 225 samples for each simulation run and we initialized the filter with $M_1=1024$ particles. At each point in time we computed estimates for the GEV parameters ${\widehat{\xi }}_k$ and ${\widehat{\beta }}_k$. We also computed the maximum likelihood (ML) estimates for the GEV parameters from (24) and used at least 25 observations ($k\ge 25$) before reporting estimation errors. For each of the 100 simulation runs, we computed a root-mean-square (RMS) error for the parameter estimates, across time, and shown in Figure 2 are the histograms for the RMS errors for both the RBPF and ML method. We can see the RBPF method outperforms with an average RMS error for $\widehat{\xi }$ of 0.11 vs. 0.17 for the ML method. Similarly, the average RMS error for $\widehat{\beta }$ was 0.32% vs. 0.41% for the ML method.

One of the improvements we made to our original algorithm was to use systematic resampling of the existing set of particles versus standard multinomial resampling which improves the filter’s performance. The other improvement was to adapt the number of particles based on an assessment of the filter’s predictive performance as in [28] although we modified their method so that new particles, if needed, are generated from the initial prior distribution. Finally, the quantile estimates, $\widehat{{\eta }_{\alpha }}={\widehat{F}}^{-1}\left(\alpha \right|y_{1:k})$ for $\alpha \in [0,1]$, based on averaging $H_{{\mathbf{x}}_{k+1}}(.)$ over the particles that approximate the predictive state distribution, were improved. In particular, we computed an estimate of the predictive distribution, via approximating Equation (8), versus using $E_{{\mathbf{x}}_{k+1}|y_{1:k}}\left[H_{{\mathbf{x}}_{k+1}}^{-1}\left(\alpha \right)\right]$. Algorithm 2 provides the details of the method used for quantile estimation, computing p-values, and to adapt the number of particles in the filter.

To test our approach, we examined the predictive performance of our method, outlined in Algorithm 2, versus the maximum likelihood. We applied a chi-square test as described in [29]. At each time k, we effectively partitioned the support of the predictive cdf estimate, $\widehat{F}\left(y_{k+1}\right|y_{1:k})$ into B equiprobable buckets. For the ML method, we did the same using the ML estimates of the GEV parameters as “truth”. We then compared the $(k+1)$st observation, $y_{k+1}$, to each bucket interval to find its location and to create an empirical frequency distribution. If the model is true, this frequency distribution will be a sample estimate of a uniform distribution and the statistic D is approximately chi-square with $B-1$ degrees of freedom (See Algorithm 2 for details).

For our test, we chose $B=20$ intervals and shown in Figure 3 are the frequency distributions for the entire 100 runs. We note that, in total, there were 20,000 sets of quantiles forecasted and each forecast was based solely on past data. More importantly, each prediction was based on a varied amount of past data. For example, in each simulation run, the first set of quantile forecasts were based on just 25 observations ($y_{1:25})$ and the last forecasts were based all that run’s observations except the last ($y_{1:224})$. We see that the RBPF method has a more uniform frequency distribution compared to the ML. The ML, similar to our previous work in [8], tended to under-estimate the 5% quantile, leading to more observations than expected. More importantly, the extreme quantiles were overestimated leading to less observations exceeding those thresholds. While one might argue that the bias to higher thresholds, for those extreme quantiles, is at least being “conservative”, it is clearly more beneficial to having a more accurate cdf and quantile forecast, particularly since these quantile forecasts are typically used to compute exceedance probabilities over multi-block periods.

Algorithm 2 Compute quantile estimates, p-values, and $M_k$.

Initialization (prior to Recursive Loop) Set ${\alpha }_j=j/B$ for $j=0,\dots ,B$ (B = # buckets) Create y-grid: $\tilde{y}=0:y_{step}:y_{max}$ and set ${[{\widehat{\xi }}_k^{-},{\widehat{\beta }}_k^{-}]}^T=exp\left({\widehat{{\mathbf{x}}_k}}^{-}\right)$. Initialize counter $C\left(j\right)=0$ for $j=1,\dots ,B$ Within Recursive Loop ( $k\ge N_{pred}$ ) Compute predictive CDF: $\widehat{F}\left(\tilde{y}\right|y_{1:k})={\sum }_{m=1}^M{w}^{\left(m\right)}{H}_{{\widehat{\xi }}_k^{-},{\widehat{\beta }}_k^{-}}\left(\tilde{y}\right)$. Invert predictive CDF: ${\eta }_{{\alpha }_j}={\widehat{F}}^{-1}\left({\alpha }_j\right)\forall j$ Find index $j^{*}$ such that $y_k\in [{\eta }_{{\alpha }_{j^{*}-1}},{\eta }_{{\alpha }_{j^{*}}})$ for $j^{*}\in \{1,\dots ,B\}$ Increment counter $C\left(j^{*}\right)$ If $k\ge 2*N_{pred}$ (compute p-values) − $N_c={\sum }_{j=1}^{B}C\left(j\right)$ − Compute ${\chi }^2$ statistic: $D=B(N_c-1){\sum }_{j=1}^B{(C\left(j\right)/N_c-1/B)}^2$ − Compute p-value using ${\chi }_{B-1}^2\left(D\right)$ test − If p-value < $p_{low}$ (double # of particles) * $M_{new}=M_k-min(2M_k,M_{max})$ * Initialize $M_{new}$ new particles (See Initialization in Algorithm 1) * $M_k=M_k+M_{new}$ − If p-value > $p_{high}$ (halve # of particles) * Randomly discard $M_k-max(M_k/2,M_{min})$ existing particles * $M_k=max(M_k/2,M_{min})$

To further illustrate the improvement in the performance of our approach, we show the histogram of the final p-values, for the 100 simulation runs, in Figure 4. The ML method has about 35 of the p-values in the range of [0,0.1), compared to an expected amount of 10, and most of the p-values are less than 0.5. This indicates that one would reject the ML method as an alternative hypothesis in most instances. In our previous work, our simulated p-values were worse than the ML method due to the reasons cited previously. Our new RBPF method, derived herein, has p-values close to the expected number for each bucket (i.e., 10), and produces quite reasonable quantile forecasts. To further see the accuracy of forecasting extreme quantiles with the RBPF, we plotted the empirical quantile estimate versus the predicted in Figure 5.

As seen in Figure 5, at the predicted 99.9%-quantile, 99.87% of the observations were below the threshold. Stated differently, 26 out of the 20,000 predictions exceeded the forecasted 99.9%-quantile threshold, which is six more than expected. We find this to be excellent given that many of the forecasts are made with limited data in hand. At the 99%-quantile, 0.15% additional observations exceeded our forecasts. Out of the 20,000 predictions, 200 exceedances were expected versus the 230 were observed. Our simulation results bear out the view that our methodology not only outperforms traditional quantile forecast methods (e.g., ML) but also our previous work in [8].

4. Discussion

From our results, illustrated in the previous section, we see that our methodology offers improvements in forecasting extreme quantiles, particularly when the time-series is non-stationary. This is true in many applications, where extreme behaviour tends to cluster in time, such as in finance, or trend over time, such as in climatology. To illustrate the application of our approach to real-world data, we applied it the data shown in Figure 1 which are the block-maximums of daily losses, or negative returns expressed in percent, for the S&P 500 stock market index, a widely used proxy for the overall stock market.

For completeness, we show in Figure 6 the likelihood surface for all of the data from 1928–2020 (93 years/block-maximums). While the ML estimates, using all of the data, are ${\widehat{\xi }}_{ML}=0.48$ and ${\widehat{\beta }}_{ML}=1.26\%$, one can see there is a large degree of uncertainty in the GEV parameters. For example, there is considerable probability that $\xi >0.5$, which would imply an underlying time-series with infinite variance and this has implications for financial models that rely on finite second moments (e.g., Black-Scholes). The Bayesian approach is to use all of the information available from the data to produce posterior densities from which to produce parameter estimates and, more importantly, to forecast. This is particularly important when forecasting quantiles since they can be quite sensitive to point estimates given the highly non-linear form of the GEV cdf.

Shown in Figure 7 are the S&P 500 block-maximums along with quantile forecasts over time. For the 90% quantile forecast, 8 out of 92 observed stock-market returns exceeded the forecast, or 8.7% of the time. This is well within what one would expect. The current, forecasted 90% quantile is 7.1% and the 99% quantile forecast is 23.8%. We should expect to see this type of market "crash," or worse, about once a century. While the range of the forecasted 90% quantile, from 5.6% in the 1970s to 8% post financial crisis, may seem minimal, we do note that more extreme quantile estimates have significant variation. For example, the 99% quantile forecast went from 18% to 25% and the 99.5% quantile went from 25.8% to 35% over the same period. This would be of use to a financial institution when constructing risk scenarios for stress testing their portfolios. Last year, in 2020, the block-maximum was a 12% loss, during the covid pandemic crisis, which was at the 96% quantile. We should expect to see at least this 12% loss about once every 25 years and it should not be considered unusual highlighting the importance of quantile forecasting-bad things do happen!

It is worth pointing out that the p-value based on the forecasted quantiles was a robust 0.75 versus 0.01 for the ML method. While one must be wary in using p-values to accept a model, we feel that the results we have shown in this paper clearly indicate an improved method of forecasting. But much work needs to be done. For one, we simply used block-maximums in this paper, which one can argue wastes data. We would like to extend our approach further to use more higher-order statistics, or exceedances above a threshold. We would also like to broaden our simulation study and apply our results to other data sets, particularly in the area of climatology. Lastly, additional computational efficiencies should be explored to allow real-time implementation to high-frequency data.