Sequential monitoring for explosive volatility regimes

Sequential monitoring for explosive volatility regimes

\address{Lajos Horv\'ath, Department of Mathematics, University of Utah, US} \address{Lorenzo Trapani, University of Leicester, UK, and Universita' di Pavia, Italy} \address{Shixuan Wang, Department of Economics, University of Reading, UK} \subjclass{Primary 62M10, 91B84; secondary 60G10, 62F12}

\doublespacing

Introduction

In recent years, developing tools for the (ex-ante or ex-post) detection of the onset or the collapse of a bubble in financial markets has been one of the most active research areas in financial econometrics; we refer the reader, in particular, to the seminal articles on ex-post detection by phillips2011, and phillips2015testing2, and also to skrobotov2023testing for a review. As far as ex-ante detection - that is, finding the onset or collapse of a bubble in real time, as new data come in - is concerned, this important issue has also been studied in numerous recent contributions; although a comprehensive literature review goes beyond the scope of this paper, we refer to the articles by homm2012testing, phillips2015testing, and whitehouse2023real, inter alia; in particular, the paper by whitehouse2023real also contains a comprehensive literature review of in-sample and online bubble detection methods. A common trait to the vast majority of the existing literature is its reliance on a linear specification, usually an AutoRegressive (AR) model, to capture regime changes in the dynamics of log prices. Whilst such a modelling choice can be justified from the theoretical point of view (see e.g. phillips2011dating), and whilst its analytical tractability offers obvious advantages from a mathematical standpoint (see e.g. phillips2007limit, and aue2007limit), a major issue is that using an AR model is fraught with difficulties when monitoring for changes from an explosive towards a stationary regime. Several promising solutions have been proposed such as the reverse regression approach by phillips2018financial. horvath2023real propose a different model, based on a Random Coefficient Autoregressive (RCA) specification, where inference is always standard normal irrespective of stationarity or the lack thereof (aue2011), and develop a family of sequential monitoring procedures based on the weighted CUSUM process, to check whether the deterministic part of the autoregressive root changes over time. Such a testing set-up also encompasses both the case of a switch from a stationary to an explosive regime (thus indicating the start of a bubble phenomenon), and a change from an explosive to a stationary regime (thus signalling the collapse of a bubble).

The theory developed in horvath2023real paves the way to a more general research question, namely developing sequential monitoring techniques which are robust to both the initial regime (i.e., which can be employed irrespective of whether the observations in the training/historical sample are stationary or explosive), and to the type of change which occurs after a changepoint (i.e., which can detect changes from stationarity to another stationary regime, or to a explosive regime; or from an explosive regime to another explosive regime or a stationary one). From a technical viewpoint, this question is nontrivial for at least three\ reasons. First, proposing a changepoint detection methodology whose asymptotics is the same irrespective of stationarity or explosivity is not easy per se, because the partial sum processes which constitute the building blocks of e.g. CUSUM-based statistics require completely different approximations depending on whether the observations are stationary or not. Second, in order to ensure timely changepoint detection, weighted versions of the CUSUM process need to be considered, with different sets of weights ensuring optimal detection delays for different changepoint locations within the monitoring horizon. Third, on account of the previous point, it is important to offer, to the applied user, a battery of results on the limiting distribution of the detection delays, so as to gauge the expected detection delay depending on the location of the changepoint and the weighing scheme employed.

Motivated by the questions above, in this paper we investigate, with an emphasis on completeness, the issue of sequential detection for regime changes in a GARCH(1,1) sequence

In particular, we develop two families of detectors based on the weighted CUSUM process of the quasi-Fisher scores associated with ((ref)): one with \textquotedblleft mild\textquotedblright\ weights, designed to detect changes that may occur \textquotedblleft not too early\textquotedblright\ after the start of the monitoring period; and one with heavy weights, designed instead to detect changes occurring \textquotedblleft very early\textquotedblright\ after the start of the monitoring period. The latter is based on applying to the CUSUM process a set of (heavy) weights, resulting in a family of test statistics known as R\'{e}nyi statistics (see horvath2021 and horvathmiller, for in-sample tests, and ghezzi2024fast, for sequential monitoring). Our methodologies can, in principle, be applied to detect any type of change in the vector $\left( \omega ,\alpha ,\beta \right) $. However, seeing as: (1) our main interest is in detecting changes between regimes (e.g., from stationarity to nonstationarity, or vice versa); the stationarity or lack thereof of ((ref)) is determined solely by $\alpha $ and $\beta $; and (3) in the presence of nonstationarity, $\omega $ is not identified ( jensen2004asymptotic), we focus on monitoring for changes only in the sub-vector $\left( \alpha ,\beta \right) $. Detecting shifts in the behaviour of the (conditional) volatility process $\sigma _{i}$ is important in general; as hillebrand2005neglecting notes, when neglecting a break inference is biased in finite samples, and the sum of the estimated autoregressive parameters $\alpha $ and $\beta $ converges to one. Furthermore, changes (and, occasionally, explosions) in the volatility of time series are often observed in practice (see e.g. bloom2007uncertainty, and jurado2015measuring). Hence, finding the start or the end of an explosive regime in $\sigma _{i}$ is of practical relevance because, as richter2023testing put it, one \textquotedblleft often sees sudden, integrated, or mildly explosive behaviour in the second moment of the process which bounces back after a while\textquotedblright\ (p. 468). Changes between stationarity and explosivity in ((ref)) can be interpreted as volatility bubbles , i.e. events in which the second moment of the data (rather than the data, e.g. prices, themselves) experiences periods of exhuberance. The link between a volatility phenomenon and a \textquotedblleft proper\textquotedblright\ bubble has not been fully explored yet (see jurado2015measuring), and, empirically, explosive regimes in volatility can be ascribed to various sources in addition to bubbles (see sornette2018can, and richter2023testing). Nevertheless, as jarrow2023explosion put it, \textquotedblleft price bubbles result from excess speculative trading decoupled from the asset's fundamentals (dividends and liquidation value), which increases the asset's price volatility to extreme levels\textquotedblright\ (p. 478). Hence, an analysis of the regimes of the volatility of financial variables is bound to contribute to a better understanding of bubble phenomena.

