Generalised least squares estimation of regularly varying space-time processes based on flexible observation schemes

TL;DR

This paper introduces a two-step estimation method for extremograms of regularly varying space-time processes, combining asymptotic theory with practical subsampling techniques, and demonstrates its effectiveness on Brown-Resnick models.

Contribution

It develops a novel two-step estimation approach for extremograms in space-time processes, including asymptotic properties and confidence interval procedures.

Findings

Method performs well for moderate sample sizes

Provides asymptotic normality results for estimators

Effective application to Brown-Resnick processes

Abstract

Regularly varying stochastic processes model extreme dependence between process values at different locations and/or time points. For such processes we propose a two-step parameter estimation of the extremogram, when some part of the domain of interest is fixed and another increasing. We provide conditions for consistency and asymptotic normality of the empirical extremogram centred by a pre-asymptotic version for such observation schemes. For max-stable processes with Fr{\'e}chet margins we provide conditions, such that the empirical extremogram (or a bias-corrected version) centred by its true version is asymptotically normal. In a second step, for a parametric extremogram model, we fit the parameters by generalised least squares estimation and prove consistency and asymptotic normality of the estimates. We propose subsampling procedures to obtain asymptotically correct confidence intervals. Finally, we apply our results to a variety of Brown-Resnick processes. A simulation study shows that the procedure works well also for moderate sample sizes.

Tables

Table 1. Table 1: True parameter values (first column) and mean, MAE, RMSE, and REL of the estimates of the parameters of model (i).

	TRUE	MEAN	MAE	RMSE	REL
${\hat{C}}_1$	0.8	0.7856	0.1353	0.1763	0.2204
${\hat{C}}_2$	0.4	0.3987	0.0785	0.0995	0.2486
${\hat{\alpha }}_1$	1.5	1.4830	0.0897	0.1131	0.0754
${\hat{\alpha }}_2$	1	0.9916	0.0625	0.0820	0.0820

Equations

$$\displaystyle n^{d}\mathbb{P}\Big{(}\frac{X_{\mathcal{I}}}{a_{n}}\in\cdot\Big{)}\stackrel{{\scriptstyle v}}{{\rightarrow}}\mu_{\mathcal{I}}(\cdot),\quad n\rightarrow\infty,$$

$${\mu }_{\mathcal{I}}(xC)=x^{-\beta }{\mu }_{\mathcal{I}}(C),x>0,$$

$${\mu }_{\{\mathbf{s}\}}(\cdot )={\mu }_{\{0\}}(\cdot )=:\mu (\cdot ),$$

$${\rho }_{AB}(\mathbf{h})=\underset{n\to \infty }{\lim }\frac{\mathbb{P}(X(0)/a_n\in A,X(\mathbf{h})/a_n\in B)}{\mathbb{P}(X(0)/a_n\in A)},\mathbf{h}\in {\mathbb{R}}^d.$$

$${\mathcal{D}}_n=\mathcal{F}\times {\mathcal{I}}_n,$$

$${\mathcal{D}}_n=\mathcal{F}\times \{1,\dots ,n\}$$

$$\mathcal{B}(0,\gamma )$$

$$\mathbf{Y}(\mathbf{s}):=X_{\mathcal{B}(\mathbf{s},\gamma )};$$

$${\mu }_{\mathcal{B}(0,\gamma )}(C)$$

$${\tau }_{\mathcal{B}(0,\gamma )\times \mathcal{B}(\mathbf{\ell },\gamma )}(C\times D)$$

$$\{\mathbf{Y}(\mathbf{s})\in D_i\}=\{X(\mathbf{s})\in A,X(\mathbf{s}+{\mathbf{h}}^{(i)})\in B\}$$

$${\mu }_{\mathcal{B}(0,\gamma )}(D_{p+1})=\underset{n\to \infty }{\lim }{n}^d\mathbb{P}(\mathbf{Y}(0)/a_n\in D_{p+1})=\underset{n\to \infty }{\lim }{n}^d\mathbb{P}(X(0)/a_n\in A)=\mu (A).$$

$$\mathcal{Z}(\mathbf{h}):=\{\mathbf{z}\in \mathcal{Z}:\mathbf{z}+\mathbf{h}\in \mathcal{Z}\},$$

$${\hat{\rho }}_{AB,m_n}(\mathbf{h}):=\frac{\frac{1}{|{\mathcal{D}}_n(\mathbf{h})|}\sum _{\mathbf{s}\in {\mathcal{D}}_n(\mathbf{h})}\text{mathbbmss}{1}_{\{X(\mathbf{s})/a_m\in A,X(\mathbf{s}+\mathbf{h})/a_m\in B\}}}{\frac{1}{|{\mathcal{D}}_n|}\sum _{\mathbf{s}\in {\mathcal{D}}_n}\text{mathbbmss}{1}_{\{X(\mathbf{s})/a_m\in A\}}},\mathbf{h}\in \mathcal{H}.$$

$${\rho }_{AB,m_n}(\mathbf{h})=\frac{\mathbb{P}\left(X(0)/a_m\in A,X(\mathbf{h})/a_m\in B\right)}{\mathbb{P}(X(0)/a_m\in A)},\mathbf{h}\in {\mathbb{R}}^d.$$

$$\mathbb{P}\Big{(}\frac{X(\boldsymbol{0})}{a_{m}}\in A,\frac{X(\boldsymbol{h})}{a_{m}}\in B\Big{)}=\mathbb{P}\Big{(}\frac{\boldsymbol{Y}(\boldsymbol{0})}{a_{m}}\in C(\boldsymbol{h})\Big{)}.$$

$${\hat{\mu }}_{\mathcal{B}(0,\gamma ),m_n}(C(\mathbf{h})):=\frac{m_n^d}{n^w}\sum _{\mathbf{i}\in {\mathcal{I}}_n}\frac{1}{|\mathcal{F}({\mathbf{h}}_{\mathcal{F}})|}\sum _{\mathbf{f}\in \mathcal{F}({\mathbf{h}}_{\mathcal{F}})}\text{mathbbmss}{1}_{\{\frac{\mathbf{Y}(\mathbf{f},\mathbf{i})}{a_m}\in C(\mathbf{h})\}}.$$

$${\mathbb{C}ov}\Big{[}\mathbbmss{1}_{\{\frac{\boldsymbol{Y}(\boldsymbol{f},\boldsymbol{i})}{a_{m}}\in D_{i}\}},\mathbbmss{1}_{\{\frac{\boldsymbol{Y}(\boldsymbol{f}^{\prime},\boldsymbol{i}^{\prime})}{a_{m}}\in D_{i}\}}\Big{]}={\mathbb{C}ov}\Big{[}\mathbbmss{1}_{\{\frac{\boldsymbol{Y}(\boldsymbol{0})}{a_{m}}\in D_{i}\}},\mathbbmss{1}_{\{\frac{\boldsymbol{Y}(\boldsymbol{\ell}_{\mathcal{F}},\boldsymbol{\ell}_{\mathcal{I}})}{a_{m}}\in D_{i}\}}\Big{]}$$

$$L_{\mathcal{F}}^{(i,j)}:=\{\mathbf{f}-{\mathbf{f}}^{\prime }:\mathbf{f}\in \mathcal{F}({\mathbf{h}}_{\mathcal{F}}^{(i)}),{\mathbf{f}}^{\prime }\in \mathcal{F}({\mathbf{h}}_{\mathcal{F}}^{(j)})\}\text{mbox}and{L}_n:=\{\mathbf{i}-{\mathbf{i}}^{\prime }:\mathbf{i},{\mathbf{i}}^{\prime }\in {\mathcal{I}}_n\}.$$

$${\text{N}}_{\mathcal{F}}^{(i,j)}({\mathbf{\ell }}_{\mathcal{F}}):=\sum _{\mathbf{f}\in \mathcal{F}({\mathbf{h}}_{\mathcal{F}}^{(i)})}\sum _{{\mathbf{f}}^{\prime }\in \mathcal{F}({\mathbf{h}}_{\mathcal{F}}^{(j)})}\text{mathbbmss}{1}_{\{\mathbf{f}-{\mathbf{f}}^{\prime }={\mathbf{\ell }}_{\mathcal{F}}\}}$$

$$\displaystyle{\mathbb{V}ar}\big{[}\widehat{\mu}_{\mathcal{B}(\boldsymbol{0},{\gamma}),m_{n}}(D_{i})\big{]}=\frac{m_{n}^{2d}}{n^{2w}|\mathcal{F}(\boldsymbol{h}_{\mathcal{F}}^{(i)})|^{2}}{\mathbb{V}ar}\big{[}\sum_{\boldsymbol{f}\in\mathcal{F}(\boldsymbol{h}_{\mathcal{F}}^{(i)})}\sum_{\boldsymbol{i}\in\mathcal{I}_{n}}\mathbbmss{1}_{\{\frac{\boldsymbol{Y}(\boldsymbol{f},\boldsymbol{i})}{a_{m}}\in D_{i}\}}\big{]}$$

$$\displaystyle=\frac{m_{n}^{2d}}{n^{2w}|\mathcal{F}(\boldsymbol{h}_{\mathcal{F}}^{(i)})|^{2}}\Big{(}|\mathcal{F}(\boldsymbol{h}_{\mathcal{F}}^{(i)})|n^{w}{\mathbb{V}ar}\big{[}\mathbbmss{1}_{\{\frac{\boldsymbol{Y}(\boldsymbol{0})}{a_{m}}\in D_{i}\}}\big{]}$$

$$\displaystyle\quad\quad\quad\quad+\sum_{\boldsymbol{f},\boldsymbol{f}^{\prime}\in\mathcal{F}(\boldsymbol{h}_{\mathcal{F}}^{(i)})}\sum_{\boldsymbol{i},\boldsymbol{i}^{\prime}\in\mathcal{I}_{n}\atop(\boldsymbol{f},\boldsymbol{i})\neq(\boldsymbol{f}^{\prime},\boldsymbol{i}^{\prime})}{\mathbb{C}ov}\big{[}\mathbbmss{1}_{\{\frac{\boldsymbol{Y}(\boldsymbol{f},\boldsymbol{i})}{a_{m}}\in D_{i}\}},\mathbbmss{1}_{\{\frac{\boldsymbol{Y}(\boldsymbol{f}^{\prime},\boldsymbol{i}^{\prime})}{a_{m}}\in D_{i}\}}\big{]}\Big{)}$$

$$\displaystyle\sim\frac{m_{n}^{d}}{n^{w}}\frac{1}{|\mathcal{F}(\boldsymbol{h}_{\mathcal{F}}^{(i)})|}\Big{(}\mu_{\mathcal{B}(\boldsymbol{0},{\gamma})}(D_{i})$$

$$\displaystyle\quad\quad\quad\quad+\sum_{\boldsymbol{\ell}_{\mathcal{I}}\in\mathbb{Z}^{w}}\frac{1}{|\mathcal{F}(\boldsymbol{h}_{\mathcal{F}}^{(i)})|}\sum_{\boldsymbol{\ell}_{\mathcal{F}}\in L_{\mathcal{F}}^{(i,i)}\atop(\boldsymbol{\ell}_{\mathcal{F}},\boldsymbol{\ell}_{\mathcal{I}})\neq\boldsymbol{0}}\textnormal{N}_{\mathcal{F}}^{(i,i)}(\boldsymbol{\ell}_{\mathcal{F}})\,\tau_{\mathcal{B}(\boldsymbol{0},{\gamma})\times\mathcal{B}((\boldsymbol{\ell}_{\mathcal{F}},\boldsymbol{\ell}_{\mathcal{I}}),{\gamma})}(D_{i}\times D_{i})\Big{)}$$

$$=:\frac{m_n^d}{n^w}{\sigma }_{\mathcal{B}(0,\gamma )}^2(D_i),n\to \infty .$$

$$\displaystyle{\mathbb{V}ar}\big{[}\widehat{\mu}_{\mathcal{B}(\boldsymbol{0},\gamma),m_{n}}(D_{i})\big{]}$$

$$\sim$$

$$=$$

$$\sqrt{\frac{n^{w}}{m_{n}^{d}}}\Big{[}\widehat{\rho}_{AB,m_{n}}(\boldsymbol{h}^{(i)})-\rho_{AB,m_{n}}(\boldsymbol{h}^{(i)})\Big{]}_{i=1,\ldots,p}\stackrel{{\scriptstyle\mathcal{D}}}{{\rightarrow}}\mathcal{N}(\boldsymbol{0},\Pi),\quad n\to\infty,$$

$${\Sigma }_{ij}=$$

$$\displaystyle+\sum_{\boldsymbol{\ell}_{\mathcal{I}}\in\mathbb{Z}^{w}}\sum_{\boldsymbol{\ell}_{\mathcal{F}}\in L_{\mathcal{F}}^{(i,j)}\atop(\boldsymbol{\ell}_{\mathcal{F}},\boldsymbol{\ell}_{\mathcal{I}})\neq\boldsymbol{0}}\textnormal{N}_{\mathcal{F}}^{(i,j)}(\boldsymbol{\ell}_{\mathcal{F}})\,\tau_{\mathcal{B}(\boldsymbol{0},{\gamma})\times\mathcal{B}((\boldsymbol{\ell}_{\mathcal{F}},\boldsymbol{\ell}_{\mathcal{I}}),{\gamma})}(D_{i}\times D_{j})\Big{)},\quad 1\leq i,j\leq p+1.$$

$$F_1=\text{diag}(\mu (A))\in {\mathbb{R}}^{p\times p},F_2=(-{\mu }_{\mathcal{B}(0,\gamma )}(D_1),\dots ,-{\mu }_{\mathcal{B}(0,\gamma )}(D_p)){}^{{}^{⊺}}.$$

$${\hat{\rho }}_{AB,m_n}({\mathbf{h}}^{(i)})\stackrel{P}{\to }{\rho }_{AB}({\mathbf{h}}^{(i)}),i=1,\dots ,p,$$

$${\hat{\rho }}_{AB,m_n}({\mathbf{h}}^{(i)})\sim \frac{{\hat{\mu }}_{\mathcal{B}(0,\gamma ),m_n}(D_i)}{{\hat{\mu }}_{\mathcal{B}(0,\gamma ),m_n}(D_{p+1})}\stackrel{P}{\to }\frac{{\mu }_{\mathcal{B}(0,\gamma )}(D_i)}{{\mu }_{\mathcal{B}(0,\gamma )}(D_{p+1})}={\rho }_{AB}({\mathbf{h}}^{(i)}),$$

Full text

AMS 2010 Subject Classifications Primary: 60F05 $\cdot$ 60G70 $\cdot$ 62F12 $\cdot$ 62G32 $\cdot$ Secondary: 37A25 $\cdot$ 62M30 $\cdot$ 62P12

1 Introduction

Max-stable processes and regularly varying processes have in recent years attracted attention as time series models, spatial processes and space-time processes. Regularly varying processes have been investigated in Hult and Lindskog (2005, 2006) and basic results for max-stable processes can be found in de Haan and Ferreira (2006). Such processes provide a useful framework for modelling and estimation of extremal events in their different settings.

Among the various regularly varying models considered in the literature, max-stable Brown-Resnick processes play a prominent role allowing for flexible fractional variogram models as often observed in environmental data. They have been introduced for time series in Brown and Resnick (1977), for spatial processes in Kabluchko et al. (2009), and in a space-time setting in Davis et al. (2013a).

For max-stable processes with parametrised dependence structure, various estimation procedures have been proposed for extremal data. Composite likelihood methods have been described in Padoan et al. (2009) and Huser and Davison (2014). Threshold-based likelihood methods have been proposed in Wadsworth and Tawn (2014) and Engelke et al. (2015). For the max-stable Brown-Resnick process asymptotic results of composite likelihood estimators have been derived in Davis et al. (2013b), Huser and Davison (2013), and Buhl and Klüppelberg (2016). In some special cases full likelihood estimation is feasible, which opens the door for frequentist or Bayesian approaches; see for example Dombry et al. (2016b); Thibaud et al. (2016).

Parameter estimation based on likelihood methods can be laborious and time consuming, and also the choice of good initial values for the optimization routine is essential. As a consequence, a semiparametric estimation procedure can be an alternative or a prerequisite for a subsequent likelihood method. Such an estimation method has been suggested and analysed for space-time processes with additively separable dependence function in Steinkohl (2013) and Buhl et al. (2017) based on the extremogram, which is a natural extremal analogue of the correlation function for stationary processes. The extremogram was introduced for time series in Davis and Mikosch (2009) and Fasen et al. (2010), and extended to a spatial and space-time settings in Steinkohl (2013) and Cho et al. (2016). Semiparametric estimation requires a parametric extremogram model. The parameter estimation is then based on the empirical extremogram, and a subsequent least squares estimation of the parameters.

The processes considered in Steinkohl (2013), Cho et al. (2016), Buhl and Klüppelberg (2017), and Buhl et al. (2017) are isotropic in space; cf. model (I) in Section 5.3 below. The central goal of this paper is to generalise the semiparametric method developed in Buhl et al. (2017) in various aspects. We list the most important extensions:

