EconBase
← Back to paper

Bubble Modeling and Tagging: A Stochastic Nonlinear Autoregression Approach

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.

61,230 characters · 10 sections · 46 citation commands

Rendered from LaTeX for readability, not typeset faithfully. Citation keys are highlighted; maths is left as source; figures, tables and equation environments are summarised rather than reproduced; unrecognised commands are greyed out so nothing is silently dropped. Email addresses are removed.

Bubble Modeling and Tagging: A Stochastic Nonlinear Autoregression Approach

\def\spacingset#1{ {#1}} \spacingset{1}

\if11 \fi

\if01 {

center[center omitted — 106 chars of source]

} \fi

abstractEconomic and financial time series can feature locally explosive behavior when a bubble is formed. The economic or financial bubble, especially its dynamics, is an intriguing topic that has been attracting longstanding attention. To illustrate the dynamics of the local explosion itself, the paper presents a novel, simple, yet useful time series model, called the stochastic nonlinear autoregressive model, which is always strictly stationary and geometrically ergodic and can create long swings or persistence observed in many macroeconomic variables. When a nonlinear autoregressive coefficient is outside of a certain range, the model has periodically explosive behaviors and can then be used to portray the bubble dynamics. Further, the quasi-maximum likelihood estimation (QMLE) of our model is considered, and its strong consistency and asymptotic normality are established under minimal assumptions on innovation. A new model diagnostic checking statistic is developed for model fitting adequacy. In addition, two methods for bubble tagging are proposed, one from the residual perspective and the other from the null-state perspective. Monte Carlo simulation studies are conducted to assess the performances of the QMLE and the two bubble tagging methods in finite samples. Finally, the usefulness of the model is illustrated by an empirical application to the monthly Hang Seng Index.

{\it Keywords:} Causal process, Financial bubble, Rational expectation, SNAR model, Speculative bubble.

\spacingset{1.9}

Introduction

Financial speculative bubbles have been attracting longstanding attention of economists and financial practitioners as an economic crisis often originates along with a burst of a bubble. In reality, however, economic or financial bubbles cannot be avoided. The presence of bubbles is partially evidenced by that many economic or financial time series possess locally explosive behavior and a subsequent burst, with such a phenomenon appearing periodically. Studying the dynamics of bubble thus becomes important and intriguing.