Based on the discussion above, in this paper we propose a battery of tests for the sequential monitoring of the volatility of financial variables, which complements the existing tests for bubbles based on the conditional mean. Specifically, we make at least three contributions. First, we study online detection of changes between stationarity and explosivity and vice versa, which - to the best of our knowledge - is a novel result in the literature, and which complements the ex-post detection statistics studied in richter2023testing. Second, we develop the full-blown theory for R \'{e}nyi statistics in the context of sequential monitoring of a GARCH(1,1) model. Third, we derive the limiting distribution of detection delays for all our monitoring schemes, including those based on R\'{e}nyi statistics; this is an entirely novel result, which complements the results in horvath2020sequential.\newline

The remainder of the paper is organised as follows. We discuss our workhorse model and the main assumptions, as well as the test statistics, in Section (ref). The theory is reported in Section (ref): in particular, the asymptotics under the null is in Section (ref), and the full-blown asymptotics of the detection delay in the presence of a changepoint is in Section (ref). We validate our theory through a comprehensive set of simulations (Section (ref)), and an empirical application to daily returns of individual stocks (Section (ref)). Section (ref) concludes.

NOTATION. We define the Euclidean norm of a vector $a$ as $\Vert a\Vert $. We denote the integer part of a number as $\left\lfloor \cdot \right\rfloor $ . We use: \textquotedblleft a.s\textquotedblright\ for \textquotedblleft almost sure(ly)\textquotedblright ; \textquotedblleft $\rightarrow $ \textquotedblright\ for the ordinary limit; \textquotedblleft $\overset{ \mathcal{D}}{\rightarrow }$\textquotedblright\ for convergence in distribution; \textquotedblleft $\overset{\mathcal{P}}{\rightarrow }$ \textquotedblright\ for convergence in probability; \textquotedblleft $ \overset{{\mathcal{D}}}{=}$\textquotedblright\ for equality in distribution. Positive, finite constants are denoted as $c_{0}$, $c_{1}$, ... and their value may change from line to line. Other, relevant notation is introduced later on in the paper.

Model, assumptions, and hypothesis testing

Model, assumptions and hypotheses of interest

The time dependent GARCH (1,1) sequence is defined by the recursion

where $\sigma _{0}^{2}$, $y_{0}^{2}$ are initial values, and $\alpha _{i}$, $ \beta _{i}$ and $\omega _{i}$ are positive parameters.

Whilst the hypothesis testing framework is spelt out below, our monitoring schemes are all based on the maintained assumption that we have $m$ observations which form a stable period (this is also known as the non-contamination assumption in CSW96), viz.

We denote the value of the common parameter in (ref) as $ \mbox{\boldmath${\theta}$}_{0}=(\alpha _{0},\beta _{0},\omega _{0})^{\top }$ . Prior to spelling out the main assumptions, we review the conditions for the stationarity of $y_{i}$. As nelson1991conditional shows (see also bougerol1992strict, and francq2012strict)

We will develop several monitoring schemes for the null hypothesis that the parameter $\mbox{\boldmath${\theta}$}_{0}$ undergoes no changes after the historical training period $1\leq i\leq m$, i.e.

Under the alternative, we assume that there is a change at time $m+k^{\ast }$ ; whilst this would correspond to having $\left( \omega _{m+k^{\ast }-j},\alpha _{m+k^{\ast }-j},\beta _{m+k^{\ast }-j}\right) $ $\neq $ $\left( \omega _{m+k^{\ast }+j},\alpha _{m+k^{\ast }+j},\beta _{m+k^{\ast }+j}\right) $ for all $j\geq 0$, it is well know that the $\omega _{i}$'s cannot be identified in explosive, nonstationary regimes ( francq2012strict). Hence, we will test for

i.e., for the possible presence of changes in $\alpha _{i}$ and $\beta _{i}$ only. Note that these are anyway the parameters of interest, since the stationarity (or lack thereof) of $y_{i}$ is not affected by $\omega _{i}$.

We require the following assumptions on $\mbox{\boldmath${\theta}$}_{0}$, and on the innovations $\{\epsilon _{i},-\infty <i<\infty \}$.

Assumptions (ref) and (ref) are standard. In particular, it is worth noting that, in Assumptions (ref), there is no requirement as to the stationarity properties of $\left\{ y_{i},1\leq i\leq m\right\} $: the observations in the training sample can belong to a stationary or explosive volatility regime.

Estimation and monitoring schemes

As stated in the Introduction, the main purpose of our analysis is to offer a detection scheme which finds changes in the parameters of a GARCH(1,1) model as soon as possible after the training period. In this section, we propose several detectors, all based on the CUSUM process of the quasi-Fisher scores.