– In Buhl et al. (2017) ordinary least squares estimation was performed separately for the spatial and the temporal dependence parameters. This was possible, since we assumed an additively separable dependence model, linear in its parameters after a suitable transformation. In the present paper we allow for a much larger class of dependence models provided they satisfy some weak regularity conditions. In particular, we allow for non-linear structures in the dependence models, and we estimate a space-time dependence model, which is not necessarily separable.

– To fit these general models to data, we develop a generalised (weighted) least squares estimation method, which estimates all dependence parameters in one go.

– We again focus on extremogram estimation, but extend the observation scheme as described below. In the context of spatial or space-time extremogram estimation based on gridded data, the observation scheme used so far in Steinkohl (2013), Cho et al. (2016), Buhl and Klüppelberg (2017), and Buhl et al. (2017) has been a regular grid in space, possibly observed at equidistant time points and assumed to expand to infinity in all spatial dimensions as well as in time. We extend this observation scheme to a more realistic setting: in practice one often observes data on a $d$-dimensional area ($d\in\mathbb{N}$), which is small with respect to some of its dimensions (for instance, the spatial dimensions) and large with respect to others (for instance, the temporal dimension). Hence, with regard to such cases, it is appropriate to assume the observed data to expand to infinity in some dimensions, but remain fixed in some others. Such observation schemes require to split up every point and every lag in its components corresponding to the fixed and increasing domains.

– For such general observation schemes we have to extend the asymptotic theory developed in Buhl and Klüppelberg (2017) considerably. The empirical extremogram estimator used in the first step of the semiparametric estimation procedure needs to be extended and asymptotic results need to be verified. For an arbitrary parametric extremogram model we then derive asymptotic results of its generalised least squares estimators, which differ considerably from those obtained when the grid increases in all dimensions.

Our paper is organised as follows. In Section 2 we introduce the theoretical framework of strictly stationary regularly varying processes. We define the extremogram, the observation scheme with its fixed and increasing dimensions as well as assumptions and asymptotic second order properties following from regular variation. Section 3 presents the empirical and the pre-asymptotic extremogram. Here we prove a CLT for the empirical extremogram centred by the pre-asymptotic version. We also specify the asymptotic covariance matrix. We prove a CLT for the empirical extremogram centred by the true extremogram under more restrictive assumptions. To formally state the asymptotic properties of the empirical extremogram, we need to quantify the dependence in a stochastic process, taking into account the different types of observation areas. For processes with Fréchet margins we prove asymptotic normality of the empirical extremogram centred by the true one. In case the required conditions are not satisfied, we provide assumptions under which a CLT for a bias corrected version of the empirical extremogram can be obtained. Section 4 is dedicated to the parameter estimation by a generalised least squares method. Under appropriate regularity conditions we prove consistency and asymptotic normality, where the rate of convergence depends on the observation scheme. We also present the covariance matrix in a semi-explicit form. In Section 5 we show our method at work for Brown-Resnick space-time processes. We state conditions for Brown-Resnick processes that imply the mixing conditions from Section 3 and are hence sufficient to obtain the corresponding CLTs for the empirical extremogram. These conditions depend highly on the model for the associated variogram. Finally, in Section 5.3 we apply these results to three different dependence models of the Brown-Resnick process, and prove the mixing conditions, which guarantee the asymptotic normality of the empirical extremogram, as well as the regularity conditions of the generalised least squares estimates. In Section 6 we examine the finite sample properties of the GLSEs in a simulation study, fitting the parametric models described in Section 5.3 to simulated Brown-Resnick processes. We apply subsampling methods to obtain asymptotically valid confidence bounds of the parameters. We examine how the sample size affects the estimates and compare with the theoretical results obtained in previous sections. Many proofs are rather technical and postponed to an Appendix.

2 Model description and the observation scheme

We consider the same theoretical framework as in Buhl and Klüppelberg (2017) and Buhl et al. (2017) of a strictly stationary regularly varying process $\{X(\boldsymbol{s}):\boldsymbol{s}\in\mathbb{R}^{d}\}$ for $d\in\mathbb{N}$, defined on a probability space $(\Omega,\mathcal{G},\mathbb{P})$. This implies that there exists some normalizing sequence $0<a_{n}\rightarrow\infty$ such that $\mathbb{P}(|X(\boldsymbol{0})|>a_{n})\sim n^{-d}$ as $n\rightarrow\infty$ and that for every finite set $\mathcal{I}\subset\mathbb{R}^{d}$ with cardinality $|\mathcal{I}|<\infty$,

[TABLE]

for some non-null Radon measure $\mu_{\mathcal{I}}$ on the Borel sets in $\overline{\mathbb{R}}^{|\mathcal{I}|}\backslash\{\boldsymbol{0}\}$, where $\overline{\mathbb{R}}=\mathbb{R}\cup\{-\infty,\infty\}$ and $X_{\mathcal{I}}$ denotes the vector $(X(\boldsymbol{s}):\boldsymbol{s}\in\mathcal{I}).$ The limit measure is homogeneous:

[TABLE]

for every Borel set $C\subset\overline{\mathbb{R}}^{|\mathcal{I}|}\backslash\{\boldsymbol{0}\}$. The notation $\stackrel{{\scriptstyle v}}{{\rightarrow}}$ stands for vague convergence, and $\beta>0$ is called the index of regular variation. Furthermore, $f(n)\sim g(n)$ as $n\to\infty$ means that $\lim_{n\to\infty}f(n)/g(n)=1$. If $\mathcal{I}$ is a singleton; i.e., $\mathcal{I}=\{\boldsymbol{s}\}$ for some $\boldsymbol{s}\in\mathbb{R}^{d}$, we set

[TABLE]

which is justified by stationarity. For more details see Buhl and Klüppelberg (2017). For background on regular variation for stochastic processes and vectors see Hult and Lindskog (2005, 2006) and Resnick (1986, 2007).

The extremogram for values in $\mathbb{R}^{d}$ is defined as follows.

Definition 1 (Extremogram)

Let $\{X(\boldsymbol{s}):\boldsymbol{s}\in\mathbb{R}^{d}\}$ be a strictly stationary regularly varying process and $a_{n}\to\infty$ a sequence satisfying (2.1). For $\mu$ as in (2.2) and two $\mu$-continuous Borel sets $A$ and $B$ in $\overline{\mathbb{R}}\backslash\{0\}$ (i.e., $\mu(\partial A)=\mu(\partial B)=0$) such that $\mu(A)>0$, the extremogram is defined as

[TABLE]

For $A=B=(1,\infty)$, the extremogram $\rho_{AB}(\boldsymbol{h})$ is the tail dependence coefficient between $X(\boldsymbol{0})$ and $X(\boldsymbol{h})$ (cf. Beirlant et al. (2004), Section 9.5.1).

For the data we allow for realistic observation schemes described in the following.

Assumption 1

The data are given in an observation area $\mathcal{D}_{n}\subset\mathbb{Z}^{d}$ that can (possibly after reordering) be decomposed into

[TABLE]

where for $q,w\in\mathbb{N}$ satisfying $w+q=d$:

1. (1)

$\mathcal{F}\subset\mathbb{Z}^{q}$* is a fixed domain independent of $n$, and* 2. (2)

$\mathcal{I}_{n}=\{1,\ldots,n\}^{w}$* is an increasing sequence of regular grids.*

This setting is similar to that used in Li et al. (2008), where asymptotic properties of space-time covariance estimators are derived. The natural extension of the regular grid $\mathcal{I}_{n}$ to grids with different side lengths only increases notational complexity, which we avoid here. Our focus is on observations schemes, which are partially fixed and partially tend to infinity.

Example 1

In the special case where the observation area is given by

[TABLE]

for $\mathcal{F}\subset\mathbb{R}^{d-1}$, we interpret the observations as generated by a space-time process $\{X(\boldsymbol{s},t):\boldsymbol{s}\in\mathbb{R}^{d-1},t\in[0,\infty)\}$ on a fixed spatial and an increasing temporal domain.

We shall need some definitions and assumptions, which we summarize as follows.

Assumption 2

$(1)$*   For some fixed ${\gamma}>0$ and $\boldsymbol{0},\boldsymbol{\ell}\in\mathbb{R}^{d}$ we define the balls*

[TABLE]

$(2)$*   The estimation of the extremogram is based on a set $\mathcal{H}=\{\boldsymbol{h}^{(1)},\ldots,\boldsymbol{h}^{(p)}\}\subset\mathcal{B}(\boldsymbol{0},{\gamma})$ of observed lag vectors.

$(3)$   We decompose points $\boldsymbol{s}\in\mathbb{R}^{d}$ with respect to the fixed and increasing domains into $\boldsymbol{s}=(\boldsymbol{f},\boldsymbol{i})\in\mathbb{R}^{q}\times\mathbb{R}^{w}$.

$(4)$   Similarly, we decompose lag vectors $\boldsymbol{h}=\boldsymbol{s}-\boldsymbol{s}^{\prime}$ or $\boldsymbol{\ell}=\boldsymbol{s}-\boldsymbol{s}^{\prime}$ for some $\boldsymbol{s},\boldsymbol{s}^{\prime}\in\mathbb{R}^{d}$ into $\boldsymbol{h}=(\boldsymbol{h}_{\mathcal{F}},\boldsymbol{h}_{\mathcal{I}})$ or $\boldsymbol{\ell}=(\boldsymbol{\ell}_{\mathcal{F}},\boldsymbol{\ell}_{\mathcal{I}})$ in $\mathbb{R}^{q}\times\mathbb{R}^{w}$. The letter $\boldsymbol{h}$ is used throughout as argument of the extremogram or its estimators.

$(5)$   We define the vectorised process $\{\boldsymbol{Y}(\boldsymbol{s}):\boldsymbol{s}\in\mathbb{R}^{d}\}$ by*

[TABLE]

*i.e., $\boldsymbol{Y}(\boldsymbol{s})$ is the vector of values of $X$ with indices in the ball $\mathcal{B}(\boldsymbol{s},{\gamma})$.

$(6)$   We shall also need the following relations, already stated in (3.3) and (3.4) of Buhl and Klüppelberg (2017). For $a_{n}\to\infty$ as in (2.1), the following limits exist by regular variation of $\{X(\boldsymbol{s}):\boldsymbol{s}\in\mathbb{R}^{d}\}$. For $\boldsymbol{\ell}\in\mathbb{R}^{d}$ and ${\gamma}>0$,*

[TABLE]

*for a $\mu_{\mathcal{B}(\boldsymbol{0},{\gamma})}$-continuous Borel set $C$ in $\overline{\mathbb{R}}^{|\mathcal{B}(\boldsymbol{0},{\gamma})|}\backslash\{\boldsymbol{0}\}$ and a $\tau_{{\mathcal{B}(\boldsymbol{0},\gamma)\times\mathcal{B}(\boldsymbol{\ell},\gamma)}}$-continuous Borel set $C\times D$ in the product space.

$(7)$   For arbitrary but fixed $\mu$-continuous Borel sets $A$ and $B$ in $\overline{\mathbb{R}}\backslash\{0\}$ such that $\mu(A)>0$, we define sets $D_{1},\ldots,D_{p},D_{p+1}$ by the identity*

[TABLE]

for $i=1,\ldots,p$, and $\{\boldsymbol{Y}(\boldsymbol{s})\in D_{p+1}\}=\{X(\boldsymbol{s})\in A\}$. Note in particular that, by the relation between $\{\boldsymbol{Y}(\boldsymbol{s}):\boldsymbol{s}\in\mathbb{R}^{d}\}$ and $\{X(\boldsymbol{s}):\boldsymbol{s}\in\mathbb{R}^{d}\}$ and regular variation,

[TABLE]

$\Box$

3 Limit theory for the empirical extremogram

We derive asymptotic properties of the empirical extremogram by formulating appropriate mixing conditions, generalising the results obtained in Buhl and Klüppelberg (2017) to the more realistic setting of this paper. The proofs are based on spatial mixing conditions, which have to be adapted to the decomposition into a fixed and an increasing observation domain. In principle, our proofs rely on general results of Ibragimov and Linnik (1971) and Bolthausen (1982).

The main theorem of this section states asymptotic normality of the empirical extremogram sampled at lag vectors $\boldsymbol{h}\in{\mathcal{H}}$ and centred by its pre-asymptotic counterpart. The empirical and the pre-asymptotic extremograms are defined in Eq. (3.2) and (3.3).

For the definition of the empirical extremogram we need the following notation: for $k\in\mathbb{N}$, an arbitrary set $\mathcal{Z}\subset\mathbb{Z}^{k}$ and a fixed vector $\boldsymbol{h}\in\mathbb{Z}^{k}$, define the sets

[TABLE]

which is the set of vectors $\boldsymbol{z}\in\mathcal{Z}$ such that with $\boldsymbol{z}$ also the lagged vector $\boldsymbol{z}+\boldsymbol{h}$ belongs to $\mathcal{Z}$.

Definition 2

Let $\{X(\boldsymbol{s}):\boldsymbol{s}\in\mathbb{R}^{d}\}$ be a strictly stationary regularly varying process, which is observed on $\mathcal{D}_{n}=\mathcal{F}\times\mathcal{I}_{n}$ as in (2.4). Let $A$ and $B$ be $\mu$-continuous Borel sets in $\overline{\mathbb{R}}\backslash\{0\}$ such that $\mu(A)>0$. For a sequence $m=m_{n}\rightarrow\infty$ and $m_{n}=o(n)$ as $n\to\infty$ define the following quantities:

1. (1)

The empirical extremogram

[TABLE]

For a fixed data set the value $a_{m}=a_{m_{n}}$ has to be specified as a large empirical quantile. 2. (2)

The pre-asymptotic extremogram

[TABLE]

Key of the proofs of consistency and asymptotic normality of the empirical extremogram below is the fact that $\widehat{\rho}_{AB,m_{n}}(\boldsymbol{h})$ is the empirical version of the pre-asymptotic extremogram $\rho_{AB,m_{n}}(\boldsymbol{h})$. This can for different $\boldsymbol{h}\in\mathcal{B}(\boldsymbol{0},{\gamma})$ in turn be viewed as a ratio of pre-asymptotic versions of $\mu_{\mathcal{B}(\boldsymbol{0},{\gamma})}(C(\boldsymbol{h}))$ (cf. Eq. (2.5)). The sets $C(\boldsymbol{h})$ are implicitly defined by $\{\boldsymbol{Y}(\boldsymbol{s})\in C(\boldsymbol{h})\}=\{X(\boldsymbol{s})\in A,X(\boldsymbol{s}+\boldsymbol{h})\in B\}$ for $\boldsymbol{s}\in\mathbb{R}^{d}$. Then in particular, for $\boldsymbol{h}\in\mathcal{B}(\boldsymbol{0},{\gamma})$,

[TABLE]

Note that, by (2.7), if $\boldsymbol{h}=\boldsymbol{h}^{(i)}\in\mathcal{H}$, then $C(\boldsymbol{h})=D_{i}$, and if $\boldsymbol{h}=\boldsymbol{0}$ and $A=B$ then $C(\boldsymbol{h})=D_{p+1}$.

In view of (2.5), $\mu_{\mathcal{B}(\boldsymbol{0},{\gamma})}(C(\boldsymbol{h}))$ can be estimated by an empirical mean, where the estimator has to cope with Assumption 1 of an observation area with fixed and increasing domain.

Definition 3

Assume the situation of Definition 2. Based on observations on $\mathcal{D}_{n}=\mathcal{F}\times\mathcal{I}_{n}$ as in (2.4) decompose the observations $\boldsymbol{s}=(\boldsymbol{f},\boldsymbol{i})\in\mathcal{F}\times\mathcal{I}_{n}$ and the lags $\boldsymbol{h}=(\boldsymbol{h}_{\mathcal{F}},\boldsymbol{h}_{{\mathcal{I}}})\in{\mathcal{H}}$ as in Assumption 2(3) and (4). For $\boldsymbol{h}_{\mathcal{F}}\in{\mathcal{H}}$ define $\mathcal{F}(\boldsymbol{h}_{\mathcal{F}})$ as in (3.1). Then an empirical version of $\mu_{\mathcal{B}(\boldsymbol{0},{\gamma})}(C(\boldsymbol{h}))$ is for $\boldsymbol{h}\in{\mathcal{H}}$ given by

[TABLE]

$\Box$

