Mixed effects models for extreme value index regression

Let ${\bm{X}}_{ij},\ i\in\mathcal{N}_{j},\ j\in\mathcal{J}$ be an independent and identically distributed (i.i.d.) random sample from some distribution. Subsequently, we assume that for each $i\in\mathcal{N}_{j}$ and $j\in\mathcal{J}$, the response $Y_{ij}$ is conditionally independently obtained from a certain conditional distribution $F_{j}(y\mid{\bm{x}})=P(Y_{ij}\leq y\mid{\bm{X}}_{ij}={\bm{x}})$ for the given ${\bm{X}}_{ij}$, where $F_{j}$ is determined for each $j\in\mathcal{J}$. In this study, we are interested in the right tail behavior of each $F_{j}$. Here, the right tail of each $F_{j}$ is modeled by the Pareto-type distribution as

$$F_{j}(y\mid{\bm{x}})=1-y^{-1/\gamma_{j}({\bm{x}})}\mathcal{L}_{j}(y;{\bm{x}}),\quad j\in\mathcal{J},$$

(2)

where $\gamma_{j}({\bm{x}})>0$ is a positive function called EVI, and $\mathcal{L}_{j}(y;{\bm{x}})$ is a conditional slowly varying function with respect to $y$ given ${\bm{x}}$, i.e., for any ${\bm{x}}$ and $s>0$, $\mathcal{L}_{j}(ys;{\bm{x}})/\mathcal{L}_{j}(y;{\bm{x}})\to 1$ as $y\to\infty$. The EVI function $\gamma_{j}({\bm{x}})$, which determines the heaviness of the right tail of $F_{j}$, is assumed here to be the classical linear model formulated as follows:

$$\gamma_{j}({\bm{x}})=\exp\left[\left({\bm{\theta}}_{j{\rm{A}}}^{0}\right)^{\top}{\bm{x}}_{\rm{A}}+\left({\bm{\theta}}_{\rm{B}}^{0}\right)^{\top}{\bm{x}}_{\rm{B}}\right],\quad j\in\mathcal{J},$$

(3)

where ${\bm{x}}\coloneqq({\bm{x}}_{\rm{A}}^{\top},{\bm{x}}_{\rm{B}}^{\top})^{\top}\in\mathbb{R}^{p_{\rm{A}}}\times\mathbb{R}^{p_{\rm{B}}}$, and ${\bm{\theta}}_{j{\rm{A}}}^{0}\in\mathbb{R}^{p_{\rm{A}}}$ and ${\bm{\theta}}_{\rm{B}}^{0}\in\mathbb{R}^{p_{\rm{B}}}$ are regression coefficient vectors. Note that ${\bm{\theta}}_{j{\rm{A}}}^{0}$ is different between areas, whereas ${\bm{\theta}}_{{\rm{B}}}^{0}$ is common to all areas. When $J=1$, the above model is reduced to that of Wang and Tsai (2009). The purpose of the model (2) with (3) is to estimate the parameter vectors ${\bm{\theta}}_{1{\rm{A}}}^{0},{\bm{\theta}}_{2{\rm{A}}}^{0},\ldots,{\bm{\theta}}_{J{\rm{A}}}^{0}$ and ${\bm{\theta}}_{\rm{B}}^{0}$. However, this model has $(J\times p_{\rm{A}}+p_{\rm{B}})$ parameters; hence, if $J$ is large, the associated estimators may be severely biased (see, Section 4 of Ruppert, Wand, and Carroll 2003; Broström and Holmberg 2011). To overcome this bias problem, we use the MEM instead of the fully parametric model (2) with (3).

For $p_{\rm{A}}\leq p$, we introduce the random effects ${\bm{U}}_{j}\in\mathbb{R}^{p_{\rm{A}}},\ j\in\mathcal{J}$ such that

$${\bm{U}}_{1},{\bm{U}}_{2},\ldots,{\bm{U}}_{J}\overset{{\rm{i.i.d.}}}{\sim}N({\bm{0}},{\bm{\Sigma}}_{0}),$$

