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.
114,477 characters · 21 sections · 83 citation commands
Binary Outcome Models with Extreme Covariates: Estimation and Prediction
\ifsubmission\else\fi
Binary outcome models are widely used in both empirical microeconomics and empirical macroeconomics research. For example, microeconomic studies may be interested in the determinants of individuals' labor force participation decisions, and macroeconomic analyses may seek to forecast the probability of recessions or country defaults. Recently, extreme events, such as the Covid-19 pandemic and its aftermath, sustained periods of high inflation, and increasingly frequent extreme weather, have highlighted the importance of studying the effect of extreme covariates on binary outcomes.
As an illustration, consider a simple example with cross-sectional data, where we have a random sample of $\{Y,X\}$ with binary outcome $Y\in \{0,1\}$ and continuous covariate $X\in \mathbb{R}$. Here and throughout, we follow the convention that uppercase letters denote random variables, while lowercase letters represent their realized values. Depending on the empirical context, one may be interested in the conditional probability
the partial effect $\partial \pi(x) /\partial x$, and the elasticity, when $x$ takes extreme values. Given a random sample with $N$ observations, we characterize the extremeness by letting $x\rightarrow\infty$ as $N\rightarrow\infty$.\footnote{Here we focus on the right tail without loss of generality. In practice, one can conduct a similar analysis for the left tail, and we allow for different tail thickness in each tail.}\textsuperscript{,}\footnote{In finite sample, our approach still shows significant improvement for moderately large $x$, roughly below the 20th or above 80th percentile of the $X$ distribution, and the corresponding $\pi(x)$ may not be very close to 0 or 1 either: see Figure (ref), for example.} In many cases, $\pi(x)\rightarrow 0$ or 1 as $x\rightarrow\infty$, so we combine both scenarios and define a unified measure of extreme elasticity as
Further details are provided in Proposition (ref) and Remark (ref). We will use sovereign debt and country default as a running example for intuitive explanations. In this context, the outcome $Y$ indicates whether the country defaults or not, the covariate $X$ is the debt-to-GDP ratio, and the conditional probability $\pi(x)$ represents the probability of default when the debt-to-GDP ratio is particularly high, which can be viewed as a counterfactual probability or a predictive probability.
Existing methods face limitations when dealing with extreme values, a common feature in economic data that manifests as heavy tails gabaix2009,Gabaix2016JEP. Parametric methods, such as Logit and Probit, assume threshold-crossing models with thin-tailed error distributions, which can lead to significant misspecification bias, particularly in the tails. Conversely, nonparametric methods, such as kernel, sieve, and spline estimators, may encounter difficulties due to limited information in the tail, resulting in highly inefficient estimators with large variances. This inefficiency can lead to imprecise forecasts, especially for extreme values.
To address these challenges, we propose a novel semiparametric approach based on Bayes' theorem and RV functions, which offer a flexible framework that encompasses a wide range of distributions. For heavy-tailed distributions, the RV condition naturally leads to a Pareto approximation in the tail (see Section (ref) for details), while making no parametric assumptions on the relationship between covariates and outcomes beyond the tail region.\footnote{To see the semiparametric nature of our tail estimator, note that traditional kernel methods focus on local observations close to a specific $x$ value of interest, leaving the model far from $x$ unspecified. Similarly, our tail estimator focuses on a specific local region, namely the extreme tail where $x\rightarrow\infty$, leaving the model for the middle part unspecified.} This is crucial because patterns in the tail and middle can differ substantially, as illustrated in the log-log plots in Figures (ref) and (ref). Specifically, by Bayes' rule, we have that
where $f_{X|Y}\left(\cdot|y\right) $ denotes the conditional density function of $X$ given $Y=y$. When $x$ takes extreme values, the RV properties allow us to approximate $f_{X|Y}\left( \cdot |y\right) $ by $x^{-\alpha^{(y)}-1}$ up to some constant, where $\alpha^{(y)}$ is the Pareto exponent for $y\in\{0,1\}$. Then, the extreme elasticity $\delta(x)$ converges to $-\left|\alpha^{(1)}-\alpha^{(0)}\right|$, as $x\rightarrow \infty $. We establish formal asymptotic results in a more general setting with additional non-extreme covariates $Z=z$, where we hold $z$ fixed and let $x\rightarrow\infty$. Note that given the limited data available in the tail for extreme $x$ values, the convergence rate of the tail estimator is slower than the standard $\sqrt N$-rate.
Furthermore, our method not only addresses the limited information in the tail but also effectively handles unobserved unit-specific tail parameters and unit-specific RV functions in panel data. The heterogeneity in the RV functions naturally accommodates unit-specific scale parameters as well. We show that, under regularity conditions, our objective function is asymptotically equivalent to that of a panel Logit regression on tail observations, using the log extreme covariate as a regressor. This equivalence simplifies the estimation procedure and makes the method more convenient for empirical studies.
Here we consider panel data with large $N$ and either small or large $T$. For small $T$, we can use the conditional MLE to eliminate unit-specific tail effects. The proof nontrivially extends WangTsai2009 to panel data with independent but not necessarily identically distributed (i.n.i.d.) observations. For large $T$, we can directly estimate the unit-specific parameters using existing methods for large-$T$ panel Logit and further provide unit-specific forecasts. Finally, we expand our discussion to dynamic panel data models, incorporating lagged outcomes into the analysis.
We assess the finite-sample performance of our tail estimator through Monte Carlo simulations under various specifications. Experiment 1 focuses on estimation accuracy in cross-sectional data, and Experiment 2 examines out-of-sample forecasts in panel data with large $N$ and large $T$. Results show that our tail estimator outperforms parametric and nonparametric alternatives, particularly in capturing heavy-tailed behavior, reducing misspecification bias, and providing more accurate forecasts.
In our empirical application, we analyze the impact of local housing price declines on the riskiness of small banks, which played an important role in the 2007--2008 financial crisis. We construct a panel dataset of loan charge-off rates of small banks, local housing prices, and unemployment rates, and again find that the tail estimator yields the most accurate pseudo out-of-sample forecasts among all methods considered. In addition, this empirical analysis also helps reveal interesting heterogeneity patterns in riskiness among small banks, varying across geographic regions and bank characteristics.
\paragraph{Related literature.} Our work draws on a wide range of econometric literature, including binary outcome models, panel data, extreme value theory, and forecasting.
First, our work builds on a large literature on the binary outcome models, especially those for panel data. However, we focus on the analysis of extreme events instead of mid-sample properties as in the existing literature.
In binary cross-sectional models, various methods have been well studied in the literature, including parametric approaches (e.g., Logit and Probit), nonparametric approaches (e.g., matzkin1992nonparametric), and semiparametric approaches, such as the single-index model klein1993efficient and the maximum score estimator manski1985semiparametric,horowitz1992smoothed. Also see horowitz2001binary for a review. Our approach falls within the semiparametric framework, approximating the heavy tails using RV functions while making no parametric assumptions in the middle. Notably, much of this existing literature assumes threshold-crossing models, whereas our approach does not require such a structure: see Section (ref) for a detailed comparison.
In contrast to common methods that directly work with $\mathbb{P}(Y=1|X=x)$, we employ Bayes' theorem to “reverse” the conditioning set. Similar transformations based on Bayes' theorem have also been applied in the discriminant analysis amemiya1985advanced with Gaussian $X|Y=y$, as well as in the semiparametric single-index model in klein1993efficient. In our context of extreme events, this transformation offers a novel way to disentangle tail behavior. Unlike those previous works focusing on mid-sample properties, we analyze tail properties with only tail observations, necessitating the development of nonstandard asymptotic theory.
In binary panel data models, the unobserved unit-specific heterogeneity can be challenging to handle due to the incidental parameter problem in the presence of nonlinear model structures. In addition to parametric, nonparametric, and semiparametric approaches, the binary panel literature can also be categorized by assumptions on unit-specific heterogeneity. See Chapter 15 in Wooldridge2010 and Chapter 7 in hsiao2022analysis for textbook discussions. In a random effects setup, the unit-specific heterogeneity is assumed to be drawn from a common underlying distribution. In a correlated random effects setup, this distribution could depend on observed covariates. In a fixed effects setup, the unit-specific heterogeneity is treated as fixed parameters unique to each unit without imposing any distributional assumptions, implicitly allowing for an arbitrary correlation between the heterogeneity and covariates. For logistic errors, Chamberlain1980, among others, employs the conditional MLE to eliminate the unit-specific intercept and estimate the common parameters; for general errors, Manski1987 extends the maximum score estimator to panel data.
Since we allow for heterogeneity in the RV functions without explicitly modeling it, our approach aligns with a fixed effects framework, where unit-specific tail thickness and RV functions can flexibly depend on the additional covariates $Z$. We show that our objective function is asymptotically equivalent to a panel Logit regression on tail observations with the log extreme covariate as a regressor. This equivalence facilitates a comparable conditional MLE analysis to estimate common parameters in small-$T$ panels. For large-$T$ panels, bias correction methods, such as those developed by fernandez2018fixed and stammann2016estimating, can be utilized to further estimate unit-specific parameters. Additionally, we extend our approach to dynamic binary panel data models, where the conditional MLE construction is related to the work of honore2000panel.
Second, our work contributes to the heavy tail and extreme value theory literature. For a comprehensive review, please refer to de2006extreme and gabaix2009,Gabaix2016JEP. The literature typically focuses on continuous outcomes with heavy tails, such as the classic hill1975 estimator, as well as estimators proposed by smith1987MLE and gabaix2011rank. See recent reviews by GomesGuillou2015review and Fedotenkov2020review for over a hundred estimators. WangTsai2009 incorporate covariates and develop a tail index regression model, and Nicolau2023 further extend this framework to time-series data under strongly mixing conditions. In contrast to the existing literature, we focus on binary outcomes and utilize Bayes' theorem to construct a novel estimator. To incorporate covariates, we build upon the tail index regression model and significantly extend the theoretical framework to handle panel data and i.n.i.d.\ random variables.
Third, our work further contributes to the recent literature on unit-specific forecasts in panel data setups. To the best of our knowledge, this paper is the first investigation into unit-specific forecasts in panel data settings with binary outcomes and extreme covariates. In linear models, liu2020forecasting and liu2023density develop empirical Bayes and full Bayesian approaches for point and density forecasts, respectively. In nonlinear models, christensen2020robust propose efficient robust forecasts for discrete outcome models, and liu2023forecasting construct set and density forecasts for Tobit models with censored outcomes. Also see Baltagi2013 for the best linear unbiased predictor, giacomini2023robust for robust forecasts, and qu2023comparing for forecasting comparison. Most existing methods focus on continuous outcomes (except for christensen2020robust) and mid-sample properties, and thus are not suitable for our purpose. Also note that panel data are particularly valuable for analyses involving extreme events, as the limited tail information makes it even more beneficial to combine information across cross-sectional units.
Finally, our work is also related to the growing interest in studying extreme events and their economic consequences, partly motivated by the 2007--2008 financial crisis and the Covid-19 pandemic. Existing literature addresses extreme events in various ways: some exclude them (e.g., schorfheide2021real), some adapt models to accommodate both extreme and ordinary observations within a unified framework (e.g., carriero2022addressing, and lenza2022estimate), and some employ quantile regressions (e.g., tobias2016covar, plagborg2020growth, and adrian2022term). Our work differs in several ways. First, their approaches are specific for continuous outcomes, not for binary ones. Second, our approach focuses on cross-sectional and panel data rather than time-series data, allowing us to exploit information across units.\footnote{We also consider dynamic panels under conditional stationarity, and our framework could be extended to accommodate strongly mixing conditions in the time-series dynamics as in Nicolau2023, though this extension is beyond the scope of this paper.} Third, we incorporate a semiparametric approach that models extreme events separately from the middle of the data, acknowledging potentially distinct patterns in normal and extreme environments. While the value-at-risk and quantile regression papers also consider the second and third points, our approach is specifically designed for binary outcomes, with a different model structure and estimation procedure.
The remainder of this paper is organized as follows. Sections (ref) and (ref) specify our methodology for cross-sectional and panel data, respectively, and derive their asymptotic properties. Section (ref) extends our estimator to various contexts, such as a dynamic panel data model. Section (ref) conducts Monte Carlo experiments to examine the finite-sample properties of our estimators. Section (ref) employs our panel data estimator to analyze how local housing price declines affected the riskiness of small banks during the 2007--2008 financial crisis. Finally, Section (ref) concludes. Appendix (ref) provides the proofs for all propositions and theorems, and Appendix (ref) contains additional tables and figures.
In this section, we continue the discussion from the introduction for cross-sectional data. Section (ref) focuses on intuition and illustrates our main idea without additional covariates. It also compares our method with the classic threshold-crossing model. Section (ref) details the estimator and establishes its asymptotic properties in a general setup with additional covariates.
Recall that in the introduction, we presented a simple cross-sectional model with binary outcome $Y$ and continuous covariate $X$. Without loss of generality, let $X\in\mathbb{R}^{+}$, and $x\rightarrow\infty $ as $N\rightarrow \infty$.\footnote{Directly incorporating multidimensional extreme covariates poses theoretical challenges. However, a practical solution is to replace $X$ with $v(X;\gamma)$, where $v(\cdot;\gamma)$ is a known function with unknown finite-dimensional parameters $\gamma$, and $\gamma$ can be consistently estimated with sufficiently fast convergence rate. For simplicity, we will treat $X$ as a scalar in the rest of the paper, but most of our discussions can be extended to the case of $v(X;\gamma)$.} We focus on three potential objects of interest: the conditional probability $\pi(x)$, the partial effect $\partial \pi(x) /\partial x$, and the extreme elasticity $\delta(x)$, as defined in (ref) and (ref).
Heavy-tailed distributions are well-documented in economic and financial data, e.g., gabaix2009,Gabaix2016JEP, and the literature has suggested various methods to assess their presence. First, a log-log plot provides a visual assessment of tail behavior by plotting the threshold $x$ against the probability of exceeding $x$, both on a log scale, highlighting the relative prevalence of large values in the data. If tail observations align around a downward-sloping line, this suggests a heavy tail pattern with the slope being approximately $-\alpha$, where $\alpha$ is the Pareto exponent or tail index. In contrast, a vertical alignment indicates a relatively thin tail. See Figures (ref) and (ref) for examples with simulated and empirical data. Second, a more rigorous evaluation involves estimating the tail index and conducting statistical tests: see for example, clauset2009power. Finally, for a data-driven evaluation, techniques like cross-validation and pseudo out-of-sample forecasting can be employed as well, as demonstrated in our Monte Carlo Experiment 2 and empirical example.
Existing methods encounter difficulties in handling extreme values. Parametric approaches, such as Logit and probit, may incur substantial misspecification bias in the tail, while nonparametric methods, such as kernel, sieve, and spline, may be highly inefficient due to a limited number of tail observations. To overcome these challenges, we propose a semiparametric approach based on Bayes' theorem and RV functions as follows.
First, we observe extreme values in the covariate $X$ and seek to understand the behavior of $Y$ under extreme $X$, so we essentially aim to analyze their comovement in the tail. To facilitate this tail analysis, we use Bayes' theorem to “reverse” the conditioning set in the conditional probability (ref), resulting in equation (ref). Note that the Bayes' theorem representation is simply an alternative characterization of the data, so we are agnostic about the causal direction between $X$ and $Y$. A similar Bayes' theorem representation has been employed in the discriminant analysis (Chapter 9.2.8 in amemiya1985advanced), where normal distributions of $X|Y=y$ leads to a quadratic Logit form of $\mathbb{P}(Y=1|X=x)$. klein1993efficient also use a similar Bayes' theorem transformation in their semiparametric single-index estimator.
Next, as $x\rightarrow \infty $, the conditional pdfs $f_{X|Y}\left( x|y\right)$ for $y\in\{0,1\}$ are dominated by their tail behaviors. If $X|Y=y$ exhibits a heavy tail, its distribution can be well approximated by a Pareto distribution smith1987MLE. Therefore, we adopt a more general concept of RV functions that leads to a Pareto approximation in the tail.\footnote{Our method can be adapted to the generalized Pareto distribution, which also encompasses thin-tailed and bounded-support distributions. However, these distributions are often easily distinguishable from heavy-tailed distributions given empirical data. Moreover, the simplicity of our Pareto-based estimator offers a practical advantage over the more involved generalized Pareto estimator.}
We first introduce the following definitions. Let $X\in \mathbb{R}^{+}$ be a generic random variable with a heavy-tailed distribution. We say the upper tail probability (or survival function) $1-F_X(x)$ is Regularly Varying (RV) at infinity with index $-\alpha$ for some $\alpha>0$, if for all $x>0$, as $\eta \rightarrow \infty $,
which is denoted as $1-F_X\in RV_{-\alpha }$. Equivalently, by Karamata's characterization theorem, we can write
where $\mathcal{L} \left( \cdot \right) $ is a slowly varying function such that for all $x>0$, $ \frac{\mathcal{L} \left( \eta x\right) }{\mathcal{L} \left( \eta \right) }\rightarrow 1, $ as $\eta \rightarrow \infty $.
Let us highlight a few key points. First, the RV condition (ref) implies that as ${\underline{x}}\rightarrow\infty$, for $x\ge{\underline{x}}$,
Thus, the tail distribution of $X$ is asymptotically proportional to a Pareto distribution as the first-order approximation. Second, the parameter $\alpha $ is also referred to as the Pareto exponent, which characterizes the tail thickness of $1-F_X\left(x \right)$. In particular, a smaller $\alpha $ indicates a heavier tail, meaning that the upper tail probability decays to zero more slowly; as $\alpha\rightarrow\infty$, the distribution becomes thin-tailed, such as logistic and normal distributions. Third, the RV condition is relatively mild and satisfied by many commonly used heavy-tailed distributions, including the Student-$t$, $F$, and Cauchy distributions. For example, in a Student-$t$ distribution, $\alpha $ is equal to the degrees of freedom. A list of distributions and their corresponding values of $\alpha $ can be found in Gabaix2016JEP. Fourth, the moments of $X$, $\mathbb{E}\left[ \left| X\right| ^{r}\right]$, are finite up to order $r<\alpha$. Note that we can handle cases with very heavy tails where $\alpha\in(0,1)$ and no moments exist. Finally, following from Proposition B.1.9(11) in de2006extreme, if the pdf $f_{X}(x)$ is non-increasing for sufficiently large $x$, we have the following approximation that as $x\rightarrow\infty$,
Back to our binary outcome model, we assume that the distribution of $X$ conditional on $Y=y$ satisfies the RV condition (ref), that is, for $y\in\{0,1\}$, \( 1-F_{X|Y}\left(\cdot|y\right)\in RV_{-\alpha^{(y)}}. \) The Pareto exponent $\alpha^{(y)}$ is indexed by $y$, allowing for different tail thickness for $y\in\{0,1\}$. If $Y$ were continuous, estimating $\alpha^{(y)}$ as a function of $y$ would be challenging without parametric assumptions. The binary outcome simplifies this, as $\alpha^{(y)}$ takes only two values, $\alpha ^{(0)}$ and $\alpha ^{(1)}$, which can be easily estimated using methods such as the classic hill1975 estimator in this simple model. Further estimation details are provided in Section (ref). Also note that our approach is semiparametric, as we impose no parametric assumptions on the relationship between $X$ and $Y$ outside the tail region.
After estimating $\alpha^{(y)}$, we can proceed with the conditional probability $\pi(x)$ as well as other objects of interest. Based on Bayes' theorem (ref) and the pdf approximation (ref), it follows that
where $f(x)\sim g(x)$ denotes $\lim_{x\rightarrow\infty} f(x)/g(x)=1$ for generic functions $f(x)$ and $g(x)$. By the RV condition (ref), given some large ${\underline{x}}^{(y)}$, the term in the denominator becomes
Then, $\pi(x)$ can be consistently estimated by plugging in $\hat\alpha^{(y)}$ and the sample analog of $\frac{\mathbb{P}\left( Y=0,X\ge{\underline{x}}^{(0)}\right) }{\mathbb{P}\left( Y=1,X\ge{\underline{x}}^{(1)}\right) }\approx\frac {N^{(0)}}{N^{(1)}}$, where $N^{(y)}=\sum_{i=1}^N \mathbf{1}\left\{ X_i^{(y)}\ge {\underline{x}}^{(y)}, Y_i=y\right\} $ is the number of tail observations in the subsample with $Y_i=y$.\footnote{The consistency result was established in a previous version of this paper and is available upon request.}
After estimating $\pi(x)$, the partial effect $\partial \pi(x) /\partial x$ can be obtained as a by-product through either analytical or numerical differentiation. It may also be of interest to consider averages of the partial effects over certain covariate values, which relate to commonly used measures, such as the average partial effects (APE) and the average marginal effects (AME) in panel data models: see AbrevayaHsu2021 for a survey of various partial effects in panels, and DaveziesDHaultfoeuilleLaage2021 for recent partial identification results of average causal effects in short-$T$ panel Logit models with fixed effects. Proposition (ref) in the Appendix establishes the existence of the tail average of partial effects.
For extreme elasticity (ref), given the Pareto approximation in the tail, the extreme elasticity is solely determined by the difference in tail indices, $\alpha^*=\alpha ^{(1)}-\alpha ^{(0)}$.
We now compare our approach with the classic threshold-crossing model, which has been extensively studied in the literature:
where $\varepsilon$ is the error term. Note that our method does not require a threshold-crossing structure. However, the following proposition reveals a direct relationship between our key parameters $\alpha ^{(y)}$ and the Pareto exponents of $X$ and $\varepsilon$ in the threshold-crossing model.
There are three points worth noting. First, the threshold-crossing structure and the unconditional RV tails provide a sufficient condition for our conditional RV tails. The comovement between $X$ and $\varepsilon$ in the tail is captured by the comovement between $X$ and $Y$ in the tail.
Second, the difference in the conditional Pareto exponents, $\alpha ^{(0)}-\alpha ^{(1)}$, equals the Pareto exponent of the unobserved error term. As the tail of the error term becomes thinner with larger $\alpha_{\varepsilon}$, this difference increases. The extreme case of a very thin-tailed error with $\alpha_{\varepsilon}\rightarrow\infty$ (e.g., close to the Logit) corresponds to $\alpha ^{(0)}\rightarrow\infty$, $\alpha ^{(1)}=\alpha_X$, and $\pi(x)\rightarrow 1$: see also the first case in Remark (ref).
Finally, directly estimating the threshold-crossing model with RV errors $\varepsilon$ could be challenging due to two related issues: first, $\varepsilon$ is not observable, making it difficult to determine which observations are in the tail; second, using the full data could lead to significant bias, as the middle part may deviate substantially from the Pareto distribution. In contrast, our method circumvents this issue by using cutoffs of observed $X$, and offers an asymptotic equivalent estimator that is relatively easy to implement and scalable to more complicated setups, such as those involving additional covariates and panel data.
In this subsection, we first extend the simple cross-sectional model in Section (ref) to incorporate additional covariates, which is more empirically relevant. We then derive the asymptotic properties of our proposed estimator.
Let $Z$ be a $d_Z$-dimensional vector of non-extreme covariates in addition to $X$. $Z$ can be either discrete or continuous. For example, consider the context of sovereign debts and country defaults, where institutional characteristics could be included as covariates given their potential impact on default risk. Now we rewrite our conditional probability and the Bayes' theorem representation by adding $Z$ to all conditioning sets:
Note that here we hold $z$ constant and let $x$ tend to infinity. Accordingly, we assume that $ 1-F_{X|Y,Z}\left(\cdot|y,z\right)\in RV_{-\alpha^{(y)}\left( z\right) }, $ where the Pareto exponent $\alpha^{(y)}\left( z\right) $ is a function of the additional covariates $z$.
While $\alpha^{(y)}(\cdot)$ could depend on $z$ in a complicated nonlinear way, it is difficult to nonparametrically estimate it due to the data limitation that only large values of $X$ are used for tail analysis. To sidestep this issue, we assume a pseudo-linear structure:
where $\theta^{(y)}$ denotes the pseudo-parameters for $y\in\{0,1\}$, and $\alpha^{(y)}\left( z\right)=\alpha^{(y)}\left( z;\theta^{(y)}_0\right)$ represents the tail index function evaluated at the pseudo-true parameter values. This simplification is also in line with the literature, such as WangTsai2009.\footnote{We adopt a pseudo-linear approximation to better align with the panel data case in Section (ref), rather than an exponential approximation in WangTsai2009.} Additionally, we require that $Z_i'\theta^{(y)}>0$ almost surely, with sufficient conditions discussed after Assumption (ref).
From the RV condition, for $y\in\{0,1\}$, given some large threshold ${\underline{x}}^{(y)}_N$, if the conditional pdf is non-increasing for sufficiently large $x$, it is asymptotically equivalent to a Pareto distribution, \[f_{X|Y,Z}\left( x\left|y,z,x\ge{\underline{x}}^{(y)}_N\right.\right)\sim\alpha^{(y)}\left( z\right)\left(\frac x {{\underline{x}}^{(y)}_N}\right)^{-\alpha^{(y)}\left( z\right)-1}.\]
Let $\Xi^{(y)}_{N,i}=\left\{ X_i^{(y)}\ge{\underline{x}}^{(y)}_N, Y_i=y\right\}$, then $N^{(y)}=\sum_{i=1}^N \mathbf{1}\left\{\Xi^{(y)}_{N,i}\right\} $ is the total number of tail observations in the subsample with $Y_i=y$. The asymptotic Pareto distribution above leads to the following pseudo-MLE
where $\Theta^{(y)}\subset\mathbb R^{d_Z}$ is a convex cone as discussed after Assumption (ref), and the pseudo-true parameter value $\theta^{(y)}_0\in {\text{int}}\left(\Theta^{(y)}\right)$, the interior of $\Theta^{(y)}$. Since the objective function is concave and the domain is convex, this MLE is a convex programming problem that can be easily solved. In the simple model without $Z$, the FOC of (ref) results in \( \hat\alpha^{(y)}=\left[\frac{1}{N^{(y)}}\sum_{i=1}^N\left( \log X_i -\log {\underline{x}}^{(y)}_N \right)\mathbf{1}\left\{\Xi^{(y)}_{N,i}\right\}\right] ^{-1}, \) which coincides with the classic hill1975 estimator, as the indicator function $\mathbf{1}\left\{\Xi^{(y)}_{N,i}\right\}$ effectively selects the largest order statistics.
To establish the asymptotic properties, we impose the following assumptions.
Condition (b) ensures non-degenerate probabilities, thus guaranteeing a non-trivial Bayes' theorem representation.
First, the expression in condition (a) implies the RV condition and resembles equation (2.2) in WangTsai2009. In the absence of additional covariates $Z$, it simplifies to the well-studied second-order condition: see for example, hall1982 and Chapter 2 in de2006extreme. The first-order term $C^{(y)}\left( z\right) x^{-\alpha^{(y)}\left( z\right) }$ is proportional to a Pareto distribution, while the second-order term $D^{(y)}\left( z\right) x^{-\beta ^{(y)}\left( z\right) }$ captures deviations from the ideal Pareto form and is essential for characterizing the asymptotic distribution. Second, to guarantee the existence of $\theta^{(y)}$ such that $\alpha^{(y)}(Z_i)=Z_i'\theta^{(y)}>0$ almost surely (i.e., the existence of the convex cone), a sufficient condition is that each element of $Z_i$ has a domain strictly above or below zero. This can be achieved through a monotonic transformation of $Z_i$, such as the exponential or probability integral transforms. For example, if $Z_i\in\mathbb R^{d_Z}_{++}$ (the set of vectors with strictly positive elements) almost surely, then we can set $\Theta^{(y)}=\mathbb R^{d_Z}_{+}\setminus{\{0\}}$ (the set of vectors with non-negative elements, excluding the zero vector). Third, the bounds on $D^{(y)}\left(z\right)$ and $\beta^{(y)}\left(z\right)$, as well as the supremum condition (b) on the remainder term, allow us to ignore higher order terms in the asymptotic analysis. Finally, the relatively mild condition (c) of a non-increasing conditional probability density function ensures that the tail behavior of the pdf can be inferred from the corresponding cdf, according to Proposition B.1.9(11) in de2006extreme.
To derive the asymptotic properties, note that $N^{(y)}$ is a random variable, so we need to carefully account for its randomness when applying the LLN and CLT. We further define
a non-random sequence representing the asymptotic proportion of tail observations for $Y_i=y$, and utilize it to characterize the asymptotic behavior of our estimator.
Conditions (a) and (b) jointly impose upper and lower bounds on the rate at which the tuning parameter ${\underline{x}}^{(y)}_N$ tends to infinity. It implies that $N^{(y)}$ goes to infinity at a slower rate than $N$. Also note that in condition (b), we select a larger ${\underline{x}}^{(y)}_N$ to eliminate the asymptotic bias, albeit at the expense of a slower convergence rate. This is close in spirit to choosing an undersmoothing bandwidth in kernel regressions. Condition (c) is a mild regularity condition that ensures the invertibility of the Hessian matrix.
The following theorem establishes the asymptotic result, building on WangTsai2009.
The convergence rate is slower than the standard parametric $\sqrt N$-rate, since the tail estimator only accounts for observations in the tail and the number of tail observations $N^{(y)}=N\xi^{(y)}_N\left(1+o_p(1)\right)$ increases at a slower rate than the total number of observations $N$. The convergence rate depends on the magnitude of the second-order term: larger $\beta^{(y)}(z)$ allows for a smaller cutoff ${\underline{x}}^{(y)}_N$ according to Assumption (ref)(b), and thus a larger $\xi^{(y)}_N$ based on equation (ref), leading to a faster convergence rate. Especially, when $\beta^{(y)}(z)\rightarrow\infty$, the distribution approaches the ideal Pareto case, and the convergence rate is close to the $\sqrt N$-rate.
For implementation, we adopt a common practice to select the threshold ${\underline{x}}^{(y)}_N$: taking empirical quantiles (e.g., 90th and 95th percentiles) as potential thresholds, and checking them through graphic diagnostics like log-log plots. While there are various threshold estimation methods in the literature,\footnote{For example, guillou2001diagnostic choose the threshold by minimizing the asymptotic MSE, clauset2009power select the threshold by maximizing the marginal likelihood or minimizing the distance between the power-law distribution and empirical data, and WangTsai2009 further consider a discrepancy measure for a tail index regression model with covariates.} our simulations and empirical analysis suggest that parameter estimates are relatively robust within a range of threshold values. Intuitively, in a log-log plot, such as Figure (ref), the threshold ${\underline{x}}^{(y)}_N$ acts as the starting point for slope estimation, and small variations in the threshold within the downward-sloping region would only minimally affect the estimated slope.
Given $\hat\theta^{(y)}$, we can obtain the conditional probability $\pi \left( x,z\right)$ similar to (ref) and (ref) in the simple model. It is more practical to estimate $\mathbb{P}\left( Y=y,X\ge{\underline{x}}^{(y)}_N|Z=z\right)$ parametrically, due to limited tail data. Analogous to Proposition (ref), if we further assume that $f_{X|Y,Z}'\left( x|y,z\right) $ is non-increasing for sufficiently large $x$, then as $x\rightarrow \infty $, the extreme elasticity is now
A consistent estimator is obtained by substituting $\hat\theta^{(y)}$ in place of $\theta^{(y)}$.
Our method is particularly useful in panel data analysis, specifically in addressing unobserved individual heterogeneity. For instance, in the context of country defaults, different countries could have various cultural and historical backgrounds that might not be fully captured by observed data. In the analysis of extreme events, this unobserved heterogeneity could manifest as unobserved unit-specific tail thickness and RV functions.\footnote{For panel data with observed heterogeneity only, a pooled estimator with observed covariates can be employed, similar to the cross-sectional case in Section (ref).}
This section focuses on panel data with large $N$ and small $T$. To build intuition, we begin in Section (ref) by considering a simple case without additional covariates $Z_i$. We then incorporate these covariates in Section (ref) and derive the asymptotic for the common parameters and extreme elasticity. Later on, we will extend our discussions to models with large $T$ in Section (ref), time-varying additional covariates in Section (ref), and dynamic panel data in Section (ref).
Suppose we observe a panel dataset $\{Y_{it},X_{it}\}$ for $i=1,...,N$ and $t=1,...,T$. For illustrative purposes, let $T=2$ (though our method is applicable to any $T \ge 2$), and $X_{it}$ be a scalar extreme covariate.
As $X_{it}$ approaches infinity, unobserved individual heterogeneity could reflect in unit-specific tail thickness and RV functions. Assume the unit-specific tail indices take the following additive form
where $\alpha^{(y)}$ represents the common component depending on $y\in\{0,1\}$, and $ \lambda_i$ denotes the unit-specific component. For example, the tail indices of stock market returns exhibit heterogeneity across regions Jondeau2003, as do Covid-19 cases and deaths across countries Einmahl2023. The additive form helps cancel out the unobserved unit-specific tail thickness, as will soon be demonstrated.
Let $\mathbb{P}_i$ be the probability measure given unit-specific quantities. Specifically, $\mathbb{P}_i(\cdot)=\mathbb{P}\left(\cdot\;;\;\lambda_i,\left\{\mathcal{L}^{(y)}_i\right\}\right)$, where $\lambda_i$ is the unit-specific tail thickness as defined above, and $\mathcal{L}^{(y)}_i(\cdot)$ is the unit-specific slowly varying function as specified below. Any functions subscripted with $i$ are implicitly conditioned on these unit-specific quantities. Note that we are working within a fixed effects framework where $\left\{\lambda_i,\left\{\mathcal{L}^{(y)}_i\right\}\right\}$ are considered as fixed for each unit $i$ and can be arbitrarily correlated with the covariates of $i$, making the panel data analysis richer and more challenging than the cross-sectional case.
Based on a simplified version of Assumption (ref) in Section (ref) on the model setup, we assume that $\{Y_{it},X_{it}\}$ are independent across units and stationary across time. The stationarity condition enables us to eliminate unit-specific tail effects, and the conditional independence condition ensures that for $y_1,y_2\in \{0,1\}$,
We rewrite the conditional probability of unit $i$ using Bayes' theorem similar to the cross-sectional case.
Here the subscript $i$ in $\pi _i(x)$ indicates potential heterogeneity across $i$, and the absence of a $t$ index reflects the stationarity condition.
To proceed, we again apply the Pareto approximation at the tail, assuming that $1-F_{i,X_{it}|Y_{it}}(\cdot|y)\in RV_{-\tilde\alpha_i^{(y)}},$ or equivalently, $1-F_{i,X_{it}|Y_{it}}(x|y) =x^{-\tilde\alpha_i^{(y)}} \mathcal{L}_i^{(y)}(x),$ for some slowly varying function $\mathcal{L}_i^{(y)}(x) $. Given the stationarity condition, $\mathcal{L}_i^{(y)}(x) $ does not change over time. From the RV condition and the pdf approximation in (ref), we have that
Let us look into each term one by one. First, note that $\frac{x^{-\tilde\alpha_i^{(0)}}}{x^{-\tilde\alpha_i^{(1)}}}=x^{\alpha ^{(1)}-\alpha ^{(0)}}$, since $\lambda_i$ enters additively into the tail indices (ref) and is thus canceled out in the ratio. Second, we have that $\frac{\mathcal{L} _i^{(0)}(x) }{\mathcal{L} _i^{(1)}(x) }\rightarrow \mathcal L^*_i$, which does not depend on $x$ asymptotically due to the slowly varying nature of $\mathcal{L}_i^{(y)}(x)$. Third, the term $\frac{\mathbb{P}_i\left(Y_{it}=0\right) }{\mathbb{P}_i\left( Y_{it}=1\right) }\frac{\tilde\alpha_i^{(0)}}{\tilde\alpha_i^{(1)}}$ could depend on $\left\{\lambda_i,\left\{\mathcal{L}^{(y)}_i\right\}\right\}$ in a complicated, possibly nonlinear manner, as they enter into the conditioning sets. Combining all three terms, let $A_i=\frac{\mathbb{P}_i\left(Y_{it}=0\right) }{\mathbb{P}_i\left( Y_{it}=1\right) }\frac{\tilde\alpha_i^{(0)}}{\tilde\alpha_i^{(1)}}\mathcal L^*_i,$ and the conditional probability admits the following asymptotic approximation
Without further assumptions, $A_i$ cannot be consistently estimated in small-$T$ panels due to the incidental parameter problem. However, we can eliminate $A_i$ by conditioning on the event $Y_{i1} + Y_{i2} = 1$. Following from the conditional independence (ref), we obtain that as $x_{1},x_2\rightarrow \infty $,
The key idea here parallels the panel data Logit model, which we will further compare in Remark (ref) below.
Let $\alpha^*=\alpha^{(1)} - \alpha^{(0)}$, then the extreme elasticity is given by \[\delta_i(x)=\frac{\partial \left(\pi_i(x)\left(1-\pi_i(x)\right)\right) }{\partial x}\frac{x}{\pi_i(x)\left(1-\pi_i(x)\right)}\rightarrow -\left|\alpha^*\right|,\] which implies that, under the RV condition together with the additive form (ref), all units share the same extreme elasticity, despite having unit-specific tail thickness and RV functions. We can estimate $\alpha^*$ via the conditional MLE on tail observations
where event $\Xi_{N,i}=\left\{Y_{i1}+Y_{i2}=1, X_{it}\ge {\underline{x}} _{N}\text{, for }t=1,2\right\},$ and the tail threshold ${\underline{x}}_{N}\rightarrow\infty$ as $N\rightarrow\infty$. While in principle the tail threshold could be specific to both the period $t$ and the value of $y$, we set it to a common ${\underline{x}}_N$ for simplicity in both theory and implementation.
In this subsection, we derive the asymptotic normality of our estimator in the panel data setup. We now extend the baseline model in Section (ref) and include additional covariates $Z_i$ to capture potential observed heterogeneity. The unit-specific tail thickness becomes
and $\tilde\alpha_i^{(y)}\left( z\right)=\tilde\alpha_i^{(y)}\left( z;\theta^{(y)}_0\right)$ denotes the tail index functions evaluated at the pseudo-true parameter values. The corresponding RV condition is $1-F_{i,X_{it}|Y_{it},Z_i}(\cdot|y,z)\in RV_{-\tilde\alpha_i^{(y)}\left(z\right) }.$ We focus on the time-invariant $Z_i$ in this subsection, and the case with time-varying $Z_{it}$ is discussed in Section (ref).
Let $\mathcal{C}_i$ denote the unit-specific collection of functions $\left\{\beta^{(y)}_i(\cdot),\,C^{(y)}_i(\cdot),\,D^{(y)}_i(\cdot)\right\}$ in the unit-specific distribution $F_{i,X_{it}|Y_{it},Z_i}$: see Assumption (ref) for more details. As $\left\{\lambda_i,\mathcal{C}_i\right\}$ are fixed for each $i$, we denote $\mathbb{P}_i$ as the probability measure given $\left\{\lambda_i,\mathcal{C}_i\right\}$, that is, $\mathbb{P}_i(\cdot) = \mathbb{P}(\cdot\;;\;\lambda_i,\mathcal{C}_i)$. Similarly, we denote $\mathbb{E}_i[\cdot] = \mathbb{E}[\cdot\;;\;\lambda_i,\mathcal{C}_i]$.
First, we adopt the following assumptions on the model setup.
In condition (a), the covariates and outcomes are i.n.i.d.\ across units due to the unobserved unit-specific heterogeneity as fixed effects. Condition (b) ensures stationarity that helps eliminate the unit-specific tail effects. Conditions (c) and (d) imply conditional independence, thus for $y_1,y_2\in \{0,1\}$,
Condition (e) guarantees non-degenerate probabilities for all units.
Second, we assume the following tail conditions.
Please refer to the discussion after Assumption (ref) for a detailed explanation. Compared to Assumption (ref), we now have all these functions as unit-specific, indexed by $i$, and condition (b) further bounds the derivative of the remainder term for the i.n.i.d.\ case. Especially, we can accommodate a unit-specific scale parameter, given by $\left[C^{(y)}_i\left( Z_i\right)\right]^{1/\tilde\alpha^{(y)}(Z_i)}$, in the Pareto tail approximation. These unit-specific functions will be absorbed into $A_i$ in (ref) below, and then differenced out for small $T$ or estimated for large $T$.
Then, $\mathcal{L}^{(y)}_i(x,z)=C^{(y)}_i\left( z\right)\left( 1+D_i^{(y)}\left( z\right) x^{-\beta^{(y)}_i\left( z\right) }+r^{(y)}_i\left(x,z\right) \right)$ is the corresponding slowly varying function, and $\frac{\mathcal{L}^{(0)}_i(x,z)}{\mathcal{L}^{(1)}_i(x,z)}\rightarrow\mathcal{L}^*_i(z)$. Define $\theta^*= \theta^{(1)}-\theta^{(0)}$, and $A_i=\frac{\mathbb{P}_i\left(Y_{it}=0|Z_i\right) }{\mathbb{P}_i\left( Y_{it}=1|Z_i\right) }\frac{\tilde\alpha_i^{(0)}}{\tilde\alpha_i^{(1)}}\mathcal L^*_i(Z_i)$. As $x\rightarrow\infty$,
Note that as $\mathbb{P}_i(\cdot) = \mathbb{P}(\cdot\;;\;\lambda_i,\mathcal{C}_i)$, $A_i$ can be viewed as the fixed effects that can potentially depend on the observed heterogeneity $Z_i$ and unobserved heterogeneity $\left\{\lambda_i,\mathcal{C}_i\right\}$ in an arbitrary way.
For small $T$, we again eliminate $A_i$ by conditioning on $Y_{i1}+Y_{i2}=1$ and construct the conditional MLE as follows:
where $\Theta^*\subset\mathbb R^{d_Z}$ is a convex cone and the pseudo-true parameter value $\theta^*_0\in {\text{int}}\left(\Theta^*\right)$. Again, this can be implemented via the conditional MLE for panel Logit models, applied to tail observations and using $Z_i\log X_{it}$ as regressors.
Finally, analogous to the cross-sectional case, the number of units contributing to the conditional likelihood, $N_{\Xi}=\sum_{i=1}^N \mathbf{1} \{\Xi_i\} $, is a random variable, and we define the following non-random sequence to capture the asymptotic proportion of these units:
Also define the Hessian terms $H_{N0,i}=\mathbb{E}_i\left[\left. \frac{\left( \frac{X_{i1}}{X_{i2}}\right) ^{Z_i'\theta^*_0}}{\left( 1+\left( \frac{X_{i1}}{X_{i2}}\right) ^{Z_i^{\prime }\theta^*_0}\right) ^2}\left( \log \frac{X_{i1}}{X_{i2}} \right) ^2Z_iZ_i'\right|\Xi_{N,i}\right]$ and \(H_{N0} = {\sum_{i=1}^N H_{N0,i}\xi_{N,i}}/{\sum_{i=1}^N\xi_{N,i}}.\) Now we assume the following conditions for estimation.
This assumption is in line with Assumption (ref) for the cross-sectional case and accommodates i.n.i.d.\ observations across $i$. Condition (a) ensures that the number of observations contributing to the likelihood increases with the cross-sectional sample size, while employing undersmoothing to eliminate asymptotic bias. Conditions (b) and (c) guarantee the validity of the LLN and Lindeberg-Feller CLT for i.n.i.d.\ observations.
Under these assumptions, we establish the asymptotic normality of our estimator for panel data models with i.n.i.d.\ random variables. As the number of tail observations $N_{\Xi}=N\xi_N\left(1+o_p(1)\right)$ increases at a slower rate than $N$, the convergence rate is slower than the parametric rate.
Moreover, the unit-specific extreme elasticity depends only on observed heterogeneity,
provided that $f_{i,X_{it}|Y_{it},Z_i}'\left( x|y,z\right)$ is non-increasing for sufficiently large $x$. In practice, we can consistently estimate this elasticity by substituting $\hat\theta^*$ for $\theta^*$.
Our previous panel data analysis has focused on small-$T$ panels. In such cases, unit-specific parameters cannot be consistently estimated, and hence we eliminate them by conditioning on the sum of $Y_{it}$ over time, which is a sufficient statistic for the individual parameters. Now, with a large $T$, we are able to estimate these unit-specific parameters.
Based on (ref) in Section (ref), let $\tilde A_i = -\log A_i$, and we have that, as $x\rightarrow\infty$,
Given this asymptotic equivalence, we can resort to panel Logit estimators for large $N$ and large $T$. For example, as shown in Example 2 in fernandez2018fixed, the MLE estimates of $\{\theta^*,\{\tilde A_i\} \}$ are jointly consistent and $\hat\theta^*-\theta^*_0$ is asymptotically distributed as $\mathcal{N}\left( \frac{B}{N},\frac{H^{-1}}{NT}\right),$ as $N,T\rightarrow \infty $. $B$ denotes the asymptotic bias due to the incidental parameter, and analytical or jackknife methods can be applied for bias correction. See also stammann2016estimating. Once $\{\theta^*,\{\tilde A_i\} \}$ are consistently estimated, we can proceed to estimate the APE and provide unit-specific forecasts of $\mathbb{P}_i\left(Y_{i,T+1}=1|X_{i,T+1}=x,Z_i\right)$ based on (ref).\footnote{Introducing time-varying thickness $\theta^{(y)}_t$ could be possible with additional parametric structure, such as combining our estimator with the autoregressive conditional density models in hansen1994autoregressive. However, this changes $A_i$ to time-varying $A_{it}=\frac{\mathbb{P}_i\left(Y_{it}=0|Z_{i}\right) }{\mathbb{P}_i\left( Y_{it}=1|Z_{i}\right) }\frac{\tilde\alpha_{it}^{(0)}(Z_{i})}{\tilde\alpha_{it}^{(1)}(Z_{i})}{\mathcal{L}_i^*(Z_{i})}$, which can no longer be differenced out even if $\left\{\lambda_i,\mathcal{C}_i\right\}$ does not change over time. To address this, we need to impose further parametric assumptions on $A_{it}$ or employ a local likelihood estimator as in Section (ref).}
In empirical studies, some additional covariates $Z$ may vary over time. For instance, in the context of country defaults, factors such as global business cycles could play a significant role. Then, the unit-specific tail indices in (ref) becomes $\tilde\alpha_i^{(y)}\left( z_{it}\right) =z_{it}'\theta^{(y)}+\lambda_i$, and $A_i$ in (ref) becomes $A_{it}=\frac{\mathbb{P}_i\left(Y_{it}=0|Z_{it}\right) }{\mathbb{P}_i\left( Y_{it}=1|Z_{it}\right) }\frac{\tilde\alpha_i^{(0)}(Z_{it})}{\tilde\alpha_i^{(1)}(Z_{it})}{\mathcal{L} _i^*(Z_{it})} $, so we cannot difference out $A_{it}$ as in (ref).
However, under regularity conditions, $\theta^*$ can be estimated using a local (conditional) likelihood estimator tibshirani1987local, honore2000panel. Intuitively, it introduces an additional kernel weight term that controls the distances across $Z_{it}$. Assuming $T=2$ for simplicity, the estimator for $\theta^*$ is given by \[ \hat\theta^*=\arg \max_{\theta^*\in\Theta^*}\sum_{i=1}^{N}k^{(d_Z)}_h\left({Z_{i1},Z_{i2}}\right)\left\{(1-Y_{i1})\log\frac{X_{i1} }{X_{i2}}\cdot Z_{i1}'\theta^*+\log\left(1+\left( \frac{ X_{i1}}{X_{i2}}\right) ^{Z_{i1}'\theta^*}\right)\right\}\mathbf{1}\{\Xi_{N,i}\} \] where $k^{(d_Z)}_h({z_{i1},z_{i2}})=\frac 1 {h^{d_Z}} \prod_{d=1}^{d_Z}k\left(\frac{z_{i1,d}-z_{i2,d}} h\right)$ is a multidimensional kernel with $k(\cdot)$ being the kernel function and $h$ being the bandwidth.
In this subsection, we extend our analysis to dynamic panel data, incorporating predetermined variables such as lagged outcomes to capture potential persistence.
For small $T$, w.l.o.g., we assume the first order Markov property: for each $i$, \[X_{it},Y_{it} \;|\;X_{i,1:t-1},Y_{i,1:t-1},Z_{i}=X_{it},Y_{it} \;|\;Y_{i,t-1},Z_i,\] given $\left\{\lambda_i,\mathcal{C}_i\right\}$, where $\mathcal{C}_i$ comprises unit-specific functions in $F_{i,X_{it}|Y_{it},Y_{i,t-1},Z_i}$, similar to the tail approximation in Assumption (ref). The first order Markov property implies the stationarity of the conditional joint distribution $Y_{it},Y_{i,t-1},X_{it}\;|\;Z_i$.
Recall that previously we applied Bayes' theorem by partitioning the data into two subsets based on $Y_{it}=0$ and 1. With dynamics, we now further partition the data according to transition dynamics. In the first order Markov case, there are four transition patterns: $0\rightarrow0$, $0\rightarrow1$, $1\rightarrow0$, and $1\rightarrow1$. Therefore, we partition the data into four subsets by $\left( Y_{it},Y_{i,t-1}\right) =(y,y_{-})$, and characterize Bayes' theorem as follows
Similar to (ref), we introduce unit-specific tail thickness as
Then, we have $1-F_{i,X_{it}|Y_{it},Y_{i,t-1},Z_i}\left(\cdot|y,y_{-},z\right) \in RV_{-\tilde\alpha_i^{(yy_{-})}(z)}.$
In the previous static panel analysis, our exercise was equivalent to normalizing $\theta^{(0)}=\mathbf0$ and estimating $\theta^*=\theta^{(1)}$. Here, we similary normalize $\theta^{(00)}=0$ to ensure identification. With four transition patterns, a minimum of five periods of data $\mathbf{Y}_i=(Y_{i1,}Y_{i2,},Y_{i3,},Y_{i4},Y_{i5})'$ are required. The following four events cover all transition patterns and help construct the conditional likelihood:
Given the stationarity of the conditional joint distribution $Y_{it},Y_{i,t-1},X_{it}\;|\;Z_i$ implied by the first order Markov assumption, we have that as $x_t\rightarrow\infty$ for $t=1,\cdots,5$,
Then, the conditional MLE for $\left(\theta^{(01)},\theta^{(10)},\theta^{(11)}\right) $ can be constructed similarly to the static panel case in Section (ref). Note that unit $i$ contributes to the conditional likelihood only if $X_{it}$ appears in the tail for at least five periods. This data requirement may be challenging, so the method would be more suitable for datasets with a larger $N$, and a fixed but slightly larger $T$.
For large $T$, the bias correction in fernandez2018fixed remains applicable to dynamic panel data models, so the estimator in Section (ref) remains valid.
We conduct two sets of Monte Carlo simulation experiments. Experiment 1 examines cross-sectional data and focuses on estimation performance, providing intuitions into when and how the proposed estimator outperforms the alternatives. Experiment 2 investigates panel data with large $N$ and large $T$,\footnote{We also conducted Monte Carlo simulations in panel data with large $N$ and small $T$ in a previous version of this paper. Results from these simulations are available upon request.} and focuses on pseudo out-of-sample forecasting performance, aligning more closely with the empirical example of bank loan charge-off rates.
We compare the proposed estimator with four alternatives: a Logit estimator using all observations (Logit, all $X$), a Logit estimator using only tail observations (Logit, tail $X$), a local linear estimator, and a local Logit estimator, where the first two are parametric estimators and the last two are nonparametric ones.
For cross-sectional data, let $\beta = (\beta_0,\,\beta_1)'$. First, the Logit estimator is defined as \[Y_i = \mathbf{1}(\beta_0+\beta_1X_i-\varepsilon_i\ge0),\quad\varepsilon_i\sim\text{ standard logistic.}\] Second, the local linear estimator is given by
where $k_h(\cdot)$ is a kernel function with bandwidth $h$. We employ a Gaussian kernel here, and choose the bandwidth based on Silverman's rule of thumb $h \approx 1.06 \hat{\sigma} N^{-1/5}$, where $\hat{\sigma}$ is the standard deviation of the $X_i$s. The results are robust with respect to a range of bandwidth choices. Finally, the local Logit estimator is a flexible nonparametric estimator that specifically accounts for binary outcomes $Y_i$. Let $p_i(\beta,x) = \frac{1}{1 + \exp\left[-(\beta_0 + \beta_1 (X_i-x))\right]}$. The local Logit maximizes the locally weighted log-likelihood function
For panel data, the Logit estimator is given by \[Y_{it} = \mathbf{1}(\beta_0+\beta_1X_{it}+C_i-\varepsilon_{it}\ge0),\quad \varepsilon_{it}\sim \text{standard logistic},\] where $C_i$ captures unobserved individual heterogeneity. In panels with large $N$ and large $T$, $\{\beta_0,\beta_1,\{C_i\}\}$ can again be jointly estimated using a fixed effects estimator with bias corrections fernandez2018fixed, stammann2016estimating. For the local linear and local Logit estimators, the setup is more flexible \[Y_{it}=g\left(X_{it},C_i,\varepsilon_{it}\right),\] where the function $g$ and the distributions of $C_i$ and $\varepsilon_{it}$ could be unknown. We incorporate a correlated random effects structure, which allows for the unobserved individual heterogeneity to be correlated with the covariates $X_{it}$ (and $Z_i)$ and thus may help enhance the performance of these alternatives. More specifically, suppose there is a sufficient statistic $V_i$ which could be multidimensional, such that $C_i|X_{i,1:T}=C_i|V_i$. One commonly used example of $V_i$ is the time sum of $X_{it}$, i.e., $V_i = \sum_t X_{it}.$ Then, $V_i$ can be included in the conditioning set for the local linear and local Logit estimators, so we essentially run a nonparametric regression to estimate the conditional mean $\mathbb{E}[Y_i|X_i,V_i]$: see details in liu2021identification for general nonlinear panel data models. Without loss of generality, let $V_i$ be a scalar for notation simplicity. Now $\beta = (\beta_0,\,\beta_1,\,\beta_2)'$, $p_{it}(\beta,x) = \frac{1}{1 + \exp\left[-(\beta_0 + \beta_1 (X_{it}-x)+\beta_2(V_i-v))\right]}$, and
Furthermore, for all these estimators, we can also incorporate additional covariates $Z_i$, and the formulas are similar to those above.
Experiment 1 is based on cross-sectional data, where we focus on comparing estimation performance across estimators.
The Monte Carlo design is summarized in Table (ref). The data are generated from a threshold-crossing model with both $X_i$ and $\varepsilon_i$ following Student-$t$ distributions.\footnote{We use the difference in medians $\left(\text{med}_X-\text{med}_{\varepsilon}\right)$ as the threshold to keep the samples more balanced between $Y_i=0$ and 1.} Here we consider a range of $\alpha$ values from 0.5 to 2. As elaborated in Proposition (ref), in the tail, $\alpha^{(0)}=\alpha_X+\alpha_{\varepsilon}$ ranges from 1 to 4, $\alpha^{(1)}=\alpha_X$ varies from 0.5 to 2, and the extreme elasticity $-\left|\alpha^{(1)}-\alpha^{(0)}\right|=-\alpha_{\varepsilon}$ spans from $-2$ to $-0.5$.\footnote{In many economic datasets, the tail index $\alpha$ typically falls between 1 and 2. For example, in a review paper, gabaix2009 mentions that “it seems that the tail exponent of wealth is rather stable, perhaps around 1.5,” referencing klass2006forbes. Also, clauset2009power remark that the tail index usually lies between 1 and 2 (note that the $\alpha$ in their notation corresponds to $\alpha-1$ in ours).} The sample size $N=10000$ is directly comparable with our empirical data sets on bank loan charge-off rates, which comprises 8538 banks in the baseline sample. For each experimental setup, we execute 1000 Monte Carlo simulations.
Our tail estimator for $\alpha^{(y)}$ is based on the “rank-1/2” estimator in gabaix2011rank, which can be viewed as a refinement of the classic Hill estimator and often performs well in finite samples. Then, we estimate the conditional probability $\pi(x)$ using the sample analog of (ref) and (ref). We also compare with the alternative estimators described in Section (ref). For both the tail estimator and the “Logit, tail $X$” estimator, we set the cutoffs ${\underline{x}}_N^{(y)}$ at the 97.5th percentile of the distributions of $X$ for $Y=0$ and 1 separately.\footnote{The tail estimator is robust with respect to a range of cutoffs ${\underline{x}}_N^{(y)}$. As evident from the log-log plot in Figure (ref), the tail exhibits a pronounced downward pattern, with the slope remaining stable across small variations in ${\underline{x}}_N^{(y)}$.} In the main text, we plot the comparisons regarding estimated parameters, conditional probability, and extreme elasticity, for the specification with $\alpha_X = 1$ and $\alpha_{\varepsilon}=1$. For detailed results across all model specifications, please refer to Tables (ref)--(ref) in the Appendix. The main messages are similar across specifications with different tail index values.
Figure (ref) depicts a log-log plot from one of the 1000 Monte Carlo simulations, where the tail observations align around a downward sloping line, contrasting with the flatter non-tail region. This distinct pattern between the tail and non-tail region is common in many empirical datasets as well, and suggests that the proposed tail estimator would be able to effectively capture the heavy tail behavior. The estimated tail indices (red solid lines) closely match their true asymptotic values (black dashed lines). Figure (ref) further shows the distributions of the parameter estimates from all 1000 Monte Carlo simulations. The distributions are both bell-shaped and centered around the true asymptotic values indicated by the blue vertical lines. Notably, $\alpha^{(1)}$ is more precisely estimated than $\alpha^{(0)}$, as the former corresponds to a thicker tail.
Figure (ref) presents the estimated conditional probabilities for $X$ at the 90th, 95th, 97.5th, and 99th percentiles of its distribution. The blue vertical lines indicate the true conditional probability \(\mathbb{P}(Y=1\mid X=x)\), which increases with the evaluation point \(x\). The proposed tail estimator dominates all other alternatives, centered around the true probabilities with lower variance. The variance of the tail estimator decreases as the evaluation point $x$ increases, reflecting reduced estimation uncertainty as the conditional probability approaches one.
In contrast, the Logit estimators are largely biased due to the model misspecification that fails to account for the heavy tail. The direction of the bias depends on the relative positions of the estimation sample and the evaluation points. When the estimation sample includes all observations, the “Logit, all $X$” estimator shows an upward bias. This occurs because the estimation sample is overweighted by non-tail observations, and thus the tail evaluation points are too extreme given the misspecified thin-tailed logistic distributions, which leads to an overestimation of the probability of $Y=1$ in the tail evaluation points. Conversely, when the estimation sample includes tail observations only, the “Logit, tail $X$” estimator tends to exhibit a downward bias, except at the 99th percentile. The reason is that given the logistic model's misspecified assumption of a thinner tail, lower evaluation points appear relatively moderate compared to the tail observations for estimation, resulting in an underestimation of the probability of $Y=1$.
The local Logit estimator performs reasonably well for less extreme evaluation points, such as when $x$ is at the 90th percentile, but it exhibits large bias and variance at more extreme points. The local linear estimator is even more flexible, and thus yields even larger variance across all evaluation points. Due to their poor estimation performance and relatively long computation times, we omit these nonparametric estimators in Experiment 2 and the empirical example below, except for the baseline case in Table (ref).
Across all comparisons, we see that our tail estimator enjoys the best of both worlds as a semiparametric estimator that puts parametric assumptions only on the tail, the region of primary interest, while remaining agnostic in the non-tail region. Furthermore, as mentioned before, our approach offers significant improvement in finite samples for moderately large $x$ where $\pi(x)$ may not be very close to 1.
Figure (ref) plots the estimated extreme elasticities evaluated at the 97.5th percentile of the distribution of $X$, and the key messages are similar---the Logit estimators induce large bias, the nonparametric estimators exhibit high variance, whereas the proposed tail estimator flexibly and efficiently capture the tail behavior and yields the most accurate estimates.
Experiment 2 examines panel data with large $N$ and large $T$, focusing on pseudo out-of-sample forecasting performance, which is closely related to the empirical analysis of bank loan charge-off rates.
The Monte Carlo design is adapted from the cross-sectional case in Experiment 1, with specifics described in Table (ref). Now we incorporate unit-specific tail thickness $\lambda_i$ for $X_{it}$, where the underlying distribution of $\lambda_i$ is bimodal. For the tail estimator, we employ bias corrections as in fernandez2018fixed and stammann2016estimating. See Sections (ref) and (ref) for additional details on the tail and alternative estimators, respectively. The cutoffs for both the tail estimator and “Logit, tail $X$” are set at the 90th percentile of the $X$ distributions. Recall that the tail estimator is robust over a range of cutoffs, as discussed in footnote (ref).
The accuracy of density forecasts is evaluated using the log predictive score (LPS), as recommended by AmisanoGiacomini2007. The LPS is calculated as $ LPS = \frac{1}{ N_{f}^{\dagger}}\sum_{i\in \mathcal N_{f}^{\dagger}}\log\hat{p}(y_{i,T+1}|D) $, where $ y_{i,T+1} $ is the outcome at time $ T+1 $, and $ \hat{p}(y_{i,T+1}|D) $ is the predictive likelihood based on the estimated model and observed data $ D $. $\mathcal N_{f}^{\dagger}$ is the set of units forecasted, characterized by following criteria: (1) $X_{it}$ exceeds the 90th percentile in the estimation sample, (2) $Y_{it}$ switches values over time among these tail observations, and (3) $X_{i,T+1}$ also falls within the tail during the forecasting period. To assess the significance of the differences in the LPS, we combine the tests from AmisanoGiacomini2007 for density forecasts and qu2023comparing for panel data.
In Table (ref), the first row in each subpanel presents extreme elasticity estimates from the tail estimator, which are close to the true values, $-\alpha_{\varepsilon}$. Subsequent rows compare forecast accuracy among estimators, and we see that the tail estimator consistently outperforms both “Logit, tail $X$” and “Logit, all $X$” across all specifications. Although “Logit, tail $X$” is better than “Logit, all $X$,” it still performs worse than the tail estimator, especially with smaller $\alpha_X$ and $\alpha_{\varepsilon}$, where heavy tails are more pronounced. Therefore, it is important to distinguish the tail from the middle of the sample as well as account for heavy tail patterns. To further demonstrate this point, Figure (ref) in the Appendix shows the scatter plots of the LPS from all Monte Carlo repetitions in the setup with $\alpha_X = 1$ and $\alpha_{\varepsilon}=1$.
Charge-off rates serve as an indicator of bank losses. A bank could be riskier for a particular type of loan if its corresponding charge-off rates exceed a certain level. In our analysis, we focus on a panel of small banks with assets less than \$1 billion, similar to liu2023forecasting. Since the banks are small, it is reasonable to assume that they operate primarily in local markets. In this empirical example, we examine the impact of substantial local housing price declines on the riskiness of small banks, considering that this channel played a pivotal role during the 2007--2008 financial crisis.
In this empirical example, we focus on the setup in (ref) in Section (ref) for panel data models with large $N$ and large $T$. The binary outcome $Y_{it}$ is a risk dummy based on the loan charge-off rate for a specific loan type of bank $i$ in quarter $t$. $Y_{it} = 1$ if the charge-off rate is greater than a level $c$. We present results for $c = 0$ in the main text and relegate robustness checks for alternative $c$ values to the Appendix. The qualitative findings are consistent across different levels of $c$. This exercise aligns with policy analysis practices, where policymakers often compare current charge-off rates to historical averages to monitor bank risk: see, for example, the Federal Reserve Board's Financial Stability Report.\footnote{\href{https://www.federalreserve.gov/publications/financial-stability-report.htm}{https://www.federalreserve.gov/publications/financial-stability-report.htm}.} For instance, the historical average charge-off rate for Residential Real Estate (RRE) loans is around 0.1%, with the average plus two standard deviations around 0.2%. Robustness checks include these levels.
The extreme regressor $X_{it}$ captures decreases in local housing prices. To convert the extreme values to the right tail, we define $X_{it}$ as the deflation rate of the local housing price in the previous quarter.\footnote{We define the tail region as observations above the 90th percentile of $X_{it}$, as in Experiment 2. Given that the 90th percentiles are positive in all our samples, $\log X_{it}$ is well-defined in the tail.} The additional covariate $Z_i$ is given by the average quarterly change in the local unemployment rate, accounting for local economic conditions.
Our data are obtained from the following sources. Bank balance sheet data at a quarterly frequency, such as loan charge-off rates, are constructed based on the Call Reports from the Federal Reserve Bank of Chicago.\footnote{Please refer to Appendix D in liu2023forecasting for details on constructing loan charge-off rates from the raw data.} The local market is defined at the county level, and the local market for each bank is determined based on the annual Summary of Deposits from the Federal Deposit Insurance Corporation.\footnote{We calculate the deposits received by each bank from every county and link the bank to the county from which it received the largest amount of deposits. We also assess the robustness of our analysis by constructing weighted averages of the covariates, using deposit proportions from each county as weights. The results are very similar, which can be attributed to the concentration of deposits across counties and the similarity in covariate values among neighboring counties. } Housing price indices at a quarterly frequency (all transactions, not seasonally adjusted) are sourced from the Federal Housing Finance Agency, and the 3-digit zip code data are converted to the county level using the HUD USPS ZIP Code Crosswalk from the Department of Housing and Urban Development.\footnote{We use the 2010Q1 Crosswalk, the earliest available version that is close to our empirical time frame. We also experimented with the county-level Zillow Home Value Index, but it unfortunately has more missing data within our study period, though the results are qualitatively similar.} The county-level not seasonally adjusted unemployment rates are obtained from the Bureau of Labor Statistics website, and we aggregate the monthly data to a quarterly frequency by time averaging. Finally, we standardize the county-level housing price deflations and unemployment rate changes using the means and standard deviations calculated from their corresponding aggregate time series.
Our baseline sample focuses on the RRE charge-off rates. The estimation sample spans from 1999Q4 to 2009Q3, comprising $N=8538$ small banks across 40 quarters.\footnote{We exclude banks with: 1) average non-missing domestic total assets exceeding \$1 billion, or 2) missing target charge-off rates for all periods in the sample. These criteria result in the removal of 583 and 162 banks, respectively, for the baseline RRE sample.} There are $N_{e}^{\dagger}=2642$ banks in the tail for more than one period and contributing to the likelihood, that is, $X_{it}$ is above the 90th percentile of the estimation sample and $Y_{it}$ switches values across time for these tail observations. We perform a pseudo out-of-sample forecast for the period of 2009Q4. There are $N_{f}^{\dagger}=2098$ banks that satisfy the conditions for $N_{e}^{\dagger}$ and additionally have $X_{i,T+1}$ fall in the tail during the forecasting period.\footnote{To account for banks' endogenous exit choices, one could extend to a panel Tobit model as in liu2023forecasting, which is left for future research.}
For robustness check, we also consider (non-farm) non-residential commercial real estate (CRE) charge-off rates, various time dimensions in the estimation samples ranging from 32 to 48 quarters, and different forecasting periods $T+1 =$ 2009Q3 and 2009Q4. The sample statistics of these samples are reported in Table (ref) in the Appendix. The main results remain consistent across different samples.
Based on the sample skewness and kurtosis in Table (ref) for all samples, as well as the log-log plot in Figure (ref) for the baseline sample,\footnote{Fitted lines, such as those in Figure (ref), are not plotted here, because these lines could be unit-specific due to observed heterogeneity $Z_i$ and unobserved heterogeneity $\{\lambda_i,\mathcal{C}_i\}$.} we see that $X_{it}|Y_{it}=y$ indeed exhibits heavy right tails, so our tail estimator would be more appealing. Furthermore, it is worth noting that in finite samples, the probability of $Y_{it}=1$ for tail $x$, $\pi(x)$, is not necessarily very close to 1. This is evident from both the empirical data (see the log-log plot in Figure (ref) and the histogram in Figure (ref) in the Appendix) and the estimated model (see the predictive probability in Figure (ref)).
Table (ref) compares forecasting performance across estimators. The estimators are similar to those in the Monte Carlo Experiment 2: see Sections (ref) and (ref) for more details. The proposed tail estimator is the overall best. “Logit, tail $X$” ranks second, yet still significantly worse than the tail estimator, indicating the importance of carefully addressing the heavy tail behavior. “Logit, all $X$” ranks third, indicating the presence of large misspecification and possible distinct pattern in the tail compared to the middle range. Local Logit yields the least accurate forecasts, suggesting that while being flexible, the nonparametric approach could be too noisy given the limited data available in the tail.
Table (ref) in the Appendix provides further details on the parameter and APE estimates based on the tail estimator. First, the coefficient on $\log X_{it}$ is always significant, with values around 0.95--1.15 for the RRE samples and around 1.2--1.4 for the CRE ones. The negative of this coefficient can be roughly viewed as the homogeneous part of the extreme elasticity, and the estimated values suggest the potential presence of heavy tails. Second, the coefficient on $Z_i\log X_{it}$ is mostly positive, being larger and significant for the RRE samples, while smaller and insignificant for the CRE ones. This difference aligns with intuition: larger increases in unemployment could directly amplify the impact of housing price drops on residential loan risk, whereas changes in unemployment may not directly affect non-residential commercial loan risk. Third, the estimated APEs are around 0.15 for the RRE samples and 0.13 for the CRE ones, which could be interpreted as that in the tail, a 1% decrease in housing prices in the previous quarter corresponds to approximately a 0.15 (0.13) increase in the probability of high risk for RRE (CRE) loans. Finally, the unobserved heterogeneity exhibits a larger dispersion than the observed heterogeneity, as the sample variances of the estimated $\tilde A_i$ are around 1 while the sample variances of $\hat\theta^*_{Z\log X}Z_i\log X_{it}$ range from 0.1--0.3.
Tables (ref) and (ref) in the Appendix also show that our results are robust with respect to the level $c=0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 1$,\footnote{As the level $c$ increases, the coefficient on $Z_i\log X_{it}$ becomes less significant, particularly becoming insignificant when $c=1$. This is because a larger $c$ reduces the occurrence of $Y_{it}=1$, leading to fewer banks in the tail with $Y_{it}$ switching values across time. Then, the effective sample size for estimation $N_{e}^{\dagger}$ substantially decreases, resulting in noisier estimates.} to both CRE and RRE loans, to the estimation sample's time dimension $T=32,36,40,44,48$, and to the forecasting period being 2009Q3 and 2009Q4.
Figure (ref) plots the average predictive probability of high risk across counties, with darker shades indicating higher risk levels in those counties. Notably, Florida and the Great Lakes regions appear as areas with elevated risk. This heterogeneity in risk may arise from observed heterogeneity, such as variations in local housing prices $X_{i,T+1}$ and local unemployment rates $Z_i$, as well as unobserved heterogeneity $\{\lambda_i,\mathcal{C}_i\}$, which is captured by $\tilde A_i$, where larger values of $\tilde A_i$ are associated with higher risk: see (ref).
Accordingly, in Table (ref), we regress the estimated $\tilde A_i$ on bank characteristics. The “Initial” column uses bank characteristics in the initial period 1999Q4, which are exogenous to subsequent dynamics. The “Average” column uses time-averaged bank characteristics over the estimation period, incorporating more recent information. The findings across both columns are in general consistent. We see that larger banks (measured by log assets), those that specialize in RRE loans (measured by RRE loans to total loans ratio), in lending activities (measured by loan to assets ratio), with more diversified earnings (measured by the share of non-interest income to total income), and with higher operational efficiency (measured by the negative of overhead costs to assets ratio), tend to have a higher $\tilde A_i$ and thus higher risk, as these characteristics could help increase banks' capacity to assume riskier RRE loans. On the other hand, higher profitability (measured by return on assets) corresponds to a lower $\tilde A_i$, potentially due to reverse causality where reduced risk boosts profitability. The credit quality (measured by ALLL to total loans ratio) and the capital-asset ratio do not have significant effects in the regression using initial bank characteristics.
This paper proposes a novel semiparametric method based on Bayes' theorem and RV functions. It models heavy tail behavior through a Pareto approximation, while making no parametric assumptions on the relationship between the outcome and covariates outside of the tail region.
The proposed method is particularly useful in panel data models, accounting for unobserved unit-specific heterogeneity. We show that under regularity conditions, our objective function asymptotically aligns with a panel Logit regression on tail observations using $\log X_{it}$ as a regressor. Then, various established econometric techniques could be applicable, which could be convenient for empirical research. Specifically, in panels with large $N$ and small $T$, the unobserved unit-specific tail thickness and RV functions can be canceled out via the conditional MLE; in panels with large $N$ and large $T$, bias correction methods fernandez2018fixed,stammann2016estimating can be employed to estimate unit-specific parameters and predict unit-specific future outcomes. Furthermore, we also extend our method to dynamic panel data models that incorporate lagged outcomes.
The potential applications of our method could encompass both microeconomics and macroeconomics studies, particularly valuable in light of recent extreme events. For example, one may be interested in analyzing the effect of large sovereign debts on country default risks, or assessing the impact of extreme weather on productivity at the regional, firm, and individual levels.
\ifsubmission
We thank Frank Diebold, Christian Haefke, Bo Honor\'{e}, Roger Klein, Ulrich M\"{u}ller, Frank Schorfheide, and seminar participants at the FRB Chicago, Monash University, University of Melbourne, Reserve Bank of Australia, University of Sydney, FRB Philadelphia, UC Irvine, UCSD, and Penn State, as well as conference participants at the Applied Time Series Econometrics Workshop at the FRB St.\ Louis, Dolomiti Macro Meetings, AiE Conference and Festschrift in Honor of Joon Y.\ Park, Midwest Econometrics Group Annual Meeting, Greater New York Econometrics Colloquium, Women in Macroeconomics Workshop II, NBER Summer Institute, and NBER-NSF Time Series Conference for helpful comments and discussions. \else\fi