Observe that for fixed $\boldsymbol{h}_{\mathcal{F}}\in\mathbb{Z}^{q}$ and observations on $\mathcal{D}_{n}=\mathcal{F}\times\mathcal{I}_{n}$ there will be points $\boldsymbol{s}=(\boldsymbol{f},\boldsymbol{i})\in\mathcal{F}(\boldsymbol{h}_{\mathcal{F}})\times\mathcal{I}_{n}$ with $\boldsymbol{i}$ near the boundary of $\mathcal{I}_{n}$, such that not all components of the vector $\boldsymbol{Y}(\boldsymbol{s})=\boldsymbol{Y}(\boldsymbol{f},\boldsymbol{i})$ are observed. However, since we investigate asymptotic properties of $\mathcal{I}_{n}$ whose boundary points are negligible, we can ignore such technical details. As will be seen in the proofs below, for every $\boldsymbol{h}\in\mathcal{H}$, the empirical extremogram $\widehat{\rho}_{AB,m_{n}}(\boldsymbol{h})$ is asymptotically equivalent to the ratio of estimates $\widehat{\mu}_{\mathcal{B}(\boldsymbol{0},{\gamma}),m_{n}}(C(\boldsymbol{h}))/\widehat{\mu}_{\mathcal{B}(\boldsymbol{0},{\gamma}),m_{n}}(D_{p+1})$.

Limit results for the empirical extremogram (3.2) involve the calculation of mean and variance of $\widehat{\mu}_{\mathcal{B}(\boldsymbol{0},{\gamma}),m_{n}}(C(\boldsymbol{h}^{(i)}))=\widehat{\mu}_{\mathcal{B}(\boldsymbol{0},{\gamma}),m_{n}}(D_{i})$ for $\boldsymbol{h}^{(i)}\in{\mathcal{H}}$. Strict stationarity and Assumption 2(6) yields immediately by a law of large numbers that

$\mathbb{E}[\widehat{\mu}_{\mathcal{B}(\boldsymbol{0},{\gamma}),m_{n}}(D_{i})]\to{\mu}_{\mathcal{B}(\boldsymbol{0},{\gamma})}(D_{i})$ as $n\to\infty$. Calculation of the variance involves the covariance structure and we decompose as in Assumption 2(4) $\boldsymbol{h}^{(i)}$ into $\boldsymbol{h}^{(i)}=(\boldsymbol{h}^{(i)}_{\mathcal{F}},\boldsymbol{h}^{(i)}_{\mathcal{I}})$ $\in\mathbb{R}^{q}\times\mathbb{R}^{w}$. We have to calculate for $\boldsymbol{f},\boldsymbol{f}^{\prime}\in\mathcal{F}(\boldsymbol{h}_{\mathcal{F}}^{(i)})$ and $\boldsymbol{i},\boldsymbol{i}^{\prime}\in\mathcal{I}_{n}$,

[TABLE]

with $\boldsymbol{\ell}_{\mathcal{F}}=\boldsymbol{f}-\boldsymbol{f}^{\prime}$ and $\boldsymbol{\ell}_{\mathcal{I}}=\boldsymbol{i}-\boldsymbol{i}^{\prime}$, where the equality holds by stationarity. The lag vectors $\boldsymbol{\ell}_{\mathcal{F}}$ and $\boldsymbol{\ell}_{\mathcal{I}}$ are contained in $L_{\mathcal{F}}^{(i,i)}$ and $L_{n}$, respectively, where for $i,j\in\{1,\ldots,p\}$,

[TABLE]

The number of appearances of the lag $\boldsymbol{\ell}_{\mathcal{F}}$ we denote for $i,j\in\{1,\ldots,p\}$ by

[TABLE]

Observe that a lag $(\boldsymbol{\ell}_{\mathcal{F}},\boldsymbol{\ell}_{\mathcal{I}})$ with $\boldsymbol{\ell}_{\mathcal{I}}=(\ell_{\mathcal{I}}^{(1)},\ldots,\ell_{\mathcal{I}}^{(w)})$ appears in $L_{\mathcal{F}}^{(i,i)}\times L_{n}$ exactly $\textnormal{N}_{\mathcal{F}}^{(i,i)}(\boldsymbol{\ell}_{\mathcal{F}})\prod_{j=1}^{w}(n-|\ell_{\mathcal{I}}^{(j)}|)$ times. We show in Lemma A.1 that

[TABLE]

Remark 1

For comparison we recall the expression in the corresponding Lemma 5.1 of Buhl and Klüppelberg (2017), where $\mathcal{F}$ is not fixed, but part of the increasing regular grid. Then $|\mathcal{F}(\boldsymbol{h}_{\mathcal{F}}^{(i)})|\sim\textnormal{N}_{\mathcal{F}}^{(i,i)}(\boldsymbol{\ell}_{\mathcal{F}})\sim n^{q}$ as $n\to\infty$, such that (3) can be approximated as follows:

[TABLE]

Thus, a difference from the setting of a partly fixed observation area $\mathcal{F}\subset{\mathcal{D}}_{n}$ is that the fixed observation terms do not disappear asymptotically, but remain as constants in the limit expression.

3.1 The extremogram for regularly varying processes

For proving asymptotic normality of the empirical extremogram we have to require appropriate mixing conditions and make use of a large/small block argument as in Buhl and Klüppelberg (2017). For simplicity we assume that $n^{w}/m_{n}^{d}$ is an integer and subdivide $\mathcal{D}_{n}$ into $n^{w}/m_{n}^{d}$ non-overlapping $d$-dimensional large blocks $\mathcal{F}\times\mathcal{B}_{i}$ for $i=1,\ldots,n^{w}/m_{n}^{d}$, where the $\mathcal{B}_{i}$ are $w$-dimensional cubes with side lengths $m_{n}^{d/w}$. From those large blocks we then cut off smaller blocks, which consist of the first $r_{n}$ elements in each of the $w$ increasing dimensions. The large blocks are then separated (by these small blocks) with at least the distance $r_{n}$ in all $w$ increasing dimensions and shown to be asymptotically independent. Such large/small block arguments are common in verifying properties of estimators in extreme value theory, in particular in a time series context, cf. for example Davis and Mikosch (2009), Section 6. For a visualization in the $2$-dimensional case $d=2$ with $w=1$ increasing dimension, see Figure 1.

In order to formulate the CLT below, in particular, the asymptotic covariance matrix, we need to compute ${\mathbb{C}ov}[\widehat{\mu}_{\mathcal{B}(\boldsymbol{0},{\gamma}),m_{n}}(D_{i}),\widehat{\mu}_{\mathcal{B}(\boldsymbol{0},{\gamma}),m_{n}}(D_{j})]$ for possibly different $i,j\in\{1,\ldots,p\}$. The asymptotic results stated in Theorem 3.1 extend those in Theorem 4.2 of Buhl and Klüppelberg (2017), where the observation area $\mathcal{D}_{n}$ is assumed to increase with $n$ in all dimensions. The decomposition (2.4) into a fixed domain $\mathcal{F}$ and an increasing domain ${\mathcal{I}}_{n}$ results in mixing conditions which focus on properties for ${\mathcal{I}}_{n}$ increasing to $\mathbb{Z}^{w}$, while ${\mathcal{F}}$ remains fix and appears in the limit, similarly as in Eq. (3).

Theorem 3.1

Let $\{X(\boldsymbol{s}):\boldsymbol{s}\in\mathbb{R}^{d}\}$ be a strictly stationary regularly varying process, which is observed on $\mathcal{D}_{n}=\mathcal{F}\times\mathcal{I}_{n}$ as in (2.4). Let ${\mathcal{H}}=\{\boldsymbol{h}^{(1)},\ldots,\boldsymbol{h}^{(p)}\}\subset\mathcal{B}(\boldsymbol{0},{\gamma})$ for some ${\gamma}>0$ be a set of observed lag vectors. Suppose that the following conditions are satisfied.

1. (M1)

$\{X(\boldsymbol{s}):\boldsymbol{s}\in\mathbb{R}^{d}\}$* is $\alpha$-mixing with respect to $\mathbb{R}^{w}$ with mixing coefficients $\alpha_{k_{1},k_{2}}(\cdot)$ defined in (A.1).*

There exist sequences $m_{n},r_{n}\to\infty$ with $m_{n}^{d}/n^{w}\to 0$ and $r_{n}^{w}/m_{n}^{d}\to 0$ as $n\to\infty$ such that:

1. (M2)

$m_{n}^{2d}r_{n}^{2w}/n^{w}\to 0$. 2. (M3)

For all $\epsilon>0$, and for all fixed $\boldsymbol{\ell}_{\mathcal{F}}\in\mathbb{R}^{q}$ with $a_{m}=a_{m_{n}}\to\infty$ as in (2.1),

$\lim\limits_{k\to\infty}\limsup\limits_{n\to\infty}\sum\limits_{\boldsymbol{\ell}_{\mathcal{I}}\in\mathbb{Z}^{w}\atop k<\|\boldsymbol{\ell}_{\mathcal{I}}\|\leq r_{n}}m_{n}^{d}\,\mathbb{P}(\max\limits_{\boldsymbol{s}\in\mathcal{B}(\boldsymbol{0},\gamma)}|X(\boldsymbol{s})|>\epsilon a_{m},\max\limits_{\boldsymbol{s}^{\prime}\in\mathcal{B}((\boldsymbol{\ell}_{\mathcal{F}},\boldsymbol{\ell}_{\mathcal{I}}),\gamma)}|X(\boldsymbol{s}^{\prime})|>\epsilon a_{m})=0$.* 3. (M4)

1. (i)

$\lim\limits_{n\to\infty}m_{n}^{d}\sum\limits_{\boldsymbol{\ell}\in\mathbb{Z}^{w}:\|{\boldsymbol{\ell}}\|>r_{n}}\alpha_{1,1}(\|\boldsymbol{\ell}\|)=0$, 2. (ii)

$\sum\limits_{\boldsymbol{\ell}\in\mathbb{Z}^{w}}\alpha_{k_{1},k_{2}}(\|{\boldsymbol{\ell}}\|)<\infty$*  ** for ** $2\leq k_{1}+k_{2}\leq 4$,* 3. (iii)

$\lim\limits_{n\to\infty}m_{n}^{d/2}n{{}^{w/2}}\ \alpha_{1,{n^{w}}}(r_{n})=0$.

Then the empirical extremogram $\widehat{\rho}_{AB,m_{n}}$ defined in (3.2), sampled at lags in ${\mathcal{H}}$ and centred by the pre-asymptotic extremogram $\rho_{AB,m_{n}}$ given in (3.3), is asymptotically normal; i.e.,

[TABLE]

where $\Pi=\mu(A)^{-4}F\Sigma F{}^{{}^{\intercal}}\in\mathbb{R}^{p\times p}$. Writing $\boldsymbol{h}^{(i)}=(\boldsymbol{h}_{\mathcal{F}}^{(i)},\boldsymbol{h}_{\mathcal{I}}^{(i)})$ for $1\leq i\leq p+1$, with the convention that $(\boldsymbol{h}_{\mathcal{F}}^{(p+1)},\boldsymbol{h}_{\mathcal{I}}^{(p+1)})=\boldsymbol{0}$, and recalling (3.5) and (3.6), the matrix $\Sigma\in\mathbb{R}^{(p+1)\times(p+1)}$ has components

[TABLE]

If $i=j$, we have $\Sigma_{ii}=\sigma^{2}_{\mathcal{B}(\boldsymbol{0},\gamma)}(D_{i})$ with $\sigma^{2}_{\mathcal{B}(\boldsymbol{0},\gamma)}(D_{i})$ specified in (3.8). The matrix $F=[F_{1},F_{2}]$ consists of a diagonal matrix $F_{1}$ and a vector $F_{2}$ in the last column:

[TABLE]

Note that condition (M3) is the analogue of condition (3.3) of Davis and Mikosch (2009) in the time series case and thus similar in spirit but weaker than the classical anti-clustering condition $D^{\prime}(\epsilon a_{n})$ as explained there.

Corollary 1

Assume the setting of Theorem 3.1 and suppose that the following conditions are satisfied.

1. (1)

$\{X(\boldsymbol{s}):\boldsymbol{s}\in\mathbb{R}^{d}\}$* is $\alpha$-mixing with respect to $\mathbb{R}^{w}$ with mixing coefficients $\alpha_{k_{1},k_{2}}(\cdot)$ defined in (A.1).* 2. (2)

There exist sequences $m:=m_{n},r:=r_{n}\to\infty$ with $m_{n}^{d}/n^{w}\to 0$ and $r_{n}^{w}/m_{n}^{d}\to 0$ as $n\to\infty$ such that (M3) and (M4i) hold.

Then, as $n\to\infty$,

[TABLE]

Proof

As in part II of the proof of Theorem 3.1 (cf. Appendix A.2), we find that for $i=1,\ldots,p$, as $n\to\infty$,

[TABLE]

where the sets $D_{i}$ and $D_{p+1}$ are defined in (2.7). Convergence in probability follows by Lemma A.1 and Slutzky’s theorem. The last identity holds by definitions (2.3) and (2.5), recalling that ${\mu}_{\mathcal{B}(\boldsymbol{0},{\gamma})}(D_{p+1})=\mu(A)>0$. $\square$

Remark 2

(i)   If the choice $m_{n}=n^{\beta_{1}}$ and $r_{n}=n^{\beta_{2}}$ with $0<\beta_{2}<\beta_{1}d/w<1$ satisfies conditions (M3) and (M4), then for $\beta_{1}\in(0,w/(2d))$ and $\beta_{2}\in(0,\min\{\beta_{1}d/w;1/2-\beta_{1}d/w\})$ the condition (M2) also holds and we obtain the CLT (3.9).

(ii)   The pre-asymptotic extremogram (3.3) in the CLT (3.9) can be replaced by the true one (2.3), if the pre-asymptotic extremogram converges to the true extremogram with the same convergence rate; i.e., if

[TABLE]

(iii)   Unfortunately, for general regularly varying processes, it is not known if the bias condition (3.11) holds, but the CLT (3.9) based on the pre-asymptotic extremogram holds. Hence, the important asymptotic interpretation of the empirical extremogram as a conditional probability of extremal events remains; cf. Cho et al. (2016), Davis and Mikosch (2009), and Drees (2015) and references therein. An important class of processes, where we know conditions such that (3.11) is satisfied or not, are the max-stable processes with finite-dimensional Fréchet marginal distributions, as defined in Section 3.2.

3.2 The extremogram of processes with Fréchet marginal distributions

We start with the definition of max-stable processes.

Definition 4 (Max-stable process)

A process $\{X(\boldsymbol{s}):\boldsymbol{s}\in\mathbb{R}^{d}\}$ is called max-stable if there exist sequences $c_{n}(\boldsymbol{s})>0$ and $d_{n}(\boldsymbol{s})$ for $\boldsymbol{s}\in\mathbb{R}^{d}$ and $n\in\mathbb{N}$ such that

[TABLE]

where $\{X_{j}(\boldsymbol{s}):\boldsymbol{s}\in\mathbb{R}^{d}\}$ are independent replicates of $\{X(\boldsymbol{s}):\boldsymbol{s}\in\mathbb{R}^{d}\}$ and the maximum is taken componentwise.

If max-stable processes have Fréchet marginal distributions, they are regularly varying. Theorem 3.2 below states a necessary and sufficient condition for such processes such that both (3.9) and (3.11) hold, yielding the CLT (3.19) for the empirical extremogram (3.2) centred by the the true one (2.3). In case this condition is not satisfied, Theorem 3.3 states conditions such that (3.19) holds for a bias corrected version of the empirical extremogram.

Theorem 3.2 (CLT for processes with Fréchet margins)

Let $\{X(\boldsymbol{s}):\boldsymbol{s}\in\mathbb{R}^{d}\}$ be a strictly stationary max-stable process with standard unit Fréchet margins, which is observed on ${\mathcal{D}}_{n}={\mathcal{F}}\times{\mathcal{I}}_{n}$ as in (2.4). Let ${\mathcal{H}}=\{\boldsymbol{h}^{(1)},\ldots,\boldsymbol{h}^{(p)}\}\subset\mathcal{B}(\boldsymbol{0},{\gamma})$ for some ${\gamma}>0$ be a set of observed lag vectors. Suppose that conditions (M1)–(M4) of Theorem 3.1 hold for appropriately chosen sequences $m_{n},r_{n}\to\infty$. Let $\rho_{AB}$ be the extremogram (2.3) and $\rho_{AB,m_{n}}$ the pre-asymptotic version (3.3) for sets $A=(\underline{A},\overline{A})$ and $B=(\underline{B},\overline{B})$ with $0<\underline{A}<\overline{A}\leq\infty$ and $0<\underline{B}<\overline{B}\leq\infty.$ Then the limit relation (3.11) holds if and only if $n^{w}/m_{n}^{3d}\to 0$ as $n\to\infty$. In this case we obtain

[TABLE]

with $\Pi$ specified in Theorem 3.1.

Proof