(4)

where $N({\bm{0}},{\bm{\Sigma}}_{0})$ refers to the multivariate normal distribution with zero mean vector and unknown covariance matrix ${\bm{\Sigma}}_{0}$. The MEM uses these random effects to express the differences in $F_{j}$ between areas as a latent structure. Let $F(y\mid{\bm{u}}_{j},{\bm{x}})\coloneqq P(Y_{ij}\leq y\mid{\bm{U}}_{j}={\bm{u}}_{j},{\bm{X}}_{ij}={\bm{x}})$ be the conditional distribution function of $Y_{ij}$ given ${\bm{U}}_{j}={\bm{u}}_{j}$ and ${\bm{X}}_{ij}={\bm{x}}$. In this expression, the information about the differences in $F_{j}$ between areas is assigned to $F$ by ${\bm{u}}_{j}$. Note that the random effects ${\bm{U}}_{j},\ j\in\mathcal{J}$ are not observed as data.

As an alternative model to (2) with (3), the Pareto-type distribution using the MEM is defined as follows:

$$F(y\mid{\bm{u}}_{j},{\bm{x}})=1-y^{-1/\gamma({\bm{u}}_{j},{\bm{x}})}\mathcal{L}(y;{\bm{u}}_{j},{\bm{x}}),\quad j\in\mathcal{J},$$

(5)

where $\mathcal{L}(y;{\bm{u}}_{j},{\bm{x}})$ conditional on ${\bm{U}}_{j}={\bm{u}}_{j}$ and ${\bm{X}}_{ij}={\bm{x}}$ is a slowly varying function with respect to $y$. Then, as an extension of (3) to the MEM, the EVI function $\gamma({\bm{u}}_{j},{\bm{x}})$ is assumed to be

$$\gamma({\bm{u}}_{j},{\bm{x}})=\exp\left[\left({\bm{\theta}}_{\rm{A}}^{0}+{\bm{u}}_{j}\right)^{\top}{\bm{x}}_{\rm{A}}+\left({\bm{\theta}}_{\rm{B}}^{0}\right)^{\top}{\bm{x}}_{\rm{B}}\right],\quad j\in\mathcal{J}.$$

(6)

Compared to (3), the area-wise differences in the slopes of log-EVI with respect to ${\bm{X}}_{{\rm{A}}ij}$ are represented by ${\bm{u}}_{j},\ j\in\mathcal{J}$ as a latent structure. Thus, the total number of parameters in the model (5) using (6) is $p+p_{\rm{A}}(p_{\rm{A}}+1)/2$, which is independent of $J$ and less than that of the fully parametric model (2) with (3) when $J$ is large. Here, $p_{\rm{A}}(p_{\rm{A}}+1)/2$ is the number of parameters included in ${\bm{\Sigma}}_{0}$.

The simplest model of (6) is the location-shifting MEM with $p_{\rm{A}}=1$ and ${\bm{X}}_{{\rm{A}}ij}\equiv 1$, denoted ${\bm{\theta}}_{\rm{A}}^{0}$ and ${\bm{u}}_{j}$ by the scalars $\theta_{\rm{A}}^{0}$ and $u_{j}$,

$$\gamma(u_{j},{\bm{x}}_{\rm{B}})=\exp\left[\theta_{\rm{A}}^{0}+u_{j}+\left({\bm{\theta}}_{\rm{B}}^{0}\right)^{\top}{\bm{x}}_{\rm{B}}\right],\quad j\in\mathcal{J},$$

(7)

which can be regarded as an EVI regression version of the nested error regression model (see, Battese, Harter, and Fuller 1988). The model (7) indicates that the intercept of $\log\gamma$ accepts the heterogeneity between areas, although each covariate has the common slope of $\log\gamma$ across all areas. The nested error regression model is useful for SAE (see, Diallo and Rao 2018; Sugasawa and Kubokawa 2020). The application in Section 4 demonstrates the analysis using the model (7) and verifies its effectiveness numerically. Alternatively, the case $p=p_{\rm{A}}$ yields the most complicated model of (6), indicating that the slope of $\log\gamma$ with respect to each covariate varies across areas.