As is typical, we estimate the parameter $\mbox{\boldmath${\theta}$}_{0}$, using the training sample, by Quasi Maximum Likelihood (QML). The $\log $ likelihood function is given by

where $\mbox{\boldmath${\theta}$}=(\omega ,\alpha ,\beta )^{\top }$, $\omega >0,\alpha >0$ and $\beta >0$, and the random functions $\bar{\sigma}_{i}^{2}( \mbox{\boldmath${\theta}$})$ given by the recursion

the recursion starts from the initial values $y_{0}$ and $\bar{\sigma} _{0}^{2}$. The QML estimator computed from the historical sample is denoted as $\hat{\mbox{\boldmath${\theta}$}}_{m}$, with

where the (compact) space $\mbox{\boldmath${\theta}$}$ is defined as

$0<\underline{\omega },\overline{\omega },\underline{\alpha },\overline{ \alpha },\underline{\beta },\overline{\beta }<\infty $.

Starting with the initial values ${y}_{m}$ and ${\sigma }_{m}^{2}$, we define the random functions $\bar{\sigma}_{m+k}^{2}(\mbox{\boldmath${ \theta}$})$ based on the observations after the historical sample by the recursion

with the log likelihood function given by

Hence, we define the CUSUM\ process of the quasi-Fisher scores as

Heuristically, under the null of no change, the scores have zero mean; hence, the partial sum process ${\mathcal{r}}_{m,k}(\mbox{\boldmath${ \theta}$})$ calculated at $\hat{\mbox{\boldmath${\theta}$}}_{m}$\ should also fluctuate around zero with increasing variance. Conversely, in the presence of a break (at, say, $k^{\ast }$), $\hat{\mbox{\boldmath${\theta}$}} _{m}$ is a biased estimator for the \textquotedblleft new\textquotedblright\ parameter $\mbox{\boldmath${\theta}$}_{m+k^{\ast }}$; thus, ${\mathcal{r}} _{m,k}(\mbox{\boldmath${\theta}$})$, calculated at $\hat{\mbox{\boldmath${ \theta}$}}_{m}$, should have a drift term. In the light of these heuristic considerations, we propose the following detector

where

Based on ((ref)), a break is flagged as soon as the detector ${ \mathcal{D}}_{m}(k)$ exceeds a threshold. We call such a threshold the boundary function. Similarly to CSW96, lajos04, horvath2020sequential, homm2012testing, we use the boundary function, designed for a closed-ended procedure - i.e., for a procedure which terminates at a certain time, say ${\mathcal{n}}$, if there is no change

On account of ((ref)) and ((ref)), a changepoint is found at a stopping time $\tau _{m}$ defined as

The (user-chosen) parameter $\eta $ in ((ref)) determines the timeliness of changepoint detection of our sequential monitoring procedure. aue2004delay and aue2008monitoring show that, as $\eta $ approaches $1$, changepoints are detected with a smaller and smaller delay depending on their location; on the other hand, different values of $\eta $ work better for different changepoint locations, as also pointed out in a recent contribution by kirch2022sequential. In particular, values of $0\leq \eta <1$ are able to offer short detection delays for breaks that do not occur \textquotedblleft too early\textquotedblright\ after $m$.

In order to detect earlier changes, ghezzi2024fast (building on previous contributions by horvath2021, and horvathmiller , developed for in-sample changepoint detection) suggest using R\'{e}nyi type statistics, with stopping rule

where $r$ is a trimming sequence specified in Assumption (ref), and

We note that ((ref)) and ((ref)) exclude the case $\eta =1$. Indeed, aue2004delay show that, as far as stopping times under the alternative are concerned, using $\eta =1$ would produce the shortest detection time. The case $\eta =1$ requires to be studied separately. Following lajos07, we modify the boundary function. Let

We use the boundary functions

and

The corresponding stopping times - denoted as $\tau _{m}^{\ast }$ and $\bar{ \tau}_{m}^{\ast }$ - are defined exactly in the same way as $\tau _{m}$ and $ \bar{\tau}_{m}$ in ((ref)) and ((ref)), respectively, using now the boundaries $g_{m}^{\ast }(k)$ and $\bar{g}_{m}^{\ast }(k)$.\newline

As mentioned above, we consider a closed-ended scheme, which terminates ${ \mathcal{n}}$ periods after $m$. The following assumptions characterise the length of the monitoring horizon and of the trimming sequence $r$ defined in ((ref)); in particular, Assumption (ref) is designed in order to consider only early changepoint detection.

Asymptotics

In this section, we investigate the asymptotic behaviour of our monitoring schemes under the null and under the alternative hypotheses.

Asymptotics under the null

Let $\mathbf{W}(t)=(W_{1}(t),W_{2}(t)),t\geq 0$ be a two dimensional standard Wiener process - i.e.,\ $\{W_{1}(t),t\geq 0\}$ and $ \{W_{2}(t),t\geq 0\}$ are two independent Gaussian processes with $ EW_{1}(t)=EW_{2}(t)=0$, and covariance kernel $EW_{1}(t)W_{1}(s)$ $=$ $ EW_{2}(t)W_{2}(s)$ $=\min (t,s)$.

Theorem (ref) offers a rule to calculate the asymptotic critical values; for a given nominal level $\alpha $, the critical value ${\mathcal{c}} _{\alpha }$ is defined as

for all $0\leq \eta <1$, and

for $\eta >1$. Using the scale transformation of the Wiener process, it immediately follows that $\{\mathbf{W}(t),t>0\}$ $\overset{{\mathcal{D}}}{=}$ $\{t\mathbf{W}(1/t),t>0\}$; hence, for all $\eta >1$

