  We introduce a class of unbiased Monte Carlo estimators for the multivariate
density of max-stable fields generated by Gaussian processes. Our estimators
take advantage of recent results on exact simulation of max-stable fields
combined with identities studied in the Malliavin calculus literature and ideas
developed in the multilevel Monte Carlo literature. Our approach allows
estimating multivariate densities of max-stable fields with precision
$\varepsilon $ at a computational cost of order $O\left( \varepsilon ^{-2}\log
\log \log \left( 1/\varepsilon \right) \right) $.