In this section, we construct estimators of the unknown parameters $\{{\bm{\theta}}_{\rm{A}}^{0},{\bm{\theta}}_{\rm{B}}^{0},{\bm{\Sigma}}_{0}\}$ included in the model (5) with (6).

Let $F_{\omega_{j}}(y\mid{\bm{u}}_{j},{\bm{x}})\coloneqq P(Y_{ij}\leq y\mid{\bm{U}}_{j}={\bm{u}}_{j},{\bm{X}}_{ij}={\bm{x}},Y_{ij}>\omega_{j})$ be the conditional distribution function given ${\bm{U}}_{j}={\bm{u}}_{j}$, ${\bm{X}}_{ij}={\bm{x}}$ and $Y_{ij}>\omega_{j}$, where $\omega_{j}\in\mathbb{R}^{+},\ j\in\mathcal{J}$ are thresholds. In this paper, we assume that $\mathcal{L}(y;{\bm{u}}_{j},{\bm{x}})$ belongs to the Hall class (see, Hall 1982), which is mentioned in (A1) of Section 3.1. From this, we have

$$F_{\omega_{j}}(y\mid{\bm{u}}_{j},{\bm{x}})\approx 1-\left(\frac{y}{\omega_{j}}\right)^{-1/\gamma({\bm{u}}_{j},{\bm{x}})},\quad j\in\mathcal{J}$$

(8)

for sufficiently large $\omega_{j}$. Using (8) instead of (5), we can remove ${\mathcal{L}}$ for the estimation of $\gamma$. Similarly, from (8) and the assumption (A1), the density of $Y_{ij}$ given ${\bm{U}}_{j}={\bm{u}}_{j}$, ${\bm{X}}_{ij}={\bm{x}}$ and $Y_{ij}>\omega_{j}$ is obtained as follows:

$$f_{w_{j}}(y\mid{\bm{u}}_{j},{\bm{x}})\approx\omega_{j}^{-1}\gamma({\bm{u}}_{j},{\bm{x}})^{-1}\left(\frac{y}{\omega_{j}}\right)^{-1/\gamma({\bm{u}}_{j},{\bm{x}})-1},\quad j\in\mathcal{J}.$$

(9)

Wang and Tsai (2009) considered the similar approximated distribution (8) and density (9) in linear extreme value index regression. Estimation using data exceeding thresholds is the so-called peak-over-threshold method (see, Hill 1975; Wang and Tsai 2009).

We assume that ${\bm{U}}_{j}$ and ${\bm{X}}_{ij}$ are independent for $i\in\mathcal{N}_{j}$ and $j\in\mathcal{J}$. Furthermore, we assume that $Y_{ij}$ given ${\bm{U}}_{j}$ and ${\bm{X}}_{ij}$ is conditionally independent for $i\in\mathcal{N}_{j}$ and $j\in\mathcal{J}$ and has the distribution function (5) (see, Jiang, Wand, and Bhaskaran 2022). Using (9), we then define the likelihood of $({\bm{\theta}}_{\rm{A}}^{0},{\bm{\theta}}_{\rm{B}}^{0},{\bm{\Sigma}}_{0})$ as

$$L({\bm{\theta}}_{\rm{A}},{\bm{\theta}}_{\rm{B}},{\bm{\Sigma}})\coloneqq\prod_{j=1}^{J}E_{{\bm{U}}_{j}}\left[\prod_{i\in\mathcal{N}_{j}:Y_{ij}>\omega_{j}}f_{w_{j}}(Y_{ij}\mid{\bm{U}}_{j},{\bm{X}}_{ij})\right],$$

