spectralmeas: Estimation of the Spectral Measure

Description

Kiriliouk et al. (2016, pp. 360--364) describe estimation of the spectral measure of bivariate data. Standardize the bivariate data as $X^\star$ and $Y^\star$ as in psepolar and select a “large” value for the pseudo-polar radius $S_f$ for nonexceedance probability $f$. Estimate the spectral measure $H(w)$, which is the limiting distribution of the pseudo-polar angle component $W$ given that the corresponding radial component $S$ is large: $$\mathrm{Pr}[W \in \cdot | S > S_f] \rightarrow H(w) \mbox{\ as\ } S_f \rightarrow \infty\mbox{.}$$ So, $H(w)$ is the cumulative distribution function of the spectral measure for angle $w \in (0,1)$. The $S_f$ can be specified by a nonexceedance probability $F$ for $S_f(F)$.

The estimation proceeds as follows:

Step 1: Convert the bivariate data $(X_i, Y_i)$ into $(\widehat{S}_i, \widehat{W}_i)$ by psepolar and set the threshold $S_f$ according to “$n/k$” (this part involving $k$ does not make sense in Kiriliouk et al. (2016)) where for present implementation in copBasic the $S_f$ given the $f$ by the user is based on the empirical distribution of the $\widehat{S}_i$. The empirical distribution is estimated by the Bernstein empirical distribution function from the lmomco package.

Step 2: Let $I_n$ denote the set of indices that correspond to the observations when $\widehat{S}_i \ge S_f$ and compute $N_n$ as the cardinality of $N_n = |I_n|$, which simply means the length of the vector $I_n$.

Step 3: Use the maximum Euclidean likelihood estimator, which is the third of three methods mentioned by Kiriliouk et al. (2016): $$\widehat{H}_3(w) = \sum_{i \in I_n} \hat{p}_{3,i} \times \mathbf{1}[\widehat{W}_i \le w ]\mbox{,}$$ where $\mathbf{1}[\cdot]$ is an indicator function that is only triggered if smooth=FALSE, and following the notation of Kiriliouk et al. (2016), the “3” represents maximum Euclidean likelihood estimation. The $\hat{p}_{3,i}$ are are the weights $$\hat{p}_{3,i} = \frac{1}{N_n}\bigl[ 1 - (\overline{W} - 1/2)S^{-2}_W(\widehat{W}_i - \overline{W}) \bigr]\mbox{,}$$ where $\overline{W}$ is the sample mean and $S^2_W$ is the sample variance of $\widehat{W}_i$ $$\overline{W} = \frac{1}{N_n} \sum_{i \in I_n} \widehat{W}_i\mbox{\quad and\quad } S^2_W = \frac{1}{N_n - 1} \sum_{i \in I_n} (\widehat{W}_i - \overline{W})^2\mbox{,}$$ where Kiriliouk et al. (2016, p. 363) do not show the $N_n - 1$ in the denominator for the variance but copBasic uses it because those authors state that the sample variance is used.

Step 4: A smoothed version of $\hat{H}_3(w)$ is optionally available by $$\tilde{H}_3(w) = \sum_{i \in I_n} \hat{p}_{3,i} \times \mathcal{B}(w;\, \widehat{W}_i\nu,\, (1-\widehat{W}_i)\nu)\mbox{,}$$ where $\mathcal{B}(x; p, q)$ is the cumulative distribution function of the Beta distribution for $p, q > 0$ and where $\nu > 0$ is a smoothing parameter that can be optimized by cross validation.

Step 5: The spectral density lastly can be computed optionally as $$\tilde{h}_3(w) = \sum_{i \in I_n} \hat{p}_{3,i} \times \beta(w;\, \widehat{W}_i\nu,\, (1-\widehat{W}_i)\nu)$$ where $\beta(x; p, q)$ is the probability density function (pdf) of the Beta distribution. Readers are alerted to the absence of the $\mathbf{1}[\cdot]$ indicator function in the definitions of $\tilde{H}_3(w)$ and $\tilde{h}_3(w)$. This is correct and matches Kiriliouk et al. (2016, eqs. 17.21 and 17.22) though this author was confused for a day or so by the indicator function in what is purported to be the core definition of $\hat{H}_l(w)$ where $l = 3$ in Kiriliouk et al. (2016, eq. 17.21 and 17.17).

Usage

spectralmeas(u, v=NULL, w=NULL, f=0.90, snv=FALSE,
                                smooth=FALSE, nu=100, pdf=FALSE, ...)

Value

An R

vector of $H_3(w)$, $\tilde{H}_3(w)$, or $\tilde{h}_3(w)$ is returned.

Arguments

u: Nonexceedance probability $u$ in the $X$ direction (actually the ranks are used so this can be a real-value argument as well);
v: Nonexceedance probability $v$ in the $Y$ direction (actually the ranks are used so this can be a real-value argument as well) and if NULL then u is treated as a two column R data.frame;
w: A vector of polar angle values $W \in [0,1]$ on which to compute the $H(w)$;
f: The nonexceedance probability $F(S_f)$ of the pseudo-polar radius in psepolar;
snv: Return the standard normal variate of the $H$ by the well-known transform through the quantile function of the standard normal, qnorm();
smooth: A logical to return $\tilde{H}_3(w)$ instead of $H_3(w)$;
nu: The $\nu > 0$ smoothing parameter;
pdf: A logical to return the smoothed probability density $\tilde{h}_3(w)$. If pdf=TRUE, then internally smooth=TRUE will be set; and
...: Additional arguments to pass to the psepolar function.

Author

William Asquith william.asquith@ttu.edu

References