Dynamic tail inference with log-Laplace volatility

TL;DR

This paper introduces a novel stochastic volatility model with log-Laplace distribution enabling dynamic estimation of extreme event probabilities in heavy-tailed, nonlinear time series, with applications in finance.

Contribution

The paper develops a new model with analytic conditional structure for tail inference, and provides a straightforward estimation method for dynamic extreme event probabilities.

Findings

Effective in simulated nonlinear series with chaotic volatility

Accurate predictive tail inference demonstrated in financial data

Computationally efficient estimation procedure

Abstract

We propose a family of models that enable predictive estimation of time-varying extreme event probabilities in heavy-tailed and nonlinearly dependent time series. The models are a white noise process with conditionally log-Laplace stochastic volatility. In contrast to other, similar stochastic volatility formalisms, this process has analytic expressions for its conditional probabilistic structure that enable straightforward estimation of dynamically changing extreme event probabilities. The process and volatility are conditionally Pareto-tailed, with tail exponent given by the reciprocal of the log-volatility's mean absolute innovation. This formalism can accommodate a wide variety of nonlinear dependence, as well as conditional power law-tail behavior ranging from weakly non-Gaussian to Cauchy-like tails. We provide a computationally straightforward estimation procedure that uses an asymptotic approximation of the process' dynamic large deviation probabilities. We demonstrate the estimator's utility with a simulation study. We then show the method's predictive capabilities on a simulated nonlinear time series where the volatility is driven by the chaotic Lorenz system. Lastly we provide an empirical application, which shows that this simple modeling method can be effectively used for dynamic and predictive tail inference in financial time series.

Tables

Table 1. Table 1: Estimation results using ( 3.5 )’s Δ ^ ​ ( Λ ) ^ Δ Λ \widehat{\Delta}\left(\Lambda\right) , averaged over 1000 simulations of T 𝑇 T samples of ( 1.1 )-( 1.4 ) with H t subscript 𝐻 𝑡 H_{t} given by ( 4.1 ) and ( 4.2 ). Model for ( 3.3 )’s H ¯ ^ t subscript ^ ¯ 𝐻 𝑡 \widehat{\bar{H}}_{t} is AR(10).

$H_t=$(4.1)
$T=625-10$
$\Delta$	$.05$	$.10$	$.15$	$.20$	$.25$	$.30$	$.35$	$.40$	$.45$	$.50$
avg. $\hat{\Delta }\left(2{\hat{\sigma }}_{\epsilon }\right)$	.28	.27	.28	.31	.35	.39	.43	.47	.50	.55
avg. $\hat{\Delta }\left(3{\hat{\sigma }}_{\epsilon }\right)$	.14	.16	.20	.26	.32	.37	.41	.47	.51	.56
avg. $\hat{\Delta }\left(4{\hat{\sigma }}_{\epsilon }\right)$	.10	.13	.19	.26	.31	.36	.42	.47	.51	.56
std. dev. $\hat{\Delta }\left(2{\hat{\sigma }}_{\epsilon }\right)$	.02	.03	.04	.05	.06	.07	.08	.09	.11	.13
std. dev. $\hat{\Delta }\left(3{\hat{\sigma }}_{\epsilon }\right)$	.02	.03	.06	.07	.09	.10	.11	.13	.14	.15
std. dev. $\hat{\Delta }\left(4{\hat{\sigma }}_{\epsilon }\right)$	.03	.05	.07	.08	.09	.11	.12	.13	.14	.16
$T=1250-10$
$\Delta$	$.05$	$.10$	$.15$	$.20$	$.25$	$.30$	$.35$	$.40$	$.45$	$.50$
avg. $\hat{\Delta }\left(2{\hat{\sigma }}_{\epsilon }\right)$	.26	.25	.27	.30	.34	.38	.41	.45	.50	.54
avg. $\hat{\Delta }\left(3{\hat{\sigma }}_{\epsilon }\right)$	.13	.15	.21	.27	.32	.37	.41	.46	.51	.55
avg. $\hat{\Delta }\left(4{\hat{\sigma }}_{\epsilon }\right)$	.09	.14	.21	.26	.32	.38	.42	.47	.52	.56
std. dev. $\hat{\Delta }\left(2{\hat{\sigma }}_{\epsilon }\right)$	.01	.02	.02	.03	.04	.05	.06	.08	.09	.11
std. dev. $\hat{\Delta }\left(3{\hat{\sigma }}_{\epsilon }\right)$	.01	.03	.06	.07	.08	.09	.10	.11	.12	.14
std. dev. $\hat{\Delta }\left(4{\hat{\sigma }}_{\epsilon }\right)$	.02	.05	.06	.07	.08	.09	.10	.12	.14	.14
$H_t=$(4.2)
$T=625-10$
$\Delta$	$.05$	$.10$	$.15$	$.20$	$.25$	$.30$	$.35$	$.40$	$.45$	$.50$
avg. $\hat{\Delta }\left(2{\hat{\sigma }}_{\epsilon }\right)$	.29	.27	.26	.26	.28	.32	.36	.41	.46	.52
avg. $\hat{\Delta }\left(3{\hat{\sigma }}_{\epsilon }\right)$	.15	.14	.17	.21	.26	.32	.38	.43	.48	.52
avg. $\hat{\Delta }\left(4{\hat{\sigma }}_{\epsilon }\right)$	.10	.11	.15	.20	.26	.32	.37	.42	.48	.53
std. dev. $\hat{\Delta }\left(2{\hat{\sigma }}_{\epsilon }\right)$	.02	.02	.03	.03	.05	.07	.09	.10	.12	.12
std. dev. $\hat{\Delta }\left(3{\hat{\sigma }}_{\epsilon }\right)$	.02	.02	.04	.06	.08	.08	.09	.10	.10	.12
std. dev. $\hat{\Delta }\left(4{\hat{\sigma }}_{\epsilon }\right)$	.02	.04	.06	.07	.08	.08	.08	.09	.10	.11
$T=1250-10$
$\Delta$	$.05$	$.10$	$.15$	$.20$	$.25$	$.30$	$.35$	$.40$	$.45$	$.50$
avg. $\hat{\Delta }\left(2{\hat{\sigma }}_{\epsilon }\right)$	.27	.26	.25	.25	.27	.32	.36	.41	.46	.52
avg. $\hat{\Delta }\left(3{\hat{\sigma }}_{\epsilon }\right)$	.13	.13	.17	.23	.28	.33	.38	.43	.48	.53
avg. $\hat{\Delta }\left(4{\hat{\sigma }}_{\epsilon }\right)$	.09	.11	.16	.22	.28	.33	.38	.43	.48	.53
std. dev. $\hat{\Delta }\left(2{\hat{\sigma }}_{\epsilon }\right)$	.01	.02	.02	.03	.05	.06	.08	.09	.10	.11
std. dev. $\hat{\Delta }\left(3{\hat{\sigma }}_{\epsilon }\right)$	.01	.02	.05	.06	.07	.07	.08	.08	.09	.10
std. dev. $\hat{\Delta }\left(4{\hat{\sigma }}_{\epsilon }\right)$	.02	.04	.05	.06	.06	.07	.07	.08	.08	.10

Table 2. Table 2: Lorenz-Driven Volatility Testing Results: Averages ( ± plus-or-minus \pm Standard Deviation)