where $E_{{\bm{U}}_{j}}$ denotes the expectation over the random effects distribution, $({\bm{\theta}}_{\rm{A}}^{\top},{\bm{\theta}}_{\rm{B}}^{\top})^{\top}\in\mathbb{R}^{p_{\rm{A}}}\times\mathbb{R}^{p_{\rm{B}}}$ is any vector corresponding to $(({\bm{\theta}}_{\rm{A}}^{0})^{\top},({\bm{\theta}}_{\rm{B}}^{0})^{\top})^{\top}$, and ${\bm{\Sigma}}\in\mathbb{R}^{p_{\rm{A}}\times p_{\rm{A}}}$ is any positive definite matrix corresponding to ${\bm{\Sigma}}_{0}$. The above $L$ is derived from the standard definition of the likelihood for the MEM, because ${\bm{U}}_{j},\ j\in\mathcal{J}$ are unobserved random variables, unlike the data (1) (see, Section 2 of Wu 2009). From (4), (6), (9) and (A1), the log-likelihood of $({\bm{\theta}}_{\rm{A}}^{0},{\bm{\theta}}_{\rm{B}}^{0},{\bm{\Sigma}}_{0})$ can be expressed as

	$\displaystyle\ell({\bm{\theta}}_{\rm{A}},{\bm{\theta}}_{\rm{B}},{\bm{\Sigma}})$	$\displaystyle\coloneqq\log L({\bm{\theta}}_{\rm{A}},{\bm{\theta}}_{\rm{B}},{\bm{\Sigma}})$
	$\displaystyle\begin{split}&\approx\sum_{j=1}^{J}\log\int_{\mathbb{R}^{p_{\rm{A}}}}\phi({\bm{u}};{\bm{0}},{\bm{\Sigma}})\exp\left(\sum_{i=1}^{n_{j}}\biggl{\{}-\left({\bm{\theta}}_{\rm{A}}+{\bm{u}}\right)^{\top}{\bm{X}}_{{\rm{A}}ij}-{\bm{\theta}}_{\rm{B}}^{\top}{\bm{X}}_{{\rm{B}}ij}\biggr{.}\right.\\ &\quad\Biggl{.}\left.-\exp\left[-\left({\bm{\theta}}_{\rm{A}}+{\bm{u}}\right)^{\top}{\bm{X}}_{{\rm{A}}ij}-{\bm{\theta}}_{\rm{B}}^{\top}{\bm{X}}_{{\rm{B}}ij}\right]\log\frac{Y_{ij}}{\omega_{j}}\right\}I(Y_{ij}>\omega_{j})\Biggr{)}d{\bm{u}}+C,\end{split}$			(10)

where $I(\cdot)$ is an indicator function that returns 1 if $Y_{ij}>\omega_{j}$ and 0 otherwise, $\phi(\cdot;{\bm{0}},{\bm{\Sigma}})$ is a density function of $N({\bm{0}},{\bm{\Sigma}})$, and $C$ is a suitable constant independent of $({\bm{\theta}}_{\rm{A}},{\bm{\theta}}_{\rm{B}},{\bm{\Sigma}})$. Again, because ${\bm{U}}_{j},\ j\in\mathcal{J}$ are not observed as data, the log-likelihood (10) includes the integral over the domain $\mathbb{R}^{p_{\rm{A}}}$ of the random effects. We denote the approximate maximum likelihood estimator of $({\bm{\theta}}_{\rm{A}}^{0},{\bm{\theta}}_{\rm{B}}^{0},{\bm{\Sigma}}_{0})$ by $(\hat{\bm{\theta}}_{\rm{A}},\hat{\bm{\theta}}_{\rm{B}},\hat{\bm{\Sigma}})$, which is the maximizer of the right-hand side of (10).

In the proposed model (5) using (6), we are not only interested in estimating the parameters $\{{\bm{\theta}}_{\rm{A}}^{0},{\bm{\theta}}_{\rm{B}}^{0},{\bm{\Sigma}}_{0}\}$, but also in predicting the random effects ${\bm{U}}_{j},\ j\in\mathcal{J}$. Here, we propose the conditional mode method to predict these random effects (see, Santner and Duffy 1989; Section 11 of Wu 2009). Now, the conditional density function of $({\bm{U}}_{1},{\bm{U}}_{2},\ldots,{\bm{U}}_{J})$ given the data (1) is proportional to