Theorem (ref) rules out the boundary case $E\log (\alpha _{0}\epsilon _{0}^{2}+\beta _{0})=0$; this is because we would need an exact (and large enough) rate of divergence for $\left\vert \sigma _{i}\right\vert $ as $ i\rightarrow \infty $, but this result is not available in the case $E\log (\alpha _{0}\epsilon _{0}^{2}+\beta _{0})=0$ (see francq2012strict , p. 823; and also Theorem 4 in HT2019, where a similar problem is encountered, and the discussion thereafter).\newline

Theorem (ref) does not consider the case $\eta =1$, which corresponds to the stopping times $\tau _{m}^{\ast }$ and $\bar{\tau}_{m}^{\ast }$ based on the boundaries defined in ((ref)) and ((ref)) respectively. The case $\eta =1$ is studied separately in the following theorem.

Theorem (ref) stipulates that the asymptotic critical values, for a given nominal level $\alpha $, can be calculated as

using ${\mathcal{c}}_{\alpha }$\ or $\overline{{\mathcal{c}}}_{\alpha }$\ according as ((ref)) or ((ref)) is employed. Although the theorem offers an explicit formula to compute asymptotic critical values, these are bound to be inaccurate due to the slow convergence to the Extreme Value distribution. In particular, simulations show that, in finite samples, asymptotic critical values overstate the true values thus leading to low power.

Asymptotics under the alternative

We now study the behaviour of our monitoring schemes under the alternative, focussing, in particular, on the limiting distribution of the detection delay. We report the limiting distribution of the detection delay when using $0\leq \eta <1$ in Section (ref); in Section (ref), we report the limiting distribution of the detection delay when using $\eta >1$.

In both cases, under the alternative $H_{A}$, of ((ref)), the parameter $ \mbox{\boldmath${\theta}$}_{0}=(\alpha _{0},\beta _{0},\omega _{0})^{\top }$ changes to $\mbox{\boldmath${\theta}$}_{A}=(\alpha _{A},\beta _{A},\omega _{A})^{\top }$ satisfying

Detection delays with $0\leq \protect\eta <1$

We begin by investigating the asymptotic behaviour of the stopping time $ \tau _{m}$ defined in ((ref)). Whilst the result in Theorem (ref) below is valid for all cases, we need to introduce some preliminary notation, separately, for the two cases: (1) when the sequence is stationary after the change and (2) when the sequence is explosive after the change.\newline

\qquad Preliminary notation\newline

We begin by introducing some preliminary notation for the former case, i.e.

Under ((ref)), after the change the observations are exponentially close to $\{\hat{x}_{i},-\infty <i<\infty \}$, a stationary GARCH (1,1) sequence given by

We also define the log likelihood function

where $\hat{h}_{i}^{2}(\mbox{\boldmath${\theta}$})=\omega +\alpha \hat{x} _{i-1}^{2}+\beta \hat{h}_{i-1}^{2}(\mbox{\boldmath${\theta}$})$. Let

and define the size of the change as

We define the covariance matrix

and

where $A_{m}=\mbox{\boldmath${\Delta}$}^{\top }\hat{\mathbf{D}}_{m}^{-1} \mbox{\boldmath${\Delta}$}$. Finally (as far as preliminary notation is concerned), we define $u_{\mathcal{n}}>0$ as the unique solution of the equation

and $u^{\ast }>0$ as the solution of

It is easy to see that $u_{\mathcal{n}}\rightarrow u^{\ast }$, and that $ u^{\ast }=1$, if $t^{\ast }=0$. We are now ready to introduce the main notation to spell out the properties of the stopping time $\tau _{m}$ when the observations change into a stationary sequence.\newline

We now introduce the preliminary notation for the case when the observations turn into an explosive sequence after the time of change, i.e.

jensen2004asymptotic proved that

exists. Similarly to $\mbox{\boldmath${\Delta}$}$ in equation ((ref) ), we define the size of the change under $H_{A}$ as $\mbox{\boldmath${ \Upsilon}$}={\mathcal{r}}^{(2)}(\mbox{\boldmath${\theta}$}_{0})\neq \mathbf{0 }$. Similarly to ((ref)), we define

where $B_{m}=\mbox{\boldmath${\Upsilon}$}^{\top }\hat{\mathbf{D}}_{m}^{-1} \mbox{\boldmath${\Upsilon}$}$. Finally, similarly to $u_{\mathcal{n}}$ and $ u^{\ast }$ we define $\tilde{u}_{\mathcal{n}}$ and $\tilde{u}^{\ast }$ as the solutions of the equations

After the change in the parameters, the gradient of the likelihood function is approximated with the sequences

and

Analogously to $\mbox{\boldmath${\Sigma}$}_{1}$ in ((ref)), we finally introduce

\qquad Main notation\newline

After spelling out the preliminary notation for the two cases of the observations begin stationary and nonstationary, we now introduce the main notation. As we will see in Theorem (ref) below, in several cases the delay $\tau _{m}-k^{\ast }$ converges to a standard normal random variable after being centered and rescaled; the centering and rescaling for $\tau _{m}-k^{\ast }$ depend on whether $t^{\ast }<\infty $ or $t^{\ast }=\infty $ ($\tilde{t}^{\ast }<\infty $ or $\tilde{t}^{\ast }=\infty $, equivalently). In particular

