fExtDepSpat: Fitting of a max-stable process

Description

This function uses the Stephenson-Tawn likelihood to estimate parameters of max-stable models.

Usage

fExtDepSpat(
  x, model, sites, hit, jw, thresh, DoF, range, smooth,
  alpha, par0, acov1, acov2, parallel, ncores, args1, args2,
  seed = 123, method = "BFGS", sandwich = TRUE,
  control = list(trace = 1, maxit = 50, REPORT = 1, reltol = 0.0001)
)
# S3 method for ExtDep_Spat
logLik(object, ...)

Value

fExtDepSpat: An object of class ExtDep_Spat in the form of a list comprising:

the vector of estimated parameters (est),
the composite likelihood order (jw),
the maximised log-likelihood value (LL),
if sandwich = TRUE, the standard errors from the sandwich information matrix (stderr.sand),
the TIC for model selection (TIC),
if the composite likelihood is considered, a matrix of all tuples with weight 1 (cmat).

logLik: The value of the maximised log-likelihood.

Arguments

x: A $(n \times d)$ matrix containing $n$ observations at $d$ locations (fExtDepSpat).
object: An object of class ExtDep_Spat (logLik).
model: A character string indicating the max-stable model, currently extremal-t ("ET") and extremal skew-t ("EST") only available. Note that the Schlather model can be obtained as a special case by specifying the extremal-t model with DoF = 1.
sites: A $(d \times 2)$ matrix corresponding to the coordinates of locations where the processes are simulated. Each row corresponds to a location.
hit: A $(n \times d)$ matrix containing the hitting scenarios for each observation.
jw: An integer between $2$ and $d$, indicating the tuples considered in the composite likelihood. If jw = d the full likelihood is considered.
thresh: A positive real indicating the threshold value for pairwise distances. See details.
DoF: A positive real indicating a fixed value of the degree of freedom of the extremal-t and extremal skew-t models.
range: A positive real indicating a fixed value of the range parameter for the power exponential correlation function (only correlation function currently available).
smooth: A positive real in $(0, 2]$ indicating a fixed value of the smoothness parameter for the power exponential correlation function (only correlation function currently available).
alpha: A vector of length $3$ indicating fixed values of the skewness parameters $\alpha_0$, $\alpha_1$, and $\alpha_2$ for the extremal skew-t model. If some components are NA, then the corresponding parameter will be estimated.
par0: A vector of initial values of the parameter vector, in order: the degree of freedom $\nu$, the range $r$, the smoothness $\eta$, and the skewness parameters $\alpha_0$, $\alpha_1$. Its length depends on the model and the number of fixed parameters.
acov1, acov2: Vectors of length $d$ representing covariates to model the skewness parameter of the extremal skew-t model.
parallel: A logical value; if TRUE parallel computing is done.
ncores: An integer indicating the number of cores considered in the parallel socket cluster of type 'PSOCK', based on the makeCluster routine from the parallel package. Required if parallel = TRUE.
args1, args2: Lists specifying details about the Monte Carlo simulation scheme to compute multivariate CDFs. See details.
seed: An integer for reproducibility in the CDF computations.
method: A character string indicating the optimisation method to be used. See optim for more details.
sandwich: A logical value; if TRUE the standard errors of the estimates are computed from the Sandwich (Godambe) information matrix.
control: A list of control parameters for the optimisation. See optim for more details.
...: For the logLik() function: can provide a digits argument, an integer indicating the number of decimal places to be reported.

Author

Simone Padoan, simone.padoan@unibocconi.it, https://faculty.unibocconi.it/simonepadoan/; Boris Beranger, borisberanger@gmail.com, https://www.borisberanger.com

Details

This routine follows the methodology developed by Beranger et al. (2021). It uses the Stephenson-Tawn likelihood which relies on the knowledge of time occurrences of each block maxima. Rather than considering all partitions of the set $\{1, \ldots, d\}$, the likelihood is computed using the observed partition. Let $\Pi = (\pi_1, \ldots, \pi_K)$ denote the observed partition, then the Stephenson-Tawn likelihood is given by $$ L(\theta; x) = \exp \left\{ - V(x; \theta) \right\} \times \prod^{K}_{k = 1} - V_{\pi_k}(x; \theta), $$ where $V_{\pi}$ represents the partial derivative(s) of $V(x; \theta)$ with respect to $\pi$.

When jw = d the full Stephenson-Tawn likelihood is considered, whereas for values lower than $d$ a composite likelihood approach is taken.

The argument thresh is required when the composite likelihood is used. A tuple of size jw is assigned a weight of one if the maximum pairwise distance between corresponding locations is less than thresh, and a weight of zero otherwise.

Arguments args1 and args2 relate to specifications of the Monte Carlo simulation scheme to compute multivariate CDF evaluations. These should take the form of lists including the minimum and maximum number of simulations used (Nmin and Nmax), the absolute error (eps), and whether the error should be controlled on the log-scale (logeps).

References

Beranger, B., Stephenson, A. G., and Sisson, S. A. (2021). High-dimensional inference using the extremal skew-t process. Extremes, 24, 653--685.

Examples

Run this code

set.seed(14342)

# Simulation of 20 locations
Ns <- 20
sites <- matrix(runif(Ns * 2) * 10 - 5, nrow = Ns, ncol = 2)
for (i in 1:2) sites[, i] <- sites[, i] - mean(sites[, i])

# Simulation of 50 replicates from the Extremal-t model
Ny <- 50
x <- rExtDepSpat(
  Ny, sites, model = "ET", cov.mod = "powexp", DoF = 1,
  range = 3, nugget = 0, smooth = 1.5,
  control = list(method = "exact")
)

# Fit the extremal-t using the full Stephenson-Tawn likelihood
args1 <- list(Nmax = 50L, Nmin = 5L, eps = 0.001, logeps = FALSE)
args2 <- list(Nmax = 500L, Nmin = 50L, eps = 0.001, logeps = TRUE)

if (FALSE) {
fit1 <- fExtDepSpat(
  x = x$vals, model = "ET", sites = sites, hit = x$hits,
  par0 = c(3, 1, 1), parallel = TRUE, ncores = 6,
  args1 = args1, args2 = args2, control = list(trace = 0)
)
stderr(fit1)
}

Run the code above in your browser using DataLab