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.
79,136 characters · 18 sections · 75 citation commands
{.26in} \thispagestyle{empty}
This paper introduces a novel method of bubble detection based on strictly stationary noncausal autoregressive processes and their tail process representation during a bubble. In this context, a simple statistic with an asymptotic normal distribution is proposed to provide an ex ante warning of a potentially emerging bubble. Ex post, this statistic is a tool for measuring the duration of the bubble and identifying its starting and ending dates. The proposed diagnostic method detects single and multiple bubbles with rates of explosion, determined by the noncausal autoregressive coefficients of the process. The bubbles can burst vertically to zero, or decline at a slower rate, depending on the presence of a causal component in the strictly stationary mixed causal-noncausal process. Our approach is also applicable to jump detection in the conventional causal autoregressive processes with heavy-tailed non-Gaussian error distributions. In addition, the proposed method yields simple prediction formulas for the time to peak of a bubble and its duration.
During the bubble growth and/or decline phases, the strictly stationary noncausal autoregressive process admits a tail process representation. The tail process was first defined by basrak2009regularly as the weak limit of finite-dimensional distributions of time series, conditional on the occurrence of an extreme value. This definition was further generalized to the spectral tail process by kulik2020heavy, which is based solely on normalized observations exceeding a high threshold, and interpreted as "a model for the clusters of exceedances". The tail process representation is valid for geometrically ergodic Markov processes, and consequently also for mixed causal-noncausal autoregressive MAR(1,1) and purely noncausal MAR(0,1) processes, where the notation MAR($r,s$) indicates a mixed model of causal order $r$ and noncausal order $s$ with locally explosive patterns including bubbles and spikes. As shown, for example, by gourieroux2017local, fries2019mixed and cavaliere2020bootstrapping, the MAR(1,1) and MAR(0,1) processes are well-suited for modeling variables such as commodity and cryptocurrency prices. We show that these processes can also represent the time series of green stock prices that display bubbles and spikes, allowing us to examine bubbles that occurred after year 2020.
The MAR processes have recently received attention in applications to bubble forecasting and testing. de2025forecasting consider the forecasts of extreme trajectories of $\alpha$-stable MAR processes, based on a measure describing the conditional distribution of a normalized path of the process, following a large value. The method is applied to predict the occurrences of the El Niño and La Niña phases of the temperatures and winds of the Pacific Ocean. In comparison, our approach is applicable to MAR models with any heavy-tailed error distribution and focuses on detecting rather than predicting the bubbles. In addition, it is computationally simple as it does not require determining the length of a future path as a tuning parameter. Blasq also consider the $\alpha$-stable distributed processes. They introduce a test for bubbles based on the detection of a large future shock in the forward-looking component of a MAR process and apply it to oil prices. Compared to their approach, our proposed test statistic has the advantage of being easy to compute and having a known asymptotic normal distribution.
Alternative methods for detecting bubbles, based on the Dickey-Fuller augmented test (e.g., SADF and GSADF), assess unit roots and explosive regimes in nonstationary autoregressive processes [phillips2011explosive, phillips2015testing]. Our approach departs from these conventional methods by employing mixed autoregressive causal-noncausal models, which capture locally explosive patterns observed in the data. Unlike those detection approaches, we assume that the series follows a strictly stationary non-Gaussian process in which bubbles are an inherent part of the dynamics, rather than distinguishing a stationary and non-stationary (unit root) regime of a time series. In addition, our approach accommodates local explosive patterns with various explosion and burst rates, which are estimable.
The rest of the paper is organized as follows. Section (ref) reviews the causal-noncausal autoregressive processes and their estimation methods, with a focus on the semi-parametric Generalized Covariance (GCov) estimator. Section 3 studies the behavior of the MAR process and its tail dynamics during a bubble period. Section 4 develops test statistics to detect the bubbles and determine their duration. In Section 5, we apply our approach to investigate the presence and duration of "green bubbles" in the green energy stock market. Section (ref) concludes. Appendices A and B contain additional technical results, and Appendix C presents the simulation tables discussed in Section 4.3 of the main paper. The following notation is used: $\{y_t, t \in \mathbb{Z}\}$, denotes the strictly stationary mixed autoregressive (MAR) process representing the green stock prices, $y$ denotes a high threshold selected among the admissible values of $y_t$ which can be exceeded at an exogenous date, and $N$ is a random variable representing distance in time to the peak of a bubble.
The causal-noncausal models represent stationary processes characterized by locally explosive patterns, such as bubbles and spikes. The univariate causal-noncausal models were examined, for example, by breid1991maximum and lanne2011noncausal, and extended to multivariate analysis by lanne2013noncausal, gourieroux2017noncausal, gourieroux2023generalized, and davis2020noncausal. In applied research, causal-noncausal models were used to study various economic and financial variables, including Bitcoin prices [hencic2015noncausal, cavaliere2020bootstrapping], stock market indices [gourieroux2017local], commodity prices [hecq2021forecasting, lof2017noncausality], and inflation rates [lanne2013noncausal, hecq2023predicting]. The main advantage of these models is their ability to capture complex nonlinear patterns, such as local trends and conditional heteroskedasticity, while still resembling traditional linear time series models in terms of the specification. However, the standard Box–Jenkins approach to identifying and estimating linear time series processes does not apply here, as it relies on the assumption of Gaussian errors. Under the assumption of Gaussian errors, the causal and noncausal dynamics cannot be distinguished (see gourieroux2015pricing). Therefore, for identification, non-Gaussian error distributions are required in the autoregressive causal–noncausal processes.
A strictly stationary univariate mixed causal-noncausal autoregressive MAR($r,s$) model is defined as:
where the error term $\epsilon_t$ is non-Gaussian, independent, identically distributed (i.i.d.) and such that $E(|\epsilon_t|^{\delta}) < \infty$ for $\delta > 0$ [Gourieroux and Zakoian (2015)]\footnote{This condition allows the errors to have infinite variance, and possibly infinite mean: for $\delta \geq 2$ the second-order moments exist, for $\delta \in [1,2)$ the errors have infinite variance, but finite first-order moment, for $\delta \in (0,1)$ the errors have no first-order moments. lanne2011noncausal assumes errors with zero mean and a finite variance. }. The polynomial $\Phi(L)$ in the lag operator $L$ is of order $r$. The polynomial $\Psi(L^{-1})$ in the leading operator $L^{-1}$ is of order $s$. Both polynomials $\Phi(L)$ and $\Psi(L^{-1})$ have roots outside the unit circle.
The MAR($r,s$) process (ref) admits a unique strictly stationary solution, which is a two-sided moving average of order infinity MA($\infty$):
in past, present and future shocks, with $c_0 = 1$ [breid1991maximum]. This MA($\infty$) representation exists and is unique, and the coefficients $c_h$ on past and future errors are distinct and uniquely defined, provided that $(\epsilon_t)$ are non-Gaussian [See, breid1991maximum and gourieroux2015uniqueness for errors with finite variance, and with infinite moments, respectively]. When $y_t$ is purely noncausal (resp. causal), the coefficients $c_h$ are zero for all $h > 0$ (resp. $h < 0$). Thus, a purely causal process is determined only by the past and present shocks, while a purely noncausal process is influenced only by the present and future shocks. For $r=s=1$, we obtain the MAR($1,1$) process:
with $|\psi|<1, |\phi|<1$, which is purely causal (resp. noncausal) if $\psi= 0$ (resp. $\phi= 0$). For each of these pure processes, the effects of a large $\epsilon_t$ are easily distinguished, as a large error leads to a (vertical) jump if $\psi = 0$ and $\phi>0$, and an explosive bubble with a (vertical) burst if $\psi > 0$ and $\phi=0$.
The MAR($1,1$) process can be decomposed into the following unobserved components [lanne2011noncausal]:
and
gourieroux2016filtering show that $u_t$ is $\epsilon$-noncausal (dependent on the future and present values of $\epsilon$) and $y$-causal (dependent on the past and present values of $y$), representing the regular dynamics of $y_t$. In contrast, $v_t$ is $\epsilon$-causal (dependent on the past and present values of $\epsilon$) and $y$-noncausal (dependent on the future and present values of $y$), representing the explosive part of the process, including bubbles and spikes. Henceforth, $u_t, v_t$ are called the unobserved causal and noncausal components of $y_t$, respectively.
The process $y_t$ has the following deterministic representation based on the above unobserved components that can be used for simulations and bootstrapping: $$ y_t = \displaystyle \frac {1}{1-\phi \psi} (\phi v_{t-1} + u_t), \; \mbox{or} \;\;y_t = \displaystyle \frac {1}{1-\phi \psi} (v_t + \psi u_{t+1}). $$
We observe that $y_t$ is a linear function of the first lag of $v_t$ and of the current value of $u_t$. Alternatively, $y_t$ can be expressed as a linear function of the current value of $v_t$ and of the first lag of $u_t$.
The autocovariances of the latent components defined in eq. (3) and (4), respectively, help us distinguish the bubble episode in the process $\{y_t\}$. We observe that the autocovariances at lag $h \geq 1$ of $u_{t+h}$ and $v_t$ conditional on $y_t=y$ are time-varying, in general, provided that $y$ is not an extreme [see Appendix B]. When $y$ is large, their behavior is different and can be examined using the tail process described in Section 3. In practice, the latent components are computed given the estimated values of the autoregressive parameters. The estimation methods are discussed in the next section.
One way to estimate and identify MAR($r,s$) models is by using a parametric non-Gaussian Maximum Likelihood (ML) approach [see, e.g. hecq2016identification]. Alternatively, the semi-parametric Generalized Covariance (GCov) estimator can be used, which does not require any distributional assumptions on the errors other than satisfying the i.i.d. and non-Gaussianity conditions, the latter one being required for identification. The GCov is a one-step estimator that is consistent, asymptotically normally distributed, and semi-parametrically efficient. It can achieve parametric efficiency in special cases [gourieroux2023generalized]. The GCov minimizes a portmanteau-type objective function involving nonlinear autocovariances, i.e., the autocovariances of nonlinear transformations of model errors, which successfully identify the causal and noncausal dynamics [chan2006note].
Let us consider the nonlinear transformations $a\left(\epsilon_t\right)=a_1(\epsilon_t),...,a_K(\epsilon_t)$ of the error process that increase its dimension from 1 to $K$. These transformations satisfying the regularity conditions given in gourieroux2023generalized are introduced to identify the noncausal and nonlinear dynamics and to ensure the existence of moments of the transformed errors. Let $\hat{\Gamma}^a\left(h; \theta\right), h=1,...,H$ denote the autocovariance matrices of the transformed errors at lags $h=0,...,H$, with $\hat{\Gamma}^a\left(0; \theta\right)$ representing their variance, and $\theta$ the vector of autoregressive parameters. The GCov estimator $\hat{\theta}_{T}$ minimizes the following objective function:
where $Tr$ denotes the trace of a matrix and $L_T( \hat{\theta}_T,H)$ is the value of the objective function at its minimum\footnote{When the number of transformations $K$ is large and the inversion of the variance matrix becomes difficult, we can replace $\hat{\Gamma}^a\left(0; \theta\right)$ in (ref) with $diag\left(\hat{\Gamma}_d(0; \theta)\right)$ containing only the diagonal elements of $\hat{\Gamma}_a(0; \theta)$. This latter version of the GCov estimator is no longer semiparametrically efficient, as it is not optimally weighted [see gourieroux2023generalized, cubadda2011testing]. Another strategy for handling problems with the inversion of matrix $\hat{\Gamma}^a\left(0; \mbox{\boldmath $\theta$}\right)$ in high-dimensional settings is to replace the GCov estimator with the regularized RGCov estimator proposed by giancaterini2025regularized, which can preserve the semiparametric efficiency under suitable conditions.}.
The choice of an informative set of transformations ($a_{k}, k=1,..., K$) depends on the specific time series under investigation. For example, in financial applications, linear and quadratic functions can be selected, such as $a_1(\epsilon_t) =\epsilon_{t}$, $a_2(\epsilon_t) = \epsilon_{t}^2$. This implies that $a_1$ is a linear function of errors in a causal-noncausal process, while $a_2$ transforms the error term by squaring it for each $t=1, \dots, T$. In application to MAR processes with heavy-tailed error distributions, such as $\alpha$-stable distributions, including Cauchy, or t-student distribution with low degrees of freedom, one can use the square root of absolute values, or other fractional powers to ensure the existence of moments of the transformed errors and the asymptotic normality of the GCov estimator.
Moreover, the objective function minimized in (ref) and evaluated at the estimated parameter $\hat{\theta}_T$ can be used to test the fit of the model. Specifically, we test the null hypothesis $H_0 : \{\Gamma^a_0 (h) =0, \; h=1,...,H\},$ using the statistic
which, under the implicit null hypothesis of serial independence of errors, has an asymptotic chi-square distribution with degrees of freedom equal to $H K^2 -dim(\theta)$ [gourieroux2023generalized].
Our approach to bubble analysis assumes that the process of interest follows a strictly stationary MAR($r,s$) process with a heavy-tailed error distribution. As mentioned earlier, among the heavy-tailed distributions are the $\alpha$-stable distributions including the Cauchy distribution and t-student distributions with low degrees of freedom. In this Section, we focus our attention on MAR processes of orders $r$ and $s$ such that the combined order $p=r+s$ is less than or equal to 2, so that the process of interest is either a MAR(0,1), i.e., a purely noncausal process of order 1, or a MAR(1,1), or a causal MAR(1,0). These processes are geometrically ergodic. gourieroux2017local and fries2019mixed also show that the MAR(0,1) and MAR(1,1) processes are Markov of order 1 and 2, respectively. The results can be generalized to higher-order MAR($r,s$) processes, which are Markov too [fries2019mixed, Proposition 3.1].
Consider the locally explosive MAR(0,1) and MAR(1,1) processes. It has been shown in the literature that during a bubble episode, and conditional on $y_t >y$, where $y$ is large, the causal-noncausal MAR(0,1) and MAR(1,1) processes with $\alpha$-stable distributed errors have a distinct dynamic [fries2022conditional], which leads to different behavior of the conditional autocovariances of its latent components. In this Section, we describe these dynamics and generalize the results in fries2022conditional using the concept of a spectral tail process, based only on observations above a high threshold $y$. Our results are more general than those in fries2022conditional or de2025forecasting as we do not need to assume a distribution such as the $\alpha$-stable for the whole distribution, but only conditions on the tail process. We consequently define a bubble in this paper as follows:
Definition: The bubble occurs when the observations exceeding a high threshold from a strictly stationary autoregressive noncausal process with i.i.d. errors and a heavy-tailed error distribution with tail index $\alpha$ increase (decline) approximately at the rate of a spectral tail process.
The spectral tail process captures and measures the extremal dependence, which determines the shape and duration of a bubble. It is characterized as follows:
Proposition 1 [kulik2020heavy, Chapter 15.3]: Let $\{y_t, t \in \mathbb{Z}\}$ be a strictly stationary process with i.i.d. errors $\epsilon_t, t=1,2,...$ and a heavy-tailed error distribution with (Pareto-type) tails of tail index $\alpha$ and with a two-sided MA representation:
$$y_t = \sum_{h=-\infty}^{+ \infty} c_h \epsilon_{t-h},$$
with nonnegative coefficients. The process $\left( \displaystyle \frac {y_{t+h}}{y_t}\right)_h, h \in \mathbb{Z}$, converges in distribution to the (spectral) tail process $(X_h)$, so that conditional on a large $y_t>y$, where $y$ is a high threshold, we have : $${\cal L} \left( \frac{y_{t+h}}{y_t},\; h=-H,...,H| y_t>y \right) \stackrel{d}{\rightarrow}{\cal L} (X_h, \; h=-H,...,H),$$
where $\stackrel{d}{\rightarrow}$ denotes convergence in distribution and ${\cal L}$ stands for "law".
It is natural to interpret this convergence as an approximation $\displaystyle \frac {y_{t+h}}{y_t} \stackrel{d}{\approx} X_h$ for $h \in \mathbb{Z}$ [drees2020peak], where
and where $X_0=1$, and $N$ is an integer-valued random variable, independent of $y$ and such that:
By setting $X_0=1$, we consider an upward (positive) bubble. Henceforth, we focus our attention on this case for ease of exposition as well as for the pattern of the bubbles observed on the "green" indicators investigated in this paper. Similar results are easily obtained for downward bubbles, conditional on $y_t<y$ where $y$ tends to $-\infty$ with $X_0=-1$.
Proposition 1 can be applied to the MAR(1,1) and MAR(0,1) processes as follows. Let us first consider the MAR(1,1) process defined in Section 2.1, assuming positive coefficients $\phi$ and $\psi$.
Proposition 2: The strictly stationary MAR(1,1) process with i.i.d. errors $\epsilon_t, t=1,2,...$ and a heavy-tailed error distribution with tail index $\alpha$: $$(1 - \phi L)(1-\psi L^{-1})y_t = \epsilon_t,$$
(i) admits a two-sided MA($\infty$) representation with the coefficients: $c_h = \displaystyle \frac {1}{1-\phi \psi} \psi^{-h}$, if $ h \leq 0$, \;\;and \;\;$c_h = \displaystyle \frac {1}{1-\phi \psi} \phi^{h}$, if $ h \geq 0$.
(ii) the tail process $ X_h$ is such that: $$X_h = \psi^{-1} X_{h-1} 1\hskip -3pt \mbox{l}_{N \leq - h} + \phi X_{h-1} 1\hskip -3pt \mbox{l}_{N > - h},$$
with $X_0=1$, and the probability $ P[N=h]$:
$P[N=h] = \displaystyle \frac {\psi^{-h \alpha}}{[\frac{1}{1-\phi^{\alpha}}+ \frac{1}{1-\psi^{\alpha}} -1]} $, if $ h \leq 0$, and $P[N=h] = \displaystyle \frac {\phi^{h \alpha}}{[\frac{1}{1-\phi^{\alpha}}+ \frac{1}{1-\psi^{\alpha}} -1]} $, if $ h \geq 0$,
with $0^0 = 1$, by convention.
Proof: See Appendices A1-A2.
The random variable $N$ determines the time to the peak of a bubble and its overall duration. We observe that the distribution of variable $N$ is a mixture of two geometric distributions (see Proposition 4).
Corollary 1: In the MAR(1,1) process, we have: $$X_h = \phi^h 1\hskip -3pt \mbox{l}_{N > Max(-h,0)} + \phi^{-h-N} \psi^{-N} 1\hskip -3pt \mbox{l}_{0<N \leq -h} + \phi^{h+N} \psi^{N} 1\hskip -3pt \mbox{l}_{-h <N \leq 0} + \psi^{-h} 1\hskip -3pt \mbox{l}_{N \leq Min(-h,0)}$$
During a bubble episode, we can distinguish the phases of growth and decline. The noncausal persistence of MAR(1,1) processes determines the rate at which a bubble keeps growing up to time $t-N$. The negative power of $\psi$ indicates that the growth is explosive while $N<0$ given that $|\psi|<1$ is assumed for strict stationarity. The above result is consistent with Proposition 4.2 in fries2022conditional which describes the dynamics of a MAR(1,1) process conditional on a large value $y_t>y$ and $\frac{y_{t+1}}{y_t} = \frac{y_{t+2}}{y_{t+1}} = 1/\psi$. In the MAR(1,1), conditional on $y_t>y$, the bubble either keeps growing to the next value $y_{t+1} = \frac{1}{\psi} y_t$, or it bursts and decreases to $y_{t+1} = \phi y_t$. This latter pattern is also consistent with formula (v) given in the proof of Corollary 1, Appendix A.2a, for $N > 0$, which additionally implies that the bubble bursts at $N=0$.
In a MAR(0,1), the bubble bursts vertically to 0, while in the causal AR(1), i.e., MAR(1,0), we observe a jump, followed by a decline determined by the autoregressive coefficient, as shown below.
Corollary 2: The strictly stationary MAR(0,1) process with i.i.d. errors $\epsilon_t, t=1,2,,..$ and a heavy-tailed error distribution with tail index $\alpha$: $$(1-\psi L^{-1}) y_t = y_t - \psi y_{t+1} = \epsilon_t,$$
admits a one sided MA($\infty$) representation with the moving average coefficients:
$$ c_h = \psi^{-h},\;\;\mbox{if}\; h \leq 0, \;\;\mbox{and} \;\; c_h = 0, \;\; \mbox{if} \; h > 0. $$
The tail process is such that: $$X_h = \psi^{-1} X_{h-1} 1\hskip -3pt \mbox{l}_{N \leq -h},$$
with $X_0=1$ and the probability $ P[N=h]$:
$P[N=h] = 0$, if $ h > 0$, and $P[N=h]= (1-\psi^{\alpha}) \psi^{-h \alpha} $, if $ h \leq 0$,
Because the bubble bursts vertically in a MAR(0,1) process, variable $N$ takes only negative values in that process. It follows from Proposition 2 that the variable $-N$ has a geometric distribution on $\mathbb{N}$ with parameter $\psi^{\alpha}$ during the growth phase of bubble.
The causal autoregressive processes with i.i.d. errors and a heavy-tailed error distribution with tail index $\alpha$ may admit jumps followed by tail process behavior.
Corollary 3: The strictly autoregressive of order 1 (causal AR(1), i.e. MAR(1,0)) process: $$y_t = \phi y_{t-1} + \epsilon_t, $$
with i.i.d. errors $\epsilon_t, t=1,2,,..$ and a heavy-tailed error distribution with tail index $\alpha$ and $0<\phi <1$ admits a one sided MA($\infty$) with coefficients:
$$ c_h = \phi^{h}\;\;\mbox{if}\; h \geq 0, \;\;\mbox{and} \;\; c_h = 0, \;\; \mbox{if} \; h < 0. $$
The tail process is such that: $$X_h = \phi X_{h-1} 1\hskip -3pt \mbox{l}_{N \geq h},$$
with $X_0=1$ and the probability $ P[N=h]$ :
$P[N=h] = 0$; if $ h < 0$, and $P[N=h] = (1-\phi^{\alpha}) \phi^{h \alpha} $; if $ h \geq 0$,
Hence, the path of $y_t$ displays jumps, with a geometric decline. The variable $N$ has a geometric distribution with parameter $\phi^{\alpha}$ during the decline following a jump. Moreover, the behavior of pure noncausal autoregressive processes of order 2 (MAR(0,2) with i.i.d. errors $\epsilon_t, t=1,2,..$ and a heavy-tailed error distribution with tail index $\alpha$ is described in Appendix B.
The latent components of a MAR(1,1) process are defined in equations (3) and (4) as: $$u_{t+1} = (1-\phi L) y_{t+1} = y_{t+1} - \phi y_{t}\; \mbox{and} \;v_t = (1-\psi L^{-1}) y_t = y_t -\psi y_{t+1}, t=1,2,...$$
We observe that the ratios of these latent components divided by $y_t$: $$\displaystyle \frac {u_{t+h}}{y_{t}} = \displaystyle \frac {y_{t+h} - \phi y_{t+h-1}}{y_{t}} \;\mbox{and}\; \displaystyle \frac {v_{t+h}}{y_{t}} = \displaystyle \frac {y_{t+h} - \psi y_{t+h+1}}{y_{t}}$$
are functions of ratios $y_{t+h}/y_t$ for $h \in -H,..,H$. From Proposition 1, it follows that at any time $t$ when $y_t>y$ where $y$ is a high threshold, the ratios of observations $y_{t+h}/y_t \stackrel{d}{\approx} X_h$ for any $h$, where $\stackrel{d}{\approx}$ denotes approximately equal to in distribution. Our goal is to replace the above ratios by the tail components: $U_h = X_{h} - \phi X_{h-1}$ and $ V_h= X_h- \psi X_{h+1}$, respectively. Then, when $y_t>y$ is large, the above ratios become linear deterministic functions of the components $X_h$ of the tail process.
Let us now examine how Proposition 2 can be used to obtain the asymptotic behavior (in distribution) of $U_h, V_h, \; h = -H,..., H$.
Proposition 3: At time $t$ such that $y_t>y$, for a high threshold $y$, we have: $$u_{t+h}/y_t \stackrel{d}{\approx} X_h - \phi X_{h-1} := U_h, \; v_{t+h}/y_t \stackrel{d}{\approx} X_h - \psi X_{h+1} := V_h,$$ and $$U_h = (\psi^{-1} - \phi) X_{h-1} 1\hskip -3pt \mbox{l}_{N \leq -h}, \; V_h = (1-\psi \phi) X_{h} 1\hskip -3pt \mbox{l}_{N > -h-1}.$$
Proof: Let us consider the causal component. It follows from Proposition 2 that: $$U_h = X_h - \phi X_{h-1} = (\psi^{-1} - \phi) X_{h-1} 1\hskip -3pt \mbox{l}_{N \leq -h}.$$ Similarly, we have: $$V_h = X_h - \psi X_{h+1} = X_h - X_h 1\hskip -3pt \mbox{l}_{N \leq -h-1} - \psi \phi X_h 1\hskip -3pt \mbox{l}_{N > -h-1} = (1-\phi \psi) X_h 1\hskip -3pt \mbox{l}_{N > -h-1}.$$
Let us now examine the values of $U_h$ and $V_h$ during a bubble episode. We can show that either one of them is zero during a bubble episode as follows:
Corollary 4: The causal and noncausal tail components $U_h, V_h$ are such that:
$U_h = 0$, if $N>-h <=> t+h > t-N$,
$V_h = 0$, if $N\leq-h-1 <=> t+h < t-N-1$,
Proposition 3 and Corollary 4 can be interpreted as follows. Let us consider an exogenous time $t$ when $y_t>y$, for large $y$. Then $t$ belongs to a bubble episode, with a peak at time $t-N$. The noncausal (resp. causal) tail component is equal to zero after (resp. before) the peak and has the behavior of a MAR(0,1) tail process before the peak (resp. MAR(1,0) tail process after the peak). This result can also be used to determine the behavior of the products of the causal and noncausal tail components. As shown below in Corollary 5, since during the bubble episode either $U_h$ or $V_h$ is zero, their product is zero during the entire bubble episode.
Corollary 5: We have $$\xi_{t,h,k} = \frac{u_{t+h}v_{t+k}}{y_{t}^2} \stackrel{d}{\approx} U_h V_k =0, \;\; \mbox{if} \;\; h-k \geq 1.$$
Proof: We see that $U_hV_k = 0 \iff (N >-h)$ or $(N \leq -k-1)$. This condition is equivalent to $N \in \mathbb{Z}$ that is always satisfied iff:
The result follows. In particular, we have: $$\xi_{t,h}= \displaystyle \frac {u_{t+h+1}v_{t+h}}{y_{t}^p} \stackrel{d} {\approx} U_{h+1} V_h = 0, \; \forall h$$
where $p=r+s=2$ is the combined autoregressive order of the process, and $$\xi_{t,0}= \displaystyle \frac {u_{t+1}v_t}{y_t^p} \stackrel{d}{\approx} U_1 V_0 = 0.$$
In particular, for the MAR(0,1) process with $p=1$, we get:
(i) at $h=0$. $$ \xi_{t,0} = \frac{v_{t}}{y_t} \stackrel{d}{\approx} 1- \psi X_1, $$ where $$X_1 = \psi^{-1} 1\hskip -3pt \mbox{l}_{N \leq -1}. $$
Then, $\xi_{t,0} = 1 - 1\hskip -3pt \mbox{l}_{N \leq -1}$ and is approximately always 0 during the bubble.
(ii) at $h=1$. $$ \xi_{t,1} = \frac{v_{t+1}}{y_t} \stackrel{d}{\approx} X_1-\psi X_2 $$ where $$X_2 = \psi^{-1} X_1 1\hskip -3pt \mbox{l}_{N \leq -2} $$
Then, $\xi_{t,1}$ is approximately equal to 0 as well.
The observed behavior of $\xi_{t,h}$ can be used for testing the hypothesis that the process is in a bubble episode at time $t+h$. This approach is pursued in Section 4. Below, we discuss the behavior of the variable $N$ and associated inference.
This section describes the stochastic properties of the random variable $N$, which determines the duration of a bubble, and of $t-N$ for $N\leq 0$, which is interpreted as the time to peak. We consider first the MAR(1,1) process.
Proposition 4: In the MAR(1,1) process, the distribution of $N$ is a mixture of two geometric distributions with $P[N \geq 0] = \displaystyle \frac {1}{1-\phi^{\alpha}}/(\frac{1}{1-\phi^{\alpha}} + \frac{\psi^{\alpha}}{1-\psi^{\alpha}})$. The distribution of $N$ given $N \geq 0$ is a geometric distribution in $\mathbb{N}$ with parameter $\phi^{\alpha}$, and the distribution of $-1-N$, given $N \leq -1$ is a geometric distribution on $\mathbb{N}$ with parameter $\psi^{\alpha}$.
Proof: From Proposition 2, it follows that: $$ P[N \geq 0] = \displaystyle \frac {1}{1-\phi^{\alpha}} / \left( \frac{1}{1-\phi^{\alpha}} + \frac{\psi^{\alpha}}{1-\psi^{\alpha}} \right).$$
Then, for $h \geq 0$:
$P(N=h| N \geq 0) = \displaystyle \frac {\phi^{h \alpha}}{1-\phi^{\alpha}}$, which is a geometric distribution on $\mathbb{N}$ with the parameter $\phi^{\alpha}$.
For $h <0$, we have:
$P(N=h| N < 0) = \displaystyle \frac {\psi^{(-h-1)\alpha}}{1-\psi^{\alpha}}$, which means that $-1-N$ has a geometric distribution on $\mathbb{N}$ with parameter $\psi^{\alpha}$.
Proposition 4 allows us to infer about the bubble duration from the moments of variable $N$. In practice, one may be interested in finding the marginal expected value of $N$ as it is informative about the shape of a bubble and can also reveal an asymmetry in the growth and burst phases.
Corollary 6: The expected value of $N$ is: $$E(N) = \displaystyle \frac {1}{\frac{1}{1-\phi^{\alpha}} + \frac{1}{1-\psi^{\alpha}} -1} \left( \displaystyle \frac {\phi^{\alpha}}{(1-\phi^{\alpha})^2} - \displaystyle \frac {\psi^{\alpha}}{(1-\psi^{\alpha})^2} \right).$$
Proof: We have
where $Z_1, Z_2$ follow geometric distributions on $\mathbb{N}$ with parameters $\phi^{\alpha}$ and $\psi^{\alpha}$, respectively. It follows that:
The marginal expected value of $N$ depends on the values of $\phi$ and $\psi$ and is symmetric in $\phi$ and $\psi$. In particular, it is equal to $0$ if $\phi=\psi$, which corresponds to a symmetric bubble. In addition, this expectation depends on $\alpha$. The effects of $\phi, \psi, \alpha$ are illustrated in Figure (ref), where we observe that, for a fixed $\psi$, the expected value $E(N)$ increases in $\phi$. For a fixed $\phi$, $E(N)$ increases in absolute values in $\psi$. In addition, the range of values of $E(N)$ depends on the tail parameter $\alpha$.
If the bubble is asymmetric, then $E(N)$ can take either a positive or negative value. Therefore, one may be interested in computing the conditional expectation of $N$ given $N<0$ to approximate the expected time to peak.
It is also interesting to compute the cumulative distribution function (cdf), or survival function of $N$, to derive the quantiles of the distribution of $N$ and obtain a prediction interval for $N$. We find the quantiles $h^*_U$ and $h^*_L$ of the distribution of $N$ at levels $(1-\gamma)$ and $\gamma$, where $\gamma$ is small, which are: $h^*_U \;\; \mbox{such that} \;\; P[N \leq h^*_U] = \gamma, \;\mbox{and} \;\; h^*_L \;\; \mbox{such that} \;\; P[N>h^*_L] = 1-\gamma$. Using small $\gamma$, we separately consider the geometric distributions of the mixture.
Corollary 7: For $h \leq 0$, sufficiently small, we have: $$ P[N \leq h] = \displaystyle \frac {\psi^{-h \alpha}}{1-\psi^{\alpha}} \displaystyle \frac {1}{[\frac{1}{1-\phi^{\alpha}} + \frac{1}{1-\psi^{\alpha}} -1]}, $$
and for $h>0$, sufficiently large, we have: $$P[N>h] = \displaystyle \frac {\phi^{h \alpha}}{1-\phi^{\alpha}} \displaystyle \frac {1}{[\frac{1}{1-\phi^{\alpha}} + \frac{1}{1-\psi^{\alpha}} -1]}. $$
Proof: see Appendix.
The expected value (median or mode) provides a point prediction of the variable $N$. The quantiles of the distribution of $N$ provide the prediction interval for $N$.
Let us now consider the pure noncausal MAR(0,1) process characterized by a vertical bubble burst. For the MAR(0,1) process, the variable $-N$ follows a geometric distribution on $\mathbb{N}$ with parameter $\psi^{\alpha}$ (by Corollary 2 or Proposition 4) during bubble growth. The cumulative probability function of $N$ is:
$$P[N \leq h] = (1-\psi^{\alpha}) \sum_{i \leq h} \psi^{-i \alpha} =(1-\psi^{\alpha}) \sum_{i \geq -h} \psi^{i \alpha} = \psi^{-h \alpha},$$
for $h \leq 0$. The mean of this geometric distribution is equal to $$E(N) = - \displaystyle \frac {1}{1-\psi^{\alpha}},$$
and can be interpreted as the average time to peak. The median is: $$med(N) =\left[ \displaystyle \frac {\log 2}{ \alpha \log (\psi)}\right],$$
where the brackets denote the integer part, and the mode is denoted by $mode(N)=0$. Since the distribution of $-N$ is skewed, the expectation, median, and mode provide three different point predictions. We illustrate the effect of the parameters $\psi$ and $\alpha$ in Figure (ref). We observe that for MAR(0,1), the absolute value of the expected time to peak increases in $\psi$. This is consistent with the dynamics of this process, with a growing phase of a bubble followed by an instantaneous crash. Moreover, the values of $E(N)$ depend on the parameter $\alpha$.
Up to the effect of the integer part, both the median and the expectation of $N$ increase in the absolute value of $\psi^{\alpha}$. To build the prediction interval, we use the quantiles of the geometric distribution, which are determined by inverting the cumulative probability function $P[N \leq h] = \psi^{-h \alpha}$. Given that the quantile at level $\gamma$ is equal to : $$h(\gamma) = - \displaystyle \frac {\log \gamma}{\alpha \log \psi},$$
the prediction interval for $N$ at level 10% is: $$(- \displaystyle \frac {\log 0.05}{\alpha \log \psi}, - \displaystyle \frac {\log 0.95}{\alpha \log \psi} ).$$
Because the variable is discrete, we can round up or down the upper and lower quantiles to integer values to make it more conservative.
In practice, the predictions and prediction intervals of $N$ can be found by replacing the unknown parameters $\phi, \psi$ by their GCov-based estimates based on a sample of $T$ observations, which are also used for computing the test statistics given in the next section. The tail parameter $\alpha$ can be approximated by the Hill estimator [see e.g. eq. (9.5.1), Section 9.5, kulik2020heavy], which suffers from the bias-variance trade-off in finite samples. kulik2020heavy show that under the tail process approach, the Hill estimator is consistent and asymptotically normally distributed.
Let us show how the results derived in Section 3.2 can be used to build diagnostic tools to test the MAR(1,1) process for bubbles. Diagnostics are based on the counterparts of $$\xi_{t,0} = \displaystyle \frac {u_{t+1} v_{t}}{y_t^2} = \left(\frac{y_{t+1}}{y_t} -\phi_0 \right) \left(1-\psi_0 \frac{y_{t+1}}{y_t} \right),$$
where $\phi_0, \psi_0$ denote the true values of the parameters for a given observed process, we wish to estimate from a sample of length $T$. The quantity $\xi_{t,0}$ depends on the observations and the unknown parameters. The unknown parameters can be replaced by consistent and asymptotically normally distributed estimators $\hat{\phi}_T, \hat{\psi}_T$ of $\phi,\psi$ obtained from a sample of $T$ observations. Then, the sample counterpart of $\xi_{t,0}$ is: $$\hat{\xi}_{t,T}(0) = \left( \frac{y_{t+1}}{y_{t}} - \hat{\phi}_T \right) \left(1- \hat{\psi}_T \frac{y_{t+1}}{y_t} \right).$$
The statistics $\hat{\xi}_{t,T}(0), t=1,...,T-1$ can be used as follows. From Corollary 5, we know that if $y_t$ is sufficiently large at time $t$, then $\xi_{t,0} \stackrel{d}{\approx} U_1 V_0$ and that $U_1 V_0 = 0$. Therefore, we expect that $\hat{\xi}_{t,T}(0)$ is close to 0. We can also consider another lag $h$ with: $$ \xi_{t,h} = \displaystyle \frac {u_{t+h+1} v_{t+h}}{y_t^2} = \left(\frac{y_{t+h+1}}{y_t} -\phi_0 \frac{y_{t+h}}{y_t} \right) \left(\frac{y_{t+h}}{y_t} -\psi_0 \frac{y_{t+h+1}}{y_t} \right). $$
It was shown in Corollary 5 that $\xi_{t,h} \stackrel{d}{\approx} U_{h+1} V_h$ and $U_{h+1} V_h = 0$. Its sample counterpart is $$ \hat{\xi}_{t,T}(h) = \displaystyle \frac {\hat{u}_{t+h+1} \hat{v}_{t+h}}{y_t^2} = \left(\frac{y_{t+h+1}}{y_t} - \hat{\phi}_T \frac{y_{t+h}}{y_t} \right) \left(\frac{y_{t+h}}{y_t} - \hat{\psi}_T \frac{y_{t+h+1}}{y_t} \right), $$
which can be used for inference at higher horizons.
Consider the statistic $\hat{\xi}_{t,T}(0)$. The zero value of the transformed tail process can be considered as the true value of some tail parameter $\theta_1 = U_1 V_0$, which is deterministic and equal to 0 if $y_t$ is sufficiently large, and it is stochastic, otherwise. Then, at each exogenous time $t$, we can consider the null hypothesis: $$H_{0,1} = \{ \theta_1 = 0\}.$$
The difficulty is that we have a double asymptotic in the level of threshold $y$ and in the number of observations $T$. To derive reliable confidence bands, we assume that during the bubble episode, the uncertainty in the approximation $\xi_{t,0} \stackrel{d}{\approx} U_1V_0$ is negligible with respect to the asymptotics in $T$. Then, during a bubble episode when $H_{0,1}$ is satisfied and $y_{t} >y$, with large $y$, we have conditional on $y_t, y_{t+1}$:
where $ \sigma_t^2$ is obtained by the delta method as follows:
$$ \sigma_t^2 = V_t \left\{ \left[ -\left( 1 - \psi_0 \frac{y_{t+1}}{y_{t}} \right), \frac{-y_{t+1}}{y_t} \left( \frac{y_{t+1}}{y_{t}} - \phi_0 \right) \right] \sqrt{T} \left(
\right) \right\}.$$
This quantity is then consistently estimated from: $$ \hat{\sigma}_{t,T}^2(0) = \left( \frac{\hat{v}_t}{y_t}, \frac{y_{t+1} \hat{u}_{t+1}}{y_{t}^2} \right) \hat{\Omega}_T \left(
\right),$$
where $\hat{\Omega}_T$ is a consistent estimator of the asymptotic variance matrix of the parameter estimators. Then, under the assumption of a MAR(1,1) process and following a large $y_t$, we have $$| \sqrt{T} \, \hat{\xi}_{t,T}(0)/\hat{\sigma}_{t,T}(0)| \leq 1.96, $$
with an asymptotic probability of 95$\%$. This leads to a functional diagnostic tool that consists of reporting for any time $t$ the quantities $\hat{\xi}_{t, T}, \; t=1,..., T-1$ along with the band
$$\left(\hat{\xi}_{t,T}(0) \pm 1.96 \frac{\hat{\sigma}_{t,T}(0)}{\sqrt{T}} \right).$$
Since the tail process does not depend on the level of $y_t$ and the values of $y_{t+1}$, the width of the band is independent of time $t$, $t=1,..., T$. Then, the times at which the statistic is inside the band are the times associated with a bubble with probability 95%.
This approach is easily extended to other lags. From the statistic $\hat{\xi}_{t,T}(h)$, we can test the null hypothesis: $$H_{0,h} = \{ \theta_h = 0\}$$
where $\theta_h = U_{h+1}V_h$ and $y_{t} >y$. Then, conditional on $y_t, y_{t+h}, y_{t+h+1}$ we have:
where $ \sigma_t^2(h)$ is: $$ \sigma_t^2(h) = V_t \left\{ \left[ \frac{-y_{t+h}}{y_{t}} \left( \frac{y_{t+h}}{y_t}- \psi_0 \frac{y_{t+h+1}}{y_t} \right), \frac{-y_{t+h+1}}{y_{t}} \left(\frac{y_{t+h+1}}{y_{t}} - \phi_0 \frac{y_{t+h}}{y_{t}} \right) \right] \sqrt{T} \left(
\right) \right\}.$$
This quantity is then consistently estimated from: $$ \hat{\sigma}_{t,T}^2(h) = \left( \frac{y_{t+h} \hat{v}_{t+h}}{y_{t}^2}, \frac{y_{t+h+1} \hat{u}_{t+h+1}}{y_{t}^2} \right) \hat{\Omega}_T \left(
\right).$$ \noindent Under the assumption of a MAR(1,1) process and following a large $y_{t}>y$: $$| \sqrt{T} \, \hat{\xi}_{t,T}(h)/\hat{\sigma}_{t,T}(h)| \leq 1.96, $$
with the asymptotic probability of 95 %. This leads to a set of functional diagnostic tools, where for any time $t$ and $h$ the quantities $\hat{\xi}_{t,T}(h), \; t=1,...,T-h-1$ are reported along with the band:
$$\left(\hat{\xi}_{t,T}(h) \pm 1.96 \frac{\hat{\sigma}_{t,T}(h)}{\sqrt{T}} \right).$$
In practice, the above procedure may not distinguish between a bubble and a short-lasting spike with the same growth rate. In addition, there can be times $t$ during the bubble when the statistic is outside the band because either (i) the MAR(1,1) model is misspecified, or (ii) the MAR(1,1) is well specified, but the value of $y_t$ is not sufficiently large.
This method easily provides the test statistics for the MAR(0,1) and MAR(1,0) processes, based on Corollary 5, with $\hat{\xi}_{t,T}(0) = \hat{v}_{t}/y_t$, and $\hat{\xi}_{t,T}(0) = \hat{u}_{t+1}/y_t$, respectively. The above results applied to MAR(1,0) processes provide a tool for jump detection in the causal autoregressive process with a heavy-tailed error distribution. In each case, the confidence bands are independent of the tail parameter $\alpha$.
$\hat{\xi}_{t,T}(0)$ can be used as a bubble detection tool in the MAR(1,1) and MAR(0,1) processes, as it is constant and close to 0 during the bubble growth and decline periods, and it is time-varying otherwise. Hence, the first difference $\Delta \hat{\xi}_{t,T}$ is also constant and close to 0 during the bubble growth and decline periods.
Consider a set of $\hat{\xi}_{t,T}(0)$, evaluated at $t=1,2,..$ following a high threshold value $y_t$. A close to zero value of $\hat{\xi}_{t+1,T}(0)$ following a large $y_t$ is a warning of an upcoming bubble. Since this statistic remains close to 0 throughout the duration of a bubble, the number of values of $\hat{\xi}_{t,T}(0)$, at $t=t+1, t+2,...$ that are not statistically significant, is a measure of the duration of that bubble. The first time $t_j$ when $\hat{\xi}_{t_j,T}(0) \approx 0$ marks the start of the bubble. The last time $t_J$, such that $\hat{\xi}_{t_J,T}(0) \approx 0$ marks the end of the bubble.
Our approach is illustrated in Figure (ref) by a bubble episode of a simulated MAR(1,1) process with $\psi=0.9, \phi=0.3$ and Cauchy distributed errors with scale coefficient 1. The top panel of Figure (ref) displays the bubble episode of the MAR(1,1) process. The first difference $\Delta \hat{\xi}_{t,T}(0)$ of the statistic, defined as:
$$\Delta \hat{\xi}_{t,T}(0) = \Delta \left( \frac{\hat{u}_{t+1} \hat{v}_t} {y_t^2} \right),\; y_t \neq 0$$
is computed and displayed graphically to detect periods when it is approximately constant and close to 0. The statistic $\Delta \hat{\xi}_{t,T}(0)$ indicates the times when the process becomes a "tail process", which can be used to approximate the start and end of a bubble. The bottom panel of Figure (ref) shows the first difference $\Delta \hat{\xi}_{t,T}(0)$. The first difference $\Delta \hat{\xi}_{t,T}(0)$ is approximately constant and close to zero during the entire bubble episode.
Figure 3 below shows an example of the path of a simulated MAR(1,1) process with $\phi=0.3$ and $\psi = 0.9$, and i.i.d. errors with a t(3) distribution.
We observe a bubble with a peak of 52.051 at $t=163$. The estimates of the GCov parameter of the MAR(1,1) model are $\hat{\phi}$=0.3085 and $\hat{\psi}=$0.908, with standard errors of 0.028 and 0.030. These estimates are based on $K=4$ power transformations and $H=3$. The test statistic $\hat{\xi}_{t,T}(0)$ evaluated during the bubble, conditional on observation $y(160)=37.951$ is -0.0220 and is within the confidence interval of $\pm 0.05646 $. Therefore, the null hypothesis $\theta_1 = 0$ at time $t=160$ is not rejected.
Next, we examine the effect of the conditioning values and perform functional diagnostics using $\hat{\xi}_{t,T}(h)$ on increasing horizons $h$=1 to 10, conditional on $y_{127}$ =-0.7522 and $y_{151}$ =21.810. The results are reported in Table (ref). The columns of Table (ref) report: the horizon (col. 1), the value of the test statistic (cols. 2 and 5), the confidence interval of the form $0.0 \pm CI$, and the result of the test of $H_{0,h}: \theta_h=0$ coded 1 for not rejected and 0 for rejected. We observe that, conditional on the large value at time $t=151$, the diagnostics do not reject the null hypothesis of a bubble over the upward-sloping sequence of the next 10 observations. Conditional on $y_t$ being close to the mean value of the process at time $t=127$, the diagnostics reject the null hypothesis of a bubble over the next 10 observations that are close to the mean.
To evaluate the performance of the test in finite samples, we perform the following experiment. We generate MAR(0,1) processes in samples of length T=400 \footnote{We generate a series of length 800 and discard the first and last 200 observations in 1000 replications}, with t(3), t(4), and t(5) distributed errors and with different values of the noncausal the coefficient $\psi$. We choose samples of 400 observations to increase the probability of occurrence of at least one bubble in a simulated path, which we do not observe. Moreover, we want to increase the probability that the order statistic-based quantile estimator used as a threshold is sufficiently close to the true quantile.
We use an upper quantile $y=q_y(0.975)$ as a conditioning large positive threshold, consistently with the statistical literature on exceedances over a high threshold [see e.g. davis2018inference]. Next, we compute the statistic at time $t$ from each simulated path, to study the size of the test. To study the power of the test, we use $y=q_y(0.525)$ as the conditioning value at time $t$ followed by a moderate decrease or increase, and preceded by a random pattern, where by moderate we mean that $y_{t+1}$ is less than or equal to $q_y(0.975)$, and $y_t$ is not zero. The empirical size and power of the test are reported in Tables 4 and 5, Appendix C. Table 6 considers the MAR(1,1) with Cauchy distributed errors. A similar exercise is performed for the MAR(1,1) processes with t-student distributed errors, estimated by the GCov with $K=2$ power transformations and lag $H=4$. The size and power are reported in Tables 7 and 8, Appendix C.
There are some challenges involved in this analysis. It is difficult to accurately estimate the quantiles of the marginal distribution of the process from 400 observations, especially for processes with high persistence and error distributions with low degrees of freedom. Hence, the estimates of $y=q_y(0.975)$ as the conditioning value, obtained from the replicated paths of the process, may vary. Moreover, without observing the trajectory, we do not know a priori whether the process admits a bubble or a jump instead, the latter occurring in the MAR(1,1) processes due to the causal component. In processes with $\psi=0.9$, bubbles are infrequent and grow slowly at a rate of about 1.1. In processes with $\psi=0.6$, the bubbles are more frequent and grow faster. We distinguish the bubble from jumps by checking if, following the conditioning value of $y=q_y(0.975)$, the process grows at a rate close to the theoretical rate. Another difficulty is the asymptotic normality of the estimators, which impacts the results in the MAR(0,1) processes with t(3) distributed errors estimated by the OLS.
In Table 4, we observe that the test over-rejects when the error distribution is t(5), especially for higher values of $\psi$. Table 5 shows that the power of the test is slightly lower when the noncausal persistence is higher. Table 6, illustrating processes with $\psi=0.9$ and Cauchy distributed error, shows that the test tends to under-reject for closer to 0.9 values of parameters $\phi$. In addition, the power of the test decreases for higher values of $\phi$. In Table 7, we observe that the size of the test is closer to the nominal value for moderate values of $\psi = 0.7$ and $\psi=0.8$. The test is conservative for $\psi=0.6$ and over-rejects for $\psi=0.9$ in processes with t(3) distributed errors. For each $\psi$, the test is the most conservative in processes with t(5) distributed errors. This error distribution is the closest to the normal among the three error distributions considered, leading to potential parameter identification problems. We find that the size is close to the nominal value for moderate values of causal persistence $\phi$, regardless of $\psi$. Table 8 shows that the test has good power, especially for $\psi = 0.7$ and $\psi=0.8$. The power deteriorates in processes with t(5) distributed errors for $\psi=0.6$ and $\psi=0.9$ and for higher values of causal persistence, making bubbles harder to distinguish from standard dynamics.
The transition to a clean energy economy has recently gained significant attention, driven by global events such as the COVID-19 pandemic and the Russia-Ukraine war. In fact, in response to these crises, governments, industries, and investors have increasingly prioritized clean energy investments, recognizing their long-term benefits for a stable and environmentally friendly energy system [mohammed2023all]. This growing focus on clean energy is also driven by the need to achieve net zero emissions by 2050, which requires substantial investment in clean energy from both developed and developing countries [khalifa2022accelerating]. In 2023, global investment in the energy sector was estimated at USD 2.8 trillion, an increase of 0.6 trillion USD from five years earlier. Almost all of this growth was directed toward clean energy and infrastructure, increasing total clean energy spending to 1.8 trillion USD, compared to 1 trillion USD for fossil fuels\footnote{IEA (2023), World Energy Outlook 2023, IEA, Paris https://www.iea.org/reports/world-energy-outlook-2023, License: CC BY 4.0 (report); CC BY NC SA 4.0 (Annex A).}. These large investments pose the risk of "green bubbles" -- rapid stock price increases followed by crashes. Specifically, such bubbles occur when overinvestment and speculative behavior drive the market value of clean energy assets beyond sustainable levels. Consequently, rising interest rates can increase financial pressure on investors, potentially triggering bubble bursts and undermining the credibility of the clean energy transition [wimmer2016green].
The literature on green energy stocks is relatively recent and focuses on return analysis. The pioneering paper by henriques2008oil shows that returns on technology stocks and oil prices are both Granger-caused by the returns of green (alternative) energy stocks, based on a vector autoregressive (VAR) model. The relationship between green stocks and oil has also been examined by sadorsky2012correlations, sadorsky2012modeling, and kumar2012stock. These articles consider carbon prices and do not find their significant impact on green energy stocks. In contrast, this impact is evidenced after the year 2007 by managi2013does. The literature has not yet examined the green stock price dynamics to detect and explain the presence of bubbles, which is done in this paper.
We consider the Renixx Index, the WHETF, and the iShare green stock ETFs. For Renixx, we use daily data from \url{https://www.renewable-energy-industry.com/stocks} and sample them at a monthly frequency by taking the last day of each month of the closing price series. When the last day of the month is a bank holiday, we consider the previous available day. The Renixx index tracks the global renewable energy market, covering sectors such as wind, solar, bioenergy, geothermal, hydropower, electronic mobility, and fuel cells. It comprises 30 companies, each of which derives more than 50% of its revenues from these sectors (see \url{www.iwr.de/renixx}). The WHETF tracks the WilderHill Clean Energy ETF, which includes companies listed in the United States that focus on developing cleaner energy and conservation efforts. In particular, WHETF allocates a minimum of 90$\%$ of its total assets to common stocks included in this ETF. It is rebalanced and reconstituted quarterly. The iShare tracks the S$\&$P Global Clean Energy ETF, which provides exposure to the top 30 largest and most liquid publicly traded companies operating worldwide in the clean energy sector, based on a modified market capitalization weighting system.
We investigate Renixx on its whole available sample from January 2002 to February 2024, and both WHETF and iShare from January 2009 to February 2024. We consequently have a total of $T=266$ monthly observations for Renixx and $T=182$ observations for WHETF and iShare, with potentially two bubble patterns for Renixx and a single bubble for both the WHETF and iShare. The Renixx index experienced significant bubbles in 2008 and 2020, coinciding with two major global events: the 2008 financial crisis and the outbreak of the COVID-19 pandemic. In 2008, the financial crisis rocked global markets, leading to widespread economic instability and investor panic. The collapse of major financial institutions, coupled with a credit crunch and falling stock markets, triggered a flight to safety among investors. This risk aversion had a significant impact on the renewable energy sector, with reduced investment in renewable energy projects and a decrease in the demand for renewable energy stocks [see giorgis2024salvation]. Similarly, in 2020, the COVID-19 pandemic caused unprecedented economic disruption around the world. Lockdown measures, supply chain disruptions, and reduced consumer spending resulted in a global economic downturn. The renewable energy sector was particularly affected by the decrease in energy demand due to reduced economic activity and travel restrictions. Figure (ref) shows the path of Renixx between January 2009 to February 2024. WETHF and iShare are displayed in Figures (ref) and 4c, respectively. However, only the COVID period bubble appears in the data, as the observations are from January 2009 to February 2024, and therefore do not include data from the global financial crisis.
Although most studies in green energy finance focus on returns [e.g., henriques2008oil; sadorsky2012correlations], our analysis is based on prices, allowing for a direct examination of bubble dynamics that would otherwise be obscured by log-differencing. In fact, in financial analysis, especially in commodity pricing, it is crucial to focus on the price process rather than returns or first differences. This approach aligns with the financial theory that underlies commodity pricing, as discussed in hull2016options. Transforming the price series into returns or first differences can obscure critical aspects of the price process, including the presence of bubbles. Moreover, differencing can eliminate noncausal components essential to understanding bubble dynamics [giancaterini2022climate]. We employ the detrended cubic spline method to address this issue and eliminate trend components while preserving significant bubble patterns, as detailed in jasiakhall. This method fits a cubic spline, a piecewise function composed of polynomial segments of degree three, to the time series data. The points where these segments connect, known as knots, allow separate cubic polynomials to be fitted within each segment. We place knots every two years to effectively detrend the series, balancing data smoothing and avoiding overfitting. By applying this approach, we successfully isolate the cyclical variations that contain the bubble patterns from the trend component. The detrended series are shown in the right panels of Figure (ref).\footnote{An alternative detrending technique that may help us preserve the bubble patterns is the HP filter method [see giancaterini2022climate, hecq2023predicting]. However, as indicated by jasiakhall, the HP filter requires you to choose a value for the smooth parameter lambda, which is typically a function that increases with the sampling frequency of the data. In our examination of monthly data, a high lambda value may result in computational inaccuracies, thus justifying our preference for using a cubic spline to detrend the data. However, we have investigated several detrending methods, including the HP filter, as well as polynomial trends of different degrees (results upon request). We only report the results with the spline detrending approach as it passes the absence of nonlinear serial dependence test [jasiak2023gcov]. Finally, note that an alternative detrending approach based on unobserved components has been developed by Blasques2023Observation}.
The detrended data is non-Gaussian, as evidenced by the Kolmogorov-Smirnov and Shapiro-Wilk's tests, which both reject the null hypothesis of normality. Table (ref) presents the summary statistics, which confirm that the distributions of the series are non-Gaussian, given the reported excess kurtosis. We test the spline-detrended data for causal and noncausal persistence by using the test introduced in jasiak2023gcov [see Section (ref)]. We choose $K=2$, including the time series and its squares as (non)linear transformations, and $H=3$ as the number of lags in the objective function of the test. The null hypothesis of the absence of nonlinear serial dependence is rejected since the test value is 514.98, while the critical value at a 5$\%$ significance level from the chi-square distribution is $\chi^2(12)$ = 21.026.
We apply the GCov estimator to obtain the parameter estimates of MAR models for the three series. Specifically, we use $H=K=2$ and $a_j(\epsilon_t)=\epsilon_t^j$, for $j=1,2$, in (ref), i.e., the residuals and their squares as transformations in the GCov estimation of Renixx [see cubadda2023optimization] and $a_j(\epsilon_t) = \log(|\epsilon|)^j$ for $j = 1, 2$, for WHETF and iShare. Table (ref) presents the estimation results. To identify the dynamics of the underlying processes, we evaluate all combinations of $r$ and $s$ such that $r+s=p$. We begin with $p=1$ and gradually increase $p$ to find the values of $r$ and $s$ that provide i.i.d. residuals based on the GCov specification test. The last row of (ref) shows the results of the GCov test. This approach allows us to identify the Renixx index and the two ETFs as MAR($1,1$) processes. The GCov specification test results indicate that these processes are correctly specified and provide a good fit to the data.
Figure 5 shows the estimated latent components $\hat{u}_t$ and $\hat{v}_t$ defined in eq. (3) and (4) for the three series of interest. We observe that the changes in the noncausal component anticipate the dynamics of the observed series, followed by the causal component. In the next section, these estimates are used to compute the statistic $\hat{\xi}_{t,T}(0)$ for each series.
The bubbles in the price series can manifest themselves as periods of high volatility. We estimate local variances to detect periods of high volatility in Renixx and the ETFs. We consider a rolling window of 5 observations and estimate the local variance for the three series. The plot of rolling estimates provides a graphical tool for preliminary analysis. Figure (ref) displays the rolling estimates of the variances of the Renixx index and two ETFs. Two bubble periods, one that occurred during the financial crisis of 2008 and one associated with the COVID pandemic, are observable in Figure (ref). Since the data on ETFs before 2009 are not available and the observations on bubbles during that period are incomplete, we focus on the bubble observed during the COVID period.
To detect and determine the dates of bubbles in the three price series, we focus specifically on the period between 12/01/2020 and 05/31/2021. Figure (ref) illustrates how the method introduced in Section 4 can be applied to Renixx and the ETFs to test for bubbles, and to determine the dates at which that bubble starts and ends. Figures (ref) and (ref) show the test statistics $\hat{\xi}_{t,T}(0)$ computed over that period and the associated confidence band.
The conditioning threshold values $q$ used in the analysis correspond to the upper tail quantiles at levels $97.5\%$ in Renixx, $98\%$ in WHETF, and $96\%$ in iShare. We use the nonparametric forward estimator of the conditional distribution function introduced in davis2018inference to estimate the probability that $y_{t+1}/y_t \geq 1$ conditional on $y_t>q$, where $q$ denotes the tail quantile from each sample given above. This provides us with the probabilities of tail process dynamics for each time series, following a threshold value $q$. For Renixx and WHETF, the estimated conditional probability is $80\%$, and for iShare it is $40\%$. Given these non-zero conditional probability estimates, we proceed with the analysis\footnote{Note that these results need to be considered with caution as the approach is valid in large samples.}.
We observe that the bubble in Renixx and WTETF started on 12/01/2020 and ended on 03/01/2021. The bubble in iShare ended on 04/30/2021. Note that in Figure (ref) the confidence intervals depend on the parameter variances reported in Table (ref) and are larger for Renixx compared to the other estimated processes.
Given a tail index value, the time to peak depends on the rate of bubble growth. It takes longer for the bubble to grow when $\hat{\psi}$ is large. As shown in Section 3.3, the time to peak is also determined by the tail parameter, which is estimated from the Hill estimator R-package. For $\hat{\alpha}=1.3$, the expected time to peak for Renixx is 1.5 months (computed as the conditional expectation of $N$ given $N<0$), and the probability that $N$ exceeds 3 months is 0.232.
Given $\hat{\alpha}=1.7$ for WHETF, the expected time to peak is 4.5 months, and the probability that it exceeds 5 months is 0.37. With $\hat{\alpha}=1.15$, the time to peak for iShare is slightly longer than 1 month, and the probability that it exceeds 3 months is 0.17.
This paper introduced a new statistical test for detecting financial bubbles, derived from the tail process representation of mixed causal-noncausal autoregressive models. Unlike traditional unit root-based methods, our approach is grounded in strictly stationary noncausal dynamics, which provide a natural framework for capturing locally explosive behavior. Classical tests such as SADF or GSADF are often sensitive to heavy tails or may misinterpret short-lived spikes as bubbles. In contrast, our test exploits the tail process properties of causal and noncausal dynamics, and has a statistic remaining approximately constant and close to zero during a bubble episode. The additional inference on time-to-peak and bubble duration allows us to separate persistent speculative bubbles from short-lived spikes and regular fluctuations.
We applied univariate MAR($r,s$) models to the Renixx index and two green energy ETFs (WHETF and iShare), estimating the underlying dynamics using the GCov approach. This combination provides a robust identification strategy: the GCov estimator ensures reliable parameter estimation under non-Gaussianity, while the proposed test captures the tail process behavior that signals bubble formation. The empirical results reveal bubbles in all three series, showing that the test performs well in detecting and dating speculative episodes in green financial markets.
For bubble detecting in real time, the proposed method can be applied sequentially when new observations become available. In that case, the model needs to be re-estimated, and the test statistic recalculated at each step. Then, one may need to adjust the test level to the nominal size using, e.g., the approaches given in virani2019sequential. Moreover, the causal-noncausal process can be forecasted out-of-sample [lanne2012optimal, gourieroux2016filtering, gourieroux2026nonlinear, hecq2021forecasting, hecq2023predicting] and the method proposed in this paper can be applied to the time series augmented by the forecast period.