$$\prod_{j=1}^{J}\left[\phi({\bm{u}}_{j};{\bm{0}},{\bm{\Sigma}}_{0})\prod_{i\in\mathcal{N}_{j}:Y_{ij}>\omega_{j}}f_{\omega_{j}}(Y_{ij}\mid{\bm{u}}_{j},{\bm{X}}_{ij})\right],$$

as a function of $({\bm{u}}_{1},{\bm{u}}_{2},\ldots,{\bm{u}}_{J})$. Then, the predictor of ${\bm{U}}_{j}$ is defined as the mode of this conditional distribution, i.e.,

$$\tilde{\bm{u}}_{j}\coloneqq\operatorname*{argmax}_{{\bm{u}}_{j}\in\mathbb{R}^{p_{\rm{A}}}}\ \phi({\bm{u}}_{j};{\bm{0}},{\bm{\Sigma}}_{0})\prod_{i\in\mathcal{N}_{j}:Y_{ij}>\omega_{j}}f_{\omega_{j}}(Y_{ij}\mid{\bm{u}}_{j},{\bm{X}}_{ij}),\quad j\in\mathcal{J},$$

(11)

where $f_{\omega_{j}}$ and $({\bm{\theta}}_{\rm{A}}^{0},{\bm{\theta}}_{\rm{B}}^{0},{\bm{\Sigma}}_{0})$ included in $f_{\omega_{j}}$ are replaced by (9) and the estimator $(\hat{\bm{\theta}}_{\rm{A}},\hat{\bm{\theta}}_{\rm{B}},\hat{\bm{\Sigma}})$, respectively.

The thresholds $\omega_{j},\ j\in\mathcal{J}$ in (10) are tuning parameters that balance between the quality of the approximation (9) and the amount of data exceeding the thresholds. By setting higher thresholds, we can generally improve the estimation bias, but the estimator becomes more unstable. Conversely, if we lower the thresholds, the estimator will behave more stably, but may be more biased. Therefore, these thresholds $\omega_{j},\ j\in\mathcal{J}$ should be chosen appropriately, considering this trade-off relationship. Here, for each area, we apply the discrepancy measure (see, Wang and Tsai 2009), which considers the goodness of fit of the model, to select the area-wise optimal threshold. In the simulation studies in Section 1.2 of our supplementary material, we verify the performance of the proposed estimator using this threshold selection method.

Let $n_{j0}\coloneqq\sum_{i=1}^{n_{j}}I(Y_{ij}>\omega_{j})$, which is the sample size exceeding the threshold $\omega_{j}$ for the $j$th area. Additionally, we define $n_{0}\coloneqq J^{-1}\sum_{j=1}^{J}n_{j0}$ as the average of the effective sample sizes of all areas. Note that $n_{j0},\ j\in\mathcal{J}$ and $n_{0}$ are random variables, not constants. In the case (iii) defined above, we assume that for each $j\in\mathcal{J}$, the threshold $\omega_{j}$ diverges to infinity in tandem with the sequence of $J$ and the $j$th within-area sample size $n_{j}$. Accordingly, we denote $\omega_{j}$ by $\omega_{(J,n_{j})}$. The asymptotic properties of the proposed estimator $(\hat{\bm{\theta}}_{\rm{A}},\hat{\bm{\theta}}_{\rm{B}},\hat{\bm{\Sigma}})$ rely on the following assumptions (A1)-(A6):