In order to present the limiting distribution for both cases, we define: the Gaussian process $\{\mbox{\boldmath${\Gamma}$}(t),t\geq 0\}$, with $E \mbox{\boldmath${\Gamma}$}(t)=\mathbf{0}$ and $E\mbox{\boldmath${\Gamma}$}(t) \mbox{\boldmath${\Gamma}$}^{\top }(s)=\min (t,s)\mathbf{D}$; the random variables

and, finally, the asymptotic variances

\newline

Let ${\mathcal{N}}$ denote a standard normal random variable, and let

In order to understand the practical implications of Theorem (ref), note that (up to some positive and finite constant)

The case ((ref)) corresponds to a \textquotedblleft very early\textquotedblright\ break; in this case, Theorem (ref) states that the expected delay is approximately $v_{1,{\mathcal{n}}}$, i.e. that it is approximately equal to $\left( m^{1-\eta }\right) ^{1/(2-\eta )}$. Clearly, as $\eta $ increases, $v_{1,{\mathcal{n}}}$ decreases; the dispersion around the expected delay, measured by $v_{2,{\mathcal{n}}}$, also decreases, indicating that the choice of $\eta $ plays a role in determining the delay in detecting (very early) changepoints, and that larger values of $\eta $ reduce such a delay. In the presence of an \textquotedblleft early, but not so early\textquotedblright\ break - corresponding to case (ii) of the theorem, where recall that $ k^{\ast }=o\left( {\mathcal{n}}\right) $, the expected delay $v_{3,{\mathcal{ n}}}$ still decreases as $\eta $ increases, as long as $k^{\ast }={\mathcal{n }}^{\gamma }$, for $\gamma >1/(2-\eta )$, but the dispersion around the expected delay - given by the standardization $v_{4,{\mathcal{n}}}$ - does not depend on $\eta $. Finally, the case of a late(r) change is studied in part (iii) of the theorem: in such a case, $\eta $ - and therefore the weight function in the definition of the detector - does not play any role.

Detection delays when $\protect\eta >1$

We now investigate the asymptotic behaviour of the stopping time $\tilde{\tau }_{m}$ defined in ((ref)) - that is, when the detector is a R\'{e} nyi type statistic with $\eta >1$. In such a case, the asymptotic behaviour of the detection delay uses the same notation irrrespective of whether $ E\log (\alpha _{A}\epsilon _{0}^{2}+\beta _{A})<0$ or $>0$.

Let

We define two independent normal random vectors $\mathbf{N}_{1}$ and $ \mathbf{N}_{2}$ such that $E\mathbf{N}_{1}=\mathbf{0},E\mathbf{N}_{2}= \mathbf{0},E\mathbf{N}_{1}\mathbf{N}_{1}^{\top }=\mathbf{D}$ and

Similarly to ${\mathcal{A}}_{1}$ and ${\mathcal{A}}_{2}(x)$ in ((ref)) and ((ref)), we define

and the centering sequence

and the rescaling sequence $v_{6,{\mathcal{n}}}=(k^{\ast })^{1/2}$.

Similarly to Theorem (ref), Theorem (ref) describes the detection delay when using R\'{e}nyi type statistics depending on the location of the break; to the best of our knowledge, this is the first time such a result has ever been derived. Part (i) of the theorem is also derived in ghezzi2024fast, and, in essence, it states that if the break occurs prior to the trimming sequence $r$ in the R\'{e}nyi type statistics, then it is identified straight at $r$ - that is, as soon as the R \'{e}nyi type statistics starts the monitoring. Parts (ii) and (iii) of the theorem refine and extend the results in ghezzi2024fast. Finally, part (iv) states that, in the case of a break occurring late - or, better, much later than $r$ - a large value of $\eta $ could even be detrimental because the centering sequence $v_{5,{ \mathcal{n}}}$ diverges with $k^{\ast }$, at a faster rate as $\tau $ increases. This confirms the common wisdom (see kirch2022asymptotic and kirch2022sequential), and the findings in ghezzi2024fast, that R\'{e}nyi type statistics are designed for the fast detection of very early occurring breaks, whereas they may yield suboptimal results for later breaks.

Simulations

In this section, we assess the finite sample performance of our monitoring procedures via Monte Carlo simulations. According to the theory in Section (ref), we can have two classes of monitoring schemes, based on $ \eta \neq 1$ (covered by Theorem (ref)) and $\eta=1$ (covered by Theorem (ref)). For the sake of brevity, here we only focus on the case $\eta \neq 1$. We consider several data generating processes (DGP). We use three lengths of the historical training sample $m=500,1000,5000$ and two lengths of the monitoring $\mathcal{n}=250,500$. The sequential procedure is performed $5,000$ times with independently generated samples, and the percentage of simulations for which the detector crosses the boundary functions is reported for several values of $\eta $. For the R\'{e}nyi type statistic based on Theorem (ref)(ii), we follow horvath2021 and set $r=\sqrt{\mathcal{n}}$. Guidelines on implementation are provided in Section (ref) of the Supplement. To obtain critical values, we simulate two independent standard Wiener processes $W_{1}(t)$ and $W_{2}(t)$ on a grid of $100,000$ equally spaced points in the unit interval $\left[ 0,1\right] $ and compute $\sup_{0<t\leq 1}\frac{1}{t^{\eta }}\left( W_{1}^{2}(t)+W_{2}^{2}(t)\right) $ and $ \sup_{0<t\leq 1}\frac{1}{t^{1-\eta }}\left( W_{1}^{2}(t)+W_{2}^{2}(t)\right) $. We repeat this by $100,000$ times and obtain the empirical $90\%$, $95\%$ , and $99\%$ percentiles of the above two quantities, corresponding to the critical values at $10\%$, $5\%$, and $1\%$ levels based on Theorem (ref) (i) and (ii). Critical values are in Table (ref).