$N_{\text{Train}}/N_{\text{Test}}$	avg. $\hat{\Delta }\left(4{\hat{\sigma }}_{\epsilon }\right)$	avg. ${\rho }_{\left|\epsilon \right|,\hat{\sigma }}$	avg. Sn.	avg. Sp.
2,994/6,987 (30/70)	.359 ($\pm .008$)	.604 ($\pm .001$)	.955 ($\pm .002$)	.745 ($\pm .008$)
3,992/5,989 (40/60)	.345 ($\pm .006$)	.614 ($\pm .000$)	.934 ($\pm .003$)	.781 ($\pm .005$)
4,990/4,991 (50/50)	.347 ($\pm .005$)	.608 ($\pm .000$)	.958 ($\pm .003$)	.736 ($\pm .006$)

Table 3. Table 3: SPX Backtesting Results: Averages ( ± plus-or-minus \pm Standard Deviation)

2014 to 2019
$N_{\text{Train}}/N_{\text{Test}}$	avg. $\hat{\Delta }\left(4{\hat{\sigma }}_{\epsilon }\right)$	avg. ${\rho }_{\left|\epsilon \right|,\hat{\sigma }}$	avg. Sn.	avg. Sp.
625/625 (50/50)	.33 ($\pm .03$)	.39 ($\pm .02$)	.78 ($\pm .00$)	.74 ($\pm .04$)
750/500 (60/40)	.28 ($\pm .02$)	.47 ($\pm .01$)	.85 ($\pm .05$)	.76 ($\pm .02$)
1990 to 2019
$N_{\text{Train}}/N_{\text{Test}}$	avg. $\hat{\Delta }\left(4{\hat{\sigma }}_{\epsilon }\right)$	avg. ${\rho }_{\left|\epsilon \right|,\hat{\sigma }}$	avg. Sn.	avg. Sp.
3,666/3,666 (50/50)	.206 ($\pm .005$)	.575 ($\pm .000$)	.893 ($\pm .010$)	.805 ($\pm .003$)
4,399/2,933 (60/40)	.199 ($\pm .003$)	.584 ($\pm .000$)	.913 ($\pm .000$)	.733 ($\pm .003$)

Equations

$${\epsilon }_t={\sigma }_t{z}_t,$$

$${\sigma }_t=\exp \left(H_t\right),$$

$$H_t=E\left\{H_t|{\mathcal{F}}_{t-1}\right\}+h_t,$$

$$p_h\left(h_t\right)=\frac{1}{2\Delta }\exp \left(-\frac{\left|h_t\right|}{\Delta }\right),$$

$$\Delta =E\left\{\left|h_t\right|\right\}.$$

$$P\left\{\left|{\epsilon }_t\right|\ge \Lambda |{\mathcal{F}}_{t-1}\right\}\sim f(\Delta )\exp \left(\frac{E\left\{H_t|{\mathcal{F}}_{t-1}\right\}}{\Delta }\right){\Lambda }^{-1/\Delta }$$

