Measuring tail risk at high-frequency: An $L_1$-regularized extreme value regression approach with unit-root predictors

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.

118,258 characters · 16 sections · 49 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.

abstractWe study tail risk dynamics in high-frequency financial markets and their connection with trading activity and market uncertainty. We introduce a dynamic extreme value regression model accommodating both stationary and local unit-root predictors to appropriately capture the time-varying behaviour of the distribution of high-frequency extreme losses. To characterize trading activity and market uncertainty, we consider several volatility and liquidity predictors, and propose a two-step adaptive $L_1$-regularized maximum likelihood estimator to select the most appropriate ones. We establish the oracle property of the proposed estimator for selecting both stationary and local unit-root predictors, and show its good finite sample properties in an extensive simulation study. Studying the high-frequency extreme losses of nine large liquid U.S. stocks using 42 liquidity and volatility predictors, we find the severity of extreme losses to be well predicted by low levels of price impact in period of high volatility of liquidity and volatility.

Introduction

Measuring tail risk at high-frequency has become of utmost importance to market players and regulators weller2018. While much efforts have been devoted to the measurement of tail risk at low-frequency nieto2016frontiers, few attempts have been made to measure risk at high-frequency, see giot2005market, dionne2009intraday and chavez2012modelling. Moreover, although these models can be very accurate, they explain the tail risk evolution in a “reduced form” manner, i.e., using autoregressive terms exploiting the persistence of the time series. They thus fail to provide a deeper structural understanding of the factors driving tail risk. As much as understanding the macroeconomic determinants of tail risk is a relevant problem at low-frequency massacci2017tail, it is important to understand how market uncertainty and trading activity impacts tail risk at high-frequency.

From a market microstructure perspective, though the intensification of high-frequency trading has improved trading costs and liquidity hendershott2011, it is also suspected to be responsible for more frequent extreme price movements over short periods of time brogaard2018. Such extreme fluctuations are often the result of an aggressive directional market making activity initiated when the market is already under stress. brogaard2018 find that market wide extreme shocks are likely to trigger the risk controls of high-frequency liquidity providers that thus withdraw from the market to reduce their risk exposure. Similarly, kirilenko2017flash find that during the market turbulence induced by the 2010 Flash Crash, many high-frequency liquidity providers withdrew from the market, thus exacerbating the price fall. Studying how market uncertainty and trading activity affect extreme losses can thus provide a deeper understanding of the evolution of tail risk at high-frequency, and this paper proposes appropriate econometric techniques to do so.

We consider a dynamic extreme value regression framework chavez2016extreme,massacci2017tail,schwaab2021modeling where the distribution of extreme losses is assumed to be well approximated by a generalized Pareto distribution (GPD) with time-varying parameters driven by exogenous preditors and autoregressive terms. To assess the impact of market uncertainty and trading activity on extreme losses, we consider several volatility predictors, proxing for market uncertainty, and liquidity predictors, characterizing trading activity. Despite extreme value regression techniques have been widely applied in finance chavez2016extreme, hambuckers2018understanding,bee2019realized, our investigation presents new challenges: (i) as the financial literature proposes several volatility and liquidity measures, we face a variable selection problem aimed at identifying predictors capturing the most relevant aspects of trading activity affecting extremes as well as improving the predictive accuracy of tail risk; (ii) volatility and liquidity measures observed at high-frequency exhibit strong persistence and seasonalities, thus violating the classical stationary assumptions required for inference with the maximum likelihood estimator (MLE). To overcome these issues, we develop a two-step adaptive $L_1$-regularized maximum likelihood estimator (ALMLE) that allows performing variable selection with both stationary and local unit-root predictors lee2021LASSO, and establish its oracle property.

We investigate the impact of 42 liquidity and volatility indicators on the distribution of high-frequency extreme losses of nine large liquid U.S. stocks observed from 2006 to 2014. We find that the severity of tail risk, as measured by the shape parameter of the GPD, is well predicted by low price impact goyenko2009liquidity during periods of high volatility of volatility and high volatility of liquidity. This finding is coherent with the evidence in brogaard2018 that market markers liquidity supply is outstripped by liquidity demand after large uncertainty shocks, and their rush to leave the market to lower their risk exposures amplify extreme price movements. Our two-step ALMLE is necessary to reveal this pattern as the standard MLE finds almost all predictors to be significant. To validate our estimating strategy, we provide an out-of-sample VaR forecast analysis and find that the estimated model performs well in the out-of-sample.

The remainder of the paper is organized as follows: Section (ref) presents the time-varying GPD model accommodating stationary and local unit-root predictors as well as autoregressive components; Section (ref) presents the MLE and shows its asymptotic non-normality when local unit-root predictors are included in the model; Section (ref) introduces the two-step ALMLE and prove the oracle property of this estimator in selecting both stationary and local unit-root predictors; Section (ref) provides an extensive simulation study comparing the performance of the two-step ALMLE to those of the MLE, showing the superiority of the former in finite samples. Section (ref) discusses the results of the empirical study whereas Section (ref) concludes. Additional results and mathematical proofs are relegated to the Appendix.

Extreme value regression

We denote the logarithmic loss and return time series of a financial asset by $\{ l_t\}_{t=1}^T$ and $\{ r_t\}_{t=1}^T$, respectively, with $l_t = - r_t$, and denote $\bm{z}_t$ a vector of exogenous predictors observed at time $t$.

assumption{M.1} $\{l_t\}_{t=1}^T$ and $\{\bm{z}_t\}_{t=1}^T$ are on a complete probability space $(\Omega, \mathcal{F}, P)$. At each time $t \in\{ 1,2,\ldots,T\}$, we have an information set $\mathcal{F}_{t-1}$ available which is the $\sigma$-algebra generated by $\{\bm{z}_{t-1}, l_{t-1}, \bm{z}_{t-2},$ $l_{t-2}, \ldots \}$.

Let assume $\{ l_t\}_{t=1}^T$ is independent and identically distributed (i.i.d.) with a cumulative distribution function (c.d.f.) $F(\cdot)$. Probabilistic results from extreme value theory show that if there exist real sequences $a_T > 0$ and $\beta_T$ such that $\lim_{T\rightarrow \infty} F^T\left( a_T\,x + \beta_T\right)$ converges to a non-degenerate distribution $G(\cdot)$, then $F(\cdot)$ belongs to the max-domain of attraction of $G(\cdot)$, i.e. $F\in \mathcal{D}(G)$, and $G(\cdot)$ must be the generalized extreme value (GEV) distribution (see Theorem 3.1.1. of coles2001introduction).

Let $\{ y_t\}_{t=1}^T$ be a censored sequence of excess losses above a high threshold $u$, such that the excess loss $y_t = l_t - u$, if $l_t > u$, and $y_t = 0$ otherwise. Define the conditional distribution of excess losses, $$ F_{|u}(y) := P\left\{ l_t - u \leq y \middle| l_t > u\right\} = P\left\{ y_t \leq y \middle| y_t > 0\right\}, \quad 0< y_t\leq L^F - u, $$ with $L^F := \sup\{x: F(x) < 1\}$ the right end point of $F(\cdot)$. pickands1975statistical and balkema1974residual show that if $F(\cdot)\in \mathcal{D}(G)$ then the limiting distribution of $F_{|u}(y)$ is a GPD, i.e.

equation[equation omitted — 131 chars of source]

where $\text{GPD}(\cdot; k,\sigma)$ denotes the GPD with shape parameter $k\in\mathbb{R}$ and scale parameter $\sigma>0$,

equation[equation omitted — 110 chars of source]

Eq. (ref) suggests that $F_{|u}(y)$ with $u$ large enough can be approximated by a $\text{GPD}(\cdot; k,\sigma_{u})$, where the scale parameter $\sigma_u$ depends on $u$. The peaks-over-threshold (POT) approach assumes this relationship holds exactly above a fixed threshold $u$ and uses the exceedances of such threshold to estimate the GPD parameters $\sigma$ and $k$ (see section 4.3 of coles2001introduction).

Time-varying peaks-over-threshold (POT) approach

The classical POT approach assumes that $\{l_t\}$ is i.i.d. However, financial data typically exhibit dependence features such as time-varying heteroscedasticity and extremal clustering that violate this assumption. To capture these aspects, we adopt a dynamic POT approach. Let $\{ y_t\}_{t=1}^T$ be a censored sequence of excess losses over a threshold time series $\{u_t\}_{t=1}^T$, we model the excess loss distribution conditional on the information set $\mathcal{F}_{t-1}$, $F_{t|u_t}(y):=P\{ y_t\leq y | y_t>0, \mathcal{F}_{t-1}\}$, using a GPD with time-varying parameters $k_t$ and $\sigma_t$. See, e.g., chavez2014extreme,massacci2017tail,bee2019realized.