All finite-dimensional distributions are max-stable distributions with standard unit Fréchet margins, hence they are multivariate regularly varying. Furthermore we can choose $a_{m}=m_{n}^{d}$ in Definition 1. Let $V_{2}(\boldsymbol{h};\cdot,\cdot)$ be the bivariate exponent measure defined by $\mathbb{P}(X(\boldsymbol{0})\leq x_{1},X(\boldsymbol{h})\leq x_{2})$ $=\exp\{-V_{2}(\boldsymbol{h};x_{1},x_{2})\}$ for $x_{1},x_{2}>0$, cf. Beirlant et al. (2004), Section 8.2.2. From Lemma A.1(b) of Buhl and Klüppelberg (2017) we know that for $\boldsymbol{h}\in\mathcal{H}$ and with $\overline{V}^{2}_{2}(\boldsymbol{h}):=\underline{A}\overline{A}/(\overline{A}-\underline{A})(V_{2}^{2}(\boldsymbol{h};\overline{A},\overline{B})-V_{2}^{2}(\boldsymbol{h};\overline{A},\underline{B})-V_{2}^{2}(\boldsymbol{h};\underline{A},\overline{B})+V_{2}^{2}(\boldsymbol{h};\underline{A},\underline{B}))$,

[TABLE]

If $\overline{A}=\infty$ and/or $\overline{B}=\infty$, appropriate adaptations need to be taken, which are described in Lemma A.1 of Buhl and Klüppelberg (2017). Hence, for $\boldsymbol{h}\in\mathcal{H}$,

[TABLE]

which converges to 0 if and only if $n^{w}/m_{n}^{3d}\to 0$.

If $n^{w}/m_{n}^{3d}\not\to 0$ in Theorem 3.2, a CLT centred by the true extremogram can still be obtained for a bias corrected empirical estimator. Eq. (3.14) is the basis for such a bias correction if the sets $A$ and $B$ are given by $A=(\underline{A},\infty)$ and $B=(\underline{B},\infty)$ with $\underline{A},\underline{B}>0$. In that case we have

[TABLE]

see Buhl and Klüppelberg (2017), Eq. (A.4). An asymptotically bias corrected estimator is given by

[TABLE]

and we set, covering both cases,

[TABLE]

Theorem 3.3 below guarantees asymptotic normality of the bias corrected extremogram for an—according to Theorem 3.1—valid sequence $m_{n}$ satisfying $n^{w}/m_{n}^{5d}\to 0$. The proof, which is given in Appendix A.3, generalises that of Theorem 4.4 of Buhl et al. (2017), which covers the special case $\underline{A}=\underline{B}=1$ for Brown-Resnick processes.

Theorem 3.3 (CLT for the bias corrected extremogram for processes with Fréchet margins)

Let $\{X(\boldsymbol{s}):\boldsymbol{s}\in\mathbb{R}^{d}\}$ be a strictly stationary max-stable process with standard unit Fréchet margins. Assume the situation of Theorem 3.2 for sets $A=(\underline{A},\infty)$ and $B=(\underline{B},\infty)$ with $\underline{A},\underline{B}>0$. Then the bias corrected extremogram (3.16) is asymptotically normal if and only if $n^{w}/m_{n}^{5d}\to 0$. In that case,

[TABLE]

where $\Pi$ is specified in Theorem 3.1.

Remark 3

From Theorems 3.2 and 3.3 in relation to Remark 2 (i) we deduce two cases:

(I) For $w/(5d)<\beta_{1}\leq w/(3d)$ we cannot replace the pre-asymptotic extremogram by the theoretical version in (3.13), but can resort to a bias correction as described in (3.16) to obtain

[TABLE]

for sets $A=(\underline{A},\infty)$ and $B=(\underline{B},\infty)$ with covariance matrix $\Pi$ specified in Theorem 3.1.

(II) For $w/(3d)<\beta_{1}<w/(2d)$ we obtain indeed

[TABLE]

with covariance matrix $\Pi$ specified in Theorem 3.1.

Observe that Remark 3 generalises Remark 4.1 of Buhl et al. (2017).

4 Generalised least squares extremogram estimates

In this section we fit parametric models to the empirical extremogram using least squares techniques for the parameter estimation. Our approach and extremogram models extend the weighted least squares estimation developed in Steinkohl (2013) and Buhl et al. (2017) considerably. In these papers isotropic space-time models such as the Brown-Resnick model (I) of Section 5.3 below have been estimated by separation of space and time, which is not possible for all models of interest. In what follows we present generalised least squares approaches to fit general parametric extremogram models taking the observation scheme ${\mathcal{D}}_{n}={\mathcal{F}}\times{\mathcal{I}}_{n}$ of a fixed and an increasing domain into account. The approach bears some similarity to the semiparametric variogram estimation in Lahiri et al. (2002).

Our setting is as follows. Let $\{\rho_{AB,\boldsymbol{\theta}}(\boldsymbol{h}):\boldsymbol{h}\in\mathbb{R}^{d},\boldsymbol{\theta}\in\Theta\}$ be some parametric valid extremogram model with parameter space $\Theta$ and continuous in $\boldsymbol{h}\in\mathbb{R}^{d}$. Assume that $\rho_{AB}(\cdot)=\rho_{AB,\boldsymbol{\theta}^{\star}}(\cdot)$ with true parameter vector $\boldsymbol{\theta}^{\star}$, which lies by assumption in the interior of $\Theta$. Denote by $\widehat{\rho}_{AB,m_{n}}(\boldsymbol{h})$ any of the estimators of Theorem 3.1, Theorem 3.2, or Theorem 3.3 for the appropriately chosen $\mu$-continuous Borel sets $A$ and $B$ such that $\mu(A)>0$ and lags $\boldsymbol{h}\in\mathcal{H}=\{\boldsymbol{h}^{(1)},\ldots,\boldsymbol{h}^{(p)}\}$.

First note that under the much weaker conditions of Corollary 1 the empirical extremogram is a consistent estimator of the extremogram such that as $n\rightarrow\infty$,

[TABLE]

Under more restrictive conditions needed for the three CLTs above,

[TABLE]

where $\Pi$ is the covariance matrix specified in Theorem 3.1.

As we shall prove below, consistency of the empirical extremogram entails consistent generalised least squares parameter estimates, whereas asymptotic normality of the empirical extremogram entails asymptotically normal generalised least squares parameter estimates.

Definition 5 (Generalised least squares extremogram estimator (GLSE))

Let $\{X(\boldsymbol{s}):\boldsymbol{s}\in\mathbb{R}^{d}\}$ be a strictly stationary regularly varying process, which is observed on $\mathcal{D}_{n}=\mathcal{F}\times\mathcal{I}_{n}$ as in (2.4). Let $A$ and $B$ be $\mu$-continuous Borel sets in $\overline{\mathbb{R}}\backslash\{0\}$ such that $\mu(A)>0$. For a sequence $m=m_{n}\rightarrow\infty$ and $m_{n}=o(n)$ as $n\to\infty$ define for $\boldsymbol{\theta}\in\Theta$ the column vector

[TABLE]

For some non-singular positive definite weight matrix $V(\boldsymbol{\theta})\in\mathbb{R}^{p\times p}$, the GLSE is defined as

[TABLE]

Assumption 3 presents a set of conditions, which imply consistency and asymptotic normality of the GLSE.

Assumption 3

Assume the situation of Definition 5. We shall require the following conditions.

1. (G1)

Consistency: $\widehat{\rho}_{AB,m_{n}}(\boldsymbol{h}^{(i)})\stackrel{{\scriptstyle P}}{{\rightarrow}}\rho_{AB,\boldsymbol{\theta}^{\star}}(\boldsymbol{h}^{(i)})$ as $n\to\infty$ for $i=1,\ldots,p.$ 2. (G2)

Asymptotic normality: $\sqrt{\dfrac{n^{w}}{m_{n}^{d}}}\widehat{\boldsymbol{g}}_{n}(\boldsymbol{\theta}^{\star})\stackrel{{\scriptstyle\mathcal{D}}}{{\rightarrow}}\mathcal{N}(\boldsymbol{0},\Pi)$ as $n\to\infty$. 3. (G3)

*Identifiability condition: For all $\epsilon>0$ there exists some $\delta>0$ such that *

\inf\Big{\{}\sum\limits_{i=1}^{p}(\rho_{AB,\boldsymbol{\theta}_{1}}(\boldsymbol{h}^{(i)})-\rho_{AB,\boldsymbol{\theta}_{2}}(\boldsymbol{h}^{(i)}))^{2}:\boldsymbol{\theta}^{(1)},\boldsymbol{\theta}^{(2)}\in\Theta,\|{\boldsymbol{\theta}^{(1)}-\boldsymbol{\theta}^{(2)}}\|\geq\epsilon\Big{\}}>\delta.* If the parameter space $\Theta$ is compact, this condition can be replaced by the weaker condition*

[TABLE] 4. (G4)

Smoothness condition 1: For all $i=1,\ldots,p$:

$\rho_{AB,\boldsymbol{\theta}}(\boldsymbol{h}^{(i)})$* has continuous partial derivatives of order $z_{1}\geq 0$ w.r.t. $\boldsymbol{\theta}$, where $z_{1}=0$ corresponds to $\rho_{AB,\boldsymbol{\theta}}(\boldsymbol{h}^{(i)})$ being continuous in $\boldsymbol{\theta}$.* 5. (G5)

Smoothness condition 2:

1. (i)

$\sup\limits_{\boldsymbol{\theta}\in\Theta}\{\|V(\boldsymbol{\theta})\|_{M}+\|V(\boldsymbol{\theta})^{-1}\|_{M}\}<\infty,$* where $\|\cdot\|_{M}$ is some arbitrary matrix norm.* 2. (ii)

The matrix valued function $V(\boldsymbol{\theta})$ has continuous derivatives of order $z_{2}\geq 0$ w.r.t. $\boldsymbol{\theta}$, where $z_{2}=0$ corresponds to $V(\boldsymbol{\theta})$ being continuous in $\boldsymbol{\theta}$. 6. (G6)

*Rank condition: For $\boldsymbol{\theta}=(\theta_{1},\ldots,\theta_{k})\in\Theta\subset\mathbb{R}^{k}$ we denote by $\mathrm{P}_{AB}(\boldsymbol{\theta})$ the Jacobian matrix of $(-\rho_{AB,\boldsymbol{\theta}}(\boldsymbol{h}^{(1)}),\ldots,-\rho_{AB,\boldsymbol{\theta}}(\boldsymbol{h}^{(p)})){}^{{}^{\intercal}}$; i.e., *

[TABLE]

*The Jacobian matrix has full rank: $\textnormal{rank}(\mathrm{P}_{AB}(\boldsymbol{\theta}^{\star}))=k$.   *$\Box$

The proof of the next theorem can be found in Appendix A.4.

Theorem 4.1 (Consistency and asymptotic normality of the GLSE)

Assume the situation of Definition 5. If Assumptions 3(G1) and (G3) hold as well as (G4) and (G5) for $z_{1}=z_{2}=0$, respectively, then the GLSE is consistent; i.e.,

[TABLE]

If Assumption 3(G2) and (G3) hold as well as (G4) and (G5) for $z_{1}=z_{2}=1$, respectively, and the rank condition (G6) holds, then the GLSE is asymptotically normal; i.e.,

[TABLE]

with asymptotic covariance matrix

[TABLE]

where B(\boldsymbol{\theta}^{\star}):=\big{(}\mathrm{P}_{AB}(\boldsymbol{\theta}^{\star}){}^{{}^{\intercal}}[V(\boldsymbol{\theta}^{\star})+V(\boldsymbol{\theta}^{\star}){}^{{}^{\intercal}}]\mathrm{P}_{AB}(\boldsymbol{\theta}^{\star})\big{)}^{-1} and $\Pi$ is the asymptotic covariance matrix in Eq. (4.2).

Remark 4

The quality of the GLSE depends on the matrix $V(\boldsymbol{\theta})$. Simple choices for the matrix $V(\boldsymbol{\theta})$ in (4.4) are the identity matrix, leading to the ordinary least squares estimator, or some general weight matrix, leading to weighted least squares estimators.

An asymptotically optimal matrix $V(\boldsymbol{\theta})$ can be obtained as follows. Let $\Pi=\Pi(\boldsymbol{\theta}^{\star})$ be the asymptotic covariance matrix of the empirical extremogram in Eq. (4.2). Assume that $\Pi(\boldsymbol{\theta}^{\star})$ has a closed form that depends on the true parameter vector $\boldsymbol{\theta}^{\star}$ which can be extended to a matrix function $\Pi(\boldsymbol{\theta})$ on the whole parameter space $\boldsymbol{\Theta}$. Assume also that the inverse $V(\boldsymbol{\theta})=\Pi^{-1}(\boldsymbol{\theta})$ exists for all $\boldsymbol{\theta}\in\boldsymbol{\Theta}$ and satisfies the Assumption 3(G5) for $z_{2}=1$. Then, as pointed out in Lahiri et al. (2002), Theorem 4.1, for spatial variogram estimators and in Einmahl et al. (2016), Corollary 2.3, for extreme parameter estimation based on iid random vector observations, the resulting asymptotic covariance matrix $\Pi_{V}=\Pi_{V}(\boldsymbol{\theta}^{\star})$ of the GLSE in (4.7) is asymptotically optimal among all valid matrices $V^{\prime}=V^{\prime}(\boldsymbol{\theta})$. This means that $\Pi_{V}$ is minimal in the sense that for all valid matrices $V^{\prime}$, the difference $\Pi_{V^{\prime}}-\Pi_{V}$ is positive semidefinite.

5 Estimation of Brown-Resnick space-time processes

5.1 Brown-Resnick processes

We consider a strictly stationary Brown-Resnick process with spectral representation

[TABLE]

where $\{\xi_{j}:j\in\mathbb{N}\}$ are points of a Poisson process on $[0,\infty)$ with intensity $\xi^{-2}d\xi$, the dependence function $\delta$ is nonnegative and conditionally negative definite, and $\{W_{j}(\boldsymbol{s}):\boldsymbol{s}\in\mathbb{R}^{d}\}$ are independent replicates of a Gaussian process $\{W(\boldsymbol{s}):\boldsymbol{s}\in\mathbb{R}^{d}\}$ with stationary increments, $W(\boldsymbol{0})=0$, $\mathbb{E}[W(\boldsymbol{s})]=0$ and covariance function

[TABLE]

Spectral representations of max-stable processes go back to de Haan (1984) and Giné, Hahn, and Vatan (1990), the specific representation (5.1) to Brown and Resnick (1977) in a time series context, to Kabluchko et al. (2009) in a spatial and to Davis et al. (2013a) in a space-time setting. The univariate margins of the process $\eta$ follow standard unit Fréchet distributions. Non-stationary Brown-Resnick models have recently been discussed and fitted to data in Asadi et al. (2015), Engelke et al. (2015), and Huser and Genton (2016).

There are various quantities to describe the dependence in (5.1), where explicit expressions can be derived:

1. $\bullet$

In geostatistics, the dependence function $\delta$ is termed the semivariogram of the process $\{W(\boldsymbol{s}):\boldsymbol{s}\in\mathbb{R}^{d}\}$ based on the fact that for $\boldsymbol{s}^{(1)},\boldsymbol{s}^{(2)}\in\mathbb{R}^{d}$,

[TABLE] 2. $\bullet$

For $\boldsymbol{h}\in\mathbb{R}^{d}$, the tail dependence coefficient is given by (see e.g. Davis, Klüppelberg, and Steinkohl (2013a), Section 3)

[TABLE]

where $\Phi$ denotes the standard normal distribution function. 3. $\bullet$

For $D=\{\boldsymbol{s}^{(1)},\ldots,\boldsymbol{s}^{(|D|)}\}$ and $\boldsymbol{y}=(y_{1},\ldots,y_{|D|})>\boldsymbol{0}$ the finite-dimensional margins are given by

[TABLE]

Here $V_{D}$ denotes the exponent measure (cf. Beirlant et al. (2004), Section 8.2.2), which is homogeneous of order -1 and depends solely on the dependence function $\delta$. For $D=\{\boldsymbol{s},\boldsymbol{s}+\boldsymbol{h}\}$ where $\boldsymbol{s}\in\mathbb{R}^{d}$ and $\boldsymbol{h}\in\mathbb{R}^{d}$ is some fixed lag vector, we get (cf. Davis et al. (2013a), Section 3)

[TABLE]

with

[TABLE] 4. $\bullet$

For $\boldsymbol{h}\in\mathbb{R}^{d}$ and sets $A=(\underline{A},\overline{A})$ and $B=(\underline{B},\overline{B})$ with $0<\underline{A}<\overline{A}\leq\infty$ and $0<\underline{B}<\overline{B}\leq\infty$, the extremogram (2.3) is given by (see Buhl and Klüppelberg (2017), Eq. (A.1))

[TABLE]

for $V_{2}$ as in (5.4). For $A=(\underline{A},\infty)$ and $B=(\underline{B},\infty)$ we get formula (31) of Cho et al. (2016):

[TABLE] 5. $\bullet$

The extremal coefficient $\xi_{D}$ (see Beirlant et al. (2004), Section 8.2.7) for any finite set $D\subset\mathbb{R}^{d}$ is defined as

[TABLE]

i.e., $\xi_{D}=V_{D}(1,\ldots,1).$ If $|D|=2$ and $\boldsymbol{h}=\boldsymbol{s}^{(1)}-\boldsymbol{s}^{(2)}$, then