The boundary functions in Section (ref) are designed for the case $ m\rightarrow \infty $. However, preliminary simulations show that the empirical sizes based on those boundary functions tend to be larger than the nominal levels in finite samples, in particular for DGPs with the Student's $ t$ errors. To make our monitoring schemes more practical under small finite samples, we suggest to \textquotedblleft tune\textquotedblright\ the boundary functions as

for ((ref)) and ((ref)) respectively. The intuition underpinning the term $\left( 1+1/\log (m)\right) ^{2}$ is to boost the boundary function in small samples. The term $\left( 1+k/m\right) ^{2}$ as is typically employed when the monitoring horizon is \textquotedblleft long\textquotedblright\ (see horvath2020sequential, horvath2021monitoring, and zenhya). Although this term is inconsequential for the asymptotic theory in our set-up, we find that it can further improve the empirical size. Both tuning terms are asymptotically negligible. and only play a role in finite samples to achieve better size control at no expense for power. The proposed tuning is tailored to DGPs with Student's $t$ errors, rather than Gaussian errors; indeed, heavy tails are a well-known stylised fact of financial returns.

Empirical size under the null

Under the null hypothesis, the realisation of GARCH(1,1) in the historical training period ($1\leqslant i\leqslant m$) and in the monitoring period ($ m+1\leqslant i\leqslant m+\mathcal{n}$) is

where $\epsilon _{i}$ follows a standard normal distribution or the Student's $t$ distribution with $7$ degrees of freedom. Since our monitoring procedure does not require the historical sample to be stationary or not, we choose the following two set of GARCH(1,1) parameters, taken from francq2012strict: (i) $(\omega _{0},\alpha _{0},\beta _{0})=(0.10,0.18,0.80)$, which represents the stationary case; (ii) $(\omega _{0},\alpha _{0},\beta _{0})=(0.10,0.30,0.80)$, corresponding to the nonstationary case since $E\log (\alpha \epsilon _{0}^{2}+\beta )>0$ under errors following either the standard normal or the Student's $t$ distributions. \newline

Table (ref) reports the empirical sizes at $5\%$ significance level for the monitoring scheme based on Theorem (ref)( i) for different values of $\eta $. A noticeable feature is that a larger $\eta $ results in a higher rejection rates, and a smaller $\eta $ is more conservative in rejection. Under the Student's $t$ errors, $\eta =0.3$ is a good choice because the monitoring procedure has reasonably good empirical sizes when $m=1,000$, and the empirical sizes for $m=5,000$ are closer to the theoretical level of $5\%$. Under Gaussian errors, the monitoring procedure is slightly under-sized, which is mainly due to the additionally tuning we imposed in (ref). For practical use, the tendency to under-reject with Gaussian errors may not necessarily be a concern, because the empirical power does not seem to be affected, as shown in Section (ref). Lastly, the simulation results show that our monitoring schemed works reasonably well for both stationary and nonstationary GARCH(1,1) models.\newline

Table (ref) in the Supplement contains the empirical sizes for R\'{e}nyi type statistics. The rejection rates are slightly higher than the $5\%$ nominal level. We note that, in principle, it would be possible to design a different tuning for R\'{e}nyi type statistics.

Empirical power under $H_{A}$

We now turn to the analysis of the empirical power. Under the alternative, the data is generated by

and

where the parameter $\boldmath{\theta }_{0}$$=(\alpha _{0},\beta _{0},\omega _{0})^{\top }$ changes to $\boldmath{\theta }_{A}$$=(\alpha _{A},\beta _{A},\omega _{A})^{\top }$ at time $m+k^{\ast }$. We consider two scenarios for the time of change: (a) $k^{\ast }=\lfloor \sqrt{\mathcal{n}} \rfloor $ corresponds to a change occurring \textquotedblleft early, but not too early\textquotedblright\ after the historical sample; (b) $ k^{\ast }=0.5{\mathcal{n}}$ indicates a change happening much later than $r$ . \newline

There are many possible ways of changes under the alternative. To keep our results clean, we set $\omega _{0}=\omega _{A}=0.1$ and $\alpha _{0}=\alpha _{A}=0.18$, and concentrate on a change in $\beta $ under the following four representative alternatives:

Tables (ref) and (ref) show the empirical power of the monitoring scheme based on Theorem (ref)( i) at $5\%$ significance level when $\mathcal{n}=500$ for a change at $k^{\ast }=\lfloor \sqrt{\mathcal{n}}\rfloor $ and $k^{\ast }=\lfloor 0.5 \mathcal{n}\rfloor $, respectively.\footnote{ The empirical power of $\mathcal{n}=250$ (not reported) is marginally lower than the empirical power of $\mathcal{n}=500$.} There are five major observations. First, our monitoring scheme is highly effective in detecting changes under $H_{A,2}$ and $H_{A,3}$ for both early and late changes. These alternatives result in a change between a stationary regime and an explosive regime, which is relatively easy to detect. Second, the monitoring scheme exhibits high power in detecting early changes under $H_{A,1}$ and $H_{A,4}$ . These alternatives represent a change within either a stationary or an explosive regime. Third, there is a deterioration in power when detecting late changes under $H_{A,1}$ and $H_{A,4}$, although satisfactory levels can be achieved by using a large(r) training sample size of $m=5,000$. Fourth, the power is relatively lower when using the Student's $t$ distribution errors compared to normal errors. Lastly, there is only a marginal decline observed in the power with a larger value of $\eta $.\newline