One classical definition of the bubble is the deviation of the market price from its fundamental value (a sum of discounted future dividends) in rational expectation price models. An important model of the rational bubble is initiated by Blanchard1982Bubbles, where the bubble process is captured via a simple stochastic autoregression (AR) with a fixed explosive rate and an absorbing state zero. Their model is then extended by Evans1991Pitfalls via adopting a stochastic rate of explosion. Primarily, the bubble is regarded as an explosive nonstationary process, which motivates to test its presence via unit root and cointegration tests Diba1988a,Diba1988b. Recently, this idea is further developed by Phillips2011Dating, Phillips2011Explosive, Phillips2015Testing, Phillips2015, Harvey2019, Harvey2020, Taoyu, Kurozumi2023Time, Esteve2023 and references therein. On the other hand, Evans1991Pitfalls also notes that periodical collapse of bubbles makes the bubble paths look more like a stationary process. Within a stationary framework, Gourieroux2017Local find that noncausal AR(1) models can characterize multiple local explosions in time series. Then this noncausal approach to bubble modelling has been extended to high-order mixed causal-noncausal\footnote{The definition of `causal' in time series can be referred to BrockwellDavis1991. Early, Lanne2011Noncausal study statistical inference on noncausal AR models with an application to U.S. inflation dynamics.} time series models, see, for example, Gouriroux2016FilteringPA, Fries2019Mixed, Cavaliere2018BootstrappingNA, Davis2020NoncausalVA, and Fries2021ConditionalMO. However, one shortcoming of the noncausal approach invites computational challenge and many resampling methods are needed. To bypass this shortcoming, Blasques2022A propose a new observation driven model with time-varying parameters and study its probabilistic properties and statistical inference. Nevertheless, their estimation heavily depends on the choice of the survival function and the asymptotics can be obtained only for a part of parameters. Motivated by all above facts, we here present a new simple time series model to describe the dynamics of bubbles.

In this paper, a first-order stochastic nonlinear autoregressive (SNAR) model $\{y_t\}$ is defined as

flaligny_t=s_t\phi_0 |y_{t-1}|+\varepsilon_t, \quad t\in\mathbb{Z}:=\{0, \pm1,\pm2, ...\},

where $\phi_0\in \mathbb{R}$, $\{\varepsilon_t: \, t\in\mathbb{Z}\}$ is a sequence of independent and identically distributed (i.i.d.) random variables on some basic probability space $(\Omega, \mathcal{F}, \mathbb{P})$, and independent of i.i.d. binary variables $\{s_t: \,t\in\mathbb{Z}\}$ with $\mathbb{P}(s_t=1)=p_0=1-\mathbb{P}(s_t=0)$, $p_0\in[0, 1)$.

Clearly, when $\phi_0>1$, $y_t$ is explosive in the periods where $s_t=1$ and creates an excursion which stops once $s_t=0$. Fig. (ref) illustrates two simulated paths of the SNAR model ((ref)) with $\varepsilon_t\stackrel{i.i.d}{\sim}\mathcal N(0,6^2)$. We can observe periodically local explosions followed by bursts. And the 'shape' of a bubble before its burst is similar to a quadratic curve, which conforms to some practitioners' pointview that the accumulation of a bubble before bursting resembles a parabola. Further, when $\mathbb{P}(s_t=1)=1$, the SNAR model ((ref)) reduces to an absolute AR model or a special threshold AR model with threshold parameter zero, which is studied in Tong1990, LiTong2020 and references therein.

figure[figure omitted — 265 chars of source]

Major contributions of our paper are as follows.

First, we introduce a simple yet useful time series model, a SNAR model, for modeling the dynamics of bubble process. We then prove that the model is always strictly stationary and geometrically ergodic under minimal assumptions on innovation and the probability $p_0$. Within a causal and stationary framework, when the parameter $\phi_0>1$, our SNAR model still displays periodically local explosions and collapses. It can create long swings or persistence and can then be used to portray the bubble dynamics. More importantly, our model is always causal in the classical sense of time series. Compared with the noncausal bubble models in the literature, our model facilitates the understanding of the dynamics of bubble, and is quite simpler and more convenient in application. Additionally, our model avoids computational burden of the noncausal approach and keeps away from the choice of the survival function in the time-varying parameter model in Blasques2022A.

It is worth mentioning that a related model to ours is a stochastic AR initiated by Blanchard1982Bubbles,\footnote{It is a simple case of a random coefficient AR model in Nicholls1982RandomCA if we regard $s_t\phi_0$ as a coefficient. Random coefficient AR models have received great attention and been widely used in econometrics, finance, engineering, among others, due to their flexibility, parsimonious representation and analytical tractability.} which is defined as $y_t=s_t\phi_0 y_{t-1}+\varepsilon_t$, $t\in\mathbb{Z}$, where $\{s_t\}$ and $\{\varepsilon_t\}$ are defined in ((ref)). It can also create long swings or persistence Johansen2013. Unfortunately, it can usually generate many negative local explosions even if $\phi_0>1$ and $\varepsilon_t\stackrel{i.i.d}{\sim}\mathcal N(0,1)$. This does not meet our common sense on financial bubble.

Second, we consider the quasi-maximum likelihood estimation (QMLE) of the SNAR model and establish its strong consistency and asymptotic normality under minimal assumptions on innovation and probability parameter $p_0$, regardless of infinite variance or heavy-tailedness of the model.

Third, we develop a new model diagnostic checking statistic since the classical portmanteau test is invalid for our model.\footnote{The reason is that the residuals cannot be obtained.}

Fourth, we consider two methods for tagging the bubbles, one from the residual point of view and the other from the null-state perspective. The problem of bubble detection has been studied in the literature; see for example Phillips2011Dating, Phillips2015Testing, Phillips2015, Blasques2022A, Kurozumi, and references therein. Existing results in this direction, however, have been mainly developed by viewing the bubble as a separate process occurring on an unknown but deterministic time interval within the observation period. The current paper, on the other hand, aims to incorporate the bubble mechanism into the data generating process to provide a stationary statistical model that can capture and interpret bubbles, including their formations and collapses. As a result, unlike existing results that typically assume the bubbles to persist for an adequate duration to achieve their consistent detection, the bubbles in the current model can be transient and thus the problem of bubble tagging can be more challenging in the current setting. For this, we consider two approaches, where the first one utilizes the nonlinear autoregressive residual from the proposed model and the second one is constructed from a hypothesis testing point of view. For both methods, we provide theoretical quantification on the finite-sample probability of correct tagging under reasonably mild conditions. Monte Carlo simulation results are provided to assess the finite-sample performance of the proposed QMLE and bubble tagging methods.

The remainder of the paper is organized as follows. Section (ref) investigates strict stationarity and geometric ergodicity of model ((ref)). Section (ref) considers the QMLE with its asymptotics. Section (ref) studies model diagnostic checking. Section (ref) considers the problem of bubble tagging, where two approaches are considered with their finite-sample probability bounds studied. Section (ref) carries out Monte Carlo simulation studies to assess the finite-sample performances of the QMLE and the two bubble tagging methods. Section (ref) gives an empirical application to illustrate the usefulness of the model. Section (ref) concludes. All technical proofs are relegated to the Supplementary Material.

Probabilistic Properties of the SNAR Model

The aim of this section is to prove the strict stationarity and geometric ergodicity of model ((ref)) under a very mild condition. We will prove the following result, using the approach developed by MeynTweedie for establishing the geometric ergodicity of Markov chains. This result is important and is a theoretical foundation of statistical inference for model ((ref)).

theoremSuppose that $\mathrm{(i)}$ $\{\varepsilon_t\}$ is i.i.d. and independent of i.i.d. binary variables $\{s_t\}$ with $0\leq p_0<1$, and $\mathrm{(ii)}$ $\varepsilon_1$ has a positive density on $\mathbb{R}$ with $\mathbb{E}(\log^+|\varepsilon_1|)<\infty$. Then there exists a strictly stationary, nonanticipative solution to $\{y_t\}$ in model $(\ref{eq.model})$ and the solution is unique and geometrically ergodic.

Next, consider the existence conditions on moments of $y_t$. Clearly, when $\mathbb{E}\varepsilon_t=0$, $\mathbb{E}\varepsilon_t^2<\infty$, and both $\varepsilon_t$ and $s_t$ are independent, then $\mathbb{E}y_t^2=\mathbb{E}\varepsilon_t^2/(1-p_0\phi_0^2)<\infty$ if $p_0\phi_0^2<1$. Fig. (ref) plots the strict stationarity region of $y_t$ with $\mathbb{E}y_t^2<\infty$.

figure[figure omitted — 244 chars of source]

Further, if $\mathbb{E} \varepsilon_t^3=0$, then a simple algebraic calculation gives the kurtosis of $y_t$:

flalign*\mathrm{kurtosis}(y_t)=\frac{\left\{6p_0\phi_0^2+\mathrm{kurtosis(\varepsilon_t)} (1-p_0\phi_0^2)\right\}(1-p_0\phi_0^2)}{1-p_0\phi_0^4},\quadif\quad p_0\phi_0^4<1.

In particular, when $\varepsilon_t\sim \mathcal{N}(0, 1)$, then

flalign*\mathrm{kurtosis}(y_t)=\frac{3(1-p_0^2\phi_0^4)}{1-p_0\phi_0^4}>3,\quadif\quad 0<p_0\phi_0^4<1,

which implies that $\{y_t\}$ must be heavy-tailed.

Quasi-Maximum Likelihood Estimation

Let $\theta_0=(\phi_0, p_0, \sigma_0^2)'$ be the true parameter with $\sigma^2_0=\mathbb{E}\varepsilon_t^2$. Denote by $\theta=(\phi, p, \sigma^2)'$ be the parameter and by $\Theta$ be the parameter space. Assume that the observations $\{y_0, y_1,..., y_n\}$ are from model ((ref)) with the true value $\theta_0$. Clearly, under Assumption (ref) below, it follows that $\mathbb{E}(y_t|y_{t-1})=p\phi |y_{t-1}|$ and $\mathrm{Var}(y_t|y_{t-1})=p(1-p)\phi^2y_{t-1}^2+\sigma^2$. Then the (conditional) log-quasi-likelihood function (omitting a constant) is

flalign*L_n(\theta)=\sum_{t=1}^n\ell_t(\theta):= \sum_{t=1}^n\left\{ \log\left[p(1-p)\phi^2y_{t-1}^2+\sigma^2\right] +\frac{(y_t-p\phi |y_{t-1}|)^2} {p(1-p)\phi^2y_{t-1}^2+\sigma^2} \right\}.

The QMLE of $\theta_0$ is defined as

flalign*\widehat{\theta}_n=\arg\min_{\theta\in\Theta}L_n(\theta).

To study the asymptotics of $\widehat{\theta}_n$, the following two assumptions are needed.

assumption$\{\varepsilon_t\}$ is i.i.d. and independent of i.i.d. binary variables $\{s_t\}$ with $p_0<1$. Further, $\varepsilon_1$ has a positive density on $\mathbb{R}$ with zero mean and finite variance.
assumptionThe parameter space $\Theta$ is a compact subset of $\{\theta=(\phi, p, \sigma^2)': \phi\neq0,\, 0< p<1, \, 0<\sigma^2<\infty\}$.

The following theorems states the strong consistency and the asymptotic normality of $\widehat{\theta}_n$.

theoremIf Assumptions $\ref{assumption}$-$\ref{para:space}$ hold, then $\widehat{\theta}_n\rightarrow \theta_0$ a.s. as $n\rightarrow\infty$.
theoremIf Assumptions $\ref{assumption}$-$\ref{para:space}$ hold, $\mathbb{E}(\varepsilon_t^4)<\infty$, and $\theta_0$ is an interior point of $\Theta$, then \begin{flalign*} \sqrt{n}\big(\widehat{\theta}_n-\theta_0\big)\stackrel{d}\longrightarrow \mathcal{N}\left(0, \,\mathcal{J}^{-1}\mathcal{I}\mathcal{J}^{-1}\right),\quad as $n\to\infty$, \end{flalign*} where `$\stackrel{d}\longrightarrow$' stands for convergence in distribution, \begin{flalign*} \mathcal{J} =\mathbb{E}\left\{\frac{1}{[p_0(1-p_0)\phi^2_0y_{t}^2+\sigma^2_0]^2}\mathbf{A}_t\right\} +\mathbb{E}\left\{\frac{2y_{t}^2}{p_0(1-p_0)\phi^2_0y_{t}^2+\sigma^2_0}\right\} \mathbf{D}, && \end{flalign*} \begin{flalign*} \mathcal{I} &= \mathbb{E}\left\{\frac{p_0(1-p_0)(1-2p_0)^2\phi^4_0y_{t}^4 +4\sigma_0^2p_0(1-p_0)\phi_0^2y_{t}^2+(\kappa_4-\sigma^4_0)}{[p_0(1-p_0)\phi^2_0y_{t}^2+\sigma^2_0]^4} \mathbf{A}_t\right\}\\ &\quad+\mathbb{E}\left\{\frac{4y^2_{t}}{p_0(1-p_0)\phi^2_0y_{t}^2+\sigma^2_0}\right\}\mathbf{D} +\mathbb{E}\left\{\frac{2p_0(1-p_0)(1-2p_0)\phi_0^3y_{t}^4+2|y_{t}|\kappa_3} {[p_0(1-p_0)\phi^2_0y_{t}^2+\sigma^2_0]^3}\mathbf{B}_t\right\}&& \end{flalign*} with $\kappa_3=\mathbb{E}(\varepsilon_t^3)$, $\kappa_4=\mathbb{E}(\varepsilon_t^4)$, and \begin{flalign*} \mathbf{A}_t&= \left( \begin{array}{ccc} 4p^2_0(1-p_0)^2\phi^2_0 y_{t}^4 & 2p_0(1-p_0)(1-2p_0)\phi^3_0 y_{t}^4 & 2p_0(1-p_0)\phi_0 y_{t}^2 \\ & (1-2p_0)^2\phi^4_0 y_{t}^4 & (1-2p_0)\phi^2_0 y_{t}^2 \\ & & 1 \end{array} \right),\\ \mathbf{B}_t&=\left( \begin{array}{ccc} 4p_0^2(1-p_0)\phi_0 y_{t}^2 & p_0(3-4p_0)\phi^2_0 y_{t}^2 & p_0 \\ & 2(1-2p_0)\phi^3_0 y_{t}^2 & \phi_0 \\ & & 0 \end{array} \right),\quad and\quad \mathbf{D}=\left( \begin{array}{ccc} p^2_0 & p_0\phi_0 & 0 \\ & \phi^2_0 & 0 \\ & & 0 \end{array} \right).&& \end{flalign*} Here, the elements in the lower triangles can be completed by symmetry.

Remark 1. From the expressions of the matrices $\mathcal{I}$ and $\mathcal{J}$ in Theorem (ref), we can see that each element of random matrices within the expectation is bounded and thus it is unnecessary to require moment conditions on $y_t$ for the asymptotics of $\widehat{\theta}_n$.

Remark 2. In practice, to make statistical inference on $\theta_0$, we must estimate the matrices $\mathcal{I}$ and $\mathcal{J}$. From the proof of Theorem (ref), they can be consistently estimated by

flalign*\widehat{\mathcal{I}}_n=\frac{1}{n}\sum_{t=1}^n\frac{\partial\ell_t(\widehat{\theta}_n)}{\partial\theta} \frac{\partial\ell_t(\widehat{\theta}_n)}{\partial\theta'}\quadand\quad \widehat{\mathcal{J}}_n=\frac{1}{n}\sum_{t=1}^n\frac{\partial^2\ell_t(\widehat{\theta}_n)} {\partial\theta\partial\theta'},

respectively. Note that the plug-in method is here invalid since both $\kappa_3$ and $\kappa_4$ in $\mathcal{I}$ cannot be estimated from the residuals. Additionally, due to the constraint $p_0\in (0, 1)$, the Delta method may be needed to construct confidence intervals of $p_0$. If necessary, for example, we can consider the transformation $g(p)=\log[(1-p)/p]$ for $p\in (0, 1)$. Note that

flalign*\sqrt{n}\big(g(\widehat{p}_n)-g(p_0)\big) \stackrel{d}{\longrightarrow}\frac{\lambda_p}{p_0(1-p_0)}\mathcal{N}(0, 1),

if $\sqrt{n}(\widehat{p}_n-p_0)\stackrel{d}{\longrightarrow} \mathcal{N}(0, \lambda_p^2)$. Then, for any fixed $\alpha\in(0, 1)$, a $100(1-\alpha)\%$ confidence interval of $p_0$ is

flalign*\left[\Big\{1+\exp\Big[g(\widehat{p}_n) -\frac{\lambda_p\,z_{\alpha/2}}{\sqrt{n}\widehat{p}_n(1-\widehat{p}_n)}\Big]\Big\}^{-1},\,\: \Big\{1+\exp\Big[g(\widehat{p}_n) +\frac{\lambda_p\,z_{\alpha/2}}{\sqrt{n}\widehat{p}_n(1-\widehat{p}_n)}\Big]\Big\}^{-1}\right],

where $z_{\alpha/2}$ is the lower $\alpha/2$-quantile of the standard normal.

Model Diagnostic Checking

Diagnostic checking is important for time series modeling. The most commonly used tool is the portmanteau test, which depends on the autocorrelation of the residuals or the squared residuals, see, e.g., McLeod1983, Li1994, LiWK2004, and Chen:Zhu2015. However, such the portmanteau test fails for the adequacy of model (ref) since the residuals cannot be obtained. In fact, the residuals should be theoretically calculated by $\widehat{\varepsilon}_t=y_t-s_t\widehat{\phi}_n|y_{t-1}|$ with the initial value $y_0$ for $i=1,...,n$. Unfortunately, the latent variables $\{s_t: 1\leq t\leq n\}$ are unknown and prevent us from getting $\{\widehat{\varepsilon}_t\}$.

To check the adequacy of model (ref), we introduce a new portmanteau test, which is constructed via a transformation of an uncorrelated sequence. Note that the sequence $\{y_t-p_0\phi_0|y_{t-1}|: t\in\mathbb{Z}\}$ is still uncorrelated when $\mathbb{E}y_t^2<\infty$ after replacing $s_t$ by its mean $p_0$ in (ref). To reduce the dependence on the moments of $y_t$, similar to Ling:2005, Ling:2007, we adopt a self-weight method and then define a new sequence $\{\eta_t\}$ by

flalign\begin{split} \eta_t:=&\eta_{t,a}=(y_t-p_0 \phi_0 |y_{t-1}|)I(|y_{t-1}|\leq a)\\ =&\phi_0(s_t-p_0)|y_{t-1}|I(|y_{t-1}|\leq a)+\varepsilon_tI(|y_{t-1}|\leq a),\quad t\in\mathbb{Z}, \end{split}

where the constant $a$ is positive and is called a tuning parameter, and $I(\cdot)$ is an indicator function. Clearly, $\{\eta_{t}\}$ is strictly stationary and ergodic since it is a measurable function of strictly stationary and ergodic sequence $(y_{t-1}, s_t, \varepsilon_t)'$. Further, by the mutually independence among $s_t$, $\varepsilon_t$, and $y_{t-1}$, a simple calculation yields that

flalign\mathbb{E}\eta_{t}=0,\,\, \sigma_\eta^2:=\mathbb{E}\eta_{t}^2=p_0(1-p_0)\phi_0^2\mathbb{E}\left\{y_1^2I(|y_{1}|\leq a)\right\} +\sigma_0^2\mathbb{P}(|y_{1}|\leq a), \,\, \mathbb{E}(\eta_{t}\eta_{t-k})=0,

for $k\geq 1$. That is, $\{\eta_t\}$ is always a white noise under Assumption (ref). Moreover, it is also a martingale difference sequence.

Let $\widehat\eta_{t}=(y_t-\widehat p_n \widehat \phi_n |y_{t-1}|)I(|y_{t-1}|\leq a)$, $1\leq t\leq n$. Intuitively, its sample autocorrelation $\widehat \rho_{nk}$ should be close to zero if model specification is correct, where

flalign*\widehat \rho_{nk}=\frac{\sum_{t=k+1}^n(\widehat \eta_{t}-\bar \eta)(\widehat \eta_{t-k}-\bar \eta)}{\sum_{t=1}^n(\widehat \eta_{t}-\bar \eta)^2}\quadwith\quad \bar \eta=n^{-1}\sum_{t=1}^n\widehat \eta_{t}.

Denote $\widehat{\bm\rho}_{n}=(\widehat \rho_{n1},\dots,\widehat \rho_{nM})'$, where $M\ge 1$ is a fixed positive integer. The following theorem gives the limiting distribution of $\widehat{\bm\rho}_n$.

theoremSuppose the conditions in Theorem $\ref{thm.asym}$ hold. If model (ref) is correctly specified, then $\sqrt{n} \widehat{\bm\rho}_n\stackrel{d}\longrightarrow\ \mathcal{N}(0,\,\mathbf{UGU}')$, where $\mathbf{G}=\mathbb{E} [v_tv_t']$ with \begin{flalign*} v_t=\left(\frac{\eta_t\eta_{t-1}}{\sigma_\eta^2},\dots, \frac{\eta_t\eta_{t-M}}{\sigma_\eta^2}, \left(-\mathcal{J}^{-1} \frac{\partial \ell_t(\theta_0)}{\partial\theta}\right)'\,\right)', \end{flalign*} and $\mathbf{U}=[\mathbf{I}_M,\sigma_\eta^{-2}(u_{1},\dots, u_{M})'(p_0,\phi_0,0)]$ with $u_k=-\mathbb{E}\{\eta_{t-k}|y_{t-1}|I(|y_{t-1}|\leq a)\}$ and $\sigma_\eta^2$ being defined in $(\ref{eta_sigma})$.

Based on Theorem (ref), our portmanteau test statistic is defined as

flalign*Q_M=n\widehat{\bm\rho}_n^{\,'} \big(\widehat{\mathbf{U}}_n \widehat{\mathbf{G}}_n \widehat{\mathbf{U}}_n'\big)^{-1} \widehat{\bm\rho}_n,

where $\widehat{\mathbf{U}}_n$ and $\widehat{\mathbf{G}}_n$ are consistently sample counterparts of $\mathbf{U}$ and $\mathbf{G}$, respectively. Under conditions of Theorem (ref), we have that $Q_M \stackrel{d}\longrightarrow\ \chi^2_M$.

Remark 3. In application, we must choose the tuning parameter $a$ when our test statistic $Q_M$ is used. According to the suggestion in the literature, see, for example, Ling:2005, Ling:2007, we can let $a$ be the 90% or 95% quantile of data $\{|y_1|,...,|y_n|\}$. Many practical experience shows that this self weight performs well, although it may not be optimal and there exist some other choices. Further, from Section (ref), we can see that $\mathbb{E} y_t^4<\infty$ if $p_0\phi_0^4<1$. Note that the condition $p_0\phi_0^4<1$ is verifiable or testable. Thus, if $\widehat{p}_n\widehat{\phi}_n^{\,4}<1$, then we can let $a=+\infty$, that is, the self-weight is redundant.

Bubble Tagging

An important problem in economic or financial data analysis is to tag the bubbles, which can help the government or financial institutions to respond timely to resolve a potential financial crisis. Being able to successfully tag a bubble can also create lucrative trading opportunities. On the other hand, many economic studies can benefit from meaningful tagging of history bubbles to help understand the economic status and explain the reasoning behind certain economic behaviors in the history. The problem of bubble tagging, however, is often not easy and requires sophisticated statistical modeling and treatment. For this, Phillips2011Dating considered decomposing the asset price process into a fundamental component determined by expected future dividends and an explosive bubble component, and proposed a recursive testing procedure. Phillips2015 modeled the null hypothesis as a random walk with asymptotically negligible drift and studied the limit theory of a dating algorithm for bubble detection; see also Phillips2015Testing. We also refer to recent papers by Blasques2022A and Kurozumi, and references therein, for additional literature. Existing results on bubble detection, nevertheless, were mainly developed in a nonstationary framework for which the bubble mechanism is not incorporated into the underlying stationary process and treated as a separated period on the timeline.

A distinguishable feature of the current paper is to incorporate the bubble mechanism into the data generating process to provide a stationary statistical model that can capture and interpret bubbles. Unlike existing results where bubbles are assumed to persist for an adequate duration to achieve consistent detection, bubble tagging in the current stationary framework can be more challenging as bubbles, especially transient bubbles or bubbles that only last for a very short time, can be easily mixed with large white noise observations. We in the following provide two different methods for bubble tagging in the current stationary framework.

A Residual-Based Method for Bubble Tagging

Given model ((ref)), we consider the difference

flalignr_t:=y_t-\phi_0|y_{t-1}| =\left \{ \begin{array}{ll} \varepsilon_t, & if $s_t=1$,\\ \varepsilon_t-\phi_0|y_{t-1}|, &if $s_t=0$. \end{array} \right.

Since $\phi_0>1$ is generally assumed in real applications, it is expected that $r_t$ will be relatively smaller for time points where $s_t = 0$. Therefore, a natural approach is to tag $s_t = 0$ if the difference $r_t \leq c_r$ for some threshold $c_r$. We in the following provide some theoretical understanding of such a tagging method. To this end, we introduce an auxiliary process $\{z_t\}$, where $z_0 = \epsilon_0$ and

flalignz_t = \phi_0|z_{t-1}| + \varepsilon_t,\quad t = 1,2,\ldots.

Unlike the full model specified in ((ref)) that is stationary for which a bubble can collapse, the process $\{z_t\}$ defined above is a pure bubble process that is explosive and nonstationary. In particular, for any $k \geq 1$, $z_k$ shares the same distribution as $y_{t+k}$ if a bubble forms at time $t+1$ and persists through time $t+k$, which we call a $k$-th cumulative bubble. This also relates to the excursive period with duration $k$ in financial applications; see for example our data analysis in Section (ref). For consistent tagging of bubbles, it is generally required that $k \to \infty$, namely the bubble has to persist for a growing horizon of time; see for example Phillips2015 and references therein.

propositionFor any time $t$, if the innovation distribution is symmetric, then the conditional probability that the collapse of a $k$-th cumulative bubble will be correctly tagged by the aforementioned method equals to $\mathbb{P}(z_k \geq -c_r)$, namely the marginal probability that the auxiliary explosive bubble process will exceed the same threshold in the other direction.

In practice, a threshold of $c_r < 0$ is typically chosen for $r_t$ defined in ((ref)), and as a result $-c_r > 0$ will be a positive threshold for $z_k$. Given the explosive nature of the bubble process $\{z_k\}$, it is expected that $\mathbb{P}(z_k > -c_r) \to 1$ as $k \to \infty$ for any chosen threshold $-c_r > 0$, and as a result the probability that the collapse of a $k$-th cumulative bubble will be correctly tagged increases to one as $k \to \infty$. This resonates the result of Phillips2015 but in very different settings. To be more specific, Phillips2015 assumed that the bubble period is a deterministic segment with an increasing number of time points within the whole observation period, while the current setting treats the bubble as an integrated part of an underlying stationary process in ((ref)).

In real applications the parameter $\phi_0$ in ((ref)) is unknown and we propose to plug in the QMLE and tag $\widehat s_t = 0$ if the residual $\widehat r_t=y_t-\widehat \phi_n|y_{t-1}| < c_r$ for some threshold $c_r$. We in the following provide some empirical reference rules for choosing the threshold $c_r$.

itemize• Rule 1 (hard threshold). Set $c_r$ as the $(1-\widehat p_n)$-th quantile of $\{\widehat r_t\}$. Such a choice is simple yet effective, and it can be seen from our simulation results in Section (ref) that it is also reasonably robust to different choices of innovation distributions.

If in particular the innovations are normal with $\varepsilon_t\sim \mathcal{N}(0, \sigma_0^2)$, then we can in addition consider the following choice of thresholds. For this, we use $(\widehat{\phi}_n,\widehat{p}_n,\widehat{\sigma}^2_n)$ to denote the QMLE of $(\phi_0,p_0,\sigma_0^2)$ as in Section (ref).

itemize• Rule 2 (conditional likelihood). From ((ref)), conditioning on $y_{t-1}$, $r_t=\varepsilon_t+\mu$, where $\mu$ can be considered as a location parameter and takes only two values, i.e., $\mu=0$ or $-\phi_0|y_{t-1}|$, corresponding to $s_t=1$ or $s_t=0$, respectively. By comparing conditional likelihood of $r_t$ given $y_{t-1}$ for each time $t$, we can determine the value $\mu=0$ or $-\phi_0|y_{t-1}|$ and then determine whether we need to tag this time $t$. Equivalently, for each time $t$, we can set $c_{r,t}=-\widehat \phi_n |y_{t-1}|/2$ and tag the time $t$ if $\widehat r_t < c_{r,t}$. • Rule 3 (time-varying quantile). Conditioning on $y_{t-1}$, we have $r_t\sim p_0 \mathcal N(0, \sigma_0^2)+(1-p_0)\mathcal N(-\phi_0|y_{t-1}|, \sigma_0^2)$. Based on such a conditional distribution, a natural choice for $c_{r,t}$ is \begin{flalign*} c_{r,t}=\inf\left\{r\in\mathbb{R}: \widehat p_n\Phi\left(r/\widehat{\sigma}_n\right)+ (1-\widehat{p}_n)\Phi\big((r+\widehat{\phi}_n|y_{t-1}|)/\widehat{\sigma}_n\big)>1-\widehat{p}_n \right\} \end{flalign*} for each time $t$, where $\Phi(\cdot)$ is the cdf of $\mathcal N(0,1)$. • Rule 4 (Bayesian). Based on the Bayes' rule, conditioning on $y_{t-1}$, we can obtain the posterior probability mass function of $s_t$ given $r_t$, that is, $\mathbb{P}(s_t=1|r_t)$ and $\mathbb{P}(s_t=0|r_t)$. Theoretically, for each time $t$, if $\mathbb{P}(s_t=1|r_t)<\mathbb{P}(s_t=0|r_t)$, then, we can tag this time $t$. Otherwise, we do not do it. Equivalently, we can tag the time $t$ if \begin{flalign*} (1-p_0)f((r_t+\phi_0|y_{t-1}|)/\sigma_0)<p_0 f(r_t/\sigma_0), \end{flalign*} where $f(\cdot)$ is the density of $\mathcal{N}(0, 1)$. Using $\{\widehat{r}_t\}$, we can tag the time $t$ if \begin{flalign*} (1-\widehat{p}_n) f\big(\big(\widehat{r}_t+\widehat{\phi}_n|y_{t-1}|\big)/\widehat{\sigma}_n\big)< \widehat{p}_n f\left(\widehat{r}_t/\widehat{\sigma}_n\right). \end{flalign*}

A Null-Based Method for Bubble Tagging

The method described in Section (ref) relies on residuals from the one-step ahead recursion specified by model ((ref)) to tag the collapse of bubbles. In essence, it treats the explosive bubble alternative as the default and aims at detecting the null of no bubble as an anomaly. We shall here consider its complement which sets the null of no bubble as the baseline and detects the formation of a bubble as an anomaly. To be more specific, when $s_t = 0$ and there is no bubble at time $t$, we have $y_t = \varepsilon_t$ which forms a stationary white noise sequence. When the bubble starts to form at time $t$, however, an explosive drift $\phi_0|y_{t-1}|$ will be cumulatively added to the otherwise white noise sequence during the whole bubble period making the observed $y_t$ to cumulatively deviate away from the baseline. Therefore, it becomes natural to tag time $t$ as a bubble if $y_t > c$ for some threshold $c$. In contrast to the approach in Section (ref) which relies exclusively on model ((ref)) to compute the residuals $\{\widehat r_t\}$, this null-based method directly works on the original observations $\{y_t\}$ and can be less model dependent. In addition, since $y_t$ is distributed as a white noise sequence under the null of no bubble, the threshold $c$ can be taken as a uniform constant, which can be a convenient feature that facilitates the decision rule visualization. It can also be more advantageous in situations when bubbles are not prevailing in the observation period. Let $\{z_t\}$ be the auxiliary process defined in ((ref)), and we in the following provide some theoretical understanding of such a null-based bubble tagging method under the fixed horizon domain.

propositionFor any time $t$, the conditional probability that a $k$-th cumulative bubble will be correctly tagged by the null-based method equals to $\mathbb{P}(z_k > c)$, namely the marginal probability that the auxiliary explosive bubble process will exceed the same threshold.

For bubbles that persist for a growing horizon of time, by the explosive nature of the auxiliary bubble process it is expected that $\mathbb{P}(z_k > c) \to 1$ as $k \to \infty$ for any given threshold $c$, and as a result the aforementioned null-based bubble tagging method can identify such a persistent bubble with probability tending to one. Phillips2015 treated the bubble period as a fixed but unknown deterministic section of the whole observation time, and provided the consistency when the length of the bubble section grows proportionally with the sample size. In contrast, the current paper treats the bubble as an intrinsic feature of a stationary data generating mechanism, which serves as an important step to provide a statistical model to understand the mechanism of an economic phenomenon. We also remark that, unlike the QMLE discussed in Section (ref), the aforementioned null-based bubble tagging method and Proposition (ref) will continue to hold for situations when the hidden state process $\{s_t\}$ exhibits dependence and forms a stationary or nonstationary process by itself. For example, it can be a stationary Markov chain or a nonstationary Markov chain with time-varying transition matrices. In addition, the proof of Proposition (ref) can be readily generalized to handle bubble mechanisms other than the one-step autoregressive recursion specified in ((ref)).

Simulation Studies

To assess the performance of the QMLE of $\theta_0$ and Rules 1-4 in finite samples, we use the sample size $n=200$, 400, and 800, each with 1000 replications for model ((ref)). The error $\varepsilon_t$ follows

itemize$\mathcal{N}(0, 1)$; • the Laplace distribution with density \begin{flalign*} h(x)=\frac{1}{\sqrt{2}}\exp\big(-\sqrt{2}|x|\big),\quad x\in\mathbb{R}; \end{flalign*} • the standardized Student's $t_5$ $(\mathrm{st}_5)$ with density \begin{flalign*} h(x)=\frac{8}{3\pi\sqrt{3}}\Big(1+\frac{x^2}{3}\Big)^{-3},\quad x\in\mathbb{R}. \end{flalign*}

Three different true values of $\theta_0=(\phi_0, p_0, \sigma^2_0)'$ are used, respectively, i.e.,

itemize• Case I: $\theta_0=(1, 0.9, 1)'$; • Case II: $\theta_0=(\sqrt{10/9}, 0.9, 1)'$; • Case III: $\theta_0=(1.2, 0.9, 1)'$.

For Case I, $y_t$ is weakly stationary since $p_0\phi_0^2<1$, while $y_t$ is an infinite-variance process in Case III since $p_0\phi_0^2>1$. For Case II, $\theta_0$ is on the boundary, i.e., $p_0\phi_0^2=1$, which is never considered in the literature.

Table (ref) reports the bias, empirical standard deviation (ESD), and asymptotic standard deviation (ASD) of the QMLE $\widehat{\theta}_n$ for Cases I-III. Here, the ASD of $\theta_0$ is simulated by extra time series of length 10,000, and 2,000 replications are used to reduce the estimated bias. From the table, we can see that the QMLE performs well irrespective of infinite variance or heavy-tailedness issues. The biases are small and all the ESDs are close to the corresponding ASDs.

table[table omitted — 4,561 chars of source]

To see the overall approximation of the QMLE $\widehat{\phi}_n$, Fig. (ref) displays the histogram of $\sqrt{n}\big(\widehat{\phi}_n-\phi_0\big)$ when the sample size $n=400$. From the figure, we can see that $\sqrt{n}\big(\widehat{\phi}_n-\phi_0\big)$ is always asymptotically normal irrespective of infinite variance or heavy-tailedness of $y_t$.

figure[figure omitted — 579 chars of source]

We shall here examine the finite-sample performance of the two bubble tagging methods described in Section (ref). For the residual-based tagging method in Section (ref) with reference rules 1--4 we denote them by RBT$_1$--RBT$_4$ respectively in our numerical study, and we abbreviate the null-based tagging method in Section (ref) as NBT hereafter. For each generated process, let $\{\widehat s_t: 1\leq t\leq n\}$ be the estimated bubble tags and $\#$ denote the set cardinality. We consider the following evaluation metrics:

itemize• P: the overall proportion of correct tagging $\#\{t:\widehat s_t=s_t\}/n$; • P0: the proportion of correctly tagged null states $\#\{t:\widehat s_t=0 , s_t=0\}/\#\{t:s_t=0\}$; • P1: the proportion of correctly tagged bubbles $\#\{t:\widehat s_t=1 , s_t=1\}/\#\{t:s_t=1\}$.

The results are summarized in Tables (ref) and (ref) based on 1000 replications for each configuration. To provide a fair comparison, we set the thresholds of different tagging methods so their estimated bubble ratios $\#\{t:\widehat s_t=1\}/n$ are controlled at the same level. From Tables (ref) and (ref), we can observe the followings.

enumerate[(i)] • For both the RBT and NBT methods, the results are reasonably close across different error types. This indicates that the bubble tagging methods considered in Sections (ref) and (ref) possess a certain degree of robustness with respect to the error distribution. • For each of the method considered, the performance in general improves when the nonlinear autoregressive coefficient $\phi_0$ increases. This is mainly because a larger value of the parameter $\phi_0$ in general leads to a stronger degree of explosiveness during the bubble period, making it relatively easier to distinguish between bubbles and null-states. • When $p_0 = 0.9$ as in Table (ref), the performance of the RBT method can vary depending on which reference rule is used to obtain the threshold. The NBT method, on the other hand, seems to deliver a performance that is between the best and worst performed RBT methods. Note that the RBT method is deigned using the residuals that are more related to the bubble alternative, it meets with our intuition that the RBT method in general outperforms the NBT for most of the threshold choices when the bubble state probability $p_0 = 0.9$ is relatively high as in Table (ref). • When the true underlying bubble state probability $p_0$ decreases to 0.5 as in Table (ref), the bubble state no longer dominates and as a result the difference between the RBT and NBT methods becomes less noticeable and all the methods considered delivered quite similar performance.
table[table omitted — 2,589 chars of source]
table[table omitted — 2,589 chars of source]

An empirical example

In this section, we analyze the monthly Hang Seng Index (HSI) from December 1986 to December 2017 with a total of 373 observations. To eliminate the effect of inflation on price, we transform nominal prices into real prices by the consumer price index, which can be obtained from the Federal Reserve Bank of St Louis. Fig (ref) (a) displays the real HSI prices, from which one can see an ascendant linear trend in the time series.

figure[figure omitted — 194 chars of source]

Thus, we first subtract such a linear trend from the series. That is, we assume that the HSI real price $x_t$ is decomposed into $$x_t=b_0+b_1 t +y_t,$$ where $b_0+b_1t$ denotes the linear trend and $y_t$ follows a SNAR model. Note that $b_i,i=0,1$ can be seen as unknown parameters and can be estimated jointly. Their estimates are $\widehat{b}_0=23.661$ and $\widehat{b}_1= 0.372$, respectively. The linear time trend is plotted in Fig. (ref) (a) by the dotted line and $\{y_t\}$ in Fig. (ref) (b). The estimates with standard deviations (SDs) of the SNAR model $\{y_t\}$ are reported in Table (ref).

table[table omitted — 409 chars of source]

All estimates are statistically significant at the $5\%$ level since their corresponding $p$-values are extremely small which are thus not reported in the table. The estimate of $\phi_0$ is large than one, and its $95\%$ confidence interval is $(1.005, 1.047)$, conforming to the locally explosive behavior of series $y_t$. For the fitting adequacy, we calculate the $p$-values of the test statistic $Q_M$ with $M=6, 12, 18$, and $24$ when the tuning parameter $a$ is the 90% or 95% quantile of $\{|y_t|,t=1\dots,n\}$, respectively. The results are summarized in Table (ref), which implies that the fitting is adequate at the $5\%$ level.

table[table omitted — 342 chars of source]

We then apply the bubble tagging methods described in Section (ref) to label each time point as either being in a bubble state or being in the null. Since the estimated bubble probability $\widehat p_0 = 0.977$ from Table (ref) which is very high, in view of the simulation results in Section (ref), we shall here consider using the residual-based method in Section (ref) to tag the collapses of bubbles for the series $\{y_t\}$. In particular, Fig. (ref) displays the selected dates of $\widehat s_t=0$ under Rules 1--4. It can be seen from Fig. (ref) that the tagging times can vary based on which Rule is used, but several important dates are identified simultaneously by at least two rules. Table (ref) summarizes such these dates,

table[table omitted — 347 chars of source]

which coinside with historical financial crises, i.e., the depression started from the Black Monday in 1987, the Asian financial crises in 1997, the global financial turmoil caused by the subprime crisis over 2007-2009, and the Hong Kong stock market plummeting in 2016.

figure[figure omitted — 178 chars of source]

Although the collapse of a bubble can be dated by $\widehat s_t=0$, the emergence and exuberance of a bubble can not be asserted by $\widehat s_t=1$ immediately. After all, a short-period deviation of the price is reasonable due to the market fluctuations. Of course, a short-period deviation might be regarded as a small bubble in some sense, which bursts quickly by the market adjustment, thus we could pay little attention and ignore them afterwards. What we really need to worry about is the bubble that can trigger tremendous harm, which emerges as the accumulation of long-lasting excursion. Specifically, if $s_t=0, s_{t+1}=s_{t+2}=\dots=s_{t+k-1}=1, s_{t+k}=0$, then we call it an excursive period that starts from $t+1$ and ends at $t+k$, and define its duration as $k$. Within an excursive period, the presence of a bubble can be supported if the duration exceeds some time span, for example, one or two years. For our application, the time span is set to be 18 months. Table (ref) summarizes the periods whose durations exceed 18 months, as well as their start and end dates.

table[table omitted — 446 chars of source]

Fig. (ref) plots those periods by gray shadows.

figure[figure omitted — 258 chars of source]

We can see explosive behaviors in most of the periods, indicating the presence and accumulation of bubbles. By Proposition (ref) in Section (ref), the RBT method is capable of detecting the collapse of an accumulated bubble when its duration $k \to \infty$; see also the same proposition for a probabilistic bound with a finite duration. Another finding is that the magnitude of a bubble is larger as the period lasts longer possibly, for example, the one reaches a value of 210 in October 2007, corresponding to the period from April 2003 to November 2007 with the duration of 56 months. Investors should be alert to such a long-time excursion along with the potential of disastrous bubbles. In the periods where the bubble lasts over 24 months (plotted by the shadow with red backslash in Fig. (ref)), one should be aware of the false boom in financial markets, and adjust asset allocation to hedge the risk of a potential bubble burst.

Conclusions

The paper has introduced a novel stochastic nonlinear autoregressive model to describe the dynamics of economic or financial bubbles within a causal and stationary framework, and discussed its strict stationarity and geometric ergodicity. The paper has further studied the quasi-maximum likelihood estimation of the model and established the asymptotics under minimal assumptions on innovation. Due to the unobservability of the latent variable $s_t$ and the resulting unavailability of the residuals, a new model diagnostic checking tool has been proposed for the adequacy of the fitting. Finally, the paper considers two approaches, one from the residual perspective and the other from the null perspective, for bubble tagging.

Although our new model is useful, the model assumption on the independence between $\{\varepsilon_t\}$ and $\{s_t\}$ seems a little bit stronger from the perspective of empirical pragmatism. To obtain more reasonable interpretation or approximation of the bubble, such an independence assumption can be relaxed. For instance, we can assume that $s_t$ depends on the history of the observed process. Specifically, we can let $\mathbb{P}(s_t=1|\mathcal{F}_{t-1})=g(\bm\beta'\mathbf{y}_{t-1})$, where $\mathcal{F}_{t-1}=\sigma(y_{t-j}: j\geq 1)$ be a sigma-field, $\mathbf{y}_{t-1}=(1, y_{t-1},...,y_{t-q})'$, and $g$ is a measurable function (e.g. a logistic function). Furthermore, we can also restrict the form of $s_t$ in macroeconomic time series analysis and let $s_t=I(\bm\beta'\mathbf{x}_t>c)$, where $\mathbf{x}_t$ may contain many exogenous macroeconomic variables or indexes and $c$ is a threshold parameter. In addition, it is possible to consider the situation when the hidden state process $\{s_t\}$ exhibits temporal dependence and forms a Markov chain. In this case, the null-based bubble tagging method in Section (ref) can be more advantageous when bubbles occur in separated but persistent clusters. Another potential topic is to study multivariate stochastic nonlinear AR models. We leave these topics for future research.