{'fraction_non_alphanumeric': 0.04811365007036119, 'fraction_numerical': 0.01839442471352945, 'mean_word_length': 4.439584472389284, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 16, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Full likelihood inference for max-stable data

We show how to perform full likelihood inference for max-stable multivariate distributions or processes based on a stochastic Expectation-Maximisation algorithm, which combines statistical and computational efficiency in high-dimensions. The good performance of this methodology is demonstrated by simulation based on the popular logistic and Brown--Resnick models, and it is shown to provide dramatic computational time improvements with respect to a direct computation of the likelihood. Strategies to further reduce the computational burden are also discussed.


Introduction
Under mild conditions, max-stable distributions and processes are useful for studying highdimensional extreme events recorded in space and time (Padoan et al., 2010;Davis et al., 2013;Davison et al., 2013;de Carvalho and Davison, 2014;Huser and Davison, 2014;. This broad but constrained class of models may, at least theoretically, be used to extrapolate into the joint tail, hence providing a justified framework for risk assessment of extreme events. The probabilistic justification is that the max-stable property arises in limiting models for suitably renormalised maxima of independent and identically distributed processes; see, e.g., Davison et al. (2012), Davison and Huser (2015) and Davison et al. (2018).
Because extremes are rare by definition, it is crucial for reliable estimation and prediction to extract as much information from the data as possible. Thus, efficient estimators play a particularly important role in statistics of extremes, and the maximum likelihood estimator is a natural choice thanks to its appealing large-sample properties. However, the likelihood function is excessively difficult to compute for high-dimensional data following a max-stable distribution. As detailed in §3, likelihood evaluations require the computation of a sum indexed by all elements of a given set P D , the cardinality of which grows more than exponentially with the dimension, D. In a thorough simulation study, Castruccio et al. (2016) stated that current technologies are limiting full likelihood inference to dimension 12 or 13, and they concluded that without meaningful methodological advances, a direct full likelihood approach will not be feasible.
To circumvent this computational bottleneck, several strategies have been advocated. Padoan et al. (2010) proposed a pairwise likelihood approach, combining the bivariate densities of carefully chosen pairs of observations. Although this method is computationally attractive and inherits many good properties from the maximum likelihood estimator, it also entails a loss in efficiency, which becomes more apparent in high dimensions . More efficient triplewise and higher-order composite likelihoods were inves-tigated by Genton et al. (2011), Huser and Davison (2013), Sang and Genton (2014) and Castruccio et al. (2016). However, they are still not fully efficient, and it is not clear how to optimally select the composite likelihood terms. Furthermore, because composite likelihoods are generally not valid likelihoods (Varin et al., 2011), the classical likelihood theory cannot be blindly applied for uncertainty assessment, testing, model validation and selection, and so forth.
Alternatively, Stephenson and Tawn (2005) suggested augmenting the componentwise block maxima data z n = (z n 1 , . . . , z n D ) T , where n is the block size, with their occurrence times.
This extra information may be summarized by an observed partition π n of the set {1, . . . , D}, which indicates whether or not these maxima occurred simultaneously. Essentially, the Stephenson-Tawn likelihood corresponds the limiting joint "density" of z n and π n , as n → ∞, and it yields drastic simplifications and improved efficiency. However, Wadsworth (2015) and  noted that this approach may be severely biased for finite n, especially in low-dependence scenarios. By fixing the limit partition, π, to the observed one, π n , a strong constraint is imposed, creating model misspecification, to which likelihood methods are sensitive.
In this paper, to mitigate the sub-asymptotic bias due to fixing the partition, we suggest returning to the original likelihood formulation, which integrates out the partition rather than treating it as known. By interpreting the limit partition π as missing data, we show how to design a stochastic Expectation-Maximisation algorithm (Dempster et al., 1977;Nielsen, 2000) for efficient inference. The quality of the stochastic approximation to the full likelihood can be controlled and set to any arbitrary precision at a computational cost. We show that higher-dimensional max-stable models may be fitted in reasonable time. Importantly, our method is based solely on max-stable data and does not require extra information about the partition or the original processes, unlike the Stephenson-Tawn or related threshold-based methods. Our approach is based on the algorithm of  for conditional simulation of the partition given the data, and it is closely related to the recent papers of Thibaud et al. (2016) and Dombry et al. (2017a), who in a Bayesian setting developed a Markov chain Monte Carlo algorithm for max-stable processes by treating the partition as a latent variable to be resampled at each iteration.
2 Max-stable processes and distributions 2.1 Definition, construction, and models Consider a sequence of independent and identically distributed processes Y 1 (s), Y 2 (s), . . . , indexed by spatial site s ∈ S ⊂ R d , and assume that there exist sequences of functions a n (s) > 0 and b n (s), such that the renormalised pointwise block maximum process (with block size n) converges in the sense of finite-distributional distributions, as n → ∞, to a process Z(s) with non-degenerate marginal distributions, i.e., Y (s) is in the max-domain of attraction of Z(s).
Then, the limit Z(s) is max-stable (see, e.g., de Haan and Ferreira, 2006, Chap. 9). That is, pointwise maxima of independent copies of the limit process Z(s) have the same dependence structure as Z(s) itself, while marginal distributions are in the same location-scale family and coincide with the generalized extreme-value distribution.

Underlying partition and extremal functions
At each site s ∈ S, the pointwise supremum, Z(s), in (2) is realised by a single profile W j (s)/P j almost surely. Such profiles are called extremal functions in Dombry et al. on the left panel of Fig. 1 indicates that the max-stable process at these five sites came from three separate independent events; in particular, the maxima at s 3 = 15, s 4 = 20, and s 5 = 25 were generated from the same profile.
Similarly, an observed partition π n of {1, . . . , D} may be defined for Z n (s) in (1) based on the original processes Y 1 (s), . . . , Y n (s). The knowledge of π n tells us if extreme events at different sites occurred simultaneously or not, so π n carries information about the strength of spatial extremal dependence. If the processes Y 1 (s), . . . , Y n (s) (suitably marginally transformed) are in the max-domain of attraction of Z(s) in (2), then the partition π n converges in distribution to π, as n → ∞, on the space of all partitions P D of {1, . . . , D} (Stephenson and Tawn, 2005).
We now describe likelihood inference for max-stable vectors: by exploiting the information on the partition π (Stephenson-Tawn likelihood), or by integrating it out (full likelihood).

Full and Stephenson-Tawn likelihoods
By differentiating the distribution (3) with respect to the variables z 1 , . . . , z D , we can deduce that the corresponding density, or full likelihood for one replicate, may be expressed as where V τ i denotes the partial derivative of the function V with respect to the variables indexed by the set τ i ⊂ {1, . . . , D} Castruccio et al., 2016). The sum in (4) is taken over all elements of P D , the size of which equals the Bell number of order D, leading to an explosion of terms, even for moderate D. Each term in (4) corresponds to a different configuration of the profiles W j (s)/P j in (2) at the sites s 1 , . . . , s D . Thus, Castruccio et al. (2016) argued that the computation of (4) is limited to dimension 12 or 13 with modern computational resources.
As detailed in the Supplementary Material using a point process argument, and originally shown by Stephenson and Tawn (2005), the joint density of the max-stable data z = (z 1 , . . . , z D ) T and the associated partition π = {τ 1 , . . . , τ |π| } ∈ P D is simply equal to hence reducing the problematic sum to a single term, making likelihood inference possible in higher dimensions and simultaneously improving statistical efficiency. Because the asymptotic partition π is not observed, Stephenson and Tawn (2005) suggested replacing it by the observed partition π n of occurrence times of maxima, which converges to π provided the asymptotic model is well specified. However, Wadsworth (2015) and  showed that lack of convergence of π n to π may result in severe estimation bias, which is especially strong in low-dependence cases. To circumvent this problem, Wadsworth (2015) proposed a bias-corrected likelihood; alternatively, we show in the next section how to design a stochastic Expectation-Maximisation algorithm to maximise (4), while taking advantage of the computationally appealing nature of (5).

Stochastic Expectation-Maximisation algorithm
It is instructive to rewrite the full likelihood (4) using (5) as g Full (z 1 , . . . , z D ) = π∈P D g ST (z 1 , . . . , z D , π) because it highlights that the full likelihood simply integrates out the latent random partition π needed for the Stephenson-Tawn likelihood. Interpreting π as a missing observation and the Stephenson-Tawn likelihood as the completed likelihood, an Expectation-Maximisation algorithm (Dempster et al., 1977) may be easily formulated. Assume that the exponent function V (z 1 , . . . , z D | θ) is parametrised by a vector θ ∈ Θ ⊂ R p . Starting from an initial guess θ 0 ∈ Θ, the Expectation-Maximisation algorithm consists of iterating the following Eand M-steps for r = 1, . . . , R: • E-step: compute the functional where the expectation is computed with respect to the discrete conditional distribution of π given the data z = (z 1 , . . . , z D ) T and the current value of the parameter θ r−1 , i.e., g(π | z, θ r−1 ) = g ST (z, π | θ r−1 )/g Full (z | θ r−1 ).
• M-step: update the parameter as θ r = arg max θ∈Θ Q(θ, θ r−1 ). Dempster et al. (1977) showed that the Expectation-Maximisation algorithm has appealing properties; in particular, the value of the log-likelihood increases at each iteration, which ensures convergence of θ r to a local maximum, as r → ∞. In our case, however, the expectation in (6) is tricky to compute: it contains again the sum over the set P D , and (7) relies on the full density g Full (z | θ r−1 ), which we try to avoid. To circumvent this issue, one solution is to approximate (6) by Monte Carlo as where the partitions π 1 , . . . , π N are conditionally independent at best, or form an ergodic sequence at least. As g(π | z, θ r−1 ) ∝ g ST (z, π | θ r−1 ), see (7), it is possible to devise a Gibbs sampler to generate approximate simulations from g(π | z, θ r−1 ) without explicitly computing the constant factor g Full (z | θ r−1 ) in the denominator of (7). Thanks to ergodicity of the resulting Markov chain, the precision of the approximation (8) may be set arbitrarily high by letting N → ∞ (and discarding some burn-in iterations). More details about the practical implementation of the Gibbs sampler are given in  and in the Supplementary Material. Although the number of iterations of the Gibbs sampler, N , will typically be much smaller than the cardinality of P D , the approximation (8) to (6) will likely be reasonably good for moderate values of N because only a few partitions π ∈ P D may be plausible or compatible with the data z = (z 1 , . . . , z D ) T .
The asymptotic properties of the stochastic Expectation-Maximisation estimator, θ SEM , were studied in details by Nielsen (2000) and compared with the classical maximum likelihood estimator, θ; see §2-3 therein, in particular Theorem 2. Dombry et al. (2017b) showed that the maximum likelihood estimator θ is consistent and asymptotically normal for the most popular max-stable models, including the logistic and Brown-Resnick models used in this paper. This suggests that these appealing asymptotic properties should also be satisfied for the estimator θ SEM , provided some additional rather technical regularity conditions detailed in Nielsen (2000) are satisfied. If so, then the asymptotic performance of θ SEM is akin to that of θ, though with a larger asymptotic variance. Finally, the inherent variability of the stochastic Expectation-Maximisation algorithm may also be a blessing: unlike the deterministic Expectation-Maximisation algorithm, it is less likely to get stuck at a local maximum of the full likelihood.

Simulation study
To assess the performance of the stochastic Expectation-Maximisation estimator θ SEM , we simulated data from the multivariate logistic max-stable distribution with exponent func- Here, the parameter θ controls the dependence strength, with θ → 0 and θ = 1 corresponding to perfect dependence and independence, respectively. This model was chosen for two main reasons: first, it is the simplest max-stable distribution, often used as a benchmark, that interpolates between perfect dependence and independence; and second, the full likelihood (4) can be efficiently computed in this case using a recursive algorithm (Shi, 1995), thus allowing us to compare θ SEM and θ in high dimensions. Further results for the popular and more realistic, but computationally demanding, spatial max-stable model known as Brown-Resnick model (Kabluchko et al., 2009), are reported in the Supplementary Material. All simulations presented below and in the Supplementary Material were performed in parallel on the KAUST Cray XC40 supercomputer Shaheen II.
We first investigated the performance of the estimator θ SEM under different scenarios.
We considered dimensions D = 2, 5, 10, 20, with 20 independent temporal replicates, and θ = 0.1, . . . , 0.9 (strong to weak dependence). Setting the initial value to θ 0 = 0.6, we Relative error (%) Figure 2: Performance of the estimator θ SEM : bias (left), standard deviation (middle) and relative error in % (right), for the logistic model with θ = 0.1, . . . , 0.9 in dimension D = 2 (thinnest black), 5 (thin red), 10 (thick green) and 20 (thickest blue), based on 20 independent temporal replicates. The number of iterations for the Expectation-Maximisation algorithm was set to R = 30, averaging the last 5 iterations, and the number of iterations of the underlying Gibbs sampler was set to 110 × D (thinned by a factor D, after a burn-in of 10 × D). The initial value was set to θ 0 = 0.6. chose R = 30 iterations for the Expectation-Maximisation algorithm, averaging the last 5 iterations, and took 110 × D iterations for the underlying Gibbs sampler. Following simulations reported in the Supplementary Material, we discarded the first 10 × D iterations as burn-in, and thinned the Markov chain by a factor D, in order to keep N = 100 roughly independent partitions π i to compute (8). We repeated the experiment 1024 times to estimate the bias, B, standard deviation, SD, root mean squared error, RMSE = (B 2 + SD 2 ) 1/2 , and relative error with respect to the maximum likelihood estimator, RE = E|( θ SEM − θ)/ θ|. Figure 2 reports the results. As expected, the bias is negligible compared to the standard deviation, and the latter decreases with increasing dimension D but increases as the data approach independence (θ → 1). The RMSE (not shown) is almost only determined by the standard deviation. The relative error is always very small (uniformly less than about 0.6%), and it decreases for the most part with D and also with N as suggested by further unreported simulations.
We now turn our attention to the computational efficiency of the stochastic Expectation-Maximisation algorithm. Considering dimensions up to D = 100 under the exact same setting as before, the leftmost panel of Fig. 3 shows that it takes on average 5-6 minutes to compute θ SEM (on a single core with 2.3GHz) when D = 100 and θ = 0.9. Recall that, according to Castruccio et al. (2016), a direct evaluation of the likelihood (4) is not possible in dimensions greater than D = 12 or 13, thus this result is a big improvement over the current existing methods. The computational time appears to be roughly linear with D, which is due to the number of iterations of the Gibbs sampler set proportional to D.
However, for more complex models such as the Brown-Resnick model, the computational time is significantly larger; see the Supplementary Material. Therefore, it makes sense to seek strategies to reduce the computational burden. One possibility is to tune the number of iterations of the stochastic Expectation-Maximisation algorithm. To investigate its speed of convergence, the two middle panels of

Discussion
To resolve the problem of inference for max-stable distributions and processes, we have proposed a stochastic Expectation-Maximisation algorithm, which does not fix the underlying partition, but instead, treats it as a missing observation, and integrates it out. The beauty of this approach is that it combines statistical and computational efficiency in high dimensions, and it does not suffer from misspecification entailed by lack of convergence of the partition.
As a proof of concept, we have validated the methodology by simulation based on the logistic model, and we have shown that in this case it is easy to make inference beyond dimension see also Genton et al. (2018). Our method is not limited to these two models and could potentially be applied to any max-stable model for which the exponent function V and its partial derivatives V τ i are available. The main computational bottleneck of our approach is that we need to generate a Gibbs sampler for each independent temporal replicate of the process. Fortunately, as this setting is embarrassingly parallel, we may thus easily take advantage of available distributed computing resources. Finally, there is a large volume of literature on the stochastic Expectation-Maximisation algorithm, and it might be possible to devise automatic stopping criteria and adaptive schemes for the Gibbs sampler to further speed up the algorithm (Booth and Hobert, 1999). the block τ i , all the points s j , j ∈ τ i , share the same extremal function that will hence be denoted by ϕ + τ i .
• By definition of the extremal functions, one gets the joint law of the hitting scenario and max-stable observations: • By integrating out the hitting scenario, one obtains the law of the max-stable observations: Equation (11) provides an alternative formula for the Stephenson-Tawn likelihood, g ST , based on the Poisson point process intensity, λ, while Equation (12) is the max-stable full likelihood, g Full . Identifying the expressions (11) and (12) above with (5) and (4) in the main paper, respectively, we can see that This relates a partial derivative of the exponent function V with a partial integral of the point process intensity λ. In particular, (13) implies that the intensity is the mixed derivative of the exponent function with respect to all arguments, i.e., Furthermore, the function V corresponds to the integrated intensity of the set A = [0, z] c , i.e.,

A.2 Computing the Poisson point process intensity
The intensity measure λ is an important feature of max-stable models and can be computed for most popular models; see  for a derivation of λ for the Brown-Resnick model (Kabluchko et al., 2009) and Ribatet (2013) for an expression of λ for the extremal-t model (Opitz, 2013). Partial integrals of λ for these models may be found in Wadsworth and Tawn (2014) and Thibaud and Opitz (2015), respectively. Using the relations (13) and (14), the intensity λ and its partial integrals can be deduced for the Reich and Shaby (2012) model from the expressions in the appendix of Castruccio et al. (2016).
Here, as a simple pedagogical illustration for many other multivariate or spatial maxstable models, we consider the multivariate logistic model, which we used in our simulation study. In this case, the function V and its partial and full derivatives can readily be obtained by direct differentiation.
Recall that the exponent function for the logistic model is It is known that the multivariate counterpart of the spectral representation (9) for the logistic model is obtained by taking W = (W 1 , . . . , W D ) T with independent and identically distributed Fréchet(β, c β ) components, where β = 1/θ and c β = 1/Γ(1 − 1/β) are shape and scale parameters, respectively; see, for example, Proposition 6 in Dombry et al. (2016). Then, and we deduce from Equation (10) that

Similar computations entail
where, for the first equality, we used

A.3 Details on the underlying Gibbs sampler
The Gibbs sampler proposed by  is designed to draw an ergodic sequence of partitions π 1 , . . . , π N conditional on the observed max-stable data z = (z 1 , . . . , z D ) T , i.e., from the discrete distribution g(π | z, θ), where θ ∈ Θ ⊂ R p is the parameter vector characterizing the max-stable dependence structure. One has The normalizing constant in the denominator of (15) is computationally demanding to compute as it involves the sum over all partitions. Nevertheless, the Gibbs sampler of  provides a way to construct a Markov chain whose stationary distribution is g(π | z, θ), while avoiding the computation of the normalizing constant. Let π t = {τ t;1 , . . . , τ t;|πt| } ∈ P D be the partition at the tth iteration of the Gibbs sampler. The idea of the Gibbs sampler is to sample the next partition π t+1 = {τ t+1;1 , . . . , τ t+1;|π t+1 | } ∈ P D by keeping all but one component fixed. Let ∈ {1, . . . , D} be the component to be updated, and let π − t and π − t+1 denote the partitions π t and π t+1 , respectively, with the th component removed. We update the partition π t by modifying the (randomly chosen) th component using the full conditional distribution The combinatorial explosion is avoided, because the number of possible updates π t+1 such that π − t+1 = π − t is at most |π t | + 1. Moreover, as we update only one component at a time, many terms in the ratio (16) cancel out, and at most four of them need to be computed, which makes it computationally attractive. However, for the same reason, the resulting partitions will also be heavily dependent, and so, intuitively, we should take the number of Gibbs iterations to be roughly proportional to the dimension D and thin the Markov chain by a factor D to get approximately independent (or weakly dependent) partitions. A suitable burn-in should also be specified to ensure that the Markov chain has appropriately converged to its stationary distribution.
In order to assess the number of iterations required for the Gibbs sampler to converge (i.e., the burn-in), we considered the logistic model defined by its exponent function V (z 1 , . . . , z D | θ) = (z −1/θ 1 + · · · + z −1/θ D ) θ , θ ∈ (0, 1]. We generated five independent copies of a logistic random vector in dimension D = 50, and considered the cases θ = 0.9, 0.6, 0.3 (weak to strong dependence). For each dataset, we ran five Gibbs samplers (one per independent replicate) for 5000 iterations. To easily visualize the resulting Markov chains and assess convergence, we display in Fig. 4  To validate such results for another max-stable model, we did the same experiment for the Brown-Resnick model (Kabluchko et al., 2009), defined by taking W j (s) = exp{ε j (s) − γ(s)} in (9), where the terms ε j (s) are independent copies of ε(s), which is an intrinsically stationary Gaussian process with zero mean and variogram 2γ(h) = var{ε(s) − ε(s + h)} such that ε(0) = 0 almost surely. Using the exact simulation algorithm of Dombry et al.
(2016), we simulated five independent replicates of the Brown-Resnick with semi-variogram γ(h) = (h/λ) ν , where λ > 0 and ν ∈ (0, 2] are the range and smoothness parameters, respectively, at D = 10 randomly generated sites s 1 , . . . , s 10 ∈ [0, 1] 2 . We considered the cases λ = 0.5, 1, 1.5 (short to long range dependence) with ν = 1.5. For each dataset, we ran five Gibbs samplers (one per independent replicate) for 1000 iterations. Figure  These results suggest to discard the first 10 × D iterations as burn-in, and to thin the resulting Markov chains by a factor D to obtain approximately (conditionally) independent partitions. With this setting, the initial partition has negligible impact on the results.
Alternatively, another natural option could be to initialize the partition randomly from its unconditional distribution, which can be easily obtained from an unconditional simulation of the max-stable distribution. This could potentially provide further computational savings by reducing the time it takes for the Gibbs sampler to converge (thus reducing the burn-in).
In order to assess the performance of the stochastic Expectation-Maximisation estimator in each scenario of Table 1 Table 1, covering short to long range dependent processes (top to bottom) and rough to smooth processes (left to right). Realizations are displayed on standard Gumbel margins. of §A.3; and (ii) a pairwise likelihood estimator θ PAIR = ( λ PAIR , ν PAIR ) T (see, e.g., Padoan et al., 2010;Huser and Davison, 2013), which maximizes the pairwise likelihood constructed by combining the likelihood contributions from all 10 2 = 45 pairs of sites together with equal weight. We repeated this experiment 1024 times to compute performance metrics, such as the root mean squared error of parameter estimates.  Table 1 based on the Brown-Resnick model with semi-variogram γ(h) = (h/λ) ν , simulated at D = 10 random sites in [0, 1] 2 , with 10 independent replicates. Left (red) and right (blue) boxplots correspond to θ PAIR = ( λ PAIR , ν PAIR ) T and θ SEM = ( λ SEM , ν SEM ) T , respectively. Each boxplot summarizes the variability of parameter estimates based on 1024 simulations. Five estimates reaching up to log( λ PAIR ) ≈ 40 were omitted in Scenario 1 for visibility purposes. Orange horizontal segments are the true values.
Expectation-Maximisation and pairwise likelihood estimators seem to work well overall with a very low bias, although the variability is in some cases very high, due to the tricky estimation exercise with only 10 replicates in dimension D = 10. Nevertheless, the inter-quartile range appears to be quite moderate in all cases. The stochastic Expectation-Maximisation estimator appears clearly superior to the pairwise likelihood estimator in the cases we have considered, as it fully utilizes the information available in the data. To investigate this more in depth, Table 2 reports relative efficiencies of the pairwise likelihood estimator with re- We now desire to check the speed of convergence of the Expectation-Maximisation algorithm, similarly to the simulations that we did for the logistic model in the main paper. The right panel of Fig. 8 also reveals that the estimated range parameter is negatively correlated with the estimated smoothness parameter. This was expected as these two parameters have an opposing effect on the dependence strength, and it suggests that alternative orthogonal parametrizations might be preferable. We leave this problem for future research. Finally, we investigate the scalability of the stochastic Expectation-Maximisation estimator for the Brown-Resnick model when the dimension D increases. To assess this, we considered the same setting as above with 10 independent replicates, same number of Gibbs iterations, true values set to λ = ν = 1.5, and we measured the computational time needed for the first Expectation-Maximisation iteration, in dimensions D = 5, 10, 15, 20. Figure 9 plots the computational time on a logarithmic scale, along with a projection for dimensions up to D = 30. Unlike the logistic model, which has an explicit exponent function V and partial derivatives V τ i , the expressions for the Brown-Resnick model involve the multivariate Gaussian distribution in dimension up to D − 1, whose computation is very demanding for large D. This significantly slows down the algorithm, whose complexity now appears to be increasing exponentially with D (by contrast with the linear increase observed for the logistic model). More precisely, a single Expectation-Maximisation iteration takes on average approximately 1.5min, 12.7min, 52.2min, and 19.8hr for dimensions D = 5, 10, 15, 20, respectively. This precludes the use of the stochastic Expectation-Maximisation algorithm beyond dimension D = 20, unless distributed computing resources are used to run all Gibbs sampler in parallel. Nevertheless, these results already represent a big advancement towards full likelihood inference for max-stable processes, as Castruccio et al. (2016) argued that a direct calculation of the likelihood was impossible beyond dimension D = 12 or 13. As shown by de Fondeville and , speed-ups can be obtained by appropriately exploiting quasi-Monte Carlo techniques for the calculation of multivariate Gaussian distributions.
A similar computational burden is expected for the extremal-t model (Opitz, 2013), which relies on the computation of multivariate Student t distributions, but a better computational efficiency should prevail for the Reich and Shaby (2012) max-stable model, for which the expressions of the exponent function V and its partial derivatives V τ i are available in explicit form; see the appendix of Castruccio et al. (2016).