Tables (ref) and (ref) in the Supplement provide the empirical power for the R\'{e}nyi type statistics based on Theorem (ref)(ii) under the same setting. When detecting early changes at $k^{\ast }=\lfloor \sqrt{\mathcal{n}}\rfloor $, similar observations as above apply; the monitoring schemes with $\eta =1.3$ and $1.5$ proves to be effective. However, one noticeable difference is that a larger value of $\eta $ is detrimental in the power. In particular, $\eta =1.7$ and $2$ suffer a remarkable loss of power under $H_{A,1}$. As far as late changes ($k^{\ast }=\lfloor 0.5\mathcal{n}\rfloor $) are concerned, the R\'{e}nyi type statistics become much less effective, as predicted by the theory. This is because R\'{e}nyi type statistics are devised for the fast detection of very early changes, whilst being suboptimal for late changes. \newline

It is also worthwhile to examine the stopping time $\tau _{m}$ and $\bar{\tau }_{m}$ in order to investigate the detection delays of our monitoring procedures. Figure (ref) shows the boxplot of the detection delays of $\tau _{m}$ and $\bar{\tau}_{m}$ for a change at $k^{\ast }=\lfloor \sqrt{\mathcal{n}}\rfloor $ under $H_{A,2}$ when $m=500$, $ \mathcal{n}=500$. For the monitoring procedure based on Theorem (ref)( i), it is consistent with our theory that larger values of $\eta $ reduce the detection delay. Considering the R\'{e}nyi type statistics based on Theorem (ref)(ii), there is only a marginal difference in using various values of $\eta $. Comparing the detection delay between the monitoring procedures based on Theorem (ref)(i) and (ii ), we can clearly see the merit of the R\'{e}nyi type statistics for the fast detection of early changes, as evidenced by shorter detection delays.

Empirical illustration

We illustrate our monitoring procedures using daily returns of individual stocks. We focus on four stocks: Apple Inc. (ticker: AAPL, Permno: 14593), Middlefield Banc Corp. (ticker: MBCN, Permno: 14932), Genetic Technologies Ltd (ticker: GENE, Permno: 90899), and NTS Realty Holdings LP (ticker: NLP, Permno: 90508). We download daily returns (without dividend) from the CRSP database.\footnote{ We choose to use the daily returns without dividend, rather than log difference of prices, to avoid the complication due to stock splits.} We consider two periods in order to showcase the detection for four types of changes (in the sense of the four different alternatives in our simulation). Depending on the specific purpose of the researcher, one can choose between the monitoring procedures based on Theorem (ref)(i) and ( ii). Based on our simulations,if the aim is to quickly detect very early changes, we suggest using the R\'{e}nyi type statistics based on Theorem (ref)(ii), with some tolerance for the compromise in size and power; conversely, if the purpose is to have good size control and high power, it is recommended to use the monitoring procedure based on Theorem (ref)(i). In this application, our preference is to have a good balance of size and power, and the procedure based on Theorem (ref)(i) (with the choice of $\eta =0.3$) delivers a good performance with sample sizes similar to the dataset used in this section. \footnote{ We relegate the results using R\'{e}nyi weights in Section (ref) of the Supplement.} Before applying our monitoring procedure, it is necessary to ensure there is no change during the historical training sample. To this end, we use the test developed by horvath2024detecting (horvath2024detecting, labeled as HW(2024) hereinafter). Their test is to detect changes in GARCH(1,1) processes without assuming stationarity, which can accommodate either stationary or nonstationary historical sample. A rejection of their test indicates there is no change of $(\alpha ,\beta )$ in the GARCH(1,1) during the historical sample. At the same time, we are keen to understand which the type of four changes may occur. Thus, we firstly examine whether our historical sample is stationary or not by employing the nonstationarity test developed by francq2012strict (francq2012strict, labeled as FZ(2012) hereinafter). At the end of our monitoring horizon, we use the FZ(2012) test again to check the stationarity of the samples after the change (if there is one).

Change from a stationary regime

To illustrate changepoint detection from a stationary regime, we choose the training period of 2016--2019 (1007 trading days) and the monitoring period of 2020--2021 (507 trading days). The training period is before the outbreak of COVID-19, while the monitoring period is in the pandemic. We apply our monitoring procedure for the stocks of AAPL and MBCN during this period. Table (ref) (Columns 1 and 2) reports the results of the sequential monitoring procedure, as well as other information, including HW(2024) test, FZ(2012) test, and parameter estimates. HW(2024) test indicates that there is no parameter change during the training sample for AAPL and MBCN. The nonstationarity test of FZ(2012) indicates that they are both stationary during the training sample. Our sequential monitoring detects a change of AAPL on July 31st, 2020 and a change of MBCN on May 8th, 2020. Based on the FZ(2012) test for the sample after the change, we can conclude that AAPL experienced a change from a stationary to another stationary regime, while MBCN shifted from a stationary regime to a nonstationary one. Figure (ref) contains returns series (upper panel) during the monitoring period and the detector versus the boundary function (lower panel).

Change from a nonstationary regime