[TABLE]

where the first identity holds in general (cf. Beirlant et al. (2004), Section 9.5.1), and the last one by (5.2).

Our aim is to fit a parametric extremogram model of a Brown-Resnick process (5.1) based on observations given in ${\mathcal{D}}_{n}={\mathcal{F}}\times{\mathcal{I}}_{n}$ as in (2.4). This approach is semiparametric in the sense that we first compute (possibly bias corrected) empirical estimates (3.16) of the extremogram ${\rho}_{AB}(\boldsymbol{h})$ for different $\boldsymbol{h}\in\mathcal{H}$, and fit a parametric model $\rho_{AB,\boldsymbol{\theta}}(\boldsymbol{h})$ by GLSE to the empirical extremogram. For sets $A=B=(\underline{A},\infty)$ with $\underline{A}>0$, this yields an estimator of the dependence function, since by (5.5) and (5.7) there is a one-to-one relation between extremogram and dependence function.

5.2 Asymptotic properties of the empirical extremogram of a Brown-Resnick process

Let $\{\eta(\boldsymbol{s}):\boldsymbol{s}\in\mathbb{R}^{d}\}$ be a strictly stationary Brown-Resnick process as in (5.1) with some valid (i.e., nonnegative and conditionally negative definite) dependence function $\delta$. Before investigating the asymptotic properties of the GLSE, we state sufficient conditions for $\delta$ so that the regularity conditions of Theorem 3.1 are satisfied.

Theorem 5.1

*Let $\{\eta(\boldsymbol{s}):\boldsymbol{s}\in\mathbb{R}^{d}\}$ be a strictly stationary Brown-Resnick process as in (5.1), observed on $\mathcal{D}_{n}=\mathcal{F}\times\mathcal{I}_{n}$ as in (2.4). Let ${\mathcal{H}}=\{\boldsymbol{h}^{(1)},\ldots,\boldsymbol{h}^{(p)}\}\subset\mathcal{B}(\boldsymbol{0},{\gamma})$ for some ${\gamma}>0$ be a set of observed lag vectors. Assume sequences *

[TABLE]

Writing $\boldsymbol{v}=(\boldsymbol{v}_{\mathcal{F}},\boldsymbol{v}_{{\mathcal{I}}})\in\mathbb{R}^{q}\times\mathbb{R}^{w}$ according to the fixed and increasing domains, assume that the dependence function $\delta$ satisfies for arbitrary fixed finite set $L\subset\mathbb{Z}^{q}:$

1. (A)

m_{n}^{d}\sum\limits_{z>r_{n}}z^{w-1}\exp\Big{\{}-\frac{1}{4}\inf\limits_{\boldsymbol{v}\in L\times\mathbb{Z}^{w}:\|\boldsymbol{v}_{\mathcal{I}}\|\geq z}\delta(\boldsymbol{v})\Big{\}}\rightarrow 0* as $n\rightarrow\infty$.* 2. (B)

m_{n}^{d/2}n^{(3w)/2}\exp\Big{\{}-\frac{1}{4}\inf\limits_{\boldsymbol{v}\in L\times\mathbb{Z}^{w}:\|\boldsymbol{v}_{\mathcal{I}}\|>r_{n}}\delta(\boldsymbol{v})\Big{\}}\rightarrow 0* as $n\rightarrow\infty.$*

Then conditions (M1)-(M4) of Theorem 3.1 are satisfied, and the empirical extremogram $\widehat{\rho}_{AB,m_{n}}$ defined in (3.2) sampled at lags in ${\mathcal{H}}$ and centred by the pre-asymptotic extremogram $\rho_{AB,m_{n}}$ given in (3.3), is asymptotically normal; i.e.,

[TABLE]

where the covariance matrix $\Pi$ is specified in Theorem 3.1.

Proof

First note that, since all finite-dimensional distributions are max-stable distributions with standard unit Fréchet margins, they are multivariate regularly varying. We first show (M3). Let $\epsilon>0$ and fix $\boldsymbol{\ell}_{\mathcal{F}}\in\mathbb{R}^{q}$. For ${\gamma}>0$ define the set

[TABLE]

Note that, writing $\boldsymbol{s}_{1}=(\boldsymbol{f}_{1},\boldsymbol{i}_{1})$ and $\boldsymbol{s}_{2}=(\boldsymbol{f}_{2},\boldsymbol{i}_{2})\in\mathbb{R}^{q}\times\mathbb{R}^{w}$ according to the fixed and increasing domains as before, it can be decomposed into $L_{{\gamma}}(\boldsymbol{\ell}_{\mathcal{F}},\boldsymbol{\ell}_{\mathcal{I}})=L_{{\gamma}}^{(1)}\times L_{{\gamma}}^{(2)}(\boldsymbol{\ell}_{{\mathcal{I}}})$ where $L_{{\gamma}}^{(1)}:=\{\boldsymbol{f}_{1}-\boldsymbol{f}_{2}:\boldsymbol{s}_{1}\in\mathcal{B}((\boldsymbol{0},\boldsymbol{0}),{\gamma}),\boldsymbol{s}_{2}\in\mathcal{B}((\boldsymbol{\ell}_{\mathcal{F}},\boldsymbol{0}),{\gamma})\}$, which is independent of $\boldsymbol{\ell}_{{\mathcal{I}}}$, and $L_{{\gamma}}^{(2)}(\boldsymbol{\ell}_{{\mathcal{I}}}):=\{\boldsymbol{i}_{1}-\boldsymbol{i}_{2}:\boldsymbol{s}_{1}\in\mathcal{B}((\boldsymbol{0},\boldsymbol{0}),{\gamma}),\boldsymbol{s}_{2}\in\mathcal{B}((\boldsymbol{\ell}_{\mathcal{F}},\boldsymbol{\ell}_{{\mathcal{I}}}),{\gamma})\}$. Then, recalling that $a_{m}=m_{n}^{d}$, and using a second order Taylor expansion as in the proof of Theorem 4.3 of Buhl et al. (2017), we have as $n\rightarrow\infty$,

[TABLE]

Therefore,

[TABLE]

Since the number of grid points $\boldsymbol{\ell}_{{\mathcal{I}}}$ in $\mathbb{Z}^{w}$ with norm $\|\boldsymbol{\ell}_{\mathcal{I}}\|=z$ is of order $\mathcal{O}(z^{w-1})$, there exists a positive constant $C$ such that the right hand side can be bounded from above by

[TABLE]

where we have used in the second last step that $1-\Phi(x)\leq\exp\{-x^{2}/2\}$ for $x>0$ and in the last step the decomposition $L_{{\gamma}}(\boldsymbol{\ell}_{\mathcal{F}},\boldsymbol{\ell}_{\mathcal{I}})=L_{{\gamma}}^{(1)}\times L_{{\gamma}}^{(2)}(\boldsymbol{\ell}_{{\mathcal{I}}})$. By condition (A), since we can neglect the constant ${\gamma}$, we have

[TABLE]

Together with $r_{n}^{w}=o(m_{n}^{d})$ as $n\rightarrow\infty$, this implies that

[TABLE]

Next we prove (M1) and (M4i)-(M4iii). To this end we bound the $\alpha$-mixing coefficients $\alpha_{k_{1},k_{2}}(\cdot)$ for $k_{1},k_{2}\in\mathbb{N}$ of $\{\eta(\boldsymbol{s}):\boldsymbol{s}\in\mathbb{R}^{d}\}$ with respect to $\mathbb{R}^{w}$, which are defined in (A.1). Observe that $d(\Lambda_{1},\Lambda_{2})$ for sets $\Lambda_{i}\subset\mathbb{Z}^{w}$ as in Definition 6 can only get large within the increasing domain. Define the set

[TABLE]

We use Eq. (5.8), as well as Dombry and Eyi-Minko (2012), Eq. (3) and Corollary 2.2 to obtain

[TABLE]

By condition (A) we have $\alpha_{k_{1},k_{2}}(z)\rightarrow 0$, since necessarily $\inf\limits_{\boldsymbol{v}\in L_{\mathcal{F}}\times\mathbb{Z}^{w}:\|\boldsymbol{v}_{\mathcal{I}}\|\geq z}\delta(\boldsymbol{v})\rightarrow\infty$ as $z\rightarrow\infty$ and, therefore, the process $\{\eta(\boldsymbol{s}):\boldsymbol{s}\in\mathbb{R}^{d}\}$ is $\alpha$-mixing; i.e., (M1) holds. We continue by estimating

[TABLE]

by condition (A). This shows (M4i). Similarly, it can be shown that (M4ii) holds, if (A) is satisfied. Finally, we show (M4iii). Using Eq. (5.11), we find

[TABLE]

as $n\rightarrow\infty$ because of condition (B).

The following is an immediate corollary of Theorem 5.1.

Corollary 2

Assume the setting of Theorem 5.1. Suppose that the dependence function $\delta$ satisfies for positive constants $C$ and $\alpha$, and for an arbitrary norm $\|\cdot\|$ on $\mathbb{R}^{w}$,

[TABLE]

for every $\boldsymbol{v}=(\boldsymbol{v}_{{\mathcal{F}}},\boldsymbol{v}_{{\mathcal{I}}})\in L\times\mathbb{Z}^{w}$, where $L\subset\mathbb{Z}^{q}$ is arbitrary, but fixed. In particular, $\delta(\boldsymbol{v})\to\infty$ if $\|\boldsymbol{v}_{{\mathcal{I}}}\|\to\infty$. With $m_{n}=n^{\beta_{1}}$ and $r_{n}=n^{\beta_{2}}$ with $\beta_{1}\in(0,w/(2d))$ and $\beta_{2}\in\min\{\beta_{1}d/w;1/2-\beta_{1}d/w\}$, the conditions of Theorem 5.1 are satisfied for $\{\eta(\boldsymbol{s}):\boldsymbol{s}\in\mathbb{R}^{d}\}$ and we conclude

[TABLE]

Proof

Due to equivalence of norms on $\mathbb{R}^{w}$ we will make no difference between the norm in (5.12) and the one used in Theorem 5.1. Clearly the sequences $m_{n}$ and $r_{n}$ satisfy the requirements $m_{n},r_{n}\rightarrow\infty$, $m_{n}^{d}/n^{w}\to 0$, $r_{n}^{w}/m_{n}^{d}\to 0$ and $m_{n}^{2d}r_{n}^{2w}/n^{w}\to 0$ as $n\to\infty.$ We have for $z>0$,

[TABLE]

Condition (B) of Theorem 5.1 is satisfied since

[TABLE]

Condition (A) holds since by Lemma A.3 of Buhl et al. (2017), there is a positive constant $K$ such that for sufficiently large $n$ the sequence $z^{w-1}\exp\{-Cz^{\alpha}/4\}$ is decreasing for $z\geq r_{n}$,

[TABLE]

With the particular choice of sequences $m_{n}=n^{\beta_{1}}$ and $r_{n}=n^{\beta_{2}}$ given in Corollary 2, we are in the setting of Remark 3. Hence, in addition to the CLT (5.13), we obtain the CLT (3.19) of the empirical extremogram centred by the true one and the CLT (3.18) corresponding to the bias corrected estimator.

Remark 5

1. (i)

Corollary 2 requires the dependence function $\delta$ of the Brown-Resnick process to be unbounded. This requirement is not satisfied, for example, by the Schlather model or extremal-$t$-models, which do not capture possible extremal independence between two process values; see for example Davison et al. (2012c), Section 6.1 and Opitz (2013), Section 4. 2. (ii)

Other prominent max-stable processes that satisfy the conditions of Theorem 3.1 are the max-moving average processes (see Example 4.6 of Buhl and Klüppelberg (2017)) or special cases of the random set model in Huser and Davison (2014).

5.3 Space-time Brown-Resnick processes: different models for the extremogram

We explore the semiparametric estimation for strictly stationary Brown-Resnick processes in their space-time form $\{\eta(\boldsymbol{s},t):\boldsymbol{s}\in\mathbb{R}^{d-1},t\in[0,\infty)\}$. For three classes of parametric models for the dependence function $\delta_{\boldsymbol{\theta}}$ we prove that the GLSE is consistent and asymptotically normal.

Note that by Eq. (5.7) every model $\{\delta_{\boldsymbol{\theta}}:\boldsymbol{\theta}\in\Theta\}$ for the dependence function yields a model $\{\rho_{AB,\boldsymbol{\theta}}:\boldsymbol{\theta}\in\Theta\}$ for its space-time extremogram. Moreover, the extremogram (5.7) is always of the same form, and only $\widetilde{\Phi}$ in (5.5) changes with the model. We consider three different model classes, which together cover a large field of environmental applications such as the modelling of extreme precipitation (cf. Davis et al. (2013a), Buhl and Klüppelberg (2016), de Fondeville and Davison (2016), Buhl et al. (2017)), extreme wind speed (cf. Engelke et al. (2015)) or extremes on river networks (cf. Asadi et al. (2015)), provided they are valid (i.e., nonnegative and conditionally negative definite) dependence functions in the considered metric.

(I)   Fractional space-time model.

Davis et al. (2013a) introduce the spatially isotropic model

[TABLE]

with parameter vector

[TABLE]

The isotropy assumption, where (5.14) depends on the norm of the spatial lag $\boldsymbol{h}$, can be relaxed in a natural way by introducing geometric anisotropy. We only discuss the case $d-1=2$, but the approach is easily transferable to higher dimensions. Let $\varphi\in[0,\pi/2)$ be a rotation angle and $R=R(\varphi)$ a rotation matrix, and $T$ a dilution matrix with $c>0$; more precisely,

[TABLE]

The geometrically anisotropic model is then given by

[TABLE]

where $A=TR$ is the transformation matrix. The parameter vector of the transformed model is

[TABLE]

For more details about geometric anisotropy see Blanchet and Davison (2011), Section 4.2, Davis et al. (2013a), Section 4.2, or Engelke et al. (2015), Section 5.2.

(II)   Spatial anisotropy along orthogonal spatial directions

Buhl and Klüppelberg (2016) generalize the fractional isotropic model (5.14) to

[TABLE]

with parameter vector

[TABLE]

It is more flexible than the isotropic model (I) as it allows for different rates of decay of extreme dependence along the axes of a $d$-dimensional spatial grid. Arbitrary principal orthogonal directions can be introduced by a rotation matrix $R$ as introduced for the isotropic model in (I), here described for the case $d-1=2$:

[TABLE]

The new parameter vector is

[TABLE]

In Buhl and Klüppelberg (2016) this model is applied to extreme precipitation in Florida and, according to a specifically developed goodness-of-fit method, performs extremely well.

(III) Time-shifted Brown-Resnick processes

With the goal to allow for some influence of the spatial dependence from previous values of the process we time-shift the Gaussian processes in the definition of the Brown-Resnick model (5.1). For $\boldsymbol{\tau}=(\tau_{1},\tau_{2})\in\mathbb{R}^{d-1}$ define

[TABLE]

Then $\{W^{(\boldsymbol{\tau})}(\boldsymbol{s},t):\boldsymbol{s}\in\mathbb{R}^{d-1},t\in[0,\infty)\}$ is also a centred Gaussian process starting in 0 with stationary increments: for $(\boldsymbol{s}^{(1)},t^{(1)}),(\boldsymbol{s}^{(2)},t^{(2)})\in\mathbb{R}^{d-1}\times[0,\infty)$, because of the stationary increments of $\{W(\boldsymbol{s},t)\}$, where $\stackrel{{\scriptstyle d}}{{=}}$ stands for equality in distribution,

[TABLE]

The corresponding time-shifted dependence function is given by

[TABLE]

which yields the covariance function

[TABLE]

By Theorem 10 of Kabluchko et al. (2009) the process

[TABLE]

defines a strictly stationary space-time Brown-Resnick process.

This method does not depend on the specific dependence function: every Brown-Resnick process $\{\eta(\boldsymbol{s},t):(\boldsymbol{s},t)\in\mathbb{R}^{d-1},t\in[0,\infty)\}$ with dependence function $\{\delta_{\boldsymbol{\theta}},{\boldsymbol{\theta}}\in\Theta\}$ results in a time-shifted Brown-Resnick process with dependence function

$\{\delta^{(\boldsymbol{\tau})}_{\boldsymbol{\theta}},{\boldsymbol{\theta}}\in\Theta,\boldsymbol{\tau}\in\mathbb{R}^{d-1}\}$. To give an example, for the Brown-Resnick process (II) without rotation, the parametrised time-shifted dependence function is given by

[TABLE]

with parameter vector

[TABLE]

This model is somewhat motivated by the time-shifted moving maxima Brown-Resnick process introduced by Embrechts et al. (2016), it is however much simpler to analyse and to estimate. As a referee has pointed out, similar models have been suggested in Section 5.3.2, models (ii)-(iv) on p. 213 in Huser (2013).

In the following we show that models (I)-(III) satisfy Assumption 3 and the conditions of Theorem 4.1 and Corollary 2.

