\documentclass{article}
\usepackage[T1]{fontenc}
\usepackage[utf8]{inputenc}
\usepackage{amsmath, amssymb, amsfonts}
\usepackage{mathpazo}
\usepackage{geometry}
\usepackage{natbib}
\usepackage{hyperref}
%\geometry{left=2cm, right=2cm, top=2.5cm, bottom=2.25cm}

\newcommand{\all}{ \; \forall \;}
\newcommand{\bs}[1]{\boldsymbol {#1}}
\newcommand{\R}{\mathbb{R}}
\renewcommand{\P}[2][]{{\mathsf P}_{#1}\left(#2\right)}
\newcommand{\E}[2][]{{\mathsf E}_{#1}\left(#2\right)}
\newcommand{\Va}[2][]{{\mathsf{Var}_{#1}}\left(#2\right)}
\newcommand{\Co}[2][]{{\mathsf{Cov}_{#1}}\left(#2\right)}
\newcommand{\code}[1]{\texttt{#1}}
\newcommand{\sumi}{\sum_{i=1}^n}
\newcommand{\pfrac}[2]{\left(\frac{#1}{#2}\right)}
\DeclareMathOperator{\diag}{diag}

\begin{document}
<<setup, include=FALSE, cache=FALSE, echo=FALSE>>=
library(knitr)
## set global chunk options
opts_chunk$set(fig.path='figure/manual-', cache.path='cache/manual-', fig.align='center', fig.show='hold', par=TRUE)
## I use = but I can replace it with <-; set code/output width to be 68
options(formatR.arrow=TRUE, width=68, digits=4)
## tune details of base graphics (http://yihui.name/knitr/hooks)
knit_hooks$set(par=function(before, options, envir){
if (before && options$fig.show!='none') par(mar=c(4,4,.1,.1),cex.lab=.95,cex.axis=.9,mgp=c(2,.7,0),tcl=-.3)
})
@
%\VignetteIndexEntry{Exact unconditional sampling from max-stable random vectors}
%\VignetteEngine{knitr::knitr_notangle}

\title{Exact unconditional sampling from max-stable random vectors}
\author{Léo Belzile}
\date{}
\maketitle
\begin{center}
{ \small
EPFL-SB-MATH-STAT \\ Station 8, Bâtiment MA \\CH-1015, Switzerland,\\ \href{leo.belzile@epfl.ch}{\texttt{leo.belzile@epfl.ch}}
}
\end{center}

Originally, the `\code{mev}' package was introduced to implement the exact unconditional sampling algorithms in \cite{Dombry:2016} to simulate simple max-stable random vectors. The implementation will work efficiently for moderate dimensions and does not currently take improvements for large spatial domains.

\section{Functions and use}

There are two main functions, \code{rmev} and \code{rmevspec}. \code{rmev} samples from simple max-stable processes, meaning it will return an $n \times d$ matrix of samples, where each of the column has a sample from a unit Frechet distribution. In constrast, \code{rmevspec} returns sample on the unit simplex from the spectral (or angular) measure. One could use this to test estimation based on spectral densities, or to construct samples from Pareto processes.

The syntax is
<<randomsim, include=TRUE, cache=FALSE, echo=TRUE>>=
library(mev)
#Sample of size 1000 from a 5-dimensional logistic model
x <- rmev(n=1000, d=5, param=0.5, model="log")
#Marginal parameters are all standard Frechet, meaning GEV(1,1,1)
apply(x, 2, function(col){ismev::gev.fit(col, show=FALSE)$mle})


#Sample from the corresponding spectral density
w <- rmevspec(n=1000, d=5, param=0.5, model="log")
#All rows sum to 1 by construction
head(rowSums(w))
#The marginal mean is 1/d
round(colMeans(w),2)
@

\section{Description of the models implemented}

The different models implemented are described in \cite{Dombry:2016}, but some other models can be found and are described here.
\begin{enumerate}

\item \textbf{logistic distribution} (\code{log})
The logistic model of \cite{Gumbel:1960} (the Gumbel Archimedean copula)
\begin{align*}
   \P{\bs{X} \leq \bs{x}}= \exp \left[ - \left(\sumi \pfrac{1}{x_i}^{\alpha}\right)^{\frac{1}{\alpha}}\right]
\end{align*}
for $\alpha>1.$ By default, \code{rmev} will transform an argument in $(0,1)$ without warning, to conform with the
implementation. The spectral measure density is
\begin{align*}
   h_{\bs{W}}(\bs{w})=\frac{1}{d}\frac{\Gamma(d-\alpha)}{\Gamma(1-\alpha)}\alpha^{d-1}\left( \prod_{j=1}^d
w_j\right)^{-(\alpha+1)}\left(\sum_{j=1}^d
   w_j^{-\alpha}\right)^{1/\alpha-d}, \qquad \bs{w} \in \mathbb{S}_d
\end{align*}

\item \textbf{asymmetric logistic distribution} (\code{alog})
This model proposed by  \cite{Tawn:1990} uses the parametrization of the \code{evd} package, merely replacing the algorithm for
the generation of logistic
variates. The distribution function is
\begin{align*}
   \P{\bs{X} \leq \bs{x}} = \exp \left[ -\sum_{b \in \mathbb{B}}\left(\sum_{i \in b} \pfrac{\theta_{i,
b}}{x_i}^{\alpha_b}\right)^{\frac{1}{\alpha_b}}\right]
\end{align*}
while the spectral density corresponds to
\begin{align*}
   h_{\bs{W}}(\bs{w})=\frac{1}{d}\frac{\prod_{j=1}^{d-1} (j\alpha_b-1)}{\prod_{i\in b} w_i} \left(\prod_{i \in
b}\frac{\theta_{i,b}}{w_i}\right)^{\alpha_b} \left(\sum_{i \in b} \pfrac{\theta_{i, b}}{w_i}^{\alpha_b}\right)^{1/\alpha_b-d}
\end{align*}

The parameters $\theta_{i, b}$ must be provided in a list and represents the asymmetry parameter.
In the $d$-variate
case, where $\mathbb{B}$ is the power set of the collection $\{1, \ldots, d\}$, without the empty set. This corresponds to
generating samples $\bs{Z}_b$ from a logistic distribution of dimension $|b|$ (or Fréchet variates if $|b|=1)$ with parameter
$\alpha_b$ (possibly recycled). Then, each marginal value corresponds to the maximum of the weighted corresponding entry. That
is, $X_{i}=\max_{b \in \mathbb{B}}\theta_{i, b}Z_{i,b}$ for all $i=1, \ldots, d$. The max-mixture is valid provided that $\sum_{b
\in \mathbb{B}} \theta_{i,b}=1$ for $i=1, \ldots, d.$ As such, empirical estimates of the spectral measure will almost surely
place mass on the inside of the simplex rather than on subfaces.
\item \textbf{negative logistic distribution} (\code{neglog})
The distribution function of the min-stable distribution due to \cite{Galambos:1975}  is
\begin{align*}
   \P{\bs{X} \leq \bs{x}} = \exp \left[ -\sum_{b \in \mathbb{B}} (-1)^{|b|}\left(\sum_{i \in b}
{x_i}^{-\alpha}\right)^{\frac{1}{\alpha}}\right]
\end{align*}
for $\alpha_b\leq 0$ and corresponds to the model in \cite{Dombry:2016}.

\begin{align*}
   h_{\bs{W}}(\bs{w}) = \frac{1}{d} \sum_{b \in B: c \subset b} (-1)^{|b|}
\left(\prod_{j=1}^{d-1}(1-j\alpha)\right)\left(\prod_{i \in c}w_j\right)^{-(\alpha+1)}\left(\sum_{i \in c}
w_i^{-\alpha}\right)^{1/\alpha-d}
\end{align*}

\item \textbf{asymmetric negative logistic distribution} (\code{aneglog})
The construction of this distribution is exactly the same as for the asymmetric logistic distribution, with distribution function
\begin{align*}
      \P{\bs{X} \leq \bs{x}} = \exp \left[ -\sum_{b \in \mathbb{B}}\sum_{c \in b} (-1)^{|c|}
\left(\sum_{i \in c}
\pfrac{x_i}{\theta_{i, b}}^{\alpha_b}\right)^{-\frac{1}{\alpha_b}}\right].
\end{align*}
In particular, it does not correspond to the ``negative logistic distribution'' given in \cite{Joe:1990}, which is not a valid
distribution function according to Stephenson.
\item \textbf{multilogistic distribution} (\code{bilog})
This multivariate extension of the logistic, proposed by \cite{Boldi:2009}, places mass on the interior of the simplex. Let $\bs{W} \in \mathbb{S}_d$ be the solution of
\begin{align*}
 \frac{W_j}{W_d}=\frac{C_jU_j^{-\alpha_j}}{C_dU_d^{-\alpha_d}}, \quad j=1, \ldots, d
\end{align*}
where $C_j=\Gamma(d-\alpha_j)/\Gamma(1-\alpha_j)$ for $j=1, \ldots, d$ and $\bs{U} \in \mathbb{S}_d$ follows a $d$-mixture of Dirichlet with the $j$th component being $\mathcal{D}(\bs{1}-\delta_{j}\alpha_j)$, so that the mixture has density function
\begin{align*}
 h_{\bs{U}}(\bs{u})=\frac{1}{d} \sum_{j=1}^d \frac{\Gamma(d-\alpha_j)}{\Gamma(1-\alpha_j)} u_j^{-\alpha_j}
\end{align*}
for $0<\alpha_j <1, j=1, \ldots, d$.
The
spectral density is given by
\begin{align*}
   h_{\bs{W}}(\bs{w}) = \frac{1}{d} \left(\sum_{j=1}^d \alpha_ju_j\right)^{-1} \left(\prod_{j=1}^d \alpha_ju_d
\right)\left(\sum_{j=1}^d \frac{\Gamma(d-\alpha_j)}{\Gamma(1-\alpha_j)}u_j^{-\alpha_j}\right)\prod_{j=1}^d w_j^{-1}
\end{align*}
for $\alpha_j \in (0,1)$
\item \textbf{Coles and Tawn Dirichlet distribution} (\code{ct})
The Dirichlet model of \cite{Coles:1991}
\begin{align*}
   h_{\bs{W}}(\bs{w}) = \frac{1}{d} \frac{\Gamma \left(1+\sum_{j=1}^d \alpha_j\right)}{\prod_{j=1}^d \alpha_jw_j}
\left(\sum_{j=1}^d \alpha_jw_j\right)^{-(d+1)}\prod_{j=1}^d \alpha_j \prod_{j=1}^d \left(\frac{\alpha_jw_j}{\sum_{k=1}^d
\alpha_kw_k}\right)^{\alpha_j-1}
\end{align*}
for $\alpha_j>0.$
\item \textbf{scaled extremal Dirichlet} (\code{dir})
The angular density of the scaled extremal Dirichlet model with parameters $\rho >0$ and $\bs{\alpha} \in \R^{d}_{+}$ is given, for all $\bs{w} \in \mathbb{S}_d$, by
\begin{align*}
	h_{\bs{W}}(\bs{w})=\frac{\Gamma(\bar{\alpha}+\rho)}{d\rho^{d-1}\prod_{i=1}^d\Gamma(\alpha_i)}
\bigl\langle\{\bs{c}(\bs{\alpha},\rho)\}^{1/\rho},\bs{w}^{1/\rho}\bigr\rangle^{-\rho-\bar{\alpha}}\prod_{i=1}^{d}
\{c(\alpha_i,\rho)\}^{\alpha_i/\rho}w_i^{\alpha_i/\rho-1}.
\end{align*}
where $\bs{c}(\bs{\alpha},\rho)$ is the $d$-vector with entries  $\Gamma(\alpha_i+\rho)/\Gamma(\alpha_i)$ for $i=1, \ldots, d.$
\item \textbf{scaled negative extremal Dirichlet} (\code{negdir})
The angular density of the scaled negative extremal Dirichlet model with parameters $0< \rho < \min(\alpha_1,\dots, \alpha_d)$ and $\bs{\alpha}$ is given, for all $\bs{w} \in \mathbb{R}_+^d$, by
\begin{align*}
h_{\bs{W}}(\bs{w}) =
\frac{\Gamma(\bar \alpha - \rho)}{d\rho^{d-1}\prod_{i=1}^d \Gamma(\alpha_i)}
\bigl\langle\{\bs{c}(\bs{\alpha},-\rho)\}^{-1/\rho},\bs{w}^{-1/\rho}\bigr\rangle^{\rho-\bar \alpha}
\prod_{i=1}^d\{c(\alpha_i,-\rho)\}^{-\alpha_i/\rho}w_i^{-\alpha_i/\rho -1}.
\end{align*}


\item \textbf{H\"usler--Reiss} (\code{hr}), due to \cite{Husler:1989}. It is a special case of the Brown--Resnick process.
While \cite{Engelke:2015} state that H\"usler--Reiss variates can be sampled following the same scheme, the spatial analog is
conditioned on a particular site ($\bs{s}_0$), which complicates the comparisons with the other methods.

Let $I_{-j}=\{1, \ldots, d\} \setminus \{j\}$ and $\lambda_{ij}^2 \geq 0$ be entries of a strictly conditionally
negative definite matrice $\bs{\Lambda}$, for which $\lambda_{ij}^2=\lambda_{ji}^2$. Then, following \cite{Nikoloulopoulos:2009}
(Remark~2.5) and \cite{Huser:2013}, we can write the distribution function as
 \begin{align*}
   \P{\bs{X} \leq \bs{x}} = \exp \left[ -\sum_{j=1}^d \frac{1}{x_j} \Phi_{d-1, \bs{\Sigma}_{-j}} \left( \lambda_{ij}-
\frac{1}{2\lambda_{ij}}  \log\pfrac{x_j}{x_i}, i \in I_{-j}\right)\right].
                  \end{align*}
                  where the partial correlation matrix $\bs{\Sigma}_{-j}$ has elements
                  \begin{align*}
                     \varrho_{i,k; j}= \frac{\lambda_{ij}^2+\lambda_{kj}^2-\lambda_{ik}^2}{2\lambda_{ij}\lambda_{kj}}
                  \end{align*}
and $\lambda_{ii}=0$ for all $i \in I_{-j}$ so that the diagonal entries $\varrho_{i,i; j}=1$.\footnote{\cite{Engelke:2015}
uses the covariance matrix with entries are $\varsigma=2(\lambda_{ij}^2+\lambda_{kj}^2-\lambda_{ik}^2)$, so the resulting
expression is evaluated at $2\bs{\lambda}_{.j}^2-\log\pfrac{x_j}{\bs{x}_{-j}}$ instead. We recover the same expression by
standardizing, since this amounts to division by the standard deviations $2\bs{\lambda}_{.j}$}



The \texttt{evd} package implementation has a bivariate implementation
of the  H\"usler--Reiss distribution with dependence parameter $r$, with $r_{ik}=1/\lambda_{ik}$ or
$2/r=\sqrt{2\gamma(\bs{h})}$ for $\bs{h}=\|\bs{s}_i-\bs{s}_i\|$ for the Brown--Resnick model. In this setting, it is particularly
easy since the only requirement is
non-negativity of the parameter. For inference in dimension $d>2$, one needs to impose the constraint  $\bs{\Lambda}=\{\lambda_{ij}^2\}_{i, j=1}^d \in
\mathcal{D}$ (cf. \cite{Engelke:2015}, p.3), where
\begin{multline*}
   \mathcal{D}=\Biggl\{\mathbf{A}\in [0, \infty)^{d\times d}: \bs{x}^\top\!\!\mathbf{A}\bs{x} <0, \all \bs{x} \in \R^{d}
\setminus\{\bs{0}\} \\ \qquad
\text{ with } \sum_{i=1}^d x_i=0, a_{ij}=a_{ji}, a_{ii}=0 \all i, j \in \{1,\ldots, d\}\Biggr\}.
\end{multline*}
An avenue to automatically satisfy these requirements is to optimize over a symmetric  positive definite  matrix parameter
$\bs{\varSigma}=\mathbf{L}^\top\mathbf{L}$, where $\mathbf{L}$ is an upper triangular matrix whose diagonal element are on the
log-scale to ensure uniqueness of the Cholesky factorization; see \cite{Pinheiro:1996}. By taking
\begin{align*}
   \bs{\Lambda}(\bs{\varSigma})= \begin{pmatrix} 0 & \diag (\bs{\varSigma})^\top \\ \diag(\bs{\varSigma}) &
\bs{1}\diag(\bs{\varSigma})^\top
+ \diag(\bs{\varSigma})\bs{1}^\top - 2 \bs{\varSigma}
\end{pmatrix}
\end{align*}
one can perform unconstrained optimization for the non-zero elements of $\mathbf{L}$ which are in one-to-one correspondence
with those of $\bs{\Lambda}$.

It easily follows that generating $\bs{Z}$ from a $d-1$ dimensional log-Gaussian distribution with covariance $\Co{Z_i,
Z_k}=2(\lambda_{ij}^2+\lambda_{kj}^2-\lambda_{ik}^2)$ for $i,
k \in I_{-j}$ with mean vector $-2\lambda_{\bullet j}^2$  gives
the finite dimensional analog of the Brown--Resnick process in the mixture representation of \cite{Dombry:2016}.

The \texttt{rmev} function checks conditional negative definiteness of the matrix. The easiest way to do so
negative definiteness of $\bs{\Lambda}$ with real entries is to form $\tilde{\bs{\Lambda}}=\mathbf{P}\bs{\Lambda}\mathbf{P}^\top$, where $\mathbf{P}$
is an $d \times d$ matrix with ones on the diagonal, $-1$ on the $(i, i+1)$ entries for $i=1, \ldots d-1$ and zeros elsewhere.
If the matrix $\bs{\Lambda} \in \mathcal{D}$, then the eigenvalues of the leading $(d-1) \times (d-1)$ submatrix of $\tilde{\bs{\Lambda}}$
will all be negative.

% using a variogram or if specifying a
% correlation matrix $\bs{\Sigma}$ with entries $\varrho_{ij}$, by taking $2/r_{ij}=\sqrt(2-2\varrho_{ij}).$
\item \textbf{Brown--Resnick} (\code{br})
The Brown--Resnick process is the extension of the H\"usler--Reiss distribution, and is a max-stable process associated with the
log-Gaussian distribution.
% One of its spectral representation is
% \begin{align*}
% \max_{i \geq 1} \zeta_i \psi_i, \qquad \psi_i(x)=   \exp(\varepsilon(x)-\gamma(x))
% \end{align*}
% where $\eps(x)$ is an intrinsically stationary Gaussian process with semivariogram $\gamma(x)$ constrained so that $\eps(o)=0$
% almost surely.

It is often in the spatial setting conditioned on a location (typically the origin). Users can provide
a variogram function or the resulting covariance function. See \cite{Engelke:2015} or \cite{Dombry:2016} for
more information.
\item \textbf{Extremal Student} (\code{extstud}) of \cite{Nikoloulopoulos:2009}, eq. 2.8, with unit Fréchet margins is
\begin{align*}
    \P{\bs{X} \leq \bs{x}} = \exp \left[-\sum_{j=1}^d \frac{1}{x_j} T_{d-1, \nu+1, \mathbf{R}_{-j}}\left(
\sqrt{\frac{\nu+1}{1-\rho_{ij}^2}}
\left[\pfrac{x_i}{x_j}^{1/\nu}\!\!\!-\rho_{ij}\right], i \in I_{-j}  \right)\right],
\end{align*}
where $T_{d-1}$ is the distribution function of the $d-1 $ dimensional Student-$t$ distribution and the partial correlation
matrix $\mathbf{R}_{-j}$ has diagonal entry \[r_{i,i;j}=1, \qquad
r_{i,k;j}=\frac{\rho_{ik}-\rho_{ij}\rho_{kj}}{\sqrt{1-\rho_{ij}^2}\sqrt{1-\rho_{kj}^2}}\] for $i\neq k, i, k \in I_{-j}$.

The user must provide a valid correlation matrix (the function checks for diagonal elements), which can be obtained from a
variogram.


\item \textbf{Dirichlet mixture} (\code{dirmix}) proposed by \cite{Boldi:2007}, see \cite{Dombry:2016} for details on the
mixture.
The spectral density of the model is
\begin{align*}
h_{\bs{W}}(\bs{w}) = \sum_{k=1}^m \pi_k \frac{\Gamma(\alpha_{1k}+ \cdots + \alpha_{dk})}{\prod_{i=1}^d \Gamma(\alpha_{ik})} \left(1-\sum_{i=1}^{d-1} w_i\right)^{\alpha_{dk}-1}\prod_{i=1}^{d-1} w_{i}^{\alpha_{ik}-1} \end{align*}
The argument \code{param} is thus a $d \times m$ matrix of coefficients, while the argument for the $m$-vector \code{weights} gives the relative contribution of each Dirichlet mixture component.

\item \textbf{Smith model} (\code{smith}), from the unpublished report of \cite{Smith:1990}. It corresponds to a moving maximum
process on a domain $\mathbb{X}$. The de Haan representation of the process is
\begin{align*}
Z(x)=\max_{i \in \mathbb{N}} \zeta_i h(x-\eta_i), \qquad \eta_i \in \mathbb{X}
\end{align*}
where $\{\zeta_i, \eta_i\}_{i \in \mathbb{N}}$ is a Poisson point process on $\R_{+} \times \mathbb{X}$ with intensity measure $\zeta^{-2}\mathrm{d} \zeta \mathrm{d} \eta$ and $h$ is the density of the multivariate Gaussian distribution. Other $h$ could be used in principle, but are not implemented.

\end{enumerate}


\bibliographystyle{apalike2}
\bibliography{bib}


\end{document}