$$p_{\sigma }({\sigma }_t|{\mathcal{F}}_{t-1})=\{\begin{array}{c}\frac{1}{2\Delta }\exp \left(-\frac{{\bar{H}}_t}{\Delta }\right){\sigma }_t^{1/\Delta -1};0<{\sigma }_t<\exp \left({\bar{H}}_t\right)\text{par}\\ \frac{1}{2\Delta }\exp \left(\frac{{\bar{H}}_t}{\Delta }\right){\sigma }_t^{-1/\Delta -1};{\sigma }_t\ge \exp \left({\bar{H}}_t\right)\text{par}\end{array}$$

$$E\left\{{\epsilon }_t^n|{\mathcal{F}}_{t-1}\right\}=\exp \left(n{\bar{H}}_t\right)\frac{(n-1)!!}{1-n^2{\Delta }^2}.$$

$$\kappa (\Delta )=3\frac{{\left(1-4{\Delta }^2\right)}^2}{1-16{\Delta }^2},$$

$$p_{\epsilon }\left({\epsilon }_t|{\mathcal{F}}_{t-1}\right)=\frac{1}{4\sqrt{\pi }\Delta }{\left(\sqrt{2}{e}^{{\bar{H}}_t}\right)}^{-1/\Delta }\Gamma \left(\frac{1-1/\Delta }{2},\frac{{\epsilon }_t^2}{2e^{2{\bar{H}}_t}}\right){\left|{\epsilon }_t\right|}^{1/\Delta -1}$$

$$+\frac{1}{4\sqrt{\pi }\Delta }{\left(\sqrt{2}{e}^{{\bar{H}}_t}\right)}^{1/\Delta }\gamma \left(\frac{1+1/\Delta }{2},\frac{{\epsilon }_t^2}{2e^{2{\bar{H}}_t}}\right){\left|{\epsilon }_t\right|}^{-1/\Delta -1}$$

$$\Gamma (a,b)={\int }_b^{\infty }{x}^{a-1}{e}^{-x}dx$$

$$\gamma (a,b)={\int }_0^b{x}^{a-1}{e}^{-x}dx$$

$$p_{\epsilon }\left({\epsilon }_t|{\mathcal{F}}_{t-1}\right)\sim \frac{1}{\sqrt{2\pi }{e}^{{\bar{H}}_t}}\frac{1}{1-{\Delta }^2}.$$

$${\tilde{\Lambda }}_t=\frac{\Lambda }{\sqrt{2}{e}^{{\bar{H}}_t}},$$

$$P\left\{\left|{\epsilon }_t\right|\ge \Lambda |{\mathcal{F}}_{t-1}\right\}=\frac{1}{2\sqrt{\pi }}\left(\gamma \left(\frac{1+1/\Delta }{2},{\tilde{\Lambda }}_t^2\right){\tilde{\Lambda }}_t^{-1/\Delta }-\Gamma \left(\frac{1-1/\Delta }{2},{\tilde{\Lambda }}_t^2\right){\tilde{\Lambda }}_t^{1/\Delta }\right)+\text{erfc}\left({\tilde{\Lambda }}_t\right),$$

$$\text{erfc}(x)=1-\text{erf}(x)=1-\frac{2}{\sqrt{\pi }}{\int }_0^x{e}^{-u^2}du=\frac{2}{\sqrt{\pi }}{\int }_x^{\infty }{e}^{-u^2}du$$

$$P\left\{\left|{\epsilon }_t\right|\ge \Lambda |{\mathcal{F}}_{t-1}\right\}=\frac{1}{2\sqrt{\pi }}\Gamma \left(\frac{1+1/\Delta }{2}\right){\tilde{\Lambda }}_t^{-1/\Delta }+O\left({\tilde{\Lambda }}_t^{-5}{e}^{-{\tilde{\Lambda }}_t^2}\right)$$

$$P\left\{{\sigma }_t\ge \Lambda |{\mathcal{F}}_{t-1}\right\}=\frac{1}{2}\exp \left(\frac{{\bar{H}}_t}{\Delta }\right){\Lambda }^{-1/\Delta }.$$

$$P\left\{\left|{\epsilon }_t\right|\ge \Lambda |{\mathcal{F}}_{t-1}\right\}\sim \frac{{\sqrt{2}}^{1/\Delta }}{\sqrt{\pi }}\Gamma \left(\frac{1+1/\Delta }{2}\right)P\left\{{\sigma }_t\ge \Lambda |{\mathcal{F}}_{t-1}\right\}.$$

$$P\left\{\left|{\epsilon }_t\right|\ge \Lambda |{\mathcal{F}}_{t-1}\right\}\sim P\left\{{\sigma }_t\ge \Lambda |{\mathcal{F}}_{t-1}\right\}$$

$$P\left\{\left|{\epsilon }_t\right|\ge \Lambda |{\mathcal{F}}_{t-1}\right\}=\frac{1}{2\sqrt{\pi }}\Gamma \left(\frac{1+1/\Delta }{2}\right){\tilde{\Lambda }}_t^{-1/\Delta }-\frac{1}{4\sqrt{\pi }{\Delta }^2}\Gamma \left(-\frac{3}{2},{\tilde{\Lambda }}_t^2\right)$$

$$+\frac{1}{2\sqrt{\pi }{\Delta }^2}\Gamma \left(-\frac{5}{2},{\tilde{\Lambda }}_t^2\right)+O\left({\tilde{\Lambda }}_t^{-9}{e}^{-{\tilde{\Lambda }}_t^2}\right).$$

$$E\left\{H_t|\log \left|{\epsilon }_t\right|\right\}=\log \left|{\epsilon }_t\right|+\frac{\log 2+\gamma }{2}\approx \log \left|{\epsilon }_t\right|+.6352,$$

$$\gamma =\underset{n\to \infty }{\lim }\left(\sum _{k=1}^n\frac{1}{k}-\log n\right)\approx .5772$$

$${\hat{H}}_t=\log \left|{\epsilon }_t\right|+\frac{\log 2+\gamma }{2}\approx \log \left|{\epsilon }_t\right|+.6352.$$

$${\hat{\bar{H}}}_t=\hat{E}\left\{H_t|{\mathcal{F}}_{t-1}\right\}=m\left({\hat{H}}_{t-1},...,{\hat{H}}_{t-p},{\vec{X}}_{t-1},...,{\vec{X}}_{t-q},...\right).$$

$$\frac{1}{T}\sum _{t=1}^T\text{mathds}1\left(\left|{\epsilon }_t\right|\ge \Lambda \right)=\frac{1}{T}\sum _{t=1}^{T}P\left\{\left|{\epsilon }_t\right|\ge \Lambda |{\hat{\bar{H}}}_t,\hat{\Delta }\right\},$$

$$\hat{\Delta }\left(\Lambda \right)={\text{argmin}}_{\Delta }\left\{\left|\sum _{t=1}^T\left(\text{mathds}1\left(\left|{\epsilon }_t\right|\ge \Lambda \right)-\frac{{\sqrt{2}}^{1/\Delta }}{2\sqrt{\pi }}\Gamma \left(\frac{1+1/\Delta }{2}\right)\exp \left(\frac{{\hat{\bar{H}}}_t}{\Delta }\right){\Lambda }^{-1/\Delta }\right)\right|\right\}$$

$$H_t=.5H_{t-1}+.4H_{t-2}+h_t$$

$$H_t=.05H_{t-1}+.05H_{t-2}+.25H_{t-3}+.2H_{t-4}+.35H_{t-5}+h_t.$$

Full text

Keywords: Stochastic volatility, Heavy tails, Extreme events, Nonlinear time series, Tail risk

AMS 2010 Subject Classification: 60G70, 62G32, 62M10, 62P05, 62P12, 91G70

1 Introduction

Financial time series data is well-known to exhibit nonlinear dependence and “fat tails”. Such time series often display trends of increasing or decreasing volatility, along with a propensity for extreme fluctuations that is far greater than what would be predicted from a Gaussian or other distribution with finite polynomial moments. The latter observation was probably most famously addressed with the early Pareto-tailed and stable models for financial time series proposed by Mandelbrot (1963) and Fama (1968), while the most well-known early approaches to the nonlinear dependence problem were given by Engle (1982) and Bollerslev (1986) with original and generalized autoregressive conditional heteroskedasticity (ARCH, GARCH) models. Since then a great deal of research has been dedicated to extreme event probability estimation and nonlinear time series modeling for financial applications (e.g., Embrechts et al. 2011 and Terasvirta et al. 2010).

Estimation of extreme event probabilities is a very important problem for risk and portfolio management. Many estimators for the tail exponent of a marginal distribution have been proposed, e.g., by Hill (1975), Pickands (1975), and deHaan and Resnick (1980), however, these estimators are very sensitive to dependence in the data (Kearns and Pagan 1997, Diebold et al. 2000). This makes them often ill-suited for application to many strongly dependent time series of interest for financial modeling, i.e., the squares and moduli of financial log-returns (Embrechts et al. 2011 p. 270, 406). Relatedly, these estimators along with much of extreme value theory (EVT) are designed for inference of stationary, rather than time-varying, tail behavior. Gardes and Girard (2008) and Gardes and Stupfler (2014) have given nonparametric estimators for time-varying tail exponents, while Kelly (2014) has given a parametric approach to dynamic power law estimation in financial time series. The parametric modeling method we propose here, however, does not assume a time-varying tail exponent. Our approach is hence closer to McNeil and Frey’s (2000) combination of stationary EVT with GARCH modeling. However, we use a novel, stochastic volatility approach to enable dynamic tail inference.

A canonical form of stochastic volatility model, first proposed by Taylor (1982, 1986), is given by

[TABLE]

where $z_{t}$ is an i.i.d. process with a mean of 0 and a variance of 1, and $\sigma_{t}$ is a non-negative process defined by

[TABLE]

where $H_{t}$ is a Gaussian process with mean $\mu_{H}$ and variance $\sigma_{H}^{2}<\infty$. An important example was the AR(1) model $H_{t}=\mu_{H}+\beta\left(H_{t-1}-\mu_{H}\right)+h_{t}$, where $\beta\in\mathbb{R}$ and $h_{t}\sim\mathcal{N}\left(0,\sigma_{h}^{2}\right)$ is i.i.d. The kind of model in (1.1) and (1.2), where $H_{t}$ is a variety of Gaussian processes, has been extensively applied and studied, e.g., to exchange rate modeling by Harvey et al. (1994) and extended to long-memory Gaussian $H_{t}$ by Breidt et al. (1998). However, the stochastic volatility $\sigma_{t}$ in (1.2) follows a log-normal distribution, which has finite polynomial moments (see Johnson et al. 1994). As a consequence, (1.1) also has bounded moments in the usual case of Gaussian $z_{t}$. This can be problematic for modeling time series with power law tails and divergent higher-order polynomial moments. Practically, such models will underestimate the probabilities of extreme events. To remedy this, $z_{t}$ has often been chosen to follow the heavier-tailed Student’s $t$-distribution, e.g., in Harvey et al. (1994), Liesenfeld and Jung (2000), and Chib et al. (2002). However, with a Gaussian or Student’s $t$ choice for $z_{t}$ and Gaussian $H_{t}$, the model (1.1) does not have convenient analytic expressions for its conditional probabilistic structure. This makes estimation of model parameters as well as outcome probabilities difficult, often requiring the use of Bayesian and Monte Carlo methods such as those described in Jacquier et al. (1994), Kim et al. (1998), Sandmann and Koopman (1998), or Chib et al. (2002). Reliable estimation of similar, conditionally Student’s t-distributed, ARCH-related models requires similar numerical procedures (see, e.g., Mousazadeh and Karimi 2007, Ardia 2008, Ardia and Hoogerheide 2010).

In this paper we propose a stochastic volatility formalism that enables straightforward and effective estimation of time-varying extreme event probabilities in time series with power law tail behavior. The models we present have the form (1.1)-(1.2), with $z_{t}\sim\mathcal{N}(0,1)$. However, instead of defining (1.2)’s $H_{t}$ as Gaussian, we make $H_{t}$ conditionally Laplace-distributed, defining $H_{t}$ as

[TABLE]

where $\mathcal{F}_{t-1}$ is the filtration up to time $t-1$ and $h_{t}$ is i.i.d. Laplace-distributed with density

[TABLE]

where

[TABLE]

This simple but important adjustment endows (1.1)-(1.2)’s $\sigma_{t}$ and $\varepsilon_{t}$ with power law-tailed conditional probability distributions, for which there are natural and convenient analytic expressions. These conditional distributions give the result

[TABLE]

as $\Lambda\rightarrow\infty$. Hence $H_{t}$’s mean absolute innovation $\Delta$ specifies the tail exponent, which allows (1.1)-(1.4) to flexibly define processes ranging from only mildly non-Gaussian with $\Delta\approx 0$ to processes with Cauchy-like tails at $\Delta=1$. The process $\varepsilon_{t}$’s tail probabilities are also strongly and explicitly dependent on the process $H_{t}$, which enables their dynamic estimation. We give a simple, probabilistic method-of-moments estimation procedure for the tail exponent, which takes advantage of the result (1.5) for $\varepsilon_{t}$’s tail probabilities. We present a simulation study to show the effectiveness of this estimator. We next demonstrate the predictive capabilities of this methodology with a simulated nonlinear time series, where the volatility is driven by the chaotic Lorenz system. We then give an empirical application to S&P 500 Index data, which shows that this modeling method can be used for dynamic and predictive estimation of volatility and extreme event probabilities in heavy-tailed financial time series. We give some concluding remarks, an appendix with proofs of the main results, and a second appendix, which shows that very similar results, including (1.5), can also be derived for the case when $z_{t}$ is Laplace-distributed.

2 Model

We will denote $H_{t}$’s conditional expectation as $\bar{H}_{t}=E\left\{H_{t}|\mathcal{F}_{t-1}\right\}$. We first give expressions for $\sigma_{t}$’s conditional probability density and $\varepsilon_{t}$’s conditional polynomial moments.

Lemma 2.1**.**

Given (1.3)-(1.4), $\sigma_{t}$ in (1.2) is conditionally distributed according to the log-Laplace probability density function

[TABLE]

Corollary 2.1**.**

Given (1.1) and (2.1) with $z_{t}\sim\mathcal{N}(0,1)$, for all even $n\geq 2$, if $\Delta\geq 1/n$, then $E\left\{\varepsilon_{t}^{n}|\mathcal{F}_{t-1}\right\}=\infty$, while if $\Delta<1/n$,

[TABLE]

Several graphs of (2.1) are given in Fig. 3. It is clear from Corollary 2.2 that the conditional kurtosis $\kappa_{t}=E\left\{\varepsilon_{t}^{4}|\mathcal{F}_{t-1}\right\}/\left(E\left\{\varepsilon_{t}^{2}|\mathcal{F}_{t-1}\right\}\right)^{2}$ diverges or is not well-defined for $\Delta\geq 1/4$. However, when (2.2) is used to calculate $\kappa_{t}$ for $\Delta<1/4$, it can be easily shown that the conditional kurtosis has the time-independent definition

[TABLE]

which is minimized at $\kappa(0)=3$, the Gaussian kurtosis. The divergent algebraic structure of (2.3) contrasts with the kurtosis’ exponential structure for Gaussian $H_{t}$, which would be written as $\kappa_{t}=3\exp\left(8\Delta^{2}\right)$. The excess kurtosis $\kappa^{\text{exc.}}(\Delta)=\kappa(\Delta)-3$ is graphed in Fig. 2 for both Laplace and Gaussian $H_{t}$. Now that we have given $\varepsilon_{t}$’s conditional moment structure, we will give its conditional probability density function.

Theorem 2.1**.**

Given (2.1) with $z_{t}\sim\mathcal{N}(0,1)$, (1.1)’s $\varepsilon_{t}$ is conditionally distributed according to the probability density function

[TABLE]

where

[TABLE]

is the upper incomplete gamma function and

[TABLE]

is the lower incomplete gamma function.

Several graphs of (2.4) are given in Fig. 3. Note that with larger $\Delta$ and smaller $\bar{H}_{t}$ the densities are more sharply peaked around the origin. This can also be seen from the following result.

Corollary 2.2**.**

For $\Delta<1$, as $\left|\varepsilon_{t}\right|\rightarrow 0$,

[TABLE]

The above result characterizes $\varepsilon_{t}$’s near-mean behavior. Next we give an analytic expression for $\left|\varepsilon_{t}\right|$’s exceedance probabilities, before characterizing $\varepsilon_{t}$’s tail behavior. We first make the notation

[TABLE]

which will make several of the next results more symbolically concise.

Theorem 2.2**.**

By (2.4)-(2.6),

[TABLE]

where

[TABLE]

is the complementary error function.

Note from the definition (2.6) that as $\widetilde{\Lambda}_{t}\rightarrow\infty$, (2.9)’s first term will tend to a power law multiplied by the gamma function and other constants. By (2.5) and (2.10) however, the second and third term of (2.9) will asymptotically become much smaller. This leads to the next result, which will be useful for straightforwardly approximating $\varepsilon_{t}$’s dynamic extreme event probabilities.

Theorem 2.3**.**

By (2.4),

[TABLE]

as $\widetilde{\Lambda}_{t}\rightarrow\infty$.

This result reduces to the asymptotic power law in (1.5). We will make extensive use of Theorem 2.3 in our estimation procedure described in the next section. The result arises from asymptotic expansions of (2.4) and its integral, which are given in (A.19)-(A.21) and (A.26) respectively. We note from integrating (2.1) with $\Lambda\geq\exp\left(\bar{H}_{t}\right)$ that

[TABLE]

Comparison of (2.12) with (2.11) after re-substituting for (2.8) gives

[TABLE]

This result shows that, asymptotically, $\varepsilon_{t}$’s large deviation probabilities are simply proportional to $\sigma_{t}$’s corresponding probabilities through a $\Delta$-dependent factor, and hence the probabilities’ dependence on $\bar{H}_{t}$ is given straightforwardly by $\exp\left(\bar{H}_{t}/\Delta\right)$. We make the additional note that $\Gamma(3/2)=\sqrt{\pi}/2$ and $\Gamma(1/2)=\sqrt{\pi}$, therefore (2.13) reduces to simply

[TABLE]

for $\Delta=1/2$ and as $\Delta\rightarrow\infty$.

It can be seen in (2.11) and (A.26) that more complex dependence on $\bar{H}_{t}$ and $\Delta$ is captured by higher-order terms with comparatively small contributions at high $\Lambda$ relative to $e^{\bar{H}_{t}}$. It is informative to write out the first few terms of (A.26)’s asymptotic expansion with (2.8), giving

[TABLE]

Note that the first order correction is negative, indicating that the approximation (2.11) has a small positive bias at large $\widetilde{\Lambda}_{t}$. Also note that the correction terms are inversely proportional to $\Delta^{2}$, meaning that (2.11)’s approximation is more accurate for heavier-tailed processes with larger $\Delta$. Now that we have presented the model’s polynomial moments and probabilistic structure, we will describe an estimation procedure using some of the above results.

3 Estimation

We begin by considering the conditional expectation of the unobservable quantity $H_{t}$, given the value of the observable $\log\left|\varepsilon_{t}\right|$.

Proposition 3.1**.**

Given (1.1) and (1.2),

[TABLE]

where

[TABLE]

is the Euler-Mascheroni constant.

By Proposition 3.1, an unbiased estimator for $H_{t}$ is given by

[TABLE]

An arbitrary regression model $m(.)$ for estimation of $E\left\{H_{t}|\mathcal{F}_{t-1}\right\}$ may then be trained using the $\widehat{H}_{t}$’s and any other relevant variables $\vec{X}_{t-1},...,\vec{X}_{t-q}$, giving

[TABLE]

We then estimate $\Delta$ by finding an approximate solution to the empirical, probabilistic moment constraint

[TABLE]

The constraint (3.4) is justified by the law of total expectation, which requires that $E\left\{\mathds{1}\left(\left|\varepsilon_{t}\right|\geq\Lambda\right)\right\}=E\left\{E\left\{\mathds{1}\left(\left|\varepsilon_{t}\right|\geq\Lambda\right)|\mathcal{F}_{t-1}\right\}\right\}$. This estimation method requires a result for the probability estimate on (3.4)’s right-hand side. Using (2.9)’s full result is less computationally straightforward because it is defined in terms of the highly nonstandard incomplete gamma functions, which often creates numerical difficulties on realistic datasets. We hence use Theorem 2.3’s result (2.11) instead, which gives an accurate approximation for sufficiently large $\Lambda/\sqrt{2}e^{\widehat{\bar{H}}_{t}}$. We therefore use the estimator

[TABLE]

In the simulation and empirical sections below, we search for (3.5)’s $\widehat{\Delta}$ over the range $[.01,1]$ with a precision of .01.

4 Simulation: $\Delta$ Estimation

In this section we present results of the estimation procedure described in the previous section for simulations of the process (1.1)-(1.4). We run 1000 simulations each, with sample sizes of 625 and 1250, for $\Delta$ ranging from .05 to .50. The model used for (3.3)’s $\widehat{\bar{H}}_{t}$ is a linear autoregressive (AR) model that uses the 10 previous values of (3.2)’s $\widehat{H}_{t}$. This AR(10) model is calibrated using the Yule-Walker method.

We present simulation results for two different, simple processes for $H_{t}$. The first is given by the AR(2) process

[TABLE]

and the second is given by the AR(5) process

[TABLE]

The results are presented in Table 1, where we give the averages and standard deviations of $\widehat{\Delta}$ over the 1000 simulations. For $\Lambda$, we choose $2\widehat{\sigma}_{\varepsilon}$, $3\widehat{\sigma}_{\varepsilon}$, and $4\widehat{\sigma}_{\varepsilon}$, where $\widehat{\sigma}_{\varepsilon}^{2}$ is $\varepsilon_{t}$’s sample variance. This gives a set of data-driven $\Lambda$ values for which (2.11) is an increasingly accurate approximation of (3.4)’s integral.

One can see in Table 1 that for $H_{t}$ given by (4.1) (resp. (4.2)) that for $\Delta\leq.30$ (resp. .25), $\widehat{\Delta}\left(k\widehat{\sigma}_{\varepsilon}\right)$’s bias appears to decrease as $k$ increases, while for $\Delta>.30$ (resp. .25) the bias seems to remain constant or slightly increase with $k$. The estimator’s variance appears to nearly always increase with $k$, which is likely due to the higher variance in probability estimates for larger deviations. It is also clear that the variance increases with $\Delta$, which could be related to the increased variance in $\widehat{\sigma}_{\varepsilon}$ for larger $\Delta$. Doubling the sample size decreased the estimators’ standard deviations by .01 to .02, and had a similar effect on the some of the biases. Although the bias actually increased in many cases by increasing the sample size.

The estimator $\widehat{\Delta}\left(2\widehat{\sigma}_{\varepsilon}\right)$’s bias is very large for small $\Delta$ and decreases greatly as $\Delta\rightarrow 1/2$. To understand this behavior, recall that (2.14)’s correction terms are inversely proportional to $\Delta^{2}$, which makes Theorem 2.3’s approximation less accurate for smaller $\Delta$. Intuitively, this is because as $\Delta$ increases, (2.4) becomes more sharply peaked around the origin. Hence a given $\Lambda$, such as $2\widehat{\sigma}_{\varepsilon}$, will become sufficiently far in the tails for (2.11) to be a good approximation to (3.4)’s integral, making (3.5) a more appropriate estimator.

Meanwhile for $k=3$ and $4$, $\widehat{\Delta}\left(k\widehat{\sigma}_{\varepsilon}\right)$’s bias appears to be much smaller than $\widehat{\Delta}\left(2\widehat{\sigma}_{\varepsilon}\right)$’s, and minimal around $.10\leq\Delta\leq.2$. The inherent upward bias in the $\widehat{\Delta}$’s is likely due to model error in $\widehat{\bar{H}}_{t}$, since decreased predictability of $H_{t}$ manifests as a higher empirical $\Delta$. This explains the increases in bias seen when doubling the sample size, because the AR(10) models for $\widehat{\bar{H}}_{t}$ overfitted the data less on the larger samples.

The processes (4.1)-(4.2) were selected for their simplicity and the superficial similarity between their short-term autocorrelation structure and those seen in financial data sets. In particular, (4.1) and (4.2)’s corresponding $\left|\varepsilon_{t}\right|$’s appear to be slightly more strongly locally correlated, but otherwise their short-term autocorrelations resemble those seen in recent daily S&P 500 and U.S. Dollar/Euro absolute log-returns. Hence these simulation results show that with a realistically sized sample and similar conditions to the financial application below, $\widehat{\Delta}\left(k\widehat{\sigma}_{\varepsilon}\right)$ with $k=3$ or $4$ can provide a computationally inexpensive and satisfactory estimate of $\Delta$, given a predictive model for $H_{t}$. Now that we have shown the effectiveness of Section 4’s $\Delta$ estimation procedure, in the next two sections we will show the predictive capabilities of our modeling method using simulated and financial time series data.

5 Simulation: Lorenz-Driven Volatility

In this section we will apply our modeling method to a simulated nonlinear time series where the volatility is driven by the Lorenz system. Note that this makes the volatility non-stochastic but instead driven by a deterministic, nonlinear system. This simulation study will show that even when the underlying process is not defined by (1.1)-(1.4), the modeling method we have presented can still deliver highly predictive, dynamic estimates of volatility and extreme event probabilities.

The Lorenz system is a well-known set of ordinary differential equations, originally formulated as a simplified model for atmospheric convection by Lorenz (1963). The system of equations, given by

[TABLE]

describes a two-dimensional fluid layer uniformly heated from below and cooled from above. The variables $x$, $y$, and $z$ are proportional to the rate of convection, and the horizontal and vertical temperature variation respectively, while $\sigma$, $\rho$, and $\beta$ are physical parameters. The Lorenz system is highly nonlinear and famously exhibits chaotic behavior, which means the system’s long-term evolution is highly sensitive to its initial conditions. Here we set the parameters $\sigma=10$, $\rho=28$, and $\beta=8/3$, which generate the system’s most well-known and chaotic behavior. We use the R package “deSolve” given in Soetaert et al. (2010) to create 10,000 samples from the discretized Lorenz system with the given parameter values and initial conditions of $x=0$, $y=1$, and $z=1$.

The Lorenz system enters into our simulated time series by driving the volatility. In particular, we set (1.2)’s $H_{t}$ as

[TABLE]

with $x_{t}$ defined by the discretized version of (5.1). Thus, $H_{t}$ is given by the centered and de-scaled $x$ output from the Lorenz system. To create the observable, heavy-tailed white noise $\varepsilon_{t}$, we then multiply (5.2)’s $\exp\left(H_{t}\right)$ by $z_{t}\sim\mathcal{N}(0,1)$ i.i.d. The Lorenz-driven $H_{t}$ with the resulting volatility $\sigma_{t}$ and time series $\varepsilon_{t}$ are shown in Fig. 4. Note the clear bursts of high volatility and extreme fluctuations. We will use our methodology to predict $\varepsilon_{t}$’s dynamically changing volatility and extreme event probabilities.

To make predictions from the observable time series $\varepsilon_{t}$, we build a sparse, linear, AR model for (3.3)’s $\widehat{\bar{H}}_{t}$ using the 20 previous values of (3.2)’s $\widehat{H}_{t}$. Hence (3.3)’s set of explanatory variables is simply $\left\{\widehat{H}_{t-1},...,\widehat{H}_{t-20}\right\}$. Because of high correlation and potential multicollinearity in this set of explanatory variables, we first apply Principal Component Analysis (PCA) to the training data to create the uncorrelated Principal Components (PCs) $\phi^{(1)}_{t},...,\phi^{(20)}_{t}$. We then use Tibshirani’s (1996) LASSO regression and Friedman et al.’s (2010) cyclic coordinate descent to give the model

[TABLE]

where

[TABLE]

and $\lambda$ minimizes $\widehat{\bar{H}}_{t}$’s mean absolute 10-fold cross-validation error. This procedure, which can be called L1-penalized PC Regression, gives a sparse AR model for $\widehat{\bar{H}}_{t}$. We will use this simple model with (3.5)’s $\widehat{\Delta}=\widehat{\Delta}\left(4\widehat{\sigma}_{\varepsilon}\right)$ to give an estimate of $\varepsilon_{t}$’s next time-step volatility $\bar{\sigma}_{\varepsilon,t}=\sqrt{E\left\{\varepsilon_{t}^{2}|\mathcal{F}_{t-1}\right\}}$ as well as a dynamic estimate of the probability that $\left|\varepsilon_{t}\right|\geq 3\widehat{\sigma}_{\varepsilon}$.

We use the moment result (2.2) to give the following model for the next time-step volatility:

[TABLE]

We also use the asymptotic approximation (2.11)-(2.8) to give the dynamic probability estimate

[TABLE]

To measure the predictive capability of (5.5) we can simply use the empirical correlation

[TABLE]

To measure the predictiveness of (5.6), we use it to create a classifier. We begin by noting that the stationary probability that $\left|\varepsilon_{t}\right|\geq 3\widehat{\sigma}_{\varepsilon}$ under a Gaussian distribution is $\approx.0027$. We then use (5.6) to define a binary classifier $\xi_{t}$ with the form

[TABLE]

Hence $\xi_{t}$ identifies times when the probability that $|\varepsilon_{t}|\geq 3\widehat{\sigma}_{\varepsilon}$ exceeds 5 times the stationary Gaussian probability. We then calculate $\xi_{t}$’s sensitivity (Sn.) and specificity (Sp.), which are defined here as

[TABLE]

and

[TABLE]

The measure (5.9) gives the percentage of 3 standard deviation events predicted by $\xi_{t}$, while (5.10) gives the percentage of non-3 standard deviation events predicted by $\xi_{t}$. We evaluate our modeling by doing three different sets of training and testing exercises for several different ratios of training sample size to testing sample size. We run 100 tests each and list averages of (3.5)’s $\widehat{\Delta}\left(4\widehat{\sigma}_{\varepsilon}\right)$, (5.7)’s $\rho_{|\varepsilon|,\widehat{\sigma}}$, and (5.9)-(5.10) in Table 2.

Note from Table 2 that the modeling is highly predictive of volatility and extreme events in $\varepsilon_{t}$. In particular, (5.6) and (5.8)’s naive classifier $\xi_{t}$ is capable of predicting 93% to 96% of 3 standard deviation events out-of-sample. The classifier does this while still maintaining a high specificity of 73% to over 78%. We also note that (5.5)’s $\widehat{\bar{\sigma}}_{\varepsilon,t}$ achieves over 60% correlation with $\left|\varepsilon_{t}\right|$ out-of-sample, which shows the modeling’s ability to predict fluctuations in volatility. This predictive performance is achieved with training data sets that are smaller than the testing data sets.

It is important to note that this simulated time series is a considerable departure from the process defined in (1.1)-(1.4), as well as from the linear autoregressive model for $H_{t}$ defined in (5.3)-(5.4). The volatility in this time series is not stochastic, but is instead driven by a deterministic, nonlinear, and chaotic system that involves unobserved variables. Nevertheless, our modeling method can produce highly predictive results. This suggests that this methodology could be useful for modeling a wide variety of nonlinear and heavy-tailed time series. Since Litimi et al. (2019) have recently demonstrated that chaos is likely to be present in financial market volatility, these simulation results are potentially directly relevant to financial applications. In the next section we will present an empirical application to financial market data where the modeling method can achieve similar predictive performance to that shown in Table 2.

6 Empirical Application

In this section we apply our modeling and estimation procedure to daily log-returns of the S&P 500 Index (SPX). We first consider a small sample of recent data from the last 5 years before next considering a much larger sample of data from the last 29 years. For both data sets we build a sparse, linear regression model for (3.3)’s $\widehat{\bar{H}}_{t}$ using two sets of covariates. The first set is the 10 previous values of (3.2)’s $\widehat{H}_{t}$, $\left\{\widehat{H}_{t-1},...,\widehat{H}_{t-10}\right\}$. The second set is the 10 previous values of $\log\left(\text{VIX}_{t}\right)$, $\left\{\log\left(\text{VIX}_{t-1}\right),...,\log\left(\text{VIX}_{t-10}\right)\right\}$, where VIX refers to the Chicago Board Options Exchange (CBOE) Volatility Index, which is an implied volatility measure computed from SPX option prices. It is perhaps the most well-known and widely disseminated measure of implied volatility for the S&P 500 Index.

Similarly to the previous section, because of potential multicollinearity in the model’s set of explanatory variables

[TABLE]

we again apply PCA to the training data to give the PCs $\phi^{(1)}_{t},...,\phi^{(20)}_{t}$. We also again use LASSO regression, giving a model of the form (5.3)-(5.4). We use the model for $\widehat{\bar{H}}_{t}$ to similarly give the dynamic volatility and probability estimates (5.5)-(5.6) as well as the classifier (5.8).

6.1 2014 to 2019

We first consider approximately 5 years of data from February 27, 2014 to February 14, 2019, giving 1,250 total samples. We evaluate our modeling using backtesting. In particular we do two sets of backtests. In the first backtesting exercise, we train the model on the first 50% of the data, or 625 samples through August 18, 2016. We then test the model’s performance on the second 50% or 625 samples of the data. In the second backtesting exercise we train the model on the first 60% of the data, or 750 samples through February 17, 2017. We test the model’s performance on the second 40% or 500 samples of the data. We run both backtests 100 times each and list averages of (3.5)’s $\widehat{\Delta}\left(4\widehat{\sigma}_{\varepsilon}\right)$, (5.7)’s $\rho_{|\varepsilon|,\widehat{\sigma}}$, and (5.9)-(5.10) in Table 3.

6.2 1990 to 2019

We now consider approximately 29 years of data from January 17, 1990 to March 1, 2019, giving 7,332 total samples. We again perform two sets of backtests. In the first backtesting exercise, we train the model on the first 50% of the data, or 3,666 samples through August 4, 2004. We then test the model’s performance on the second 50% or 3,666 samples of the data. In the second backtesting exercise we train the model on the first 60% of the data, or 4,399 samples through July 3, 2007. We test the model’s performance on the second 40% or 2,933 samples of the data. We again run both backtests 100 times each and list averages of (3.5)’s $\widehat{\Delta}\left(4\widehat{\sigma}_{\varepsilon}\right)$, (5.7)’s $\rho_{|\varepsilon|,\widehat{\sigma}}$, and (5.9)-(5.10) in Table 3.

6.3 Discussion

It is clear from Table 3 that the modeling is highly predictive of extreme events in daily fluctuations of the S&P 500 Index. (5.6) and (5.8)’s naive classifier $\xi_{t}$ can predict 85% of 3 standard deviation events out-of-sample, and was observed to predict as much as 91% of such events. The classifier does this while maintaining a specificity of 76% on average, with an observed range of 71% to 80%. We also note that (5.5)’s $\widehat{\bar{\sigma}}_{\varepsilon,t}$ achieves nearly 50% correlation with $\left|\varepsilon_{t}\right|$ out-of-sample, showing the model’s ability to predict fluctuations in market volatility. These results are achieved using very small amounts of training data.

On the larger data set, $\xi_{t}$ shows a still better balance of sensitivity and specificity, predicting 89% to 91% of 3 standard deviation events on average, while maintaining a specificity of 73% to 81%. Additionally, the correlation between $\widehat{\bar{\sigma}}_{\varepsilon,t}$ and $\left|\varepsilon_{t}\right|$ reaches 58%, indicating even higher capability for predicting fluctuations in volatility. Note that the larger-sample testing data included the exceptionally high volatility periods of the financial crisis.

We note that, by Theorem 2.3, the approximation (5.6) diverges beyond unity for high values of $\widehat{\bar{H}}_{t}$. This exceedance beyond unity occasionally occurs in these empirical applications as well as in the previous section’s application to Lorenz-driven volatility. In such extreme environments, where $\widehat{\bar{H}}_{t}$ is very large, (5.6) cannot be relied upon to give accurate probability estimates. Rather, (2.11)-(2.8) and corresponding estimates with the form of (5.6) can be used as early warning tools for upcoming extreme events.

We also note that since the model for $\widehat{\bar{H}}_{t}$ uses the 10 previous values of $\widehat{H}_{t}$ and $\log\left(\text{VIX}_{t}\right)$, it has a memory length of about two trading weeks. This relatively naive, short-memory modeling of $H_{t}$ could potentially be improved for this application, since many researchers have found empirical evidence of long-memory in SPX volatility (see Ding et al. 1993, Bollerslev and Mikkelsen 1996, Lobato and Savin 1998, Ray and Tsay 2000, Grau-Carles 2000). The modeling could also be adjusted to accomodate the asymmetric tails observed in stock market data, often called the “leverage effect” (see Black 1976, Christie 1982, Nelson 1991, Engle and Ng 1993). Overall, though, the study here shows that a simple, short-memory model for $H_{t}$ that uses only recent realized and implied volatility measures, with the parametrization (1.1)-(1.5) enables straightforward and predictive tail inference in financial time series. Plots of the SPX daily log-returns from 1990 to 2019 along with corresponding model outputs are given in Fig.5.

7 Conclusion

We have presented a family of stochastic volatility models that enable dynamic tail inference in heavy-tailed time series. The family’s conditional probabilistic structure allows for straightforward, effective, and computationally inexpensive estimation of model parameters and outcome probabilities. We have shown that this modeling formalism can be useful for predictive inference of dynamically changing extreme event probabilities in nonlinear and financial time series data. Current and future directions of research involve seeking formal results on model estimation, long-memory modeling for conditionally Laplace processes, accomodation of asymmetric volatility and tails, as well as time aggregation.

Acknowledgements

The author thanks Prof. Gennady Samorodnitsky, Cornell University, and Prof. Richard Davis, Columbia University, for their comments on a very early version of this work in 2016, as well as Dr. Dobrislav Dobrev, Federal Reserve Board, for helpful comments on an early version of this manuscript at the 2017 10th Annual Society for Financial Econometrics (SoFiE) Conference. The author is also grateful to Dr. Altan Allawala, JPMorgan Chase & Co., for his insightful comments. The author also thanks Prof. Richard Kleeman, Courant Institute of Mathematical Sciences, and Prof. Clifford Hurvich, NYU Stern, for many hours of motivating discussions, along with Prof. Fahad Saled, McGill University, and Franz Hinzen, NYU Stern, for encouragement and conversation on these topics. Finally, the author is very grateful to Prof. David Jablons, University of California San Francisco, for his generous encouragement and support.

Appendix A Appendix 1: Proofs

A.1 Proof of Lemma 2.1

Proof.

We first make the substitution $y=\log\left(y^{\prime}/e^{\bar{H}_{t}}\right)$. Then by (1.2) and (1.4), for $h_{t}\geq 0$,

[TABLE]

We then differentiate w.r.t. $y^{\prime}$ and let $y^{\prime}=\sigma_{t}$, which gives (2.1)’s result for $\sigma_{t}\geq\exp\left(\bar{H}_{t}\right)$. Next, by (1.2), (1.4), for $h_{t}<0$,

[TABLE]

We again differentiate w.r.t. $y^{\prime}$ and let $y^{\prime}=\sigma_{t}$ to give (2.1)’s result for $0\leq\sigma_{t}<\exp\left(\bar{H}_{t}\right)$. ∎

A.2 Proof of Corollary 2.2

Proof.

We note from (2.1) that for any $n\geq 1$, the conditional expected value of $\sigma_{t}^{n}$ is equal to the sum of integrals

[TABLE]

the second term of which diverges for $\Delta\geq 1/n$, while for $\Delta<1/n$ we have after integration

[TABLE]

We next note that since $z_{t}$ is i.i.d. and $z_{t}$ and $\sigma_{t}$ are independent, the conditional expectation $E\left\{\varepsilon_{t}^{n}|\mathcal{F}_{t-1}\right\}=Ez_{t}^{n}E\left\{\sigma_{t}^{n}|\mathcal{F}_{t-1}\right\}$. Since $z_{t}$ is standard normal, its $n$th moment is simply $(n-1)!!=(n-1)(n-3)(n-5)...1$. Therefore $E\left\{\varepsilon_{t}^{n}|\mathcal{F}_{t-1}\right\}=(n-1)!!E\left\{\sigma_{t}^{n}|\mathcal{F}_{t-1}\right\}$. Using (A.1) then gives (2.2)’s result. ∎

A.3 Proof of Theorem 2.1

Proof.

We note from Rohatgi (1976 p. 141) that since (1.1)’s $z_{t}$ and $\sigma_{t}$ are independent, the distribution of their product $\sigma_{t}z_{t}=\varepsilon_{t}$ is given by the formula

[TABLE]

We then substitute the standard normal density for $p_{z}$ and (2.1) for $p_{\sigma}$ into (A.2) to give

[TABLE]

We note from (2.5) that the first integral in (A.3) can be written in the following form after simplifying

[TABLE]

To simplify (A.4) further, we note from Abramowitz and Stegun (1965 p. 263) that the upper incomplete gamma function satisfies

[TABLE]

as $b\rightarrow\infty$. Letting $(1-1/\Delta)/2=a$ and $\varepsilon_{t}^{2}/2\sigma_{t}^{2}=b$ in (A.5) shows that the second term in (A.4) vanishes. Final simplification of (A.4) gives the result

[TABLE]

The definition (2.5) can be used again to write the second integral in (A.3) as

[TABLE]

after simplifying. The first upper incomplete gamma function term in (A.7) is simply $\Gamma\left((1+1/\Delta)/2\right)$. We then note that

[TABLE]

Using this identity with (A.7) and simplifying gives

[TABLE]

Adding (A.6) and (A.9) gives the result in (2.4). ∎

A.4 Proof of Corollary 2.7

Proof.

We first note the identity shown in Jameson (2016, 2017)

[TABLE]

as $b\rightarrow 0$ for $a<0$. For $\Delta<1$, $(1-1/\Delta)/2<0$, and hence using (A.10) and rewriting gives

[TABLE]

as $\left|\varepsilon_{t}\right|\rightarrow 0$. We next note the identity

[TABLE]

as $b\rightarrow 0$. Using (A.12) and rewriting then gives

[TABLE]

as $\left|\varepsilon_{t}\right|\rightarrow 0$. Substituting (A.11) and (A.13) into (2.4) and cancelling terms gives

[TABLE]

as $\left|\varepsilon_{t}\right|\rightarrow 0$, which can be simplified to give the result (2.7). ∎

A.5 Proof of Theorem 2.2

Proof.

We first make the substitution

[TABLE]

in (2.4)’s integral from $\Lambda$ to $\infty$. Multiplying by 2, cancelling terms, and using the notation (2.8) then gives

[TABLE]

We begin with the first term of (A.15). We apply integration by parts with

[TABLE]

Then using the definition (2.5) for differentiation of $u$ and cancelling terms gives

[TABLE]

Recall from (A.5) that (A.16)’s first term vanishes. We then proceed to (A.15)’s second term, again applying integration by parts with

[TABLE]

Then using (2.6) for differentiation of $u$ and cancelling terms gives

[TABLE]

Note that (A.17)’s first term also vanishes. Combining (A.16) and (A.17) and noting the definition (2.10) gives the result (2.9). ∎

A.6 Proof of Theorem 2.3

Proof.

We will use the following asymptotic expansion of (2.5), which can be derived from repeated integration by parts (see Digital Library of Mathematical Functions 8.11).

[TABLE]

as $b\rightarrow\infty$, where

[TABLE]

and $(.)_{k}$ denotes the rising factorial. We then make the following substitutions in (2.4),

[TABLE]

We then recall (A.8) and substitute (A.18) into (2.4)’s two terms (A.6) and (A.9). Simplifying then gives the series representation of (2.4) as $b_{t}\rightarrow\infty$,

[TABLE]

We next note that, after re-substituting for $b_{t}$ and $a_{2}$, the first term in (A.21) is equal to

[TABLE]

which is the limit of the term (A.9) as $\left|\varepsilon_{t}\right|\rightarrow\infty$. Multiplying (A.22) by 2 and integrating from $\Lambda$ to $\infty$ gives the first term of (2.11) with (2.8). We then proceed to (A.21)’s second term, re-substituting for $b_{t}$ and cancelling factors to give

[TABLE]

We next note the following integral,

[TABLE]

So integrating (A.23) from $\Lambda$ to $\infty$, applying (A.24), and cancelling factors gives

[TABLE]

Because (A.21)’s remainder at any $n$ does not change sign over all $b_{t}$, we can similarly apply (A.24) to give the remainder term’s integral. Finally, multiplying (A.21) by 2 and integrating from $\Lambda$ to $\infty$, applying the integral of (A.22), the result (A.25), using the clearer notation for $u_{k}(a)$, and fully re-substituting gives the asymptotic expansion

[TABLE]

Re-applying (A.18) to (A.26)’s remainder term with $n=1$ gives the remainder term in (2.11) with (2.8). ∎

The series in (A.21) and (A.26) are generally not convergent as $n\rightarrow\infty$. However, they can be truncated at low order $n$ to closely approximate (2.4) and its integral. This is because for $n\ll\infty$ the series’ terms rapidly become very small for values of $\Lambda$ that are relatively large compared to $e^{\bar{H}_{t}}$. Visualizations of (A.26)’s asymptotic series are given in Fig. 6.

A.7 Proof of Proposition 3.1

Proof.

By (1.1) and (1.2), $\left|\varepsilon_{t}\right|=\exp\left(H_{t}\right)\left|z_{t}\right|$. Taking the logarithm gives

[TABLE]

Since $z_{t}$ and $\varepsilon_{t}$ are independent, $E\left\{\log\left|z_{t}\right||\log\left|\varepsilon_{t}\right|\right\}=E\log\left|z_{t}\right|$. Then since $z_{t}$ is standard normal, $z_{t}^{2}\sim\chi^{2}(1)$. It can then be noted, e.g., from Breidt et al. (1998) that

[TABLE]

where $\psi(x)$ is the digamma function

[TABLE]

and $\psi(1/2)=-2\log 2-\gamma$ (see Abramowitz and Stegun 1965 p. 258-259). Rearranging (A.27), taking the expected value $E\left\{H_{t}|\log\left|\varepsilon_{t}\right|\right\}$, and substituting (A.28) gives the result (3.1). ∎

Appendix B Appendix 2: The $z_{t}\sim Lap(0,1)$ Case

In this section we briefly show that very similar results to those in Appendix 1 for (1.1)-(1.4) with $z_{t}\sim\mathcal{N}(0,1)$ can be derived for the case where $z_{t}$ is instead distributed according to the standard Laplace density

[TABLE]

which gives $Ez_{t}=0$ and $E\left|z_{t}\right|=1$. We first note that an equivalent argument to that in A.2 can be made. We use the fact that (B.1)’s $Ez_{t}^{n}=n!$ with the result (A.1) to give

[TABLE]

for even $n\geq 2$. The right hand side of (B.2) with $n=1$ additionally gives $E\left\{\left|\varepsilon_{t}\right||\mathcal{F}_{t-1}\right\}$.

We next note that an equivalent argument to the one in A.3 can be given, using (2.1) and (B.1) with (A.2) to give an expression composed of two integrals, each of which can be written very similarly to (A.4) and (A.7). Then (A.5) and (A.8) can be used in the same way as above to simplify these terms. Adding these two terms gives the result

[TABLE]

The identities (A.10) and (A.12) in A.4 can then be used to show that as $\left|\varepsilon_{t}\right|\rightarrow 0$,

[TABLE]

Making the substitution $b=\left|\varepsilon_{t}\right|/e^{\bar{H}_{t}}$ in (B.3) and applying integration by parts as in A.5 gives the result

[TABLE]

where $\widetilde{\Lambda}_{t}=\Lambda/e^{\bar{H}_{t}}$. Lastly, the argument in A.6 can also be repeated, using (A.18)-(A.19) with the substitutions

[TABLE]

to write (B.3) as an asymptotic expansion that is equivalent in form to (A.21), with the only change being that the factor $1/\left(4\Delta\sqrt{2\pi}e^{\bar{H}_{t}}\right)$ is replaced with $1/\left(4\Delta e^{\bar{H}_{t}}\right)$. Then re-substituting for $b_{t}$ and $a_{2}$ in the first term gives

[TABLE]

which is the limit of (B.3)’s second term as $\left|\varepsilon_{t}\right|\rightarrow\infty$. Multiplying (B.6) by 2 and integrating from $\Lambda$ to $\infty$ will give the first term of our result. Next, fully re-substituting in the second term and remainder term of (B.3)’s asymptotic expansion, applying the integral

[TABLE]

multiplying by 2, and cancelling factors then gives the full result

[TABLE]

Re-applying (A.18) to (B.8)’s remainder term with $n=1$ shows that as $\widetilde{\Lambda}_{t}\rightarrow\infty$,

[TABLE]

Then comparing (B.9) and (2.12) gives

[TABLE]

Since $\Gamma(1)=\Gamma(2)=1$, (B.10)’s second line reduces to $P\left\{\left|\varepsilon_{t}\right|\geq\Lambda|\mathcal{F}_{t-1}\right\}\sim P\left\{\sigma_{t}\geq\Lambda|\mathcal{F}_{t-1}\right\}$ for $\Delta=1$ and as $\Delta\rightarrow\infty$. Several graphs of (B.3) are given in Fig. 7.

SUPPLEMENTARY MATERIAL

R scripts:

R scripts containing code to perform the simulation studies in Section 4 and Section 5 as well as the empirical application in Section 6. (.R files)

Data sets:

Data sets used in the simulation study in Section 5 and the empirical application in Section 6. (.xlsx files)

For supplementary materials email gchavez@novocure.com.