Consider the vector-valued time series of $p\in \mathbb{N}$ explanatory variables $\{\bm{z}_t := [z_{1,t},\ldots,z_{p,t}]'\}^T_{t=1}$. Given the information set $\mathcal{F}_{t-1}$, we consider the following specification for $\{ (k_t, \sigma_{t})\}$,

align[align omitted — 292 chars of source]

We impose that $0<k_t<0.5$ and $\sigma_t>0$ (see hosking1987parameter) to ensure a finite conditional variance of $y_t$ and numerical stability in the estimation. As the scale parameter $\sigma_t$ can be associated with the variance of the underlying distribution $F_{t}(\cdot)$, we accommodate an autoregressive term in $\log(\sigma_t)$ in the spirit of GARCH models engle2001garch. We allow for both stationary and unit-root explanatory variables in (ref) and (ref), such that persistent predictors can be accommodated.

Maximum likelihood estimation

Let $\bm{\beta}: = [\beta_{1,0}, \beta_{1,1}, \ldots, \beta_{1,p}, \beta_{2,0}, \beta_{2,1}, \ldots, \beta_{2,p+1}]'$ denote the vector of the model coefficients in (ref)-(ref), and define the coefficient space $\Theta$ of $\bm{\beta}$ as a subspace of $\mathbb{R}^{2p + 2}\times(-1,1)$ accomodating all permissible coefficient vectors $\bm{\beta}$. We present the MLE of the model coefficients in (ref)-(ref) and show it is consistent but asymptotically non-normal when local unit-root explanatory variables are included in the model.

Maximum likelihood estimator

assumption{M.2} We assume that for a given $\{u_t\}$, the conditional c.d.f. $F_{t}(\cdot)$ of $l_t$ given $\mathcal{F}_{t-1}$ exists for $t\in \{1,2,\ldots, T\}$ and $y_t:=l_t-u_t>0$ follows a time-varying GPD, i.e. \begin{equation} F_{t|u_t}(y_t) = GPD(y_t;k_t,\sigma_{t}) = 1 - \left( 1 + k_t\frac{y_t}{\sigma_t} \right)^{-\frac{1}{k_t}}, \end{equation} where $\{k_t\}$ and $\{\sigma_t\}$ are specified by (ref)-(ref) with the true coefficient vector $\bm{\beta}^o\in \Theta\subset \mathbb{R}^{2p + 2}\times(-1,1)$. Moreover, $\{u_t\}$ returns a constant unconditional exceedance rate, i.e., $P\{y_t>0\}=\tau$ for all $t\in\{1,\ldots,T\}$ with a constant $\tau$ close to zero.

assumption{M.3} Among the explanatory variables in Model (ref)-(ref), we assume that $\{ z_{i,t}, i = 1, \ldots, p_0\} \in I(0)$ and $\{ z_{j,t}, j = p_0+1, \ldots, p\} \in I(1)$ with $\epsilon_{j,t} := z_{j,t} - z_{j,t-1} $ and $\{ \epsilon_{j,t}, j = p_0+1, \ldots, p\} \in I(0)$. We denote by $I(0)$ and $I(1)$ the set of stationary and unit-root predictors, respectively.

Under Assumption (ref), the conditional probability density function (p.d.f.) of $y_t |\{y_t>0, \mathcal{F}_{t-1}\}$ is

equation[equation omitted — 137 chars of source]

and the log-likelihood function $L(\cdot)$ of $\{y_t| y_t>0, \mathcal{F}_{t-1}\}$ can be defined as schwaab2021modeling,

equation[equation omitted — 343 chars of source]

where

equation[equation omitted — 358 chars of source]

for $t = 1,2,\ldots,T$, with $\mathbbm{1}\{\cdot\}$ the indicator function taking value one if the input is true and zero otherwise.

We consider standardized predictors $\{\bm{z}^*_t\}$ in the estimation to get stochastically bounded variables, i.e., for each $t\in\{1,\ldots,T\}$, we standardize $\bm{z}_t$ as follows,

equation[equation omitted — 133 chars of source]

Replacing $\{\bm{z}_t\}$ with $\{\bm{z}^*_t\}$ into the likelihood function in (ref) and maximizing we obtain

equation[equation omitted — 206 chars of source]

where $\Theta \subset \mathbb{R}^{2p + 2}\times(-1,1)$. We denote the corresponding vector of true coefficients $\bm{\beta}^{o*}$.

Remark.

Assumption (ref) assumes a constant unconditional probability for the exceedance $\mathbbm{1}\{y_t>0\}$ for $t\in\{1,\ldots,T\}$, which is more general than assuming a constant conditional probability for $\mathbbm{1}\{y_t>0|\mathcal{F}_{t-1}\}$. This causes us no extra burden to obtain the limiting behaviour of the MLE because $\mathbbm{1}\{y_t>0\}$ is bounded and not a function of the model coefficients. Assumption (ref) allows for both stationary and unit-root predictors among $\bm{z}_t$.

Asymptotic properties of the MLE

smith1985maximum establishes the asymptotic properties of the MLE of a GPD with constant $k$ and $\sigma$ in an i.i.d. setting. We extend smith1985maximum establishing the consistency and limiting distribution of the MLE of the dynamic GPD with stationary and unit-root predictors in (ref). In what follows, we list the assumptions required to derive the asymptotic behaviour of the MLE, and establish the consistency and limiting distribution of $\widehat{\bm{\beta}}$.

assumption{M.4} We assume that, \begin{equation} \left\{ \begin{aligned} & \beta_{s,i}^{o*}:= \beta_{s,i}^o = O(1) , \quadfor\quad i = 0, 1, \ldots, p_0, \;and\; s = 1,2 ; \\ & \beta_{s,j}^{o*}:=\sqrt{T}\beta_{s,j}^o = O(1), \quadfor\quad j = p_0+1,\ldots, p_0 + p \;and\; s = 1,2 ; \\ & \beta_{2,p+1}^{o*}:=\beta_{2,p+1}^o \in (-1,1), \end{aligned} \right. \end{equation} and $\bm{\beta}^{o*}:= [\beta_{1,0}^{o*}, \ldots, \beta_{1,p}^{o*}, \beta_{2,0}^{o*}, \ldots, \beta_{2,p+1}^{o*}]'\in\mathbb{R}^{2p+2}\times(-1,1)$.

assumption{M.5} $\{\bm{\epsilon}_t := \left[ z_{1,t}, \ldots, z_{p_0,t}, \epsilon_{p_0+1,t}, \ldots, \epsilon_{p,t} \right]'\}_{t=1}^T$ is assumed $\text{i.i.d.}\,(\bm{0}, \Sigma^{(0)}) $ with mean $\bm{0}$ and positive definite covariance matrix $\Sigma^{(0)}$. With $ \left\{\bm{z}_t^* := [\; z_{1,t},\ldots, z_{p_0,t}, \frac{z_{p_0+1,t}}{\sqrt{T}}, \ldots, \frac{z_{p,t}}{\sqrt{T}} \;]'\right\} $, we assume that as $T\rightarrow\infty$, we have that \begin{equation*} \begin{aligned} & (1) \left\{ \begin{aligned} & \frac{1}{\sqrt{T}}\sum_{t=1}^{T} z_{i,t}^* = O_p(1), \quad \frac{1}{T}\sum_{t=1}^{T} (z^*_{i,t})^2 = O_p(1), \quad i = 1,2,\ldots,p_0; \\ & \frac{1}{T}\sum_{t=1}^{T} z^*_{j,t} = O_p(1), \quad \frac{1}{T}\sum_{t=1}^{T} (z^*_{j,t})^2 = O_p(1), \quad j = (p_0+1), \ldots,p; \\ & \frac{1}{T} \sum_{t=1}^{T} \bm{z}_t^* \bm{z}_t^{*'} is positive definite in probability one ; \end{aligned} \right. \\ & and there exists a positive definite matrix $\Sigma:=[\Sigma_{i,j}]_{i,j = 1,\dots, p}$ such that \\ & (2) \left\{ \begin{aligned} & \lim_{T\rightarrow\infty}\frac{1}{T}\sum_{t=1}^{T} z^*_{i,t} = 0, \quad \lim_{T\rightarrow\infty}\frac{1}{T}\sum_{t=1}^{T} (z^*_{i,t})^2 = 1, \quad \frac{1}{\sqrt{T}}\sum_{t=1}^{t_0} z^*_{i,t} \overset{D}{\sim} \,W_i( \frac{t_0}{T}), \quad i = 1,2,\ldots,p_0; \\ & \frac{1}{T}\sum_{t=1}^{T} z^*_{j,t} \overset{D}{\sim} \int^1_0 \Sigma_{j,j}^{1/2}\, W_j(t)\,dt, \quad \frac{1}{T}\sum_{t=1}^{T} (z^*_{j,t})^2 \overset{D}{\sim} \int^1_0 \Sigma_{j,j} W_j^2( t)\,dt, \qquad j = (p_0+1), \ldots, p; \\ & \frac{1}{T}\sum_{t=1}^{T} \bm{z}_t^* \bm{z}_t^{*'} \overset{D}{\sim} \int^1_0 \left(\Sigma^{1/2}\,\bm{W_{z*}}(t)\right)\,\left(\Sigma^{1/2}\,\bm{W_{z*}}(t)\right)'\,dt, \\ & \frac{1}{T} \sum_{t=1}^{T} \bm{z}_t^* \bm{z}_t^{*'} is positive definite in probability one; \end{aligned} \right. \end{aligned} \end{equation*} where $1\leq t_0 \leq T$, and $W_i(\cdot), W_j(\cdot)$ are independent Brownian motions. And denote \lq\, $\overset{D}{\sim}$ \rq \,for convergence in distribution and $\bm{W_{z*}}(t) : = [\frac{d W_1(t)}{\sqrt{dt}}, \ldots, \frac{d W_{p_0}(t)}{\sqrt{dt}}, W_{p_0+1}(t), \ldots, W_{p}(t)]' $.

assumption{M.6} $\Theta$ is a compact subspace of $\mathbb{R}^{2p+2}\times(-1,1)$ containing the true coefficient vector $\bm{\beta}^{o*}$ such that $ H_{\mathcal{L}}(\cdot)$ is positive definite in $\Theta$ almost surely.

assumption{M.7} Given Assumptions (ref) and (ref), we further assume that as $T \rightarrow \infty$, it holds that \begin{equation} \frac{1}{\sqrt{T}} \sum^T_{t=1}\, \bm{\psi}_t(\bm{\beta}^{o*}) \overset{\mathcal{D}}{\sim} S_{\psi}, \end{equation} where $S_{\psi}$ is a non-degenerate distribution: \begin{equation} S_{\psi} = \Sigma^{1/2}_{\psi} \int_0^1\; \begin{bmatrix} dW_{g_k}(t) \\ \bm{W_{z*}}(t) dW_{g_k}(t) \\ dW_{g_{\sigma}}(t) \\ \bm{W_{z*}}(t) dW_{g_{\sigma}}(t) \\ d W_{ar}(t)\,dW_{g_{\sigma}}(t) \end{bmatrix} \end{equation} where $\Sigma_{\psi}$ is a positive definite $(2p+3)\times(2p+3)$ matrix and $\{W_{g_k}(t)\}, \{W_{g_{\sigma}}(t)\}, \{W_{ar}(t)\}$ are Brownian motions.

assumption{M.8} $\frac{H_L(\bm{\beta}^{o*}; \{l_t\}, \{\bm{z}^*_t\})}{T}$ at $\bm{\beta}^{o*}$ is assumed to weakly converge to a stochastic integral $\Omega_{H}$, i.e., \begin{equation} \frac{H_{\mathcal{L}}( \bm{\beta}^{o*}; \{l_t\}, \{\bm{z}_t^*\}}{T} \overset{D}{\sim} \Omega_{H}, \qquad when T\rightarrow\infty, \end{equation} where $\Omega_{H}^{-1}$ exists upon the limiting behaviours of ${\bm{z}^*_t}$ in Assumption (ref).

Remark.

Assumption (ref) imposes the orders of magnitude of $\bm{\beta}^{o}$ to ensure that the unit-root explanatory variables $z_{p_0+1,t}, \ldots, z_{p,t} $ have coefficients of local-to-zero rate being $\frac{1}{2}$, see phillips2013predictive and lee2021LASSO. Assumption (ref) ensures that the partial sums of $\{\bm{z}_t^*\}$ and $\{\bm{z}_t^*\,\bm{z}_t^{*'}\}$ converge at specific rates. Assumption (ref) is also used by saikkonen1993continuous, saikkonen1995problems,lee2021LASSO, and was shown to hold for time series with a moderate degree of temporal dependence and heteroscedasticity of $\{\bm{\epsilon}_t\}$. See, e.g.,Theorem 18.2 of billingsley2013convergence, phillips1986multiple, phillips1991optimal. Assumption (ref) restricts the permissible parameter space $\Theta$ for the ML estimation, especially maintaining the positive definiteness of $H_{\mathcal{L}}(\bm{\beta})$ in analogy to Assumption 9 of smith1985maximum for settling the uniqueness of the estimator. Assumption (ref) assumes the limiting distribution of the likelihood gradient function at $\bm{\beta}^{o*}$, see Lemma A.1.1 of lee2016predictive. Assumption (ref) assumes the existence of the limiting distribution of the likelihood Hessian matrix at $\bm{\beta}^{o*}$, see Lemma A.1.2 of lee2016predictive.

theorem[MLE consistency]\ \\ Under Assumptions (ref), (ref), (ref), (ref)(1), (ref) and (ref), and for any $\epsilon >0$, \begin{equation} \lim_{T\rightarrow\infty} P\left\{\left\lVert\widehat{\bm{\beta}}^{mle} - \bm{\beta}^{o*} \right\rVert>\epsilon \right\} = 0, \end{equation}

proofSee Appendix (ref).

theorem[MLE asymptotics]\ \\ Under Assumptions (ref) to (ref), we have \begin{equation} \sqrt{T}\left(\widehat{\bm{\beta}}^{mle}-\bm{\beta}^{o*}\right) \overset{D}{\sim} \Omega_{H}^{-1} \; S_{\psi}, \qquad as T\rightarrow\infty. \end{equation}

proofSee Appendix (ref).

Adaptive $L_1$-regularized maximum likelihood estimation

Variable selection facilitates interpretation of a regression model and solves the trade-off issue between bias and efficiency so as to achieve predictive accuracy, see james2013introduction. Although variable selection performed via inferential tests based on the asymptotic normality of the MLE might seem a viable solution, it is not appropriate in our setting because of the following three issues: (i) the inability to control type I error for multiple predictor selection; (ii) severe size distortion for selecting unit-root predictors because of the non-normal limiting distribution; (iii) low power in selecting predictors for the shape parameter due to high standard errors of coefficients, see simulations in Section (ref).

To circumvent these issues, we adopt $L_1$-regularized MLE for automatic variable selection tibshirani1996regression. Due to the constraining nature of $L_1$-regularization, this estimator sets some coefficients exactly to zero so as to perform variable selection. zou2006adaptive explore the advantages of using weighted $L_1$-regularization on model coefficients and proposed the adaptive LASSO. With proper adaptive weights, the adaptive LASSO exhibits the oracle property, which produces an asymptotic efficient estimator of variable selection consistency as if the true underlying model were given from the outset. medeiros2016 prove the oracle property for the adaptive LASSO in high-dimensional time series with non-Gaussian and heteroscedastic errors as well as with highly correlated regressors. kock2016consistent show that the adaptive LASSO is oracle efficient in stationary and non-stationary autoregressions. lee2021LASSO prove the oracle property of the adaptive LASSO with stationary and local unit-root predictors, and propose a novel post-selection adaptive LASSO for selecting mixed-root predictors i.e. stationary, local unit root, and cointegrated predictors.

Drawing on this literature, we extend the adaptive LASSO to the MLE in (ref) to estimate and select stationary and local unit-root predictors in (ref)-(ref). A general form of adaptive $L_1$-regularized maximum likelihood estimator (ALMLE) can be drawn directly from zou2006adaptive and is formulated as follows:

equation[equation omitted — 330 chars of source]

where $\mathcal{L}(\bm{\beta}; \{ y_t \}, \{ \bm{z}_t^* \})$ is the log-likelihood function specified in (ref); $\lambda_{k, T}, \lambda_{\sigma, T}>0$ are tuning parameters; and $w_{k,i}, w_{\sigma,j}$ are adaptive weights for penalizing coefficients differently. We consider two tuning parameters instead of one to be less restrictive on tuning parameter selection, thereby stabilizing the variable selection for both the shape and scale models in (ref)-(ref). To set the tuning parameters, we start off with large enough values of $\lambda_{k, T}$ and $\lambda_{\sigma, T}$ such that no predictors are selected by $\widehat{\bm{\beta}}^{\text{al}}$, and denote these two values as $\lambda_{k, T, \max}$ and $\lambda_{\sigma, T, \max}$, respectively. We then search for the optimal tuning parameters using an information criterion (IC) over equally-spaced grids of $n_{\lambda_k}$ and $n_{\lambda_{\sigma}}$ nodes\footnote{We use $n_{\lambda_k}=50$ and $n_{\lambda_{\sigma}}=30$ across this paper unless stated otherwise. We also have tried $n_{\lambda_k}=100$ and $n_{\lambda_{\sigma}}=100$ to check the sufficiency of $n_{\lambda_k}=50$ and $n_{\lambda_{\sigma}}=30$, and found that differences in the results are small.} defined on the intervals $[\lambda_{k, T, \max},10^{-6}]$ and $[\lambda_{\sigma, T, \max},10^{-6}]$. Formally, the grids for the shape and scale parameters are defined as $S_{\lambda_{k,T}}:=\{\exp(\log(\lambda_{k, T, \max}) - j\,\frac{\log(\lambda_{k, T, \max}) - \log(10^{-6})}{n_{\lambda_k} - 1} ), j=0,1,\ldots,(n_{\lambda_k}-1)\}$ and $S_{\lambda_{\sigma,T}}:=\{\exp(\log(\lambda_{\sigma, T, \max}) - j\,\frac{\log(\lambda_{\sigma, T, \max}) - \log(10^{-6})}{n_{\lambda_{\sigma}} - 1} ), j=0,1,\ldots,(n_{\lambda_{\sigma}}-1)\}$, respectively. We consider different information criteria, namely the Bayesian Information Criterion (BIC), the Hannan–Quinn information criterion (HQ) and the Akaike Information Criterion (AIC), and thus select the optimal tuning parameters $(\widehat{\lambda}_{k,T}, \widehat{\lambda}_{\sigma,T})$ according to the following rules,

align[align omitted — 1,659 chars of source]

The sequential strong rules of tibshirani2012strong is typically employed for computing LASSO-type problems. However, when $k_t$ presents persistent dynamics the sequential strong rules for $\widehat{\bm{\beta}}^{\text{al}}$ fails to screen among truly active and inactive predictors due to estimation bias when the tuning parameters are not small enough, and, as a byproduct, favors the boundary solution $k_t=0.5$. To reach variable selection consistency, it is necessary to enforce the optimizer to stay away from the boundary of the parameter space. Theorem (ref) illustrates the restriction on the permissible coefficient space $\Theta$ in order to achieve the model selection consistency of $\widehat{\bm{\beta}}^{\text{al}}$, i.e.,

equation[equation omitted — 103 chars of source]

where $ \mathcal{A}_T^{\text{al}} := \mathcal{A}_{k,T}^{\text{al}} \cup \mathcal{A}_{\sigma,T}^{\text{al}}$ with $ \mathcal{A}_{k,T}^{\text{al}} := \left\{ (1,i): i\geq 1, \widehat{\beta}_{1,i}^{\text{al}} \neq 0\right\}$ and $\mathcal{A}_{\sigma,T}^{\text{al}} := \left\{ (2,j): j\geq 1, \widehat{\beta}_{2,j}^{\text{al}} \neq 0\right\}$, and $ \mathcal{A} := \mathcal{A}_{k} \cup \mathcal{A}_{\sigma}$ with $ \mathcal{A}_{k} := \left\{ (1,i): i\geq 1, \beta_{1,i}^{o*} \neq 0\right\}$ and $\mathcal{A}_{\sigma} := \left\{ (2,j): j\geq 1, \beta_{2,j}^{o*} \neq 0\right\}$.

theoremUnder the assumptions in Theorem (ref), if there is no $\widehat{\bm{\beta}}^{\text{al}}(\lambda_{k, T}, \lambda_{\sigma, T})$ with $\lambda_{k, T}, \lambda_{\sigma, T}\in O(T^{\frac{1}{2}})$ such that \begin{equation} \det\left( \frac{\partial^2 \mathcal{L}(\bm{\beta}) }{\partial [\bm{\beta}_{\mathcal{A}_k}', \bm{\beta}_{\mathcal{A}_{\sigma}}']' \partial [\bm{\beta}_{\mathcal{A}_k}', \bm{\beta}_{\mathcal{A}_{\sigma}}'] } \middle|_{\bm{\beta} = \widehat{\bm{\beta}}^{al}(\lambda_{k, T}, \lambda_{\sigma, T})} \right) \neq 0, \end{equation} then $\lim_{T\rightarrow\infty} P\left\{ \mathcal{A}_T^{\text{al}} = \mathcal{A}\right\} \neq 1$, where $\text{det}(\cdot)$ is the matrix determinant operator; $w_{k,i}$ and $w_{\sigma,j}$ are set using the MLE in section (ref) such that $\sqrt{T}(\frac{1}{w_{k,i}} - \beta_{1,i}^{*o})=O_p(1) $ and $\sqrt{T}(\frac{1}{w_{\sigma,j}} - \beta_{2,j}^{*o})=O_p(1) $, for $i=1,\ldots, p$, $j=1,\ldots, p+1$.

proofSee Appendix (ref).

Theorem (ref) shows that if not all the truly active predictors are able to enter the regression model with $\lambda_{k, T}, \lambda_{\sigma, T}\in O(T^{\frac{1}{2}})$, then truly inactive predictors start to be selected for compensating for the missing ones since $\lambda_{k, T}, \lambda_{\sigma, T}\in O(T^{\frac{1}{2}})$ and thereby fail $ \widehat{\bm{\beta}}^{\text{al}}$ in the variable selection. The necessary condition in Theorem (ref) tends to be broken when the underlying $\{k_t(\bm{\beta}^{o*})\}$ involves local unit-root predictors. To solve this issue we propose a two-step ALMLE and prove its oracle property.

Two-Step ALMLE

From the previous discussion, we know that ALMLE can be improved if we ensure the estimation to stay away from $\{k_t(\cdot)=0.5\}$ for every $\lambda_{k, T}$. Therefore, we propose a two-step ALMLE, denoted as $\widehat{\bm{\beta}}^{\text{tal}}$, to avoid the local minimizer issue of $\widehat{\bm{\beta}}^{\text{al}}$ by selecting predictors for the shape at the first step and running the ALMLE in (ref) at the second step with the selected $\widehat{\lambda}_{k, T}$ in the first step. Specifically, the two-step ALMLE $\widehat{\bm{\beta}}^{\text{tal}}$ is obtained using the following procedure:

description• {Step 1}: Select the optimal tuning parameter $\widehat{\lambda}_{k, T}\in S_{\lambda_{k,T}}$ using an IC as follows, \begin{align*} AIC:\qquad & \widehat{\lambda}_{k,T} = \underset{\lambda_{k,T}\in S_{\lambda_{k,T}}}{\operatorname{arg}\operatorname{min}}\;\, -2\log(\mathcal{L}(\widehat{\bm{\beta}}^{k,al}(\lambda_{k, T}); \{ y_t \}, \{ \bm{z}_t^* \})) + 2\,\sum_{i=1,\ldots,p} \mathbbm{1}\{\widehat{\beta}^{k,al}_{1,i}\neq 0\} \\ HQ:\qquad & \widehat{\lambda}_{k,T} = \underset{\lambda_{k,T}\in S_{\lambda_{k,T}}}{\operatorname{arg}\operatorname{min}}\;\, -2\log(\mathcal{L}(\widehat{\bm{\beta}}^{k,al}(\lambda_{k, T}); \{ y_t \}, \{ \bm{z}_t^* \})) + 2\,\log(\log(T)) \sum_{i=1,\ldots,p} \mathbbm{1}\{\widehat{\beta}^{\text{k,al}}_{1,i}\neq 0\} \\ \text{BIC:}\qquad & \widehat{\lambda}_{k,T} = \underset{\lambda_{k,T}\in S_{\lambda_{k,T}}}{\operatorname{arg}\operatorname{min}}\;\, -2\log(\mathcal{L}(\widehat{\bm{\beta}}^{\text{k,al}}(\lambda_{k, T}); \{ y_t \}, \{ \bm{z}_t^* \})) + \log(T) \sum_{i=1,\ldots,p} \mathbbm{1}\{\widehat{\beta}^{\text{k,al}}_{1,i}\neq 0\} \,, \end{align*} where $\widehat{\bm{\beta}}^{\text{k,al}}(\lambda_{k, T}):= [\widehat{\beta}^{\text{k,al}}_{1,0},\ldots,\widehat{\beta}^{\text{k,al}}_{1,p},\widehat{\beta}^{\text{k,al}}_{2,0}, 0,\ldots,0]'$ restricts $\widehat{\beta}^{\text{k,al}}_{2,1},\dots,\widehat{\beta}^{\text{k,al}}_{2,p+1}$ to zero and define \begin{equation} \left[\widehat{\beta}^{\text{k,al}}_{1,0},\ldots,\widehat{\beta}^{\text{k,al}}_{1,p}, \widehat{\beta}^{\text{k,al}}_{2,0}\right] = \underset{\beta_{1,0},\ldots,\beta_{1,p},\beta_{2,0}}{\operatorname{arg}\operatorname{min}}\; - \mathcal{L}(\bm{\beta}; \{ y_t \}, \{ \bm{z}_t^* \}) + \lambda_{k, T}\sum_{i=1}^{p} \widetilde{w}_{k,i} |\beta_{1,i}|. \end{equation} • {\textbf{Step 2}}: Select the optimal tuning parameter $\widehat{\lambda}_{\sigma, T} \in S_{\lambda_{\sigma,T}}$ using the IC and $\widehat{\lambda}_{k,T}$ from Step 1 as follows, \begin{align*} \text{AIC:}\qquad & \resizebox{.9\hsize}{!}{$ \widehat{\lambda}_{\sigma,T} = \underset{ \lambda_{\sigma,T}\in S_{\lambda_{\sigma,T}}}{\operatorname{arg}\operatorname{min}}\;\, -2\log(\mathcal{L}(\widehat{\bm{\beta}}^{\text{tal}}( \lambda_{\sigma, T}); \{ y_t \}, \{ \bm{z}_t^* \})) + 2\left( \sum_{i=1,\ldots,p} \mathbbm{1}\{\widehat{\beta}^{\text{tal}}_{1,i}\neq 0\} + \sum_{j=1,\ldots,(p+1)} \mathbbm{1}\{\widehat{\beta}^{\text{tal}}_{2,j}\neq 0\} \right) $} \\ \text{HQ:}\qquad & \resizebox{.9\hsize}{!}{$ \widehat{\lambda}_{\sigma,T} = \underset{ \lambda_{\sigma,T}\in S_{\lambda_{\sigma,T}}}{\operatorname{arg}\operatorname{min}}\;\, -2\log(\mathcal{L}(\widehat{\bm{\beta}}^{\text{tal}}(\lambda_{\sigma, T}); \{ y_t \}, \{ \bm{z}_t^* \})) + 2\log(\log(T))\left( \sum_{i=1,\ldots,p} \mathbbm{1}\{\widehat{\beta}^{\text{tal}}_{1,i}\neq 0\} + \sum_{j=1,\ldots,(p+1)} \mathbbm{1}\{\widehat{\beta}^{\text{tal}}_{2,j}\neq 0\} \right) $} \\ \text{BIC:}\qquad & \resizebox{.9\hsize}{!}{$ \widehat{\lambda}_{\sigma,T} = \underset{ \lambda_{\sigma,T}\in S_{\lambda_{\sigma,T}}}{\operatorname{arg}\operatorname{min}}\;\, -2\log(\mathcal{L}(\widehat{\bm{\beta}}^{\text{al}}(\lambda_{\sigma, T}); \{ y_t \}, \{ \bm{z}_t^* \})) + \log(T)\left( \sum_{i=1,\ldots,p} \mathbbm{1}\{\widehat{\beta}^{\text{tal}}_{1,i}\neq 0\} + \sum_{j=1,\ldots,(p+1)} \mathbbm{1}\{\widehat{\beta}^{\text{tal}}_{2,j}\neq 0\} \right)\,, $} \end{align*} where $\widehat{\bm{\beta}}^{\text{tal}}(\lambda_{\sigma, T}):= [\widehat{\boldsymbol{\beta}}^{\text{tal}'}_{1\cdot},\widehat{\boldsymbol{\beta}}^{\text{tal}'}_{2\cdot}]'=[\widehat{\beta}^{\text{tal}}_{1,0},\ldots,\widehat{\beta}^{\text{tal}}_{1,p},\widehat{\beta}^{\text{tal}}_{2,0},\ldots,\widehat{\beta}^{\text{tal}}_{2,(p+1)}]'$, with $\widehat{\beta}^{\text{tal}}_{1,i}=\widehat{\beta}^{\text{al}}_{1,i}=0, \, \forall (1,i) \not\in \mathcal{A}_T^{k,al}$ and \begin{equation} \resizebox{.95\hsize}{!}{$ \begin{aligned} \left[\left[\widehat{\beta}^{\text{tal}}_{1,i}\right]_{(1,i) \in \{(1,0)\}\cup \mathcal{A}_T^{k,al}}, \widehat{\boldsymbol{\beta}}^{\text{tal}'}_{2\cdot} \right] & = \underset{\left\{\beta_{1,i}|(1,i) \in \{(1,0)\}\cup \mathcal{A}_T^{k,al}\right\}, \boldsymbol{\beta}_{2\cdot}}{\operatorname{arg}\operatorname{min}}\; - \mathcal{L}(\bm{\beta}; \{ y_t \}, \{ \bm{z}_t^* \}) + \widehat{\lambda}_{k, T} \sum_{i=1}^{p} \widetilde{w}_{k,i}|\beta_{1,i}| + \lambda_{\sigma, T}\sum_{j=1}^{p+1} \widetilde{w}_{\sigma,j} |\beta_{2,j}|, \end{aligned} $} \end{equation} where $\mathcal{A}_T^{k,al} := \left\{(1,i):\,i\geq 1, \widehat{\beta}_{1,i}^{k,al} \neq 0\right\}$.

The final two-step ALMLE $\widehat{\bm{\beta}}^{tal}$ is obtained using the optimal tuning parameters $\widehat{\lambda}_{k,T}$ and $\widehat{\lambda}_{\sigma,T}$.

We use two MLEs to set up $\widetilde{w}_{k,i}$ and $\widetilde{w}_{\sigma,j}$ as the two-step ALMLE involves two different likelihood functions in each step. Specifically, we set

equation[equation omitted — 366 chars of source]

where $\widehat{\bm{\beta}}^{\text{mle}} := [\widehat{\beta}^{\text{mle}}_{1,0}, \ldots,\widehat{\beta}^{\text{mle}}_{1,p},\widehat{\beta}^{\text{mle}}_{2,0}, \ldots,\widehat{\beta}^{\text{mle}}_{2,p+1}]$ is the full-model MLE (ref) and $ \widehat{\bm{\beta}}^{\text{k,mle}} := [\widehat{\beta}^{\text{k,mle}}_{1,0}, \ldots,\widehat{\beta}^{\text{k,mle}}_{1,p},\widehat{\beta}^{\text{k,mle}}_{2,0},$ $ 0,\ldots,0]$ is the partial-model MLE defined below

equation[equation omitted — 268 chars of source]

In this way, we choose $\widetilde{w}_{k,i}$ and $\widetilde{w}_{\sigma,j}$ such that truly active predictors are ensured to be selected efficiently with $S_{\lambda_{k,T}}$ and $ S_{\lambda_{\sigma,T}}$ before the truly inactive ones in both Step 1 and Step 2. Therefore, we achieve the oracle property of $\widehat{\bm{\beta}}^{tal}$ as shown in Theorem (ref).

assumption{L1} There exist $\lambda_{k,T} = O(T^{\frac{1}{2} - \gamma_1})$ and $\lambda_{\sigma,T} = O(T^{\frac{1}{2} - \gamma_2})$ with $0<\gamma_1<\frac{1}{2}$ and $0<\gamma_2<\frac{1}{2}$.

assumption{L2} We assume that there exists $\bm{\beta}^{\text{k,o}}:=[\beta^{\text{k,o}}_{1,0},\beta^{\text{k,o}}_{1,1},\ldots,\beta^{\text{k,o}}_{1,p},\beta^{\text{k,o}}_{2,0},0,\ldots,0 ]'\in \{\bm{\beta}\in \mathbb{R}^{2p+3} | \beta_{2,j} = 0, j=1,\ldots,p+1 \}$ such that for any $\epsilon>0$ \begin{equation} \lim_{T\rightarrow\infty} P\left\{ \left| \widehat{\beta}^{k,mle}_{1,i} - \beta^{k,o}_{1,i} \right| > \epsilon \right\} = 0,\quad i=1,\ldots,p; \end{equation} and $\beta^{\text{k,o}}_{1,i}\neq 0$ for any $(1,i)\in\mathcal{A}_k$.

theorem[Oracle Property of $\widehat{\bm{\beta}}^{tal}$] \\ Under Assumptions (ref), (ref) and the assumptions in Theorem (ref), we have that \\ (a) Model selection consistency: \\ \begin{equation} \lim_{T\rightarrow\infty}P\left\{ \mathcal{A}_{T}^{tal} = \mathcal{A}\right\} = 1, \end{equation} where $\mathcal{A}_{T}^{tal}:= \mathcal{A}_{k,T}^{tal}\cup\mathcal{A}_{\sigma,T}^{tal}$ with $\mathcal{A}_{k,T}^{tal}:=\left\{ (1,i): \widehat{\beta}_{1,i}^{tal} \neq 0, i=1,\ldots, p.\right\} $ and $\mathcal{A}_{\sigma,T}^{tal}:=\{ (2,j): \widehat{\beta}_{2,j}^{tal} \neq 0,$ $ j=1,\ldots,p+1.\} $. \\ (a) Limiting distribution of $\widehat{\bm{\beta}}^{tal}$: \\ \begin{equation} \begin{aligned} & \sqrt{T}\left(\widehat{\bm{\beta}}^{tal}_{\mathcal{A}} - \bm{\beta}^{o*}_{\mathcal{A}}\right) \overset{D}{\sim} \Omega^{-1}_{H_\mathcal{A}} \; S_{\psi_{\mathcal{A}}}, \\ & \sqrt{T}\left(\widehat{\bm{\beta}}^{tal}_{\mathcal{A}^c} - \bm{\beta}^{o*}_{\mathcal{A}^c}\right)\rightarrow 0, \end{aligned} \end{equation} as $T\rightarrow\infty$, where $S_{\psi_\mathcal{A}}$ and $\Omega_{H_\mathcal{A}}$ are defined in Assumption (ref) and (ref) under the model specification with only the truly active predictors involved and ordered according to $\mathcal{A}$.

proofSee Appendix (ref).

The superiority of the proposed two-step ALMLE to the ALMLE (ref) is not just in the oracle property when local unit-root predictors are included in the regression model but also in the computing cost. The ALMLE (ref) is computed over a two-dimensional tuning parameter grid in order to select an optimal pair of $(\lambda_{k,T}, \lambda_{\sigma,T}) \in S_{\lambda_{k,T}}\times S_{\lambda_{\sigma,T}}$, while the two-step ALMLE is computed over two separate one-dimensional tuning parameter grids in order to select the optimal $\lambda_{k,T}\in S_{\lambda_{k,T}}$ first and $\lambda_{\sigma,T}\in S_{\lambda_{\sigma,T}}$ after.

Simulation study

We assess the finite sample properties of $ \widehat{\bm{\beta}}^{mle}$ and $\widehat{\bm{\beta}}^{tal}$ from the perspectives of their biases, mean square errors (MSEs) and model selection using four data generating processes (DGPs). These four DGPs are designed to reflect the characteristics of the high-frequency financial data used in Section (ref). First, DGPs are heteroscedastic and the conditional exceedance rates can change over time. Second, DGPs involve predictors which are functions of lagged loss rates characterizing the serial dependence structure in $\{(k_t,\sigma_t)\}$. Third, we consider either stationary or local unit-root predictors or both.

We simulate $\{ l_t\}$ from the following conditional distribution,

equation[equation omitted — 439 chars of source]

where $\{\tau_t\}$ is i.i.d. standard uniform distributed, $F_{t\left(\nu\right)}(\cdot)$ and $F^{-1}_{t\left(\nu\right)}(\cdot)$ denote the distribution and quantile functions of a Student's t distribution with $\nu$ degrees of freedom. The processes of $\{k_t\}$ and $\{ \sigma_t\}$ are specified according to the following specifications:

equation[equation omitted — 1,106 chars of source]

We set $u = F^{-1}_{t(3)}(0.8)$ and $ r_m = 0.05$, but use different $\bm{\phi}$ and $\bm{\beta^o}$ to obtain different degrees of serial dependence.

description• There are five truly active stationary predictors for both $\{k_t\}$ and $\{\sigma_t\}$, namely $\log(|l_{t-1}| + 1 - r_m)$, $z_{1,t-1},\ldots, z_{4,t-1}$. Among truly inactive predictors $z_{5,t-1},\ldots, z_{14,t-1}$, two of them are local unit-root, i.e. $z_{13,t-1}$ and $z_{14,t-1}$, and the others are stationary. \begin{equation} \left\{ \begin{aligned} & \bm{\phi} = [0,0,0,0,0\ldots,0,1,1], \\ & \bm{\beta^o_{1\cdot}} = [-1, 0.3, -0.4, 0.2, 0.6, 0.6, 0, \ldots,0 ]'\,, \\ & \bm{\beta^o_{2\cdot}} = [-1, 0, 0.7, 0.4, 0.3, 0.5, 0.6, 0, \ldots,0 ]'\,, \end{aligned} \right. \end{equation} • As DGP 1 but with the difference that $\beta_{2,1}^o$ is changed to nonzero, and hence $\log(\sigma_{t-1})$ is now truly active. We set $\beta_{2,1}^o=0.7$ and keep the true values of the other coefficients unchanged. \begin{equation} \left\{ \begin{aligned} & \bm{\phi} = [0,0,0,0,0\ldots,0,1,1], \\ & \bm{\beta^o_{1\cdot}} = [-1, 0.3, -0.4, 0.2, 0.6, 0.6, 0, \ldots,0 ]'\,, \\ & \bm{\beta^o_{2\cdot}} = [-1, 0.7, 0.7, 0.4, 0.3, 0.5, 0.6, 0, \ldots,0 ]'\,, \end{aligned} \right. \end{equation} • As DGP 1 but with the difference that $\phi_{4} = 1$ and $(\beta_{1,5}^o,\beta_{2,6}^o) = (\frac{0.6}{\sqrt{T}}, \frac{0.6}{\sqrt{T}})$. \begin{equation} \left\{ \begin{aligned} & \bm{\phi} = [0,0,0,1,0\ldots,0,1,1], \\ & \bm{\beta^o_{1\cdot}} = [-1, 0.3, -0.4, 0.2, 0.6, \frac{0.6}{\sqrt{T}}, 0, \ldots,0 ]'\,, \\ & \bm{\beta^o_{2\cdot}} = [-1, 0, 0.7, 0.4, 0.3, 0.5, \frac{0.6}{\sqrt{T}}, 0, \ldots,0 ]'\,, \end{aligned} \right. \end{equation} • As DGP 3 but with the difference that $\log(\sigma_{t-1})$ is truly active. We set $\beta_{2,1}^o=0.7$ and and keep the true values of the other coefficients unchanged. \begin{equation} \left\{ \begin{aligned} & \bm{\phi} = [0,0,0,1,0\ldots,0,1,1], \\ & \bm{\beta^o_{1\cdot}} = [-1, 0.3, -0.4, 0.2, 0.6, \frac{0.6}{\sqrt{T}}, 0, \ldots,0 ]'\,, \\ & \bm{\beta^o_{2\cdot}} = [-1, 0.7, 0.7, 0.4, 0.3, 0.5, \frac{0.6}{\sqrt{T}}, 0, \ldots,0 ]'\,, \end{aligned} \right. \end{equation}

In each simulation, we obtain a sample $\{ l_t\}_{t=1}^T$ of $T$ observations, and extract the excess time series $\{ y_t=\max(l_t - u, 0)\}$ using the true threshold $u$. We standardize the predictors using their empirical standard deviations. We then fit the full model specification (ref) to $\{y_t\}$ using standardized predictors, estimating the model parameter by $\widehat{\bm{\beta}}^{mle}$ and $ \widehat{\bm{\beta}}^{tal}$. Bias and mean squared error (MSE) are then computed as

equation[equation omitted — 306 chars of source]

equation[equation omitted — 294 chars of source]

where $\#\bm{\beta}^o$ denotes the number of parameters in $\bm{\beta}^o$, $\widehat{s}_{i,j}$ denotes the empirical standard deviation of the $(i,j)$-th predictor, and $s_{i,j}=1$ for $I(0)$ predictors and $s_{i,j}=\sqrt{T}$ for $I(1)$ predictors.

{Table (ref) presents the average absolute bias and average MSE of the coefficient estimates obtained over 100 replications. These results show that $\widehat{\bm{\beta}}^{mle}$ and $\widehat{\bm{\beta}}^{tal}$ have decreasing biases and MSEs when $T$ increases, coherently with the theoretical results presented in Sections (ref) and (ref). Moreover $\widehat{\bm{\beta}}^{tal}$ under BIC always has the lowest bias and MSE across the DGPs, supporting the use of $\widehat{\bm{\beta}}^{tal}$ with BIC in the empirical section. Boxplots for the bias in Figure (ref) support these conclusions.

Table (ref) presents the variable selection results for both $\widehat{\bm{\beta}}^{tal}$ and $\widehat{\bm{\beta}}^{mle}$. Note that for the latter, we perform variable selection based on the significance of the t-statistics associated to the candidate predictors. To measure the ability to select the correct predictors, we assess the average selection rates of truly active and inactive predictors for both the shape and scale parameters. Moreover, we compute the correct classification rate (CCR) of each estimator, i.e. the proportion of selected truly active and unselected truly inactive predictors on the total candidate predictors. Results in Table (ref) show that variable selection improves as $T$ increases for each estimator. For $\widehat{\bm{\beta}}^{mle}$ the average selection rates of truly inactive stationary predictors approach the significance level $\alpha=0.05$, whereas the average selection rates of truly inactive local unit-root predictors are much higher than $\alpha=0.05$, for both $k$ and $\sigma$, and regardless of the DGP. These results are coherent with the asymptotic results derived in Section (ref), and echo the size distortion concerns of using t-tests to select non-stationary predictors discussed in Section (ref). Remarkably, the average selection rates of truly active predictors for $\widehat{\bm{\beta}}^{mle}$ are much lower than those for $\widehat{\bm{\beta}}^{tal}$. Moreover, we see that the power of t-tests performed with $\widehat{\bm{\beta}}^{mle}_{1\cdot}$ is lower than the one for $\widehat{\bm{\beta}}^{mle}_{2\cdot}$ due to the uncertainty in the estimation of $\widehat{\bm{\beta}}^{mle}_{1\cdot}$. Finally, Table (ref) shows that $\widehat{\bm{\beta}}^{tal}$ with BIC always has the highest CCR and produces the most accurate selection regardless the DGP, supporting the use of $\widehat{\bm{\beta}}^{tal}$ with BIC for the empirical application.

table[table omitted — 2,414 chars of source]

figure[figure omitted — 1,007 chars of source]

sidewaystable[htbp] \caption{Average selection rate across 100 replication for truly active (t.p.) and truly inactive (f.p.) stationary (I(0)) and local unit-root (I(1)) predictors in the shape ($k$) and scale ($\sigma$) parameters, correct classification rate (CCR), and selection rates of $\log(\sigma_{t-1})$.} \resizebox{0.9\textwidth}{!}{ \begin{tabular}{cc|l|l|cccccccccc} \toprule \toprule DGPs& T & estimators & selection criteria & \multicolumn{1}{l}{ t.p.(k) of $I(0)s$\;\;\;\;} & \multicolumn{1}{l}{f.p.(k) of $I(0)s$\;\;\;\;} & \multicolumn{1}{l}{t.p.(k) of $I(1)s$\;\;\;\;} & \multicolumn{1}{l}{f.p.(k) of $I(1)s$\;\;\;\;} & \multicolumn{1}{l}{t.p.($\sigma$) of $I(0)s$\;\;\;\;} & \multicolumn{1}{l}{f.p.($\sigma$) of $I(0)s$\;\;\;\;} & \multicolumn{1}{l}{t.p.($\sigma$) of $I(1)s$\;\;\;\;} & \multicolumn{1}{l}{f.p.($\sigma$) of $I(1)s$\;\;\;\;} & \multicolumn{1}{l}{ CCR\;\;\;\;} & \multicolumn{1}{l}{ $\log(\sigma_{t-1})$} \\ \midrule \multirow{12}[6]{*}{DGP 3 } & \multirow{4}[2]{*}{25,000} & $\widehat{\bm{\beta}}^{mle}$ & t test ($\alpha=0.05$) & 0.546 & 0.063 & $-$ & 0.110 & 1.000 & 0.110 & $-$ & 0.215 & 0.858 & 0.060 \\ & & \multirow{3}[1]{*}{$\widehat{\bm{\beta}}^{tal}$} & AIC & 0.998 & 0.160 & $-$ & 0.465 & 1.000 & 0.156 & $-$ & 0.220 & 0.869 & 0.150 \\ & & & HQ & 0.996 & 0.123 & $-$ & 0.420 & 1.000 & 0.024 & $-$ & 0.060 & 0.930 & 0.020 \\ & & & BIC & 0.984 & 0.076 & $-$ & 0.350 & 1.000 & 0.006 & $-$ & 0.010 & 0.953 & 0.000 \\ \cmidrule{2-14} & \multirow{4}[2]{*}{50,000} & $\widehat{\bm{\beta}}^{mle}$ & t test ($\alpha=0.05$) & 0.718 & 0.050 & $-$ & 0.085 & 1.000 & 0.112 & $-$ & 0.185 & 0.892 & 0.050 \\ & & \multirow{3}[1]{*}{$\widehat{\bm{\beta}}^{tal}$} & AIC & 0.998 & 0.100 & $-$ & 0.450 & 1.000 & 0.147 & $-$ & 0.160 & 0.892 & 0.110 \\ & & & HQ & 0.998 & 0.075 & $-$ & 0.410 & 1.000 & 0.032 & $-$ & 0.040 & 0.942 & 0.020 \\ & & & BIC & 0.998 & 0.053 & $-$ & 0.330 & 1.000 & 0.004 & $-$ & 0.015 & 0.963 & 0.000 \\ \cmidrule{2-14} & \multirow{4}[2]{*}{100,000} & $\widehat{\bm{\beta}}^{mle}$ & t test ($\alpha=0.05$) & 0.788 & 0.045 & $-$ & 0.130 & 1.000 & 0.109 & $-$ & 0.220 & 0.900 & 0.040 \\ & & \multirow{3}[1]{*}{$\widehat{\bm{\beta}}^{tal}$} & AIC & 0.988 & 0.045 & $-$ & 0.435 & 1.000 & 0.154 & $-$ & 0.200 & 0.901 & 0.160 \\ & & & HQ & 0.988 & 0.034 & $-$ & 0.385 & 1.000 & 0.030 & $-$ & 0.040 & 0.953 & 0.020 \\ & & & BIC & 0.988 & 0.024 & $-$ & 0.345 & 1.000 & 0.001 & $-$ & 0.010 & 0.969 & 0.000 \\ \midrule \multirow{12}[6]{*}{DGP 4} & \multirow{4}[2]{*}{25,000} & $\widehat{\bm{\beta}}^{mle}$ & t test ($\alpha=0.05$) & 0.588 & 0.071 & $-$ & 0.100 & 1.000 & 0.101 & $-$ & 0.195 & 0.870 & 1.000 \\ & & \multirow{3}[1]{*}{$\widehat{\bm{\beta}}^{tal}$} & AIC & 0.982 & 0.414 & $-$ & 0.665 & 1.000 & 0.176 & $-$ & 0.155 & 0.792 & 1.000 \\ & & & HQ & 0.966 & 0.213 & $-$ & 0.525 & 1.000 & 0.044 & $-$ & 0.040 & 0.892 & 1.000 \\ & & & BIC & 0.916 & 0.086 & $-$ & 0.450 & 1.000 & 0.004 & $-$ & 0.005 & 0.934 & 1.000 \\ \cmidrule{2-14} & \multirow{4}[2]{*}{50,000} & $\widehat{\bm{\beta}}^{mle}$ & t test ($\alpha=0.05$) & 0.718 & 0.055 & $-$ & 0.110 & 1.000 & 0.085 & $-$ & 0.170 & 0.900 & 1.000 \\ & & \multirow{3}[1]{*}{$\widehat{\bm{\beta}}^{tal}$} & AIC & 0.996 & 0.310 & $-$ & 0.630 & 1.000 & 0.145 & $-$ & 0.140 & 0.832 & 1.000 \\ & & & HQ & 0.992 & 0.178 & $-$ & 0.535 & 1.000 & 0.038 & $-$ & 0.025 & 0.907 & 1.000 \\ & & & BIC & 0.982 & 0.109 & $-$ & 0.500 & 1.000 & 0.001 & $-$ & 0.005 & 0.936 & 1.000 \\ \cmidrule{2-14} & \multirow{4}[2]{*}{100,000} & $\widehat{\bm{\beta}}^{mle}$ & t test ($\alpha=0.05$) & 0.828 & 0.039 & $-$ & 0.090 & 1.000 & 0.094 & $-$ & 0.190 & 0.920 & 1.000 \\ & & \multirow{3}[1]{*}{$\widehat{\bm{\beta}}^{tal}$} & AIC & 1.000 & 0.248 & $-$ & 0.560 & 1.000 & 0.134 & $-$ & 0.100 & 0.859 & 1.000 \\ & & & HQ & 1.000 & 0.111 & $-$ & 0.435 & 1.000 & 0.036 & $-$ & 0.040 & 0.931 & 1.000 \\ & & & BIC & 0.998 & 0.060 & $-$ & 0.345 & 1.000 & 0.000 & $-$ & 0.000 & 0.962 & 1.000 \\ \midrule \multirow{12}[6]{*}{DGP 5} & \multirow{4}[2]{*}{25,000} & $\widehat{\bm{\beta}}^{mle}$ & t test ($\alpha=0.05$) & 0.458 & 0.065 & 0.280 & 0.140 & 1.000 & 0.120 & 1.000 & 0.225 & 0.832 & 0.080 \\ & & \multirow{3}[1]{*}{$\widehat{\bm{\beta}}^{tal}$} & AIC & 0.998 & 0.161 & 0.930 & 0.410 & 1.000 & 0.152 & 1.000 & 0.175 & 0.874 & 0.190 \\ & & & HQ & 0.998 & 0.098 & 0.900 & 0.365 & 1.000 & 0.039 & 1.000 & 0.045 & 0.934 & 0.050 \\ & & & BIC & 0.998 & 0.083 & 0.890 & 0.350 & 1.000 & 0.006 & 1.000 & 0.005 & 0.950 & 0.000 \\ \cmidrule{2-14} & \multirow{4}[2]{*}{50,000} & $\widehat{\bm{\beta}}^{mle}$ & t test ($\alpha=0.05$) & 0.620 & 0.058 & 0.390 & 0.135 & 1.000 & 0.111 & 1.000 & 0.215 & 0.862 & 0.060 \\ & & \multirow{3}[1]{*}{$\widehat{\bm{\beta}}^{tal}$} & AIC & 0.993 & 0.085 & 0.990 & 0.355 & 1.000 & 0.151 & 1.000 & 0.175 & 0.899 & 0.150 \\ & & & HQ & 0.993 & 0.058 & 0.980 & 0.285 & 1.000 & 0.036 & 1.000 & 0.020 & 0.954 & 0.030 \\ & & & BIC & 0.993 & 0.044 & 0.970 & 0.260 & 1.000 & 0.003 & 1.000 & 0.005 & 0.969 & 0.000 \\ \cmidrule{2-14} & \multirow{4}[2]{*}{100,000} & $\widehat{\bm{\beta}}^{mle}$ & t test ($\alpha=0.05$) & 0.755 & 0.036 & 0.510 & 0.095 & 1.000 & 0.092 & 1.000 & 0.180 & 0.899 & 0.040 \\ & & \multirow{3}[1]{*}{$\widehat{\bm{\beta}}^{tal}$} & AIC & 0.993 & 0.033 & 1.000 & 0.300 & 1.000 & 0.124 & 1.000 & 0.170 & 0.924 & 0.140 \\ & & & HQ & 0.993 & 0.025 & 1.000 & 0.220 & 1.000 & 0.028 & 1.000 & 0.040 & 0.968 & 0.010 \\ & & & BIC & 0.993 & 0.020 & 1.000 & 0.195 & 1.000 & 0.000 & 1.000 & 0.005 & 0.981 & 0.000 \\ \midrule \multirow{12}[5]{*}{DGP 6} & \multirow{4}[2]{*}{25,000} & $\widehat{\bm{\beta}}^{mle}$ & t test ($\alpha=0.05$) & 0.480 & 0.058 & 0.290 & 0.095 & 1.000 & 0.091 & 1.000 & 0.215 & 0.852 & 1.000 \\ & & \multirow{3}[1]{*}{$\widehat{\bm{\beta}}^{tal}$} & AIC & 0.890 & 0.304 & 0.970 & 0.640 & 1.000 & 0.148 & 1.000 & 0.135 & 0.818 & 1.000 \\ & & & HQ & 0.873 & 0.219 & 0.970 & 0.535 & 1.000 & 0.039 & 1.000 & 0.020 & 0.880 & 1.000 \\ & & & BIC & 0.828 & 0.144 & 0.950 & 0.445 & 1.000 & 0.004 & 1.000 & 0.000 & 0.909 & 1.000 \\ \cmidrule{2-14} & \multirow{4}[2]{*}{50,000} & $\widehat{\bm{\beta}}^{mle}$ & t test ($\alpha=0.05$) & 0.683 & 0.070 & 0.370 & 0.095 & 1.000 & 0.098 & 1.000 & 0.195 & 0.877 & 1.000 \\ & & \multirow{3}[1]{*}{$\widehat{\bm{\beta}}^{tal}$} & AIC & 0.940 & 0.241 & 0.990 & 0.600 & 1.000 & 0.175 & 1.000 & 0.180 & 0.834 & 1.000 \\ & & & HQ & 0.930 & 0.140 & 0.990 & 0.500 & 1.000 & 0.034 & 1.000 & 0.045 & 0.911 & 1.000 \\ & & & BIC & 0.903 & 0.076 & 0.990 & 0.425 & 1.000 & 0.011 & 1.000 & 0.015 & 0.936 & 1.000 \\ \cmidrule{2-14} & \multirow{4}[1]{*}{100,000} & $\widehat{\bm{\beta}}^{mle}$ & t test ($\alpha=0.05$) & 0.813 & 0.048 & 0.450 & 0.070 & 1.000 & 0.069 & 1.000 & 0.155 & 0.914 & 1.000 \\ & & \multirow{3}[0]{*}{$\widehat{\bm{\beta}}^{tal}$} & AIC & 0.918 & 0.095 & 0.990 & 0.510 & 1.000 & 0.125 & 1.000 & 0.090 & 0.894 & 1.000 \\ & & & HQ & 0.918 & 0.036 & 0.990 & 0.385 & 1.000 & 0.033 & 1.000 & 0.025 & 0.945 & 1.000 \\ & & & BIC & 0.913 & 0.023 & 0.990 & 0.340 & 1.000 & 0.001 & 1.000 & 0.005 & 0.960 & 1.000 \\ \bottomrule \bottomrule \end{tabular} }

Empirical Study

We study the high-frequency excess loss distributions of nine large liquid U.S. stocks: American Express (AXP), Boeing (BA), General Electric (GE), Home Depot (HD), IBM, Johnson and Johnson (JNJ), JPMorgan Chase (JPM), Coca-Cola (KO), and ExxonMobil (XOM). Our data covers all transactions observed from January 2006 to December 2014. Market uncertainty and liquidity being elusive concepts, we study their impact on the excess loss distribution using as predictors several high-frequency volatility and liquidity indicators, and select the most appropriate ones with the two-step ALMLE developed in Section (ref). We perform an in-sample analysis providing an economic interpretation for the impact of the selected predictors on the excess loss distribution, and an out-of-sample VaR forecast analysis to assess the goodness of fit of the predicted excess loss distribution.

Variables description

The raw intraday data of the studied stocks contain transaction timestamps in milliseconds, transaction prices per share, and transaction volume in shares for each trade. We cleaned the raw data according to standard procedures in brownlees2006financial and barndorff2009realized. Since transaction data are irregularly-spaced, we need to define an equally-spaced grid at a fixed frequency to analyse losses with our model. We choose to analyse losses at the five minute frequency. Let $P_{t,i}$ be the transaction price of the $i$-th trade in the $t$-th five minute interval, and let $V_{t,i}$ be the corresponding quantity of traded shares, with $0\leq i \leq n_t$ where $n_t$ is the number of trades in the $t$-th five minute interval and $0<t\leq T$. We define 5-min prices, $P_t$, as the median transaction price in the $t$-th five minute interval, and compute 5-min losses as the negative $t$-th return, $R_{t}:= \log( P_{t} ) - \log(P_{t-1} )$. To obtain the time series of excess losses we consider a dynamic threshold accounting for the time-varying behavior of losses at high-frequency. Specifically, the threshold $u_t$ at time $t$ is defined as the $90\%$-quantile of the losses observed over the period $\left(t-1,t-h\right)$, with $h>1$ the moving window size. We consider 12 possible values of $h$ ranging from one week to twelve weeks.

Liquidity refers to the ability to trade large volume of a financial instrument with low price impact, cost and postponement. As liquidity can be decomposed into different dimensions (Harris et al., 1990), we consider several liquidity indicators as possible predictors. Similarly, to characterize market uncertainty we consider several indicators for the observed dispersion of transaction prices. Moreover, to disentangle the impact of trading activity at different frequencies, we build our set of candidate predictors considering both information within the $t$-th five minute interval and across neighbourhoods of the $t$-th five minute interval. Let $P_{t,BU}:=[P_{t,1},\ldots,P_{t,n_t}]'$ and $R_{t,BU} :=[R_{t,1},\ldots,R_{t,n_t}]'$ be the vectors of traded prices and trade returns observed within the $t$-th five minute interval, with $R_{t,i}:= \log(P_{t,i}) - \log(P_{t,i-1})$. Let $T_w$ be a neighborhood size, and define $P_{t,T_w}:= [P_t,\ldots, P_{t-T_w+1}]'$ the vector of 5-min prices within a neighborhood of size $T_w$ and $ R_{t,T_w}: = \log(P_{t,T_w}) - \log(P_{t-1,T_w})$ the corresponding vector of returns. Let $\text{dur}_{t,i}$ denote the execution duration of the $i$-th transaction in the $t$-th five minute interval, i.e. the time difference between the order executed time and order placed time. Table (ref) lists the liquidity predictors we consider in the analysis. They are classified according to their frequency, i.e. within or across the five minute interval, and by their nature of price impact or spread proxies goyenko2009liquidity or volatility of liquidity measures. Table (ref) lists the volatility predictors we consider in the analysis and are classified according to the frequency at which they are computed, i.e., within or across the five minute interval.

table[table omitted — 3,628 chars of source]

table[table omitted — 1,080 chars of source]

In-sample estimates

We divide each time series into an in-sample period covering the first 90% of the observations and an out-of-sample period spanning the last 10% of the sample. We model the excess losses $\{y_t\}$ of each stock with the time-varying GPD regression model in (ref)-(ref), using the variables defined in Tables (ref) and (ref), with $T_w\in\{2,6,12\}$, as possible predictors in both scale and shape parameters. Coefficient estimates obtained with the two-step ALMLE are presented in Tables (ref) and (ref) for the shape and scale parameters, respectively.

Results for the shape parameter in Table (ref) show that estimated coefficients have almost always the same sign across the stocks. As to liquidity predictors, we find that price impact proxies are selected for almost all the stocks, suggesting that they better capture liquidity effects on extreme losses. In particular, TV and TQ display positive coefficients while AM and EAM display negative coefficients, entailing that larger extreme losses are associated with high levels of liquidity in the last five minutes. Although counter-intuitive at first, this result is very interesting when read together with the other selected variables. As to the volatility of liquidity, we notice that RTVV$(T_w=6)$ and RTQV$(T_w=6)$ are selected across most of the stocks and display large and positive coefficients, indicating that extreme losses tend to be larger during periods of high volatility of liquidity. Almost for every stock, we select the ratio MRV2RV($T_w=12$), essentially capturing the impact of the volatility of volatility or jump risk on extreme losses, and associate a positive coefficient to it, conveying the idea that extreme losses tend to be larger during periods of high uncertainty. Altogether these results are coherent with the findings in brogaard2018, i.e. that market markers amplify extreme price movements while withdrawing from the market after large uncertainty shocks that caused their liquidity supply to be outstripped by liquidity demand.

Table (ref) shows that more variables are selected for the scale parameter but their pattern is less stable across stocks. In general, we notice that the autoregressive component contributes to the dynamics, and that the realized volatility predictor computed within the five-minute interval is always selected and displays positive coefficient. This is coherent with the fact that the scale parameter captures the time-varying heteroscedasticity in the data.

For comparison purposes, we report estimated regression coefficients for the shape and scale parameters obtained with MLE in Tables (ref)-(ref). All of the estimated coefficients are nonzero and we cannot compute the corresponding standard errors because the obtained Fisher information matrix of the MLE is not positive definitive. This makes variable interpretation very difficult if not impossible.

table[table omitted — 3,598 chars of source]

table[table omitted — 3,858 chars of source]

table[table omitted — 13,581 chars of source]

table[table omitted — 13,923 chars of source]

Out-of-sample VaR forecast

The coefficient estimates $\widehat{\bm{\beta}}$ obtained on the in-sample period are used to compute a one-step ahead VaR prediction in the out-of-sample period. Specifically, the VaR of each stock at a risk level $\alpha$ at time $t$ given $\bm{x}_{t-1}$ and $\widehat{\bm{\beta}}$ is obtained as

equation[equation omitted — 327 chars of source]

where $F_t(\widehat{u}_t)$ is the probability of exceeding the threshold $\widehat{u}_t$ and is fixed to 90%. The coverage rate of $\{\widehat{\text{VaR}_t}(\alpha)\}_{t=T_{is} + 1}^{T_{is} + T_{os}}$ for the out-of sample period is obtained as follows,

equation[equation omitted — 155 chars of source]

Table (ref) shows the coverage rate of $\{\widehat{\text{VaR}_t}(\alpha)\}_{t=T_{is} + 1}^{T_{is} + T_{os}}$ at the risk level $\alpha$ for various $\alpha\in[90\%,100\%)$. We resort to the Kolmogorov–Smirnov (K-S) test to test the goodness of fit of the predicted GPD over the out-of-sample period, i.e., we test whether $\{\widehat{F}(y_t | y_t >0 ) = \text{GPD}(y_t ;\bm{x}_{t-1},\widehat{\bm{\beta}})\}$ follows a standard uniform distribution. The p-values of the K-S tests in Table (ref) indicate that we reject the regression model on three stocks out of nine at the 1% significance level.

table[table omitted — 1,681 chars of source]

Conclusion

This paper proposes a novel extreme value regression framework to study the dynamics of high-frequency tail risk. The proposed model allows for both stationary and local unit-root predictors to capture the persistence of high-frequency extreme losses. We propose a two-step regularized approach to perform automatic variable selection, and establish the oracle property of the corresponding ALMLE in selecting stationary and local unit-root predictors. We use the proposed approach to investigate the predictive content of 42 liquidity and volatility indicators on the distribution of extreme losses for nine large liquid U.S. stocks. Our variable selection procedure reveals that the severity of tail risk is strongly associated to low price impact in periods of high volatility of liquidity and volatility of volatility. These findings can contribute to timely alert high-frequency traders of rising risk levels and facilitate improvements of their algorithmic trading practices for financial risk management. Moreover, it provides incentives for market markers to absorb liquidity demand in periods of instability. Finally, it suggests a set of predictors to regulators investigating trading activities that can help defining proper regulation guidelines and safeguard financial stability.