1. (A1)

   $\mathcal{L}(y;{\bm{u}},{\bm{x}})$ in (5) belongs to the Hall class (see, Hall 1982), that is,

   $$\mathcal{L}(y;{\bm{u}},{\bm{x}})=c_{0}({\bm{u}},{\bm{x}})+c_{1}({\bm{u}},{\bm{x}})y^{-\beta({\bm{u}},{\bm{x}})}+\lambda(y;{\bm{u}},{\bm{x}}),$$

   (12)

   where $c_{0}({\bm{u}},{\bm{x}})>0$, $c_{1}({\bm{u}},{\bm{x}})$, $\beta({\bm{u}},{\bm{x}})>0$ and $\lambda(y;{\bm{u}},{\bm{x}})$ are continuous and bounded. Furthermore, $\lambda(y;{\bm{u}},{\bm{x}})$ satisfies

   $$\sup_{{\bm{u}}\in\mathbb{R}^{p_{\rm{A}}},\ {\bm{x}}\in\mathbb{R}^{p}}\left[y^{\beta({\bm{u}},{\bm{x}})}\lambda(y;{\bm{u}},{\bm{x}})\right]\to 0\quad{\rm{as}}\quad y\to\infty.$$

2. (A2)

   There exists a bounded and continuous function $\delta:\mathbb{R}^{p_{\rm{A}}}\times\mathbb{R}^{p}\to\mathbb{R}^{+}$ such that

   $$\sup_{{\bm{u}}\in\mathbb{R}^{p_{\rm{A}}},\ {\bm{x}}\in\mathbb{R}^{p}}\left\lvert\frac{P(Y_{ij}>y\mid{\bm{U}}_{j}={\bm{u}},{\bm{X}}_{ij}={\bm{x}})}{P(Y_{ij}>y\mid{\bm{U}}_{j}={\bm{u}})}-\delta({\bm{u}},{\bm{x}})\right\lvert\to 0\quad{\rm{as}}\quad y\to\infty.$$

3. (A3)

   As $n_{j}\to\infty,\ j\in\mathcal{J}$ and $J\to\infty$,

   $$\inf_{j\in\mathcal{J},\ {\bm{u}}\in\mathbb{R}^{p_{\rm{A}}}}n_{j}P(Y_{ij}>\omega_{(J,n_{j})}\mid{\bm{U}}_{j}={\bm{u}})\to\infty.$$

4. (A4)

   There exist some bounded and continuous functions $d_{j}:\mathbb{R}^{p_{\rm{A}}}\to\mathbb{R}^{+},\ j\in\mathcal{J}$ such that under given ${\bm{U}}_{j}={\bm{u}}$, $n_{j0}/n_{0}\to^{P}d_{j}({\bm{u}})$ uniformly for all $j\in\mathcal{J}$ and ${\bm{u}}\in\mathbb{R}^{p_{\rm{A}}}$ as $n_{j}\to\infty,\ j\in\mathcal{J}$ and $J\to\infty$, where the symbol “$\to^{P}$” represents convergence in probability.

5. (A5)

   $n_{0}/J\to^{P}0$ as $n_{j}\to\infty,\ j\in\mathcal{J}$ and $J\to\infty$.

6. (A6)

   There exist some bounded and continuous functions ${\bm{b}}_{{\rm{K}}j}:\mathbb{R}^{p_{\rm{A}}}\to\mathbb{R},\ j\in\mathcal{J},\ {\rm{K}}\in\{{\rm{A}},{\rm{B}}\}$ such that

   $$\sup_{j\in\mathcal{J},\ {\bm{u}}\in\mathbb{R}^{p_{\rm{A}}}}\left\lVert\frac{J^{1/2}n_{0}^{1/2}{E_{{\bm{X}}_{ij}}}\left[{\bm{X}}_{{\rm{K}}ij}\zeta_{j}({\bm{u}},{\bm{X}}_{ij})\right]}{P(Y_{ij}>\omega_{(J,n_{j})}\mid{\bm{U}}_{j}={\bm{u}})^{1/2}}-{\bm{b}}_{{\rm{K}}j}({\bm{u}})\right\rVert\to 0\\$$

   as $n_{j}\to\infty,\ j\in\mathcal{J}$ and $J\to\infty$, where

   $$\zeta_{j}({\bm{u}},{\bm{x}})\coloneqq\frac{c_{1}({\bm{u}},{\bm{x}})\gamma({\bm{u}},{\bm{x}})\beta({\bm{u}},{\bm{x}})}{1+\gamma({\bm{u}},{\bm{x}})\beta({\bm{u}},{\bm{x}})}\omega_{(J,n_{j})}^{-1/\gamma({\bm{u}},{\bm{x}})-\beta({\bm{u}},{\bm{x}})}$$

   and $\lVert\cdot\rVert$ refers to the Euclidean norm.

