Max-stable processes are the extension of the multivariate extreme value theory to the infinite dimensional setting. More precisely consider a stochastic process , having continuous sample paths. Then the limiting process

where are independent replications of , and are sequences of continuous functions and the limiting process is assumed to be non degenerate. de Haan [1994] shows that the class of the limiting processes corresponds to that of max-stable processes and hence emphasizes on their use to model spatial extremes.
Interestingly there are two different ways of charactaresing max-stable processes: these are known as the spectral characterisations. An especially useful special case of the characterisation of de Haan [1984], is

where are the points of a Poisson point process on with intensity measure , is a probability density function on . The process defined above is a max-stable process with unit Frechet margins. Taking as the multivariate Normal density with zero mean and covariance matrix gives the Smith model [Smith, 1990] for which the bivariate distribution function is

where is the standard normal distribution function and ,
Another very useful spectral characterisation for unit Frechet max-stable processes is [Schlather, 2002]

where are the points of a Poisson point process on with intensity measure and are independent replications of a positive stochastic process having continuous sample paths such that for all .
Currently there are several useful models based on Schlather's characterisation. The first model, the Schlather's model, is to use , where is a standard Gaussian process. This leads to the bivariate distribution function

Another possibility is to take , where . This is known as the Geometric gaussian model for which the bivariate distribution is similar to the Smith model with . A last possibility, which is a generalisation of the geometric Gaussian model, is to take , where is a zero mean Gaussian process having stationary increments and (semi)variogram such that almost surely. Its bivariate distribution function is again similar to the Smith model with .

Because max-stable processes are an extension of the multivariate extreme value theory to the infinite dimensional case when trying to fit such processes to we are faced to the same problem has in the finite dimensional setting.
Suppose we have observed the process at , then the multivariate distribution is of the form

where is a function having some homogeneity property called the exponent measure. Consequently if one wants to use the maximum likelihood estimator this will yield to a combinatorial explosion and the likelihood will be intractable. Instead one can have resort to composite likelihood estimator and an especially convenient choice is to maximize the pairwise log-likelihood

where is a set of appropriate weights and is the bivariate density of a max-stable process --- for which closed forms exist. Under some regularity conditions and in particular if is identifiable from the bivariate densities, then

where and .

Well we can simulate a max-stable process using its limiting characterisation i.e. find the normalizing functions and and simulate indendant replications of a stochastic process for which the normalizing function are based. This naive simulation technique has the disadvantage of requiring a large number of independent replicates but also pose the question of how many replicates should be simulated? A better strategy consists in using one of the spectral representation. In most cases this will be more efficient since one can simulate the points of a Poisson point process with intensity can be generated as follows

Now as standard exponential random variables are positive, as . If we further assume that the stochastic process used in the spectral characterisation is bounded above, we will require only a finite number of independent replications. If the process is not bounded then Schlather [2002] shows that it is however possible to get accurate simulations.
The methodology introduced above was the first approach to generate (approximate) realizations from a max-stable process. It is now possible to get exact simulation using the methodology developped by C. Dombry, S. Engelke and M. Oesting. The package uses such a strategy whenever it is possible (or implemented ;-) )