  Statistical modeling of a nonstationary spatial extremal dependence structure
is challenging. Max-stable processes are common choices for modeling
spatially-indexed block maxima, where an assumption of stationarity is usual to
make inference feasible. However, this assumption is often unrealistic for data
observed over a large or complex domain. We propose a computationally-efficient
method for estimating extremal dependence using a globally nonstationary, but
locally-stationary, max-stable process by exploiting nonstationary kernel
convolutions. We divide the spatial domain into a fine grid of subregions,
assign each of them its own dependence parameters, and use LASSO ($L_1$) or
ridge ($L_2$) penalties to obtain spatially-smooth parameter estimates. We then
develop a novel data-driven algorithm to merge homogeneous neighboring
subregions. The algorithm facilitates model parsimony and interpretability. To
make our model suitable for high-dimensional data, we exploit a pairwise
likelihood to draw inferences and discuss computational and statistical
efficiency. An extensive simulation study demonstrates the superior performance
of our proposed model and the subregion-merging algorithm over the approaches
that either do not model nonstationarity or do not update the domain partition.
We apply our proposed method to model monthly maximum temperatures at over 1400
sites in Nepal and the surrounding Himalayan and sub-Himalayan regions; we
again observe significant improvements in model fit compared to a stationary
process and a nonstationary process without subregion-merging. Furthermore, we
demonstrate that the estimated merged partition is interpretable from a
geographic perspective and leads to better model diagnostics by adequately
reducing the number of subregion-specific parameters.