(A1) and (A2) regularize the tail behavior of the conditional response distribution (5) (see, Wang and Tsai 2009; Ma, Jiang, and Huang 2019). (A3)-(A6) impose the constraints on the divergence rates of the thresholds $\omega_{(J,n_{j})},\ j\in\mathcal{J}$. (A3) implies that for each $j\in\mathcal{J}$, the effective sample size $n_{j0}$ asymptotically diverges to infinity. Under (A4), $n_{j0},\ j\in\mathcal{J}$ are not critically different. Furthermore, (A5) means that the number of areas $J$ is relatively larger than the effective sample sizes $n_{j0},\ j\in\mathcal{J}$. (A5) mathematically links the divergence rates of $\omega_{(J,n_{j})},\ j\in\mathcal{J}$ and $J$. (A6) is related to the asymptotic bias of the proposed estimator. If (A6) fails, the consistency of the proposed estimator may not be guaranteed.

Let ${\bm{M}}$ be a matrix of zeros and ones such that ${\bm{M}}{\rm{vech}}({\bm{A}})={\rm{vec}}({\bm{A}})$ for all symmetric matrices ${\bm{A}}\in\mathbb{R}^{p_{\rm{A}}\times p_{\rm{A}}}$, where ${\rm{vec}}(\cdot)$ is a vector operator, and ${\rm{vech}}(\cdot)$ is a vector half operator that stacks the lower triangular half of a given $d\times d$ square matrix into the single vector of length $d(d+1)/2$ (see, Magnus and Neudecker 1988). The Moore-Penrose inverse of ${\bm{M}}$ is ${\bm{M}}_{*}\coloneqq({\bm{M}}^{\top}{\bm{M}})^{-1}{\bm{M}}^{\top}$.

For the maximum likelihood estimator $(\hat{\bm{\theta}}_{\rm{A}},\hat{\bm{\theta}}_{\rm{B}},\hat{\bm{\Sigma}})$, we obtain the following Theorem 1.

As described in Section 2.1, an important example of (6) is the type of nested error regression model (7). For the model (7), Theorem 1 can be simplified to the following Corollary 1. Let define $\sigma_{0}^{2}\coloneqq{\rm{var}}[U_{j}]$ for the random effects $U_{j},\ j\in\mathcal{J}$ and denote its proposed estimator as $\hat{\sigma}^{2}$.