We now consider detection from a nonstationary regime, and use 2007--2010 (1011 trading days) as the training period and 2011-2022 (505 trading days) as the monitoring period. The training period covers the global financial crisis (GFC), while the monitoring period follows the GFC but includes the European debt crisis. In this period, we monitor GENE and NLP. The results of the sequential monitoring procedure, alongside other supporting information, are displayed in Columns 3 and 4 of Table (ref). Based on the HW(2024) test, we cannot reject that the return series of GENE and NLP have change in the training period. As evidenced by the nonstationarity test of FZ(2012), both stocks are in the nonstationary regime during the training period. Our sequential monitoring procedure reveals a change of GENE on April 27th, 2011 and a change of NLP on September 4th, 2012. After applying FZ(2012) nonstationarity test on the sample after the change, it is found that the change of GENE is from a nonstationary regime to a stationary regime, whilst the change of NLP is from a nonstationary to another nonstationary regime. It is also interesting to note that GENE after the change is in a strict stationary regime, but not in a second-order stationary regime. Figure (ref) shows their returns series (upper panel) during the monitoring period and the detector versus the boundary function (lower panel).

Conclusions and discussions

In this paper, we complement the existing literature on (ex-ante) testing for bubble phenomena by proposing a family of weighted, CUSUM-based statistics to detect changes in the parameters of a GARCH(1,1) process. Our monitoring procedure can be applied irrespective of whether, in the training sample, the observations are stationary or explosive, and it is able to detect all types of changes: (a) from a stationary to another stationary regime (which is helpful in order to avoid the issues concerning the consistent estimation of a GARCH(1,1) process spelt out in hillebrand2005neglecting); (b) from a stationary to an explosive regime (which contains information on the possible inception of a bubble); (c) from an explosive to a stationary regime (which, in the light of the previous point, could shed light on the cooling off the turbulence associated with a bubble on a financial market); and (d) from an explosive to another explosive regime (which, depending on the direction of the change - towards a more or a less explosive regime - could indicate whether exuberant volatility is heating up or cooling down). Technically, we propose two families of statistics, both based on weighted versions of the CUSUM process of the quasi-Fisher scores: one family uses lighter weights, and it is designed to detect, optimally, changes occurring not immediately after the start of the monitoring horizon; the other family uses heavier, R\'{e} nyi-type weights, which make it more sensitive to changepoints occurring immediately after the end of the training period. For both cases, we study the limiting distribution of the detection delays; to the best of our knowledge, no such results exist for the case of a GARCH(1,1) models, and no results in general exist for the case of R\'{e}nyi statistics. Given the interest in the detection of bubble phenomena, and the scant amount of contributions in the context of detection of changes in the volatility, we believe that our paper should be a useful addition to the toolbox of the financial econometrician.

{ {\ } }

\setcounter{section}{0} \setcounter{subsection}{-1} \setcounter{subsubsection}{-1} \setcounter{equation}{0} \setcounter{lemma}{0} \setcounter{theorem}{0}

Implementation guidelines, further Monte Carlo evidence, and further empirical results

Practical guidance on the implementation

In this section, we provide a detailed step-by-step guidance on implementing our sequential monitoring procedures (based on Theorem (ref)), which could be useful for researchers who want to apply it without studying the underlying theory. The steps are as follows:

Additional simulation results

Further results - empirical application using R\'{e}nyi weights

Figure (ref) presents the detector $\mathcal{D}_m(k)$ (with the choice of $\eta=1.5$) versus the boundary function $\bar{\mathfrak{g}} _m(k)$ based on Theorem 3.1 (ii) for the four stocks during the same periods analysed in Section 5. The monitoring procedure based on the R\'enyi type statistics detects a change of AAPL on March 2nd, 2020, a change of MBCN on March 13th, 2020, a change of GENE on April 27th, 2011, and a change of NLP on September 4th, 2012. Such results illustrate the merit of R\'enyi type statistics for the fast detection of very early changes, in particular for AAPL (nearly 5 months earlier) and MBCN (about 2 months earlier), compared to the procedure based on Theorem 3.1(i).

\setcounter{subsection}{-1} \setcounter{subsubsection}{-1} \setcounter{equation}{0} \setcounter{lemma}{0} \setcounter{theorem}{0}

Technical lemmas

Henceforth, we use $\nabla $ and $\nabla ^{2}$ for the gradient vector and the Hessian matrix. We begin with some facts and notation which we will use throughout the proofs.\newline

If (ref) holds, then there is a unique, stationary, non anticipative sequence $\{x_{i},-\infty <i<\infty \}$ satisfying

(cf. berkes2003garch, and francq2019garch). Next we define

$\mbox{\boldmath${\theta}$}=(\omega ,\alpha ,\beta )^{\top }$. The stationary version of the likelihood function is

and we define

Similarly, the stationary version of ${\mathcal{r}}_{m,k}( \mbox{\boldmath${\theta}$})$ is

berkes2003garch and francq2019garch showed that the matrix

exists and it is non singular.

After the change the observations are given by

Next we assume that (ref) holds. The observations are generated by (ref) but in contrast to the previous case, the sequence is explosive and cannot be approximated with a stationary sequence. However, we can approximate the gradient of the log likelihood function with a stationary sequence, we need only minor modifications of the arguments used before. We use the following sequences to approximate the gradient vector of the log likelihood function when (ref) holds for an explosive sequence:

and

Let

\setcounter{subsection}{-1} \setcounter{subsubsection}{-1} \setcounter{equation}{0} \setcounter{lemma}{0} \setcounter{theorem}{0}

Proofs