Asymptotic properties of models (I)-(III)

As before, we assume space-time observations on $\mathcal{D}_{n}={\mathcal{S}}\times\mathcal{T}=({\mathcal{S}}\times\mathcal{T})(n)$, where ${\mathcal{S}}\subset\mathbb{Z}^{d-1}$ are the spatial and $\mathcal{T}\subset\mathbb{Z}$ the time series observations. Moreover, we assume that they decompose into $\mathcal{D}_{n}=\mathcal{F}\times\mathcal{I}_{n}$, where $\mathcal{F}\subset\mathbb{Z}^{q}$ is some fixed domain and $\mathcal{I}_{n}=\{1,\ldots,n\}^{w}$ is a sequence of regular grids, and $q+w=d$.

For two points $(\boldsymbol{s}^{(1)},t^{(1)})$ and $(\boldsymbol{s}^{(2)},t^{(2)})\in\mathbb{R}^{d-1}\times[0,\infty)$, we denote by $(\boldsymbol{h},u)=(\boldsymbol{s}^{(1)},t^{(1)})-(\boldsymbol{s}^{(2)},t^{(2)})\in\mathbb{R}^{d}$ their space-time lag vector. Furthermore, we choose Borel sets $A=B=(\underline{A},\infty)$ for some $\underline{A}>0$. We denote by $\widehat{\rho}_{AB,m_{n}}(\boldsymbol{h},u)$ the (possibly bias-corrected) empirical space-time extremogram (3.16), sampled at lags in ${\mathcal{H}}\subset\mathbb{R}^{d}$, and by $\widehat{\boldsymbol{\theta}}_{n,V}$ the GLSE (4.4), referring to some positive definite weight matrix $V$.

To show consistency and asymptotic normality of the corresponding GLSE, we need to verify the assumptions required in Theorem 4.1; i.e. the relevant parts of Assumption 3. Note that Corollary 2 applies for all models, since they all satisfy $\delta_{\boldsymbol{\theta}}(\boldsymbol{h},u)\geq C|u|^{\alpha}$ for $C>0$ and $\alpha\in(0,2]$. Thus we obtain the CLTs of the empirical extremogram centred by the pre-asymptotic extremogram (5.13), centred by the true one (3.13) and of the bias corrected empirical extremogram centred by the true one (3.18). Hence (G1) and (G2) hold for the empirical extremogram. Furthermore, we assume that the parameter space $\Theta\subset\mathbb{R}^{k}$, which contains the true parameter $\boldsymbol{\theta}^{\star}$ as an interior point, is a compact subset of the spaces introduced above for the corresponding models.

The following requirements concern the model-independent assumptions.

1. $\bullet$

In order to determine the GLSE we need to choose a positive definite matrix $V(\boldsymbol{\theta})$ for $\boldsymbol{\theta}\in\Theta$, and we take one, which satisfies condition (G5ii) with $z_{2}=1$. Due to compactness of the parameter space $\Theta$, condition (G5i) is therefore automatically satisfied. 2. $\bullet$

We require that $|\mathcal{H}|\geq k$, such that the rank condition (G6) can be satisfied.

Next we discuss the model-dependent assumptions. First note that the smoothness condition (G4) is satisfied for $z_{1}=0$ for all models $\{\rho_{AB,\boldsymbol{\theta}}(\cdot)\}$ (equivalently $\{\delta_{\boldsymbol{\theta}}(\cdot)\}$). Furthermore, due to compactness of the parameter space, it suffices to show condition (G3’) in order to verify identifiability of the models. Condition (G3’) is satisfied for models (I)-(III) if for two distinct parameter vectors $\boldsymbol{\theta}^{(1)}\neq\boldsymbol{\theta}^{(2)}$ there is at least one $(\boldsymbol{h},u)\in\mathcal{H}$ such that $\rho_{AB,\boldsymbol{\theta}^{(1)}}(\boldsymbol{h},u)\neq\rho_{AB,\boldsymbol{\theta}^{(2)}}(\boldsymbol{h},u)$ or, equivalently, $\delta_{\boldsymbol{\theta}^{(1)}}(\boldsymbol{h},u)\neq\delta_{\boldsymbol{\theta}^{(2)}}(\boldsymbol{h},u)$. This holds due to the power function structure of the models. For the geometric anisotropic model in (I) we need to exclude $c=1$ to ensure identifiability of the angle $\varphi$; however, if $c=1$ then $\varphi$ has no influence on the dependence function and can be neglected. Thus, the GLSEs are consistent according to Theorem 4.1.

We now turn to the CLT (4.7), where it remains to show (G4) for $z_{1}=1$. Difficulties arise due to norms and absolute values of certain parameters in the model equations:

1. $\bullet$

In their basic forms without rotation or dilution, models (I) and (II) are infinitely often continuously partially differentiable in the model parameters. Hence asymptotic normality of the GLSEs follows by Theorem 4.1. 2. $\bullet$

If rotation and/or dilution parameters are included, continuous partial differentiability still holds under the following restrictions: Let $\alpha_{1}$ (for model (I)) or

$\alpha_{1},\ldots,\alpha_{d-1}$ (for model (II)) be the spatial smoothness parameters. Since they are the powers of some norm or absolute value, restricting them to values in $[1,2]$ makes the models continuously partially differentiable; otherwise, they are partially differentiable everywhere but not in 0. As to model (II), in the case $d-1=2$, one of the parameters $\alpha_{1}$ and $\alpha_{2}$ being larger than 1 is already sufficient. To see this, recall that the spatial part of the dependence function is given by

[TABLE]

Assume w.l.o.g that $\alpha_{2}>1$. Then critical values of $\varphi\in[0,\pi/2)$ are the roots of $h_{1}\cos\varphi-h_{2}\sin\varphi$. Given a value $h_{2}\in\mathbb{R}$ we need to choose $h_{1}\in\mathbb{R}$ such that $h_{1}\neq h_{2}\tan\varphi$ for all $\varphi\in[0,\pi/2)$. Since $\tan\varphi>0$ for $\varphi\in[0,\pi_{2})$, we can choose $h_{1}$ such that $\text{sgn}(h_{1})=-\text{sgn}(h_{2})$. If all lags $(h_{1},h_{2},u)\in\mathcal{H}$ are chosen such that $(h_{1},h_{2})$ have opposite signs (or, trivially, are equal to $(0,0)$) and if ${\rm rank}(\mathrm{P}_{AB}(\boldsymbol{\theta}^{\star}))=k$, then the GLSE is asymptotically normal. 3. $\bullet$

Model (III) is continuous partially differentiable, if the spatial smoothness parameters $\alpha_{i}$ for $i=1,\ldots,d-1$ are all larger than 1. If $\alpha_{i}\leq 1$ for some $i$, then the term $C_{i}|h_{i}-u\tau_{i}|^{\alpha_{i}}$ is, as a function of $\tau_{i}$, not differentiable at $\tau_{i}=h_{i}/u\in\mathbb{R}$. However, it is possible to restrict the parameter space such that such equalities do not occur.

6 Simulation study

Specifications

Consider the framework of Section 5.3. In particular, let $\{\eta(\boldsymbol{s},t):\boldsymbol{s}\in\mathbb{R}^{2},t\in[0,\infty)\}$ be a strictly stationary space-time Brown-Resnick process (5.1) observed on ${\mathcal{D}}_{n}={\mathcal{F}}\times{\mathcal{I}}_{n}$. Denote by $\widehat{\rho}_{AB,m_{n}}(\boldsymbol{h},u)$ the space-time version of the (possibly bias corrected) empirical extremogram given in (3.16), sampled at lags in ${\mathcal{H}}\subset\mathbb{R}^{d}$, where ${\mathcal{H}}$ is specified below and we choose the sets $A=B=(1,\infty)$. As already indicated in its Definition 2(1), the computation involves the practical issue of choosing the value $a_{m_{n}}=m_{n}=:q$ as a large quantile, where the first equality is due to the standard unit Fréchet distribution of the marginals of the Brown-Resnick model, so that $q$ should be chosen as a large quantile of the standard unit Fréchet distribution. In a data example it should be chosen from a set $Q$ of large empirical quantiles of $\{\eta(\boldsymbol{s},t):(\boldsymbol{s},t)\in{\mathcal{D}}_{n}\}$ for which the empirical extremograms $\widehat{\rho}_{AB,q}(\boldsymbol{h},u)$, are robust. For a practical guideline see Davis and Mikosch (2009), Section 3.4 and the upper left panel of their Figure 1, and also Davis et al. (2013c) after their Theorem 2.1. In the following simulation scenarios we choose the lowest quantile of a given level of the sets $Q$. Note that due to the variability of the large empirical quantiles, this might involve (as below) the choice of different quantiles in different data examples.

In order to test the small sample performance of the GLSE $\widehat{\boldsymbol{\theta}}_{n,V}$ defined in (4.4), we consider some of the models (I)-(III) for the dependence function $\delta_{\boldsymbol{\theta}}$. For the simulations we use the R-package RandomFields (Schlather ) and the exact method via extremal functions proposed in Dombry et al. (2016), Section 2. In this simulation study we use standardised univariate margins. If this in not the case (as for instance in the data example treated in Section 5 of Buhl and Klüppelberg (2016)), they need to be estimated and standardised first, which naturally might lead to inferior estimation results.

**(i) Spatially isotropic fractional space-time model

**We generate 100 realisations from the model (5.14) on a grid of size 15x15x300. This corresponds to the situation of a fixed spatial and an increasing temporal observation area; i.e., it is given by ${\mathcal{D}}_{n}={\mathcal{F}}\times{\mathcal{I}}_{n}$ with ${\mathcal{F}}=\{1,\ldots,15\}^{2}$ and ${\mathcal{I}}_{n}=\{1,\ldots,300\}$. We simulate the model with the true parameter vector

[TABLE]

which we assume to lie in a compact subset of

[TABLE]

As the large empirical quantile $q$ we take the $96\%$-quantile of $\{\eta(\boldsymbol{s},t):(\boldsymbol{s},t)\in{\mathcal{D}}_{n}\}$.

**(ii) Geometrically anisotropic fractional space-time model

**We generate 100 realisations from model (5.15) on a grid of size 15x15x300. This corresponds to the same situation as in (i). We simulate the model with the true parameter vector

[TABLE]

which we assume to lie in a compact subset of

[TABLE]

where we choose $\alpha_{1}\geq 1$ to ensure differentiability of the model, cf. the discussion in Section 5.3. As the large empirical quantile $q$ we take the $97\%$-quantile of $\{\eta(\boldsymbol{s},t):(\boldsymbol{s},t)\in{\mathcal{D}}_{n}\}$.

**(iii) Spatially anisotropic time-shifted model

**We generate 100 realisations from model (5.19) on a grid of size 40x40x40, and consider this as a situation where the observation area increases in all dimensions; i.e., it is given by ${\mathcal{D}}_{n}={\mathcal{I}}_{n}$ with ${\mathcal{I}}_{n}=\{1,\ldots,40\}^{3}$. We simulate the model with the true parameter vector

[TABLE]

which we assume to lie in a compact subset of

[TABLE]

where we choose $\alpha_{1},\alpha_{2}\geq 1$ to ensure differentiability of the model, cf. the discussion in Section 5.3. As the large empirical quantile $q$ we take the $95\%$-quantile of $\{\eta(\boldsymbol{s},t):(\boldsymbol{s},t)\in{\mathcal{D}}_{n}\}$.   $\Box$

In all three settings we base the estimation on the set $\mathcal{H}$ of lags given by

[TABLE]

With this choice we ensure that the lag vectors vary in all three dimensions so that we obtain reliable estimates. Generally one should choose $\mathcal{H}$ such that the whole range of clear extremal dependence is covered. However, beyond that, no lags should be included for the estimation, since independence effects can introduce a bias in the least squares estimates, similarly as in pairwise likelihood estimation; cf. Buhl and Klüppelberg (2016), Section 5.3. One way to determine the range of extremal dependence are permutation tests, which are described in Buhl et al. (2017), Section 6. From those tests we know that our choice of lags satisfies this requirement for all three models.

For the weight matrix $V$ in (4.4) we propose two choices, which yield equally good results in our statistical analysis. The first choice is $V_{1}=\operatorname*{diag}\{\exp(-\|(\boldsymbol{h},u)\|^{2}):\boldsymbol{(}\boldsymbol{h},u)\in\mathcal{H}\},$ which reflects the exponential decay of the tail dependence coefficients $\rho_{(1,\infty)(1,\infty)}(\boldsymbol{h},u)$ of Brown-Resnick processes given by tail probabilities of the standard normal distribution. The second choice is to include the (possibly bias corrected) empirical extremogram estimates as in $V_{2}=\operatorname*{diag}\{\widehat{\rho}_{(1,\infty)(1,\infty),q}(\boldsymbol{h},u):(\boldsymbol{h},u)\in\mathcal{H}\}$ (provided this is a valid choice; i.e., $V_{2}$ has only positive diagonal entries). Since the so defined weight matrix is random, what follows is conditional on its realisation. It is in practice not possible to incorporate the asymptotic covariance matrix $\Pi$ of the empirical extremogram estimates $(\widehat{\rho}_{(1,\infty)(1,\infty),q}(\boldsymbol{h},u):(\boldsymbol{h},u)\in\mathcal{H})$ (cf. Remark 4) to obtain a weight matrix that is optimal in theory. As can be seen from its specification in Theorem 3.1, it contains infinite sums and is, hence, numerically hardly tractable.

Results

For each of the scenarios (i)-(iii) we report the mean, the mean absolute error (MAE), the root mean squared error (RMSE), and a relative root mean squared error (REL) of the resulting GLSEs for the 100 simulations. Exemplary for the parameter $C_{1}$, the REL is defined as

[TABLE]

where $C_{1}^{\star}$ denotes the true parameter value and $\widehat{C}_{1,j}$ the $j$th parameter estimate.

As weight matrix we choose $V_{2}=\operatorname*{diag}\{\widehat{\rho}_{(1,\infty)(1,\infty),q}(\boldsymbol{h},u):(\boldsymbol{h},u)\in\mathcal{H}\}$ defined above. The average computing time per simulation depends on the complexity of the model (i.e., the number of parameters to be estimated) and more crucially on the chosen set $\mathcal{H}$ and on the grid size. We report an average time of 14.51 seconds for scenario (i), 14.95 seconds for scenario (ii) and 14.63 seconds for scenario (iii). The estimation results are summarised in Tables 1-3. Furthermore, in Figures 2-4 we plot the parameter estimates and add $95\%$-confidence bounds found by subsampling; cf. Politis et al. (1999), Chapter 5. We use subsampling methods, since the asymptotic covariance matrix $\Pi_{V}$ specified in Theorem 4.1 contains the matrix $\Pi$ as specified in Theorem 3.1, which is, as explained above, hardly tractable. The fact that subsampling yields asymptotically valid confidence intervals for the true parameter vectors $\boldsymbol{\theta}_{i}^{\star}$ for $i=1,2,3$ can be proved analogously to the proof of Theorem 3.5 in Buhl et al. (2017) based on Corollary 5.3.4 of Politis et al. (1999). It requires mainly the existence of continuous limit distributions of $\sqrt{n^{w}/m_{n}^{d}}\|(\widehat{\boldsymbol{\theta}}_{n,V}-\boldsymbol{\theta}_{i}^{\star})\|$, which are guaranteed by Theorem 4.4, and some conditions on the $\alpha$-mixing coefficients, which can be shown similarly as those required in Theorem 3.1.

Summarising our results, we find that the GLSE estimates the model parameters accurately. Bias and variance are largest for the parameter estimates of model (ii). There are two main reasons for this. Compared to model (i), for model (ii) we estimate two more parameters based on the same observation scheme. However, one is a direction, which is non-trivial to estimate and decreases the overall quality of the estimates. For the estimation of model (iii) the observation scheme is different; in particular, there is a relatively large number of both spatial and temporal observations available. In contrast, in the setting of models (i) and (ii) only the number of temporal observations is large.

From Tables 1 and 2 we conclude that bias and REL of the spatial parameter estimates $\widehat{C}_{1}$ and $\widehat{\alpha}_{1}$ are comparable with those of the temporal parameter estimates $\widehat{C}_{2}$ and $\widehat{\alpha}_{2}$. Bias of the spatial estimates is slightly larger than bias of the temporal estimates, which might be due to the fact that only the number of temporal observations is large.

From Table 3 we read off that the RELs of the estimates $\widehat{C}_{1}$ and $\widehat{\alpha}_{1}$, which correspond to the first spatial dimension, are slightly smaller than those of $\widehat{C}_{2}$ and $\widehat{\alpha}_{2}$. A reason for this might be the choice of the lag vectors which we included in the set $\mathcal{H}$ and which show more variation with respect to the first dimension than with respect to the second.