In this section, we analyze the returns of a large number of cryptocurrencies. Over the past decade, the cryptocurrency market has grown rapidly, led by Bitcoin. As of January 2024, the number of active cryptocurrency stocks is estimated to be over 9,000 (see, https://www.statista.com/statistics/863917/number-crypto-coins-tokens/). However, despite the large number of stocks, most studies are limited to a few major cryptocurrencies such as Bitcoin. Gkillas and Katsiampa (2018) studied the extreme value analysis of the returns of five cryptocurrencies, namely Bitcoin, Ethereum, Ripple, Litecoin and Bitcoin Cash. We extend their work to a larger number of stocks including these five cryptocurrencies. Our goal is to demonstrate the advantage of simultaneously analyzing many cryptocurrency stocks using the MEM. In Section 4.3, we report the interpretation of our model. In Section 4.4, we compare the performance of our method with that of stock-by-stock analysis using the method of Wang and Tsai (2009).

The cryptocurrency dataset is available on CoinMarketCap (see, https://coinmarketcap.com/). We use this dataset from the last ten years, from January 1, 2014 to December 31, 2023. Then, we cover the 413 stocks that are in the top 500 on CoinMarketCap as of February 2, 2024 and have returns of 364 days or more. From the daily closing prices for each stock, the returns can be calculated as $\log P_{tj}-\log P_{(t-1)j}$, where $P_{tj}$ refers to the closing price on the $t$th day in the $j$th cryptocurrency. Let the non-missing returns of the $j$th cryptocurrency be denoted as $\{Y_{ij},\ i=1,2,\ldots,n_{j}\}$ for $j=1,2,\ldots,413$. The sample size $n_{j}$ of the $j$th cryptocurrency is primarily determined by the launch date. Note that the index $i$ does not represent the same date for all stocks.

We examine the tail behavior of both the negative and positive returns for each stock. In many applications of EVT, the sample kurtosis has been used to check the effectiveness of using the Pareto-type distribution (see, Wang and Tsai 2009; Ma, Jiang, and Huang 2019). In this study, the sample kurtosis of $\{Y_{ij}\}_{i\in\mathcal{N}_{j}}$ was greater than zero for each $j=1,2,\ldots,413$, and the minimum sample kurtosis for the 413 stocks was 1.833, implying that the returns of each cryptocurrency are heavily distributed. Therefore, we use the Pareto-type distribution instead of the GPD to analyze high threshold exceedances for the both negative and positive returns. The following Sections 4.3 and 4.4 describe our method only for positive returns. However, using the same method by replacing $\{Y_{ij}\}_{i\in\mathcal{N}_{j},j\in\mathcal{J}}$ with $\{-Y_{ij}\}_{i\in\mathcal{N}_{j},j\in\mathcal{J}}$, we can also implement the analysis of extreme negative returns.

For each pair of the 413 stocks, we estimated the tail dependence parameter from returns above the 95th percentile for the both negative and positive returns (see, Reiss and Thomas 2007). The results showed that the percentage of combinations with a tail dependence over 0.7 was only 0.074% for negative returns and only 0.032% for positive returns. Therefore, in this analysis, we do not consider the dependence between each pair of stocks.

In many cryptocurrency applications, dummy variables such as year, month, and day of the week have often been utilized as covariates to account for the non-stationarity of returns. In Longin and Pagliardi (2016) and Gkillas and Katsiampa (2018), the returns were adjusted in terms of variance using some dummy variables and were analyzed under the assumption that the EVI of each cryptocurrency remained constant over the period, despite the dramatic growth of the market. However, in their approach, it is unclear which variable influences returns to what extent. Therefore, in our application, we employ EVI regression to examine the effects of dummy variables with year (9 variables), month (11 variables) and day of the week (6 variables) on the tail behavior of the distribution of returns. The second column of Table 1 shows the assignment of 26 dummy variables. We denote the vector of these dummy variables as ${\bm{X}}_{ij}\in\prod_{k=1}^{26}\{0,1\}\subset\mathbb{R}^{26}$ for $i=1,2,\ldots,n_{j}$ and $j=1,2,\ldots,413$. For our model in Section 2.1, we set ${\bm{X}}_{{\rm{A}}ij}\equiv 1$ and ${\bm{X}}_{{\rm{B}}ij}={\bm{X}}_{ij}$ and denote the random effects as $U_{1},U_{2},\ldots,U_{413}\stackrel{{\scriptstyle\rm{i.i.d.}}}{{\sim}}N(0,\sigma^{2})$, where $\sigma^{2}>0$ is an unknown variance. Then, for the underlying Pareto-type distribution (5), the EVI (6) conditional on $U_{j}=u_{j}$ and ${\bm{X}}_{ij}={\bm{x}}$ is modeled as

$$\gamma(u_{j},{\bm{x}})=\exp\left(\theta_{\rm{A}}+u_{j}+{\bm{\theta}}_{\rm{B}}^{\top}{\bm{x}}\right),\quad j=1,2,\ldots,413,$$

(14)

where $\theta_{\rm{A}}\in\mathbb{R}$ and ${\bm{\theta}}_{\rm{B}}\in\mathbb{R}^{26}$ are unknown coefficients. Below, we analyze the returns of the 413 cryptocurrencies simultaneously using (14). For the implementation, we then use the glmer() function in the package lme4 (see, Bates et al. 2015) within the R computing environment (see, R Core Team 2021). A more detailed explanation is given in Section 1 of our supplementary material.