Here are some notes illustrating how to fit the spectral density ratio model in de Carvalho and Davison (2014). The methods can be implemented with the `alabama`

package or with a heuristic two-step approach described in Section 3.2 of the paper. The file *aux.R* includes auxiliary routines which I will be invoking below, such as simulation from the beta-spectral density ratio model (using `rBeta`

), the log-likelihood function (`LL`

), etc. Some notation in the codes may not match that used in the final version of the paper.

```
rm(list = ls())
require(emplik); require(nnet); gamma <- 0
source("aux.R")
```

The goal of the example below will be on illustrating the spectral density ratio model and to compare it with the empirical likelihood spectral measure Einmahl and Segers (2009); the data are simulated according to the following settings:

```
n0 <- 25 ; N <- c(50, 25) ; n <- n0 + sum(N)
rho <- N / n0 ; K <- length(N) + 1 ; K1 <- length(N)
D <- 2
## set true values of the tilt parameters
a0 <- 0.1 ; a <- c(1, 5)
## fix seed for reproducibility
set.seed(1)
```

We simulate the data using the routine from , and use the distortion function \(c(w) = \log\{w (1 - w)\}\):

```
v <- z <- rBeta(c(n0, N), c(a0, a))
V <- log(v * (1 - v))
Y <- c(rep(0, n0), rep(1, N[1]),rep(2, N[2]))
```

The heuristic two-stage approach can be executed using:

```
## nested approach: outer
fit <- multinom(Y ~ V)
```

```
## # weights: 9 (4 variable)
## initial value 109.861229
## iter 10 value 74.784268
## final value 72.762709
## converged
```

```
gamma0 <- c(coefficients(fit)[1, 1] + log(n0 / N[1]),
coefficients(fit)[2, 1] + log(n0 / N[2]),
coefficients(fit)[1, 2], coefficients(fit)[2, 2])
## nested approach: inner
delta0 <- inner(v, gamma0)
```

Simulation results in the paper suggest that this heuristic provides a good approximation of the tilting parameters. Yet, such inferences do not obey the normalization and moment constraints:

`check.constraints(v, gamma0, delta0)`

```
## $normalization
## [1] 1.0196937 0.9917950 0.9967164
##
## $moments
## [1] 0.5098469 0.4958975 0.4983582
```

Full optimization can be obtained using the `alabama`

package; the `auglag`

routine may yield some warnings. These are due to the fact that the algorithm attempts to search outside the domain of the objective function, and thus can be safely ignored.

```
require(alabama)
param <- try(auglag(par = c(gamma0, delta0), fn = LL,
heq = constraints)$par)
# if(is.na(param[2]) == TRUE)
# param <- try(auglag(par = c(gamma0,rep(0, K)), fn = LL,
# heq = constraints)$par, silent = TRUE)
alpha <- as.numeric(param[1:(K - 1)])
beta <- as.numeric(param[K:(2 * (K - 1))])
gamma <- as.numeric(c(alpha, beta))
delta <- as.numeric(param[(2 * (K - 1) + 1):(2 * (K - 1) + K)])
```

The inferences arising from the full optimization would obey moment constraints:

`check.constraints(v, c(alpha, beta), delta)`

```
## $normalization
## [1] 1 1 1
##
## $moments
## [1] 0.5 0.5 0.5
```

To compare with the maximum empirical likelihood estimator, we compute the Lagrange multipliers:

```
delta0 <- el.test(v[1:n0], 1 / 2)$lambda
delta1 <- el.test(v[(n0 + 1):(n0 + N[1])], 1 / 2)$lambda
delta2 <- el.test(v[(n0 + N[1] + 1):n], 1 / 2)$lambda
```

The results are summarized in Figure 1. On a final note, an alternative approach has been recently proposed by Castro and de Carvalho (2017).

Castro, D., and de Carvalho, M. (2017), “Spectral Density Regression for Bivariate Extremes.” *Stochastic Environmental Research and Risk Assessment*, 31, 1603–1613.

de Carvalho, M., and Davison, A. C. (2014), “Spectral Density Ratio Models for Multivariate Extremes.” *Journal of the American Statistical Association*, 109, 764–776.

Einmahl, J. H. J., and Segers, J. (2009), “Maximum Empirical Likelihood Estimation of the Spectral Measure of an Extreme-Value Distribution.” *The Annals of Statistics*, 37, 2953–89.