In her PhD thesis, Steinkohl (2013) compares computing times of the commonly applied pairwise likelihood estimation with the semiparametric method described in Buhl et al. (2017), which can be regarded as a special case of the method described in this paper. She reports in Table 6.4 a reduction of computing time by about a factor 15. Furthermore, in Section 5 of Buhl et al. (2017) we show that the semiparametric methods are more robust against small deviations from the model assumptions such as measurement errors.

Further insight

(a) Influence of the choice of lags

In order to understand how the choice of lags in $\mathcal{H}$ influences computing times and the quality of the estimates, we repeat simulation scenario (i) for different sets $\mathcal{H}_{\ell}$ where $\ell=1,\ldots,5$. These are given by

[TABLE]

From Table 4 we read off roughly stable results across all choices. As to the computational burden inherent with the choice of lags we observe from Table 5 that computing times increase roughly linearly with $|\mathcal{H}|$; more precisely, computing times approximately double when $|\mathcal{H}|$ doubles. Hence, it is advisable to choose $\mathcal{H}$ such that its cardinality is minimal across a selection of valid choices.

[TABLE]

[TABLE]

(b) Effect of the sample size

We extend the simulation scenario (i) by repeating the procedure with an increased sample size. Since the number of spatial points is considered as fixed, this involves an increase of the number of time points. In a first run, the observation area is now given by ${\mathcal{D}}_{n}={\mathcal{F}}\times{\mathcal{I}}_{n}$ with ${\mathcal{F}}=\{1,\ldots,15\}^{2}$ and ${\mathcal{I}}_{n}=\{1,\ldots,500\}$; i.e., the process is observed at 500 time points (instead of 300 as before). In a second run, the time points are extended to ${\mathcal{I}}_{n}=\{1,\ldots,1000\}$. Compared to the original scenario, everything else remains unchanged; in particular, as the large quantile $q$ we choose as before the $96\%$-quantile of $\{\eta(\boldsymbol{s},t):(\boldsymbol{s},t)\in{\mathcal{D}}_{n}\}$.

With regard to the results summarised in Tables 6 and 7, we notice that there is no significant change in mean; the confidence bounds (cf. Figure 2) are too wide to support such a hypothesis. However, the RMSE and the MAE (and thus the empirical standard deviation) of the estimates decrease considerably. This is not an unexpected behaviour: since we do not change $q$, we increase the number of observed points used for the estimation of the empirical extremogram and thus decrease its variance without introducing additional bias. In theory, we expect from Theorem 4.1 and Remark 3 that an increase of the number of time points by a factor $k$ leads to a decrease of the standard deviation of the estimates by a factor $f_{k}(\beta_{1})=(1/k)^{(w-\beta_{1}d)/2}=(1/k)^{(1-3\beta_{1})/2}$ for $\beta_{1}\in(w/(5d),w/(2d))=(1/15,1/6)$, possibly after a bias correction. The extensions from 300 to 500 and that from 300 to 1000 time points correspond to $k=5/3$ and $k=10/3$, respectively. The theoretical factors $f_{k}(\cdot)$ for $k=5/3$ and $k=10/3$ therefore lie in the intervals $(0.81,0.88)$ and $(0.62,0.74)$, respectively. This behaviour should be confirmed by the empirical standard deviation and related measures. Indeed, dividing the RMSE of the individual estimates of the four parameters based on 500 and 1000 time points by the RMSE based on 300 time points, we obtain factors $0.80$, 0.82, 0.89, 0.75 (mean value $0.82$) and $0.70$, 0.65, 0.70, 0.55 (mean value $0.65$), which all lie in the corresponding theoretical intervals or are close to them. Reasons for slight deviations from theory are of course sampling variability and the fact that in practice, the sequence $m_{n}$ is, as explained above, chosen as a large empirical quantile of the observations. Our findings are visualised in Figure 5.

[TABLE]

[TABLE]

Acknowledgements

Sven Buhl acknowledges support by the Deutsche Forschungsgemeinschaft (DFG) through the TUM International Graduate School of Science and Engineering (IGSSE).

Appendix A Appendix

A.1 $\alpha$-mixing with respect to the increasing dimensions

We need the concept of $\alpha$-mixing for the process $\{X(\boldsymbol{s}):\boldsymbol{s}\in\mathbb{R}^{d}\}$ with respect to $\mathbb{R}^{w}$. In a space-time setting with fixed spatial setting and increasing time series this is called temporal $\alpha$-mixing.

Definition 6 ($\alpha$-mixing and $\alpha$-mixing coefficients)

Consider a strictly stationary process $\left\{X(\boldsymbol{s}):\boldsymbol{s}\in\mathbb{R}^{d}\right\}$ and let $\|\cdot\|$ be some norm on $\mathbb{R}^{d}$. For $\Lambda_{1},\Lambda_{2}\subset\mathbb{Z}^{w}$ define

[TABLE]

Further, for $i=1,2$ denote by $\sigma_{{\mathcal{F}}\times\Lambda_{i}}=\sigma\left\{X(\boldsymbol{s}):\boldsymbol{s}\in\mathcal{F}\times\Lambda_{i}\right\}$ the $\sigma$-algebra generated by $\{X(\boldsymbol{s}):\ \boldsymbol{s}\in\mathcal{F}\times\Lambda_{i}\}$.

1. (i)

We define the $\alpha$-mixing coefficients with respect to $\mathbb{R}^{w}$ for $k_{1},k_{2}\in\mathbb{N}$ and $z\geq 0$ as

[TABLE] 2. (ii)

We call $\{X(\boldsymbol{s}):\boldsymbol{s}\in\mathbb{R}^{d}\}$ $\alpha$-mixing with respect to $\mathbb{R}^{w}$, if $\alpha_{k_{1},k_{2}}(z)\to 0$ as $z\to\infty$ for all $k_{1},k_{2}\in\mathbb{N}$.

We have to control the dependence between vector processes $\{\boldsymbol{Y}(\boldsymbol{s})=X_{\mathcal{B}(\boldsymbol{s},{\gamma})}:\boldsymbol{s}\in\Lambda_{1}^{\prime}\}$ and $\{\boldsymbol{Y}(\boldsymbol{s})=X_{\mathcal{B}(\boldsymbol{s},{\gamma})}:\boldsymbol{s}\in\Lambda_{2}^{\prime}\}$ for subsets $\Lambda_{i}^{\prime}\subset\mathbb{Z}^{w}$ with cardinalities $|\Lambda_{1}^{\prime}|\leq k_{1}$ and $|\Lambda_{2}^{\prime}|\leq k_{2}.$. This entails dealing with unions of balls $\Lambda_{i}=\cup_{\boldsymbol{s}\in\mathcal{F}\times\Lambda_{i}^{\prime}}\mathcal{B}(\boldsymbol{s},{\gamma})$. Since ${\gamma}>0$ is some predetermined finite constant independent of $n$, we keep notation simple by redefining the $\alpha$-mixing coefficients corresponding to the vector processes for $k_{1},k_{2}\in\mathbb{N}$ and $z\geq 0$ as

[TABLE]

A.2 Proof of Theorem 3.1

The proof of Theorem 3.1 is divided into two parts. In the first part we prove a LLN and a CLT in Lemmas A.1 and A.2 for the estimators $\widehat{\mu}_{\mathcal{B}(\boldsymbol{0},{\gamma}),m_{n}}$ in (3.4). In the second part of the proof we derive the CLT for the empirical extremogram $\widehat{\rho}_{AB,m_{n}}$ in (3.2), and compute the asymptotic covariance matrix $\Pi$. The proof generalizes corresponding proofs in Buhl and Klüppelberg (2017) (where the observation area increases in all dimensions) in a non-trivial way. We recall the separation of every point and every lag in its components corresponding to the fixed domain, indicated by the sub index ${\mathcal{F}}$, and the remaining components, indicated by ${\mathcal{I}}$, from Assumption 2. In particular, we decompose $\boldsymbol{h}^{(i)}=(\boldsymbol{h}_{\mathcal{F}}^{(i)},\boldsymbol{h}_{{\mathcal{I}}}^{(i)})\in{\mathcal{H}}$.

The separation of the observation space with its fixed domain has to be introduced into the proofs given in Buhl and Klüppelberg (2017), which is even in the regular grid situation highly non-trivial. We will give detailed references to those proofs, whenever possible, to support the understanding. On the other hand, if arguments just follow a previous proof line by line we avoid the details.

**Part I: ** LLN and CLT for $\widehat{\mu}_{\mathcal{B}(\boldsymbol{0},{\gamma}),m_{n}}$

As in Buhl and Klüppelberg (2017), Section 5, we make use of a large/small block argument. For simplicity we assume that $n^{w}/m_{n}^{d}$ is an integer and subdivide $\mathcal{D}_{n}$ into $n^{w}/m_{n}^{d}$ non-overlapping $d$-dimensional large blocks $\mathcal{F}\times\mathcal{B}_{i}$ for $i=1,\ldots,n^{w}/m_{n}^{d}$, where the $\mathcal{B}_{i}$ are $w$-dimensional cubes with side lengths $m_{n}^{d/w}$. From those large blocks we then cut off smaller blocks, which consist of the first $r_{n}$ elements in each of the $w$ increasing dimensions. The large blocks are then separated (by these small blocks) with at least the distance $r_{n}$ in all $w$ increasing dimensions and shown to be asymptotically independent.

We divide the lags in $L_{n}$ into different sets according to the large and small blocks. Recall the notation of (3.5) and around. Observe that a lag $(\boldsymbol{\ell}_{\mathcal{F}},\boldsymbol{\ell}_{\mathcal{I}})$ with $\boldsymbol{\ell}_{\mathcal{I}}=(\ell_{\mathcal{I}}^{(1)},\ldots,\ell_{\mathcal{I}}^{(w)})$ appears in $L_{\mathcal{F}}^{(i,i)}\times L_{n}$ exactly $\textnormal{N}_{\mathcal{F}}^{(i,i)}(\boldsymbol{\ell}_{\mathcal{F}})\prod_{j=1}^{w}(n-|\ell_{\mathcal{I}}^{(j)}|)$ times, where $\textnormal{N}_{\mathcal{F}}^{(i,i)}(\boldsymbol{\ell}_{\mathcal{F}})$ is defined in (3.6). This term will replace $\prod_{j=1}^{d}(n-|h_{j}|)$ in the proofs of Buhl and Klüppelberg (2017).

Lemma A.1

Let $\{X(\boldsymbol{s}):\boldsymbol{s}\in\mathbb{R}^{d}\}$ be a strictly stationary regularly varying process observed on ${\mathcal{D}}_{n}=\mathcal{F}\times\mathcal{I}_{n}$ as in (2.4). For $i\in\{1,\ldots,p\}$, let $\boldsymbol{h}^{(i)}=(\boldsymbol{h}_{\mathcal{F}}^{(i)},\boldsymbol{h}_{\mathcal{I}}^{(i)})\in\mathcal{H}\subseteq\mathcal{B}(\boldsymbol{0},{\gamma})$ for some ${\gamma}>0$ be a fixed lag vector and use as before the convention that $(\boldsymbol{h}_{\mathcal{F}}^{(p+1)},\boldsymbol{h}_{\mathcal{I}}^{(p+1)})=\boldsymbol{0}$. Suppose that the following mixing conditions are satisfied.

1. (1)

$\{X(\boldsymbol{s}):\boldsymbol{s}\in\mathbb{R}^{d}\}$* is $\alpha$-mixing with respect to $\mathbb{R}^{w}$ with mixing coefficients $\alpha_{k_{1},k_{2}}(\cdot)$ defined in (A.1).* 2. (2)

There exist sequences $m:=m_{n},r:=r_{n}\to\infty$ with $m_{n}^{d}/n^{w}\to 0$ and $r_{n}^{w}/m_{n}^{d}\to 0$ as $n\to\infty$ such that (M3) and (M4i) hold.

Then for every fixed $i=1,\ldots,p+1$, as $n\to\infty$,

[TABLE]

with $\sigma_{\mathcal{B}(\boldsymbol{0},\gamma)}^{2}(D_{i})$ specified in (3.8). If $\mu_{\mathcal{B}(\boldsymbol{0},{\gamma})}(D_{i})=0$, then (A.4) is interpreted as {\mathbb{V}ar}\big{[}\widehat{\mu}_{\mathcal{B}(\boldsymbol{0},\gamma),m_{n}}(D_{i})\big{]}=o(m_{n}^{d}/n^{w}). In particular,

[TABLE]

Proof (Proof of Lemma A.1.)

We suppress the superscript $(i)$ of $\boldsymbol{h}^{(i)}$ (respectively $\boldsymbol{h}_{\mathcal{F}}^{(i)}$) for notational ease. Strict stationarity and relation (2.5) imply that

[TABLE]

As to the asymptotic variance, we start from (3), where it has been calculated that

[TABLE]

By (2.5) and since $\mathbb{P}(\boldsymbol{Y}(\boldsymbol{0})/a_{m}\in D_{i})\to 0$,

[TABLE]

Counting the lags as explained above this proof, for fixed $k\in\mathbb{N}$ we have by stationarity the analogy of (5.6) in Buhl and Klüppelberg (2017)

[TABLE]

Concerning $A_{21}$ we have,

[TABLE]

With (2.5) and (2.6) we obtain by dominated convergence,

[TABLE]

As to $A_{22}$, observe that for all $n\geq 0$ we have $\prod\limits_{j=1}^{w}(1-\frac{|\ell_{\mathcal{I}}^{(j)}|}{n})\leq 1$ for $\boldsymbol{\ell}_{\mathcal{I}}\in L_{n}$. Furthermore, since $D_{i}$ is bounded away from $\boldsymbol{0}$, there exists $\epsilon>0$ such that $D_{i}\subset\{\boldsymbol{x}\in\overline{\mathbb{R}}^{|\mathcal{B}(\boldsymbol{0},{\gamma})|}:\|\boldsymbol{x}\|>\epsilon\}.$ Hence, we obtain

[TABLE]

which differs from the corresponding expression in Buhl and Klüppelberg (2017) only by finite factors. Thus by an obvious modification of the arguments in that paper it follows that, using $r_{n}^{w}/m_{n}^{d}\to 0$ and condition (M3),

[TABLE]

Using the definition (A.1) of $\alpha$-mixing for $A_{1}=\{\boldsymbol{Y}(\boldsymbol{0})/a_{m}\in D_{i}\}$ and $A_{2}=\{\boldsymbol{Y}(\boldsymbol{\ell}_{\mathcal{F}},\boldsymbol{\ell}_{\mathcal{I}})/a_{m}\in D_{i}\}$, we obtain by (M4i),

[TABLE]

Summarising these computations, we conclude from (A.7) and (A.8) that for $n\to\infty$,

[TABLE]

and, therefore, (A.6) implies (A.4). Since ${m_{n}^{d}/n^{w}\to 0}$ as $n\to\infty$, equations (A.3) and (A.4) imply (A.5). $\square$

Lemma A.2

Let $\{X(\boldsymbol{s}):\boldsymbol{s}\in\mathbb{R}^{d}\}$ be a strictly stationary regularly varying process observed on $\mathcal{D}_{n}=\mathcal{F}\times\mathcal{I}_{n}.$ For $i\in\{1,\ldots,p\}$, let $\boldsymbol{h}^{(i)}=(\boldsymbol{h}_{\mathcal{F}}^{(i)},\boldsymbol{h}_{\mathcal{I}}^{(i)})\in\mathcal{H}\subseteq\mathcal{B}(\boldsymbol{0},{\gamma})$ for some ${\gamma}>0$ be a fixed lag vector and take as before the convention that $(\boldsymbol{h}_{\mathcal{F}}^{(p+1)},\boldsymbol{h}_{\mathcal{I}}^{(p+1)})=\boldsymbol{0}$. Let the assumptions of Theorem 3.1 hold. Then for every fixed $i=1,\ldots,p+1$,

[TABLE]

with $\widehat{\mu}_{\mathcal{B}(\boldsymbol{0},{\gamma}),m_{n}}(D_{i})$ as in (3.4), $\mu_{\mathcal{B}(\boldsymbol{0},{\gamma}),m_{n}}(D_{i})):=m_{n}^{d}\mathbb{P}(\boldsymbol{Y}(\boldsymbol{0})/a_{m}\in D_{i})$ and $\sigma_{\mathcal{B}(\boldsymbol{0},\gamma)}^{2}(D_{i})$ given in (3.8).

Proof

Again we suppress the superscript $(i)$ of $\boldsymbol{h}^{(i)}$ and $\boldsymbol{h}_{\mathcal{F}}^{(i)}$. As for the proof of consistency above, we generalise the proof of the CLT in Buhl and Klüppelberg (2017) (based on Bolthausen (1982)) to the new setting. We consider the process

[TABLE]

observed on the $w$-dimensional regular grid $\mathcal{I}_{n}$. In analogy to (5.11) in Buhl and Klüppelberg (2017) define

[TABLE]

and note that by stationarity,

[TABLE]

The boundary condition required in Eq. (1) in Bolthausen (1982) is satisfied for the regular grid ${\mathcal{I}}_{n}$. By the same arguments as in Buhl and Klüppelberg (2017),

[TABLE]

such that $\sum_{\boldsymbol{i},\boldsymbol{i}^{\prime}\in\mathbb{Z}^{w}}{\mathbb{C}ov}[I(\boldsymbol{i}),I(\boldsymbol{i}^{\prime})]>0$. Replacing ${\mathcal{S}}_{n}$ in Buhl and Klüppelberg (2017) by $\mathcal{I}_{n}$ and $n^{d}$ by $n^{w}$, we define

[TABLE]

and obtain by the same arguments that

[TABLE]

Now note that

[TABLE]

as in (A.9), with mixing coefficients defined in (A.1). Therefore,

[TABLE]

The standardized quantities are again as in Buhl and Klüppelberg (2017), with ${\mathcal{S}}_{n}$ replaced by ${\mathcal{I}}_{n}$ and $n^{d}$ by $n^{w}$, by

[TABLE]

The proof continues in Buhl and Klüppelberg (2017), with $n^{d}$ replaced by $n^{w}$, by estimating the quantities $B_{1}$, $B_{2}$ and $B_{3}$. The estimation of $B_{1}$ follows the same lines of the proof, resulting in

[TABLE]

We use definition (A.1) of the $\alpha$-mixing coefficients for

[TABLE]

then $|\Lambda_{1}^{\prime}|,|\Lambda_{2}^{\prime}|\leq 2$ and for $d(\Lambda_{1}^{\prime},\Lambda_{2}^{\prime})$ we consider the following two cases:

1. (1)

$\|\boldsymbol{i}-\boldsymbol{j}\|\geq 3r_{n}.$ Then $2r_{n}\leq(2/3)\|\boldsymbol{i}-\boldsymbol{j}\|$ and $d(\Lambda_{1}^{\prime},\Lambda_{2}^{\prime})\geq\|\boldsymbol{i}-\boldsymbol{j}\|-2r_{n}$. Since indicator variables are bounded and $\alpha_{2,2}$ is a decreasing function,

[TABLE] 2. (2)

$\|\boldsymbol{i}-\boldsymbol{j}\|{<3r_{n}.}$ Set $z:=\min\{\|\boldsymbol{i}-\boldsymbol{j}\|,\|\boldsymbol{i}-\boldsymbol{j}^{\prime}\|,\|\boldsymbol{i}^{\prime}-\boldsymbol{j}\|,\|\boldsymbol{i}^{\prime}-\boldsymbol{j}^{\prime}\|\}$, then $d(\Lambda_{1}^{\prime},\Lambda_{2}^{\prime})\geq z$ and, hence,

[TABLE]

Therefore,

[TABLE]

The analogous argument as in Buhl and Klüppelberg (2017) yields

[TABLE]

Next, $\mathbb{E}[|B_{2}|]\to 0$ as $n\to\infty$ by the same arguments as in Buhl and Klüppelberg (2017) replacing ${\mathcal{S}}_{n}$ by ${\mathcal{I}}_{n}$ and $n^{d}$ by $n^{w}$. Then we find for $B_{3}$ with the same replacements

[TABLE]

We use definition (A.1) of the $\alpha$-mixing coefficients for

[TABLE]

such that $|\Lambda_{1}^{\prime}|=1$, $|\Lambda_{2}^{\prime}|\leq n^{w}$ and $d(\Lambda_{1}^{\prime},\Lambda_{2}^{\prime})>r_{n}$. Abbreviate

[TABLE]

then $I(\boldsymbol{0})$ and $\eta(r_{n})$ are measurable with respect to ${\sigma}_{\Lambda_{1}}$ and ${\sigma}_{\Lambda_{2}}$, respectively, where $\Lambda_{i}=\cup_{\boldsymbol{s}\in\mathcal{F}\times\Lambda_{i}^{\prime}}\mathcal{B}(\boldsymbol{s},{\gamma})$. Now we apply Theorem 17.2.1 of Ibragimov and Linnik to obtain

[TABLE]

where convergence to 0 is guaranteed by condition (M4iii).

**Part II: ** CLT for $\widehat{\rho}_{AB,m_{n}}$ and limit covariance matrix

Recall the definition of $\mathcal{H}=\{\boldsymbol{h}^{(1)},\ldots,\boldsymbol{h}^{(p)}\}$. For $i\in\{1,\ldots,p\}$, write $\boldsymbol{h}^{(i)}=(\boldsymbol{h}_{\mathcal{F}}^{(i)},\boldsymbol{h}_{\mathcal{I}}^{(i)})$ with respect to the fixed and increasing domains $\mathcal{F}$ and $\mathcal{I}_{n}$. Write further $\boldsymbol{h}_{\mathcal{F}}^{(i)}=(h_{\mathcal{F}}^{(i,1)},\ldots,h_{\mathcal{F}}^{(i,q)})$ and $\boldsymbol{h}_{\mathcal{I}}^{(i)}=(h_{\mathcal{I}}^{(i,1)},\ldots,h_{\mathcal{I}}^{(i,w)})$. Now we define the ratio

[TABLE]

and the corresponding empirical estimator

[TABLE]

using that ${\mathcal{F}}(\boldsymbol{0})=\mathcal{F}.$ Observe that

[TABLE]

Then the empirical extremogram as defined in (3.2) for $\mu$-continuous Borel sets $A,B$ in $\overline{\mathbb{R}}\backslash\{0\}$ satisfies as $n\to\infty$,

[TABLE]

by definition (2.7) of the sets $D_{i}$ for $i=1,\ldots,p$. The remaining proof follows exactly as that of Theorem 4.2 in Buhl and Klüppelberg (2017), where in the last part the decomposition into a fixed and increasing grid has to be taken into account.   $\Box$

A.3 Proof of Theorem 3.3

Throughout this proof, we suppress the sub index $m_{n}$ of $\widehat{{\rho}}_{AB,m_{n}}$ and $\widehat{{\rho}}_{AB,m_{n}}$ for notational ease. The case, where $n^{w}/m_{n}^{3d}\to 0$ as $n\to\infty$, is covered by Theorem 3.2, so we assume that $n^{w}/m_{n}^{3d}\not\to 0$. Hence, by definition (3.16) we have to consider

[TABLE]

Observe that for $\boldsymbol{h}\in\mathcal{H}=\{\boldsymbol{h}^{(1)},\ldots,\boldsymbol{h}^{(p)}\}$, as $n\rightarrow\infty,$

[TABLE]

Since the conditions of Theorem 3.1 are satisfied we have that

[TABLE]

and thus, by the continuous mapping theorem, it remains to show that for $\boldsymbol{h}\in{\mathcal{H}}$,

[TABLE]

We rewrite the latter as

[TABLE]

As to $A_{1}$, we calculate

[TABLE]

By Theorem 3.1, the first term converges weakly to a normal distribution. Since $\widehat{\rho}_{AB}(\boldsymbol{h})\stackrel{{\scriptstyle P}}{{\rightarrow}}\rho_{AB}(\boldsymbol{h})$ and $\rho_{AB,m_{n}}(\boldsymbol{h})\to\rho_{AB}(\boldsymbol{h})$ as $n\to\infty$, the second term converges to $1$ in probability. Slutzky’s theorem hence yields that $A_{1}\stackrel{{\scriptstyle P}}{{\rightarrow}}0$. As to $A_{2}$, observe that

[TABLE]

Therefore $A_{2}$ converges to [math] if and only if $\sqrt{n^{w}/m_{n}^{3d}}m_{n}^{-d}=\sqrt{n^{w}/m_{n}^{5d}}$ converges to [math].   $\Box$

A.4 Proof of Theorem 4.1

We start with the proof of consistency and use a subsequence argument. Let $n^{\prime}=n^{\prime}(n)$ be some arbitrary subsequence of $n$. We show that there exists a further subsequence $n^{\prime\prime}=n^{\prime\prime}(n^{\prime})$ such that $\widehat{\boldsymbol{\theta}}_{n^{\prime\prime},V}\stackrel{{\scriptstyle\rm a.s.}}{{\rightarrow}}\boldsymbol{\theta}^{\star}$ as $n\rightarrow\infty$, which in turn implies (4.6).

By (G1) we have for $i=1,\ldots,p$ that $\widehat{\rho}_{AB,m_{n}}(\boldsymbol{h}^{(i)})\stackrel{{\scriptstyle P}}{{\rightarrow}}\rho_{AB,\boldsymbol{\theta}^{\star}}(\boldsymbol{h}^{(i)})$ as $n\rightarrow\infty.$ Hence, there exists a subsequence $n^{\prime\prime}$ of $n^{\prime}$ such that

[TABLE]

as $n\rightarrow\infty.$ For $\boldsymbol{\theta}\in\Theta$, we define the column vector and the quadratic forms

[TABLE]

where we recall from (4.3) that \widehat{\boldsymbol{g}}_{n}(\boldsymbol{\theta})=\big{[}\widehat{\rho}_{AB,m_{n}}(\boldsymbol{h}^{(i)})-\rho_{AB,\boldsymbol{\theta}}(\boldsymbol{h}^{(i)})\big{]}{}^{{}^{\intercal}}_{i=1,\ldots,p}. Assumptions (G1) and (G3) imply that $Q(\boldsymbol{\theta})>0$ for $\boldsymbol{\theta}^{\star}\neq\boldsymbol{\theta}\in\Theta$ and that $Q(\boldsymbol{\theta}^{\star})=0$, so $\boldsymbol{\theta}^{\star}$ is the unique minimizer of $Q$. Smoothness and continuity of the functions $\rho_{AB,\boldsymbol{\theta}}(\boldsymbol{h}^{(i)})$ and $V(\boldsymbol{\theta})$ (Assumptions (G4) and (G5) with $z_{1}=z_{2}=0$) and (A.16) yield

[TABLE]

Now assume that there exists some $\omega\in\Omega$ such that (A.17) holds, but $\widehat{\boldsymbol{\theta}}_{n^{\prime\prime},V}(\omega)\not\rightarrow\boldsymbol{\theta}^{\star}$. Then there exist $\epsilon>0$ and a subsequence $n^{\prime\prime\prime}=n^{\prime\prime\prime}(n^{\prime\prime})$ such that for all $n\geq 1$,

[TABLE]

Thus,

[TABLE]

for all $n\geq n_{0}$ for some $n_{0}\geq 1.$ But this contradicts the definition of $\widehat{\boldsymbol{\theta}}_{n^{\prime\prime\prime},V}$ as the minimizer of $\widehat{Q}_{n^{\prime\prime\prime}}(\boldsymbol{\theta}),$ $\boldsymbol{\theta}\in\Theta.$ Hence $\widehat{\boldsymbol{\theta}}_{n^{\prime\prime},V}\stackrel{{\scriptstyle\rm a.s.}}{{\rightarrow}}\boldsymbol{\theta}^{\star}$ as $n\rightarrow\infty$ and this shows (4.6).

To prove the CLT (4.7), we introduce the following notation:

1. $\bullet$

We set $\rho_{AB,\boldsymbol{\theta}}^{(\ell)}(\boldsymbol{h}^{(i)}):=\frac{\partial}{\partial\theta_{\ell}}\rho_{AB,\boldsymbol{\theta}}(\boldsymbol{h}^{(i)})$ for $1\leq i\leq p,1\leq\ell\leq k$ and 2. $\bullet$

$\boldsymbol{\rho}_{AB}^{(\ell)}(\boldsymbol{\theta}):=(\rho_{AB,\boldsymbol{\theta}}^{(\ell)}(\boldsymbol{h}^{(i)}):i=1,\ldots,p){}^{{}^{\intercal}}$ for $1\leq\ell\leq k.$ The Jacobian matrix $\mathrm{P}_{AB}(\boldsymbol{\theta})$ (4.5) can then be written as

[TABLE] 3. $\bullet$

We denote by $\boldsymbol{e}_{\ell}\in\mathbb{R}^{k}$ the $\ell$th unit vector. 4. $\bullet$

For $1\leq i,j\leq p$, let $v_{ij}(\boldsymbol{\theta}):=(V(\boldsymbol{\theta}))_{ij}$ be the entry in the $i$th row and $j$th column of $V(\boldsymbol{\theta})$. 5. $\bullet$

Set $v_{ij}^{(\ell)}(\boldsymbol{\theta}):=\frac{\partial}{\partial\theta_{\ell}}v_{ij}(\boldsymbol{\theta})$ and $V^{(\ell)}(\boldsymbol{\theta}):=(v_{ij}^{(\ell)}(\boldsymbol{\theta}))_{1\leq i,j\leq p},\quad 1\leq\ell\leq k.$

As $\widehat{\boldsymbol{\theta}}_{n,V}$ minimizes $\widehat{\boldsymbol{g}}_{n}(\boldsymbol{\theta}){}^{{}^{\intercal}}V(\boldsymbol{\theta})\widehat{\boldsymbol{g}}_{n}(\boldsymbol{\theta})$ w.r.t. $\boldsymbol{\theta}$, we obtain for $1\leq\ell\leq k$,

[TABLE]

Now define the $p\times k$-matrix $\widehat{\mathrm{P}}_{AB,n}:=\int_{0}^{1}\mathrm{P}_{AB}(u\boldsymbol{\theta}^{\star}+(1-u)\widehat{\boldsymbol{\theta}}_{n,V})\,\mathrm{d}u$, where the integral is taken componentwise. Assumptions (G4) and (G5) with $z_{1}=z_{2}=1$ allow for a multivariate Taylor expansion of order 0 with integral remainder term of $\widehat{\boldsymbol{g}}_{n}(\widehat{\boldsymbol{\theta}}_{n,V})$ around the true parameter vector $\boldsymbol{\theta}^{\star}$, which yields

[TABLE]

Plugging this into (A.18) and rearranging terms, we find

[TABLE]

for $1\leq\ell\leq k.$ Defining $\widehat{R}_{n,V}$ as the $k\times k$-matrix whose $\ell$th row is given by

[TABLE]

the system of equations (A.19) can be written in compact matrix form as

[TABLE]

Hence, multiplying (A.20) by $\sqrt{n^{w}/m_{n}^{d}}$ and rearranging terms, we have,

[TABLE]

Observe that the smoothness conditions (G4) and (G5) and the rank condition (G6) ensure invertibility of the terms in curly brackets and boundedness of its inverse. For the remainder of the proof, we can hence use Slutsky’s theorem; to this end note that, as $n\rightarrow\infty$:

• •

By conditions (G4) and (G5ii) with $z_{1}=z_{2}=1$, the matrices $V(\boldsymbol{\theta})$ and $\mathrm{P}_{AB}(\boldsymbol{\theta})$ are continuous in $\boldsymbol{\theta}$, hence $V(\widehat{\boldsymbol{\theta}}_{n,V})\stackrel{{\scriptstyle P}}{{\rightarrow}}V(\boldsymbol{\theta}^{\star})$ and $\mathrm{P}_{AB}(\widehat{\boldsymbol{\theta}}_{n,V})\stackrel{{\scriptstyle P}}{{\rightarrow}}\mathrm{P}_{AB}(\boldsymbol{\theta}^{\star})$ by continuous mapping.

• •

Using (4.6), we find that $(\widehat{\boldsymbol{\theta}}_{n,V}-\boldsymbol{\theta}^{\star})\stackrel{{\scriptstyle P}}{{\rightarrow}}\boldsymbol{0}$, $\widehat{R}_{n,V}\stackrel{{\scriptstyle P}}{{\rightarrow}}(\boldsymbol{0},\ldots,\boldsymbol{0})$ and $\widehat{\mathrm{P}}_{AB,n}\stackrel{{\scriptstyle P}}{{\rightarrow}}\mathrm{P}_{AB}(\boldsymbol{\theta}^{\star})$.

• •

The previous bullet point directly implies that $C\stackrel{{\scriptstyle P}}{{\rightarrow}}\boldsymbol{0}.$

• •

As to $A$, condition (G2) directly yields $\sqrt{\frac{n^{w}}{m_{n}^{d}}}\widehat{\boldsymbol{g}}_{n}(\boldsymbol{\theta}^{\star})\stackrel{{\scriptstyle\mathcal{D}}}{{\rightarrow}}\mathcal{N}(\boldsymbol{0},\Pi)$.

• •

Furthermore, $\widehat{\boldsymbol{g}}_{n}(\boldsymbol{\theta}^{\star})\stackrel{{\scriptstyle P}}{{\rightarrow}}\boldsymbol{0}$ by (G1) and therefore $B\stackrel{{\scriptstyle P}}{{\rightarrow}}\boldsymbol{0}$.

Finally, summarising those results, with B(\boldsymbol{\theta}^{\star})=\big{(}\mathrm{P}_{AB}(\boldsymbol{\theta}^{\star}){}^{{}^{\intercal}}[V(\boldsymbol{\theta}^{\star})+V(\boldsymbol{\theta}^{\star}){}^{{}^{\intercal}}]\mathrm{P}_{AB}(\boldsymbol{\theta}^{\star})\big{)}^{-1}, we obtain (4.7).   $\Box$