# Spatio-temporal generation of precipitation using a Hidden Markov Model, extended extreme distributions and conditional Gaussian fields

Caroline Cognot, Liliane Bel, David Métivier, Sylvie Parey

► To cite this version:

Caroline Cognot, Liliane Bel, David Métivier, Sylvie Parey. Spatio-temporal generation of precipitation using a Hidden Markov Model, extended extreme distributions and conditional Gaussian fields. 2026. hal-05523768

HAL Id: hal-05523768

https://hal.science/hal-05523768v1

Preprint submitted on 23 Feb 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Spatio-temporal generation of precipitation using a Hidden Markov Model, extended extreme distributions and conditional Gaussian fields

Caroline Cognot¹,²* Liliane Bel² David Métivier³ Sylvie Parey¹

¹ EDF R&D, EDF, 7 Bd Gaspard Monge, 91120 Palaiseau, France

² UMR MIA Paris-Saclay, Université Paris-Saclay, AgroParisTech, INRAE, Place de l'Agronomie, 91120 Palaiseau, France

³ UMR MISTEA, Université de Montpellier, INRAE, Institut Agro, 2 place Pierre Viala, 34000 Montpellier, France

## Abstract

Current weather generators struggle to capture large-scale spatio-temporal variability in rainfall occurrence and intensity. To address this challenge, this study develops a spatio-temporal weather generator for daily precipitation across multiple stations by extending the classical multivariate hidden Markov model to a spatial framework. Both occurrence and intensity processes are driven by spatio-temporal latent fields with distinct covariance structures, allowing for different dynamical behaviours. The rainfall intensity process is simulated conditionally on the occurrence process, ensuring continuity between dry and wet conditions in both time and space. Model inference relies on a novel hybrid Expectation-Maximisation algorithm designed to address the intractability of the complete likelihood, while achieving substantial computational gains compared to composite-likelihood EM approaches. The generator is applied to 51 years of daily precipitation data from 37 weather stations across France. Four interpretable hidden states are identified, whose most probable temporal sequence is consistent with known meteorological patterns. The generator is assessed through several model-agnostic metrics, targeting marginal, pairwise, spatial, and spatio-temporal properties. Results demonstrate that the proposed generator accurately reproduces rainfall statistics across a wide range of spatial and temporal scales, including droughts and extreme precipitation events, validating its relevance for climate impact assessment.

**Keywords:** Composite likelihood, Extended Generalised Pareto Distribution, Gaussian random field, Hidden Markov model, Rainfall simulation, Spatio-temporal modelling

## 1 Introduction

Climate change amplifies both the variability and the frequency of extreme events, making accurate modelling essential for impact assessment. Many sectors, including electricity production and consumption, hydrology, and agriculture, are increasingly affected by these changes. Evaluating their vulnerability requires models that accurately represent key climate variables, e.g., temperature and precipitation together produce river flow (Garçon 1996), climatic variables determine crop growth and yields (Brisson et al. 2003). While historical records provide valuable data, they represent only a single realisation of the climate system and hardly reflect its variability and extremes. They need to be complemented by simulations that better sample the full range of possible climatic conditions. In particular, they must explore rare but high-impact events (for example, cold spells or droughts) and their large-scale spatial co-occurrences (for instance, a dry episode affecting several European river basins at once); as these situations often exert the greatest stress on environmental and energy systems.

Physical climate models, such as the ones supporting the Intergovernmental Panel on Climate change (IPCC) reports, can generate multiple realisations of future climate conditions. The high

<footnote>*Corresponding author</footnote>
computational cost of these models means that uncertainty is typically reproduced using only a limited number of simulations. For instance, in its latest assessment, the French weather service Météo-France relied on only 17 simulations (Soubeyroux et al. 2024). Moreover, these climate models are focused on solving all the atmospheric circulation equations together with physical parameterisation, and generate a wide array of 3D variables that may not be necessary for the desired applications. They also give the evolution of weather variables on a coarse grid, requiring further statistical processing like downscaling and bias correction to adapt to finer scales. In this context, stochastic weather generators have emerged as powerful tools for simulating weather data that capture both spatial and temporal variability (see Intergovernmental Panel On Climate Change (IPCC) 2023, section 10.3.3.7). These generators consider observations of meteorological variables as one of the infinite possible trajectories of a stochastic model. The primary objective is to simulate numerous plausible sequences of the desired variables that share common statistical properties with the observations, focusing on the characteristics that are deemed important for the applications. These stochastic generators are lightweight and fast, allowing a more comprehensive sampling of climate variables than physical models. Since the very first generator (Richardson 1981) that models temperature and rain occurrence at a single location, many others have been proposed, relying on state-of-the-art statistical models and methods. They can be designed for univariate, multivariate, single-site, multisite, or even spatially, depending on the desired applications.

The precipitation variable is inherently difficult to model, as the presence of a mass in zero pre-
vents modelling by a continuous variable. There are several different ways to represent its intermit-
tent nature in both space and time. Point process models can be used to generate the occurrences
(Rodriguez-Iturbe and Eagleson 1987), but they are not easily generalisable to intensity and potentially
multivariate models. Other authors use non-parameterised or semi-parameterised resampling from the
observations (Rajagopalan and Lall 1999; Yiou 2014; X. Li and Babovic 2019). In that case, the spatial
and temporal structure in the simulation is inherited from the structure of the observations. It can be
argued that these methods may be unable to produce possible unobserved events and underestimate
the variability due to constrained spatial dependence.

To overcome the complex modelling of precipitation, it can be split into different weather types,
with the most simple being used by Richardson 1981: dry or wet day at a single location. This idea
of weather type or state has been used in many studies, see Ailliot, Allard, et al. 2015 for a review.
They represent typical situations, and the weather variables are modelled conditionally to a weather
type. This helps break the complexity of the model into different sub-models, and also represents the
temporal non-stationarity. The sequence of weather states can be considered as known and fixed, but
more usually is considered to be a Markov Chain. The weather types can be external to the model or
observed : Vaittinada Ayar et al. 2020; Obakrim, Benoit, and Allard 2025 define the weather types
using circulation variables, first estimate the Markov chain parameters, and then consider them known
as a discrete covariate. Other approaches use a Hidden Markov Model (HMM), where the weather
states are directly learnt by the model as a latent variable and condition the observed distributions.
Defining a spatially coherent HMM is not straightforward and can be done using different model
structures. Dawkins et al. 2022 sets one chain by location, then uses copulas to correlate them, which
makes the model inference have a high computational cost even for a low number of stations. On the
contrary, approaches in which all stations are simultaneously in a shared state have been proposed.
Gobet, Métivier, and Parey 2025 generates rain occurrences using spatial independence of all locations
conditionally to a shared weather state, then generates intensity independently of the occurrence in
each state using Gaussian copulas. Hughes, Guttorp, and Charles 1999; Kirshner 2005 propose models
that do not rely on conditional dependence but introduce many parameters to represent pairwise
interactions, which can cause computational issues when the number of locations to model increases.
Ailliot, Thompson, and Thomson 2009 uses conditional dependence in the form of a latent Gaussian
process through its covariance matrix, but is limited to a few locations (7 in their application) and
uses simple distributions for the rain intensity. For more flexible, possibly non-Gaussian dependence
in the rain occurrence, Dawkins et al. 2022 uses vine copulas , which rely on pairs of locations, making
it intractable for a high spatial dimension.

Another popular choice is to model precipitation using a latent meta-Gaussian field, censored into
zero and non-zero precipitations depending on a threshold (Durban and Glasbey 2001; Hundecha,
Pahlow, and Schumann 2009; Bourotte, Allard, and Porcu 2016; Baxevani and Lennartsson 2015;
Benoit, Allard, and Mariethoz 2018; Sparks et al. 2018; Zhou, Y. Meng, and Abbaspour 2019; Vait-
tinada Ayar et al. 2020; Obakrim, Benoit, and Allard 2025). On the one hand, this approach is interesting in that it represents well the spatio-temporal continuity of the precipitations : locations close in time or space will tend to have close precipitations. On the other hand, other authors point out that the underlying spatial or temporal dependence may be different in the occurrence and the intensity process and use a different latent Gaussian field for each variable, with the occurrence process acting as a mask on the intensity process (Kleiber, R. W. Katz, and Rajagopalan 2012; Evin, Favre, and Hingray 2018). This makes it difficult to represent the spatial continuity of the observed rain intensity.

In addition to its complex spatio-temporal structure, the model for the rain intensity must represent well the low precipitations (linked to droughts) and the high precipitations (floods). This has been achieved in the literature using various marginals: exponential mixtures for the low/moderate/high precipitations (D. S. Wilks 1998; Daniel S. Wilks 2009; Z. Li 2014; Touron 2019; Gobet, Métivier, and Parey 2025), Gamma and Gamma mixtures(Furrer and R. Katz 2007; Kleiber, R. W. Katz, and Rajagopalan 2012; Verdin et al. 2019; Chandler 2020), empirical quantile fitting (Sparks et al. 2018; Obakrim, Benoit, and Allard 2025), mixture of exponential and Generalised Pareto (Dawkins et al. 2022).

To the best of our knowledge, no generator has yet been proposed to tackle a high number of corre-
lated locations using the HMM framework for both the occurrence and the appropriate rain intensity
process. In this study, a large-scale spatio-temporal stochastic weather generator for precipitation is
proposed. Several ideas from the literature are combined to jointly model precipitation occurrence
and intensity, while ensuring scalability to a high number of locations. Arbitrarily closely spaced
stations can be accommodated due to conditional dependence across all locations driven by a single
Markov chain in the occurrence process. Compared to multi-chain HMM approaches, this formulation
is more scalable to a large number of stations. Explicit modelling of the occurrence process through
an HMM enhances interpretability, compared to approaches based on externally defined weather types
derived from covariates, and yields more realistic temporal dynamics. The hidden states are naturally
interpretable as weather regimes. The rain intensity distributions at each location are accurately rep-
resented through the use of Extended Generalised Pareto Distributions (EGPDs), as in Evin, Favre,
and Hingray 2018. Seasonality is fully accounted for through smoothly varying seasonal parameters, as
in Touron 2019, which represents an improvement over much of the existing literature, such as Ailliot,
Thompson, and Thomson 2009; Evin, Favre, and Hingray 2018; Obakrim, Benoit, and Allard 2025,
where parameters are typically assumed to be constant within each month or season. Two methodolog-
ical novelties are introduced. Efficient inference for a hidden Markov model (HMM) with conditionally
dependent occurrence distributions across multiple locations is achieved via an EM algorithm with a
modified M-step, in which the objective function involving high-dimensional computations is replaced
by an objective that is easier to maximise. This modification transforms a computationally infeasible
optimisation into a fast optimisation step. In addition, a novel conditional simulation method is pro-
posed to couple the occurrence process, generated by the HMM, with the intensity process, through a
truncated latent Gaussian field.

The generator is fitted to data from 37 stations in France, and its performance is evaluated by
comparing a range of spatio-temporal metrics between observations and simulations, at both local and
large scales relevant to risk assessment.

The paper is organised as follows. The statistical model, its inference, and the conditional simu-
lation procedure, along with relevant indicators, are presented in Sect. 2. The application to a set of
French station data is described in Sect. 3. Finally, the paper is concluded with a discussion in Sect. 4.

2 Methodology

2.1 Model definition

Unlike other meteorological variables such as temperature or moisture, precipitation cannot be mod-
elled as a spatio-temporal continuous variable, due to having a mass in zero. Two procedures are
usually considered to handle this specificity. The first one is to consider a continuous variable for
the rain intensity and threshold it such that values below the threshold are considered as zeroes. A
drawback of this strategy is that the occurrence process and the intensity process share the same
spatio-temporal structure. A second option is to consider the product of a binary variable (presence-
absence) and a continuous variable representing the intensity of precipitation, but in this case the two processes are assumed independent which is not realistic. A third strategy is adopted here. First, the binary occurrence model is estimated, and the continuous intensity model is then estimated conditionally on the results of the occurrence model in order to preserve the dependence structure. This strategy allows for distinct spatial and temporal dynamics in the occurrence and intensity processes while maintaining their dependence.

### 2.1.1 Occurrence process

The rain occurrence at D locations and N time steps $((Y_s^{(n)})_{s \in \{1,...,D\}})_{n \in \{1,...,N\}}$ is represented as a non-homogeneous Hidden Markov Model with seasonal parameters.

**The general Hidden Markov Model (HMM)** The Hidden Markov Model provides a framework for representing complex dynamics by linking observed weather variables to an underlying hidden discrete variable that evolves over time according to a Markov process. Let Z be the hidden variable, its values can be viewed as weather states, shared by all locations $s \in \{1,...,D\}$. It is assumed to take values $k \in \{1,...,K\}$ and its temporal evolution $\{Z^{(n)}, n \in \{1,...,N\}\}$ is assumed to follow a Markov chain of order 1. The Markov chain is parameterised by its initial and transition probabilities:

$$ \eta = (\eta_1, ..., \eta_K) = \mathbb{P}(Z^{(1)} = k) \qquad \text{Initial probabilities} \qquad (1) $$

$$ Q^{(n)}(k, \ell) = \mathbb{P}(Z^{(n+1)} = \ell | Z^{(n)} = k) \qquad \text{Transition probabilities.} \qquad (2) $$

Conditionally to the unobserved weather state $Z^{(n)}$ at time n, the observed multisite rain occurrence vector $\mathbf{Y}^{(n)} = (Y_s^{(n)})_{s \in \{1,...,D\}} | Z^{(n)}$, is modelled using the conditional distribution $L_{Y|Z^{(n)}}^{(n)} = \mathbb{P}(Y_1^{(n)}, ..., Y_D^{(n)} | Z^{(n)})$.

**A correlated binary-valued distribution** Previous approaches to multisite rain occurrence modelling with HMMs introduced in Zucchini and Guttorp 1991 rely on a conditional independence assumption, which states that rain occurrences at different stations are independent Bernoulli variables given the hidden state, i.e. $\mathbb{P}(Y_1^{(n)}, ..., Y_D^{(n)} | Z^{(n)}) = \prod_{s \in D} \mathbb{P}(Y_s^{(n)} | Z^{(n)})$. This assumption is justified when stations are far apart but is no longer valid for stations at short distances and hence exhibit dependence. Introducing multisite correlation via several dependent Markov chains is possible but does not scale efficiently. To introduce spatial dependence among Bernoulli variables, previous studies (Kirshner 2005) have used Chow-Liu trees, approximating the joint distribution with a tree of pairwise dependencies that maximizes mutual information. D. S. Wilks 1998; Evin, Favre, and Hingray 2018 used a latent Gaussian process $\mathbf{X}_Y \in \mathbb{R}^D$ with an empirical covariance matrix, without hidden states. In the present approach, a spatial covariance function is employed, similarly as in spatial logistic regression models (Cessie and Houwelingen 1994; Heagerty and Lele 1998). Let $\lambda = (\lambda_1, ..., \lambda_D) = (\mathbb{P}(Y_1 = 1), ..., \mathbb{P}(Y_D = 1))$ be the probability of rain occurrence on the D locations and $C_Y$ a correlation function with parameters $\theta_{CY}$. The correlated binary occurrence vector $\mathbf{Y} = (Y_1, ..., Y_D)$ is generated in two steps by

$$ (X_{Y_1}, ..., X_{Y_D}) \sim \mathcal{N}(0, C_Y), \quad \forall s, Y_s = \begin{cases} 1 & \text{if } X_{Y_s} \le \Phi^{-1}(\lambda_s), \\ 0 & \text{else} \end{cases}, \qquad (3) $$

with $\Phi$ the cumulative distribution function (CDF) of the standard normal distribution. This construction, sometimes referred to as a dichotomised Normal (Cox and N. Reid 2004) or Gaussian copula for Bernoulli margins (Fan et al. 2017), ensures that the marginal distribution of $Y_s$ is Bernoulli($\lambda_s$) while introducing spatial dependence through $C_Y$ with a minimal number of additional parameters. More general multivariate Bernoulli constructions exist (Dai, Ding, and Wahba 2013) but typically require parameterising $2^D - D + 1$ independent probabilities, which is not scalable.

Let $f_{\theta_{CY}}$ be the probability density function of the multivariate Gaussian field with covariance $C_Y$ and parameters $\theta_{CY}$. The likelihood of the model defined in Eq. (3) for a sample $\mathbf{y} = (y_1, ..., y_D)$ is

$$ L_Y(\mathbf{y}; \lambda, \theta_{CY}) = \int_{a_1}^{b_1} \cdots \int_{a_D}^{b_D} f_{\theta_{CY}}((x_1, ..., x_D)) d(x_1, ..., x_D) \qquad (4) $$

$$ a_i = -\infty \text{ if } y_i = 1, \Phi^{-1}(\lambda_i) \text{ else, } \quad b_i = \infty \text{ if } y_i = 0, \Phi^{-1}(\lambda_i) \text{ else.} $$
<figure><img src="image_2.png" /><figcaption>Figure 1: Hidden Markov Model with time-varying coefficients and marginal memory.</figcaption></figure>

**Marginal memory** As noted by Ailliot, Thompson, and Thomson 2009, standard HMMs fail to capture 1-day lagged correlations and often underestimate the length of wet or dry spells. To produce more realistic sequences, an additional first-order temporal dependence is introduced, as in Gobet, Métivier, and Parey 2025, in the marginal distributions of $\mathbf{Y}$ through the parameter $\lambda_{k,s,y_s^{(n-1)}}^{(n)} = \mathbb{P}(Y_s^{(n)} = 1 | Z^{(n)} = k, Y_s^{(n-1)} = y_s^{(n-1)})$. This formulation accounts for the previous day's rainfall at each location, effectively doubling the number of rain probability parameters. Introducing higher-order temporal dependencies would further increase the number of parameters and the complexity of the model. The model with marginal memory is illustrated in Figure 1.

### 2.1.2 Intensity process

The rain intensity process $(\mathbf{R}^{(n)})_{n \in \{1, \dots, N\}}$ representing the non-zero precipitation amount is modelled as a truncated Gaussian spatial process, after an adequate transformation of the marginals according to the weather states identified in the occurrence step.

**Marginal parameterisation** The rain intensity marginals are strictly positive distributions that must capture both the bulk and the tail of precipitation. Standard models, such as double exponential or Gamma distributions, are easy to fit on the bulk but are not appropriate for the tail. High precipitation values beyond a threshold are usually represented using the Generalised Pareto Distribution (GPD, Pickands 1975; Coles 2001).

To represent the whole distribution, some authors (Carreau, P. Naveau, and Sauquet 2009) resort to a mixture of the GPD with some other distribution intended to model the bulk below the threshold. In this work, the selection of the threshold is avoided by using Extended Generalised Pareto Distributions (EGPD) (Philippe Naveau et al. 2016), where only a very small threshold $t_{EGPD}$ must be chosen (in this study, $t_{EGPD} = 1$ mm of rainfall). The EGPD is designed to be suitable for both the bulk and the tail of the distribution by applying a monotone transformation $G$ to the GPD CDF. This transformation is chosen to ensure that the bulk of the distribution near 0 behaves like $x^{\kappa}$ for $\kappa$ a positive parameter, while the tail behaves like a GPD. A random variable $X$ with an EGPD distribution can be written as

$$X = \sigma_s H_{\xi_s}^{-1}[G^{-1}(U)] \quad (5)$$

where $U \sim \mathcal{U}([0, 1])$, and $H_{\xi_s}(x) : H(x; 1, \xi_s)$ is the CDF of the GPD distribution of unit variance and shape $\xi_s$.

Here the truncated-Beta transform $G_{\kappa}$ proposed by Gamet and Jalbert 2022 is chosen as it is more flexible in the behaviour of the distribution near the origin than those proposed by Philippe Naveau et al. 2016, and therefore is easier to fit to the precipitation data.

To complete the distribution on the very small rain values $0 < r \le t_{EGPD}$, a simple uniform distribution is used. It is intended to take into account that rain gauges typically measure rainfall in 0.1 mm increments, and observations in this range are already affected by measurement noise. The resulting CDF is

$$F(r) = pF_1(r) \mathbf{1}_{\{0<r \le t_{EGPD}\}} + [p + (1-p)F_2(r - t_{EGPD})] \mathbf{1}_{\{r > t_{EGPD}\}}, \quad (6)$$

where $p = \mathbb{P}(0 < r \le t_{EGPD})$, $F_1$ is the uniform CDF on $[0, t_{EGPD}]$, and $F_2$ is the EGPD CDF defined by $F_2(r) = G_{\kappa}(H_{\sigma, \xi}(r))$. For one location $s$ in each class $k$ at time $n$, 4 parameters are to be estimated: $\kappa_{k,s}^{(n)}, \sigma_{k,s}^{(n)}, \xi_{k,s}^{(n)}$ for the EGPD distribution, and $p_{k,s}^{(n)}$ for the probability of low rain.
**Spatio-temporal structure** A latent Gaussian field $X_R$ is introduced to represent the spatio-temporal dependence of precipitations, it is truncated according to a rain/no-rain threshold. Let

$$T(s, n) = \Phi^{-1}(1 - \lambda_{z^{(n)}, s, y_{s^{(n-1)}}}^{(n)}) \quad (7)$$

a location-time specific threshold from which are derived lower and upper bounds given $Y_s^{(n)}$:

$$l(s, n) = -\infty; \quad u(s, n) = T(s, n) \quad \text{if } Y_s^{(n)} = 0, \quad (8)$$

$$l(s, n) = T(s, n); \quad u(s, n) = +\infty \quad \text{if } Y_s^{(n)} = 1. \quad (9)$$

and

$$X_R(., n) | \mathbf{Y}^{(n)} \sim \text{TMVN}(0, C_R, l(., n), u(., n)), \quad (10)$$

a truncated Gaussian field with correlation function $C_R$, parameters $\theta_{CR}$, and bounds $(\ell, u)$. The threshold definition ensures that marginally $X_R(s, \cdot)$ is standard Gaussian distributed. The value of $X_R(s, n)$ is below the threshold if $Y_s^{(n)} = 0$, and above if $Y_s^{(n)} = 1$.

The rain intensity is recovered by back-transformation of the adequate CDF:

$$R_s^{(n)} = \begin{cases} (F_{z^{(n)}, s}^{(t_n)})^{(-1)} (\Phi_{T(s,n)}(X_R(s,n))) & \text{if } Y(s,n) = 1 \\ 0 & \text{if } Y(s,n) = 0, \end{cases} \quad (11)$$

where $\Phi_b$ is the truncated Gaussian CDF, $\Phi_b(z) = \frac{\Phi(z) - \Phi(b)}{1 - \Phi(b)} 1_{z>b}$.

The spatial-temporal correlation function of $X_R$ is taken in the non-separable class Gneiting 2002 as

$$C_R(h, u) = \frac{1}{\left(\left(\frac{u}{\tau}\right)^{2\alpha} + 1\right)^{\beta+\delta}} \mathcal{M} \left( \frac{h}{\sqrt{\left(\left(\frac{u}{\tau}\right)^{2\alpha} + 1\right)^{\beta} - r; \nu}} \right) \quad (12)$$

where $\mathcal{M}(x; r; \nu)$ is the value of the Matérn covariance function at distance $x$, with smoothness parameter $\nu$ and scale parameter $r$. As in Bourotte, Allard, and Porcu 2016; Obakrim, Benoit, and Allard 2025, this choice is made because it is flexible enough to capture the spatio-temporal dependence in meteorological variables, where separable functions are not sufficient. The separability parameter $\beta$ is equal to 0 when the covariance reduces to a separable function, a product of spatial and temporal functions. When $\beta > 0$, the model is said to be non-separable and denotes a complex interaction between space and time. $\delta$ is a temporal covariance parameter for additional flexibility.

### 2.1.3 Overview of the model parameters

**Periodic parameterisation** All parameters are described using a periodic function, ensuring a smooth representation of seasonal variability throughout the year. Each time step $n$ corresponds to a calendar date $t_n$, scaled within the yearly cycle $t \in [1, T]$, with $T = 366$. $t = 60$ is skipped 3 times out of 4 as it corresponds to 29 February. More precisely, for $\theta$ a parameter, let $P_\theta(t)$ be a trigonometric polynomial of degree $d$:

$$P_{\theta}(t) = c_{0}^{\theta} + \sum_{j=1}^{d} (c_{2j-1}^{\theta} \cos(2\pi j t/T) + c_{2j}^{\theta} \sin(2\pi j t/T)) \quad (13)$$

The parameters are then obtained through an appropriate transformation $\theta^{(n)} = \psi_{\theta}(P_{\theta}(t_n)) = \psi_{\theta}(P_{\theta}(n))$ to enforce constraints, the exponential function if the parameter has to be positive, the sigmoid function $\zeta(x) = \frac{e^x}{1+e^x}$ if it is a probability and has to lie in $[0, 1]$.

**Rain occurrence parameters** The transition matrices coefficients of the Markov chain in Eq. (2) are written:

$$Q^{(n)}(k, \ell) = \frac{\exp(P_{Q(k,l)}(n))}{1 + \sum_{l=1}^{K-1} \exp(P_{Q(k,l)}(n))} \quad \text{for } l < K, \qquad Q^{(n)}(k, K) = \frac{1}{1 + \sum_{l=1}^{K-1} \exp(P_{Q(k,l)}(n))}, \quad (14)$$
and the rain probabilities

$$ \lambda_{k,s,y}^{(n)} = \zeta(P_{\lambda_{k,s,y}}(n)). \qquad (15) $$

Preliminary estimations of the rain occurrence covariance $C_Y$ with Matérn covariance led to a smoothness parameter close to 1/2, which motivated the use of the simpler exponential covariance function with periodic range parameter

$$ C_{Y,k}^{(n)}(h; \rho) = \exp(-h/\rho_{CY,k}^{(n)}), \quad \rho_{CY,k}^{(n)} = \exp(P_{\rho_{CY,k}}(n)). \qquad (16) $$

**Rain intensity parameters** Similarly, the parameters $p_{k,s}^{(n)} \in [0, 1]$, probability of low rain (below the threshold $t_{EGPD} = 1.0$ mm, $\sigma_{k,s}^{(n)} > 0$, $\kappa_{k,s}^{(n)} > 0$ the EGPD parameters are written:

$$ p_{k,s}^{(n)} = \zeta(P_{p_{k,s}}(n)), \quad \sigma_{k,s}^{(n)} = \exp(P_{\sigma_{k,s}}(n)), \quad \kappa_{k,s}^{(n)} = \exp(P_{\kappa_{k,s}}(n)) \qquad (17) $$

It is commonly accepted in the statistician and hydrologist communities (Yiou et al. 2008; Ragulina and Reitan 2017) that the shape parameter of the extreme distribution $\xi_{k,s}^{(n)}$ when dealing with precipitation data is positive and less than 1. Accordingly, the constraint $\xi < 1$ is imposed, and this parameter is modelled as

$$ \xi_{k,s}^{(n)} = \zeta(P_{\xi_{k,s}}(n)). \qquad (18) $$

## 2.2 Model inference

### 2.2.1 Rain occurrence

The parameters to be estimated in the rain occurrence model (section 2.1.1), are

$$ \theta = ((Q(k, \ell))_{k, \ell \in \{1, \dots, K\}}, (\eta_k)_{k \in \{1, \dots, K\}}, (\lambda_{k, \mathbf{y}})_{k \in \{1, \dots, K\}}, \mathbf{y} \in \{0, 1\}^D, (\rho_{CY,k})_{k \in \{1, \dots, K\}}). $$

**EM algorithm** Models with hidden variables are classically inferred using the Expectation-Maximisation (EM) algorithm Dempster, Laird, and D. B. Rubin 1977, an iterative two-step procedure that is guaranteed to converge to a local maximum of the likelihood. An EM iteration can be written as

$$ \theta^{(q+1)} \in \arg \max_{\theta} \left( \mathbb{E}_{(Z|\mathbf{Y}, \theta^{(q)})} \log(L(\mathbf{Y}, Z | \theta)) \right), \qquad (19) $$

where $L$ denotes the likelihood function. The first step(E) is devoted to the computation of the expectation $\mathbb{E}_{(Z|\mathbf{Y}, \theta^{(q)})} \log(L(\mathbf{Y}, Z | \theta))$ and the second one to its maximisation (M) in $\theta$. In the hidden Markov model (HMM) setting, the computation of the expectation is based on the smoothing probabilities Baum et al. 1970,

$$ \pi_{n|N}(k) = \mathbb{P}(Z^{(n)} = k | \mathbf{y}^{(1)}, \dots, \mathbf{y}^{(N)}), \quad \pi_{n,n+1|N}(k, \ell) = \mathbb{P}(Z^{(n)} = k, Z^{(n+1)} = \ell | \mathbf{y}^{(1)}, \dots, \mathbf{y}^{(N)}). $$

These quantities are computed using the Forward-Backward algorithm, see Gobet, Métivier, and Parey 2025, Appendix E for a detailed description in a similar setting. The maximisation step reduces to maximising

$$ \begin{aligned} \mathcal{R}(\theta, \theta^{(q)}) = & \sum_{k=1}^{K} \pi_{1|N}^{(q)}(k) \log(\eta_k) + \sum_{k, \ell=1}^{K} \sum_{n=1}^{N-1} \pi_{n,n+1|N}^{(q)}(k, \ell) \log(Q^{(n)}(k, \ell)) \quad (20) \\ & + \sum_{k=1}^{K} \sum_{n=1}^{N} \pi_{n|N}^{(q)}(k) \log L_Y(\mathbf{y}^{(n)} | \lambda_{k, \mathbf{y}^{(n-1)}}^{(n)}, \rho_{CY,k}^{(n)}). \end{aligned} $$

The three terms of Eq.(20) involve respectively the initial probabilities, the transition probabilities and the emission parameters. They can be maximised separately, as the sets of parameters do not overlap. The last term involves the weighted log-likelihood of samples from the spatial Bernoulli model in Eq. (3). Updating the corresponding parameter requires maximising for each $k$ a sum $N$ integrals of the $D$-variate normal distribution, which is highly intractable.
**Pairwise-likelihood EM** A suitable alternative is provided by Composite-Likelihood EM (CLEM) algorithms Varin, Høst, and Skare 2005; Gao and Song 2011; Stoehr and Robin 2024. The composite likelihood (Varin, Nancy Reid, and Firth 2011) is defined as

$$ \log CL(Y | \theta) = \sum_{a \in \mathcal{A}} w_a \log p(\mathbf{Y}_a | \theta), \qquad (21) $$

where $\mathcal{A}$ denotes a set of index subsets of $Y$ on which lower-dimensional likelihoods are computed, $w_a$ are positive weights, and $Y_a$ is the subvector of $Y$ indexed by $a$. When the subsets are pairs of locations ($a = \{i, j\}$), this leads to a pairwise EM (PEM) formulation. One iteration of the PEM algorithm is given by

$$ \theta^{(q+1)} \in \arg \max_{\theta} \left( \sum_{i,j} w_{ij} \mathbb{E}_{(Z|\mathbf{Y}_{ij}, \theta^{(q)})} \log(L(\mathbf{Y}_{ij}, Z | \theta)) \right). \qquad (22) $$

Each PEM iteration increases the pairwise likelihood and the procedure converges to a local maximum of this objective. When a large number of independent replicates of the spatio-temporal process are available, the estimator given by the pairwise likelihood is normally distributed, similarly to the maximum likelihood estimator (Varin, Nancy Reid, and Firth 2011, Sec. 2.3). When the sample consists of a large number of correlated observations in both time and space, the estimator is still asymptotically unbiased, and Cox and N. Reid 2004 argue that the asymptotic normality may still hold when the sample consists of a weakly dependent process, provided the dependence is not too strong and the chosen low-dimensional likelihoods retain identifiability of the parameters. Notably, one of their positive examples is a simplified version of the Spatial Bernoulli in Eq.(3), with constant covariance between all pairs $i \neq j$ and fixed success probability $\lambda_s$.

In the HMM case, the criterion to maximise in the M-step of the PEM algorithm is

$$ \mathcal{R}_c(\theta, \theta^{(q)}) = \sum_{i,j} w_{ij} \sum_{k=1}^{K} \pi_{ij,1|N}^{(q)}(k) \log(\eta_k) + \sum_{i,j} w_{ij} \sum_{k,\ell=1}^{K} \sum_{n=1}^{N-1} \pi_{ij,n,n+1|N}^{(q)}(k, \ell) \log(Q^{(n)}(k, \ell)) \quad (23) \\ + \sum_{k=1}^{K} \sum_{n=1}^{N} \sum_{i,j} w_{ij} \pi_{ij,n|N}^{(q)}(k) \log L_{Y2}(\mathbf{y}_{ij}^{(n)} | \mathbf{\lambda}_{k,ij,\mathbf{y}^{(n-1)}}^{(n)}, \rho_{CY,k}^{(n)}), $$

and the E-step is devoted to the calculations of

$$ \pi_{ij,n|N}(k) = \mathbb{P}(Z^{(n)} = k | \mathbf{y}_{ij}^{(1)}, \dots, \mathbf{y}_{ij}^{(N)}), \quad \pi_{ij,n,n+1|N}(k, \ell) = \mathbb{P}(Z^{(n)} = k, Z^{(n+1)} = \ell | \mathbf{y}_{ij}^{(1)}, \dots, \mathbf{y}_{ij}^{(N)}). $$

This E-step can be interpreted as a collection of local Forward-Backward algorithms applied to low-dimensional subsets of the data indexed by $\{i, j\}$, while the M-step aggregates the resulting local posterior information.

**Hybrid EM algorithm** Computing the posterior distributions $\pi_{ij,n|N}$ in the PEM algorithm becomes computationally demanding when the number of locations is large, leading to practical difficulties in parameter estimation. To address this issue, an hybrid EM algorithm that exploits the structure of the classical EM algorithm is proposed. The key idea is to replace the pairwise posterior probabilities $\pi_{ij,n|N}$ in the composite criterion $\mathcal{R}_c$ by the classical EM smoothing probabilities $\pi_{n|N}$, which are computed efficiently using the standard Forward-Backward algorithm. In the M-step, this modified criterion $\mathcal{R}'_c$ is maximised. These modifications, constituting the hybrid M-step, are not guaranteed to yield a monotonic increase of either the full log-likelihood or the composite log-likelihood, and the objective function must therefore be monitored during optimisation. Nevertheless, the parameter values proposed by the hybrid M-step are expected to be close to those obtained with a classical M-step. In preliminary experiments, the hybrid EM algorithm converged substantially faster to the true parameter than the PEM algorithm while producing comparable results, if not better.

The hybrid M-step consists, as its EM (see Eq.(20)) and PEM (see Eq.(23)) counterparts, in the independent maximisation of three sums. The maximisation of the first two sums, with respect to the transition parameters $(Q(k, \ell))_{k,\ell \in \{1, \dots, K\}}$, $(\eta_k)_{k \in \{1, \dots, K\}}$, is identical to that in classical EM. Updating the parameters $\lambda_k$ and $\rho_{CY,k}$ requires maximising, for each $k$, a criterion of the form

$$ J_k = \sum_{n=1}^{N} \sum_{i,j} w_{ij} \pi_{n|N}^{(q)}(k) \log L_{Y2}(\mathbf{y}_{ij}^{(n)} | \mathbf{\lambda}_{k,ij,\mathbf{y}^{(n-1)}}^{(n)}, \rho_{CY,k}^{(n)}), \qquad (24) $$
with $L_{Y2}$ given in Equation (4) for two variables. The maximisation can be split into sequential tasks following the Expectation-Conditional-Maximisation algorithm X.-L. Meng and Donald B. Rubin 1993, with rain probabilities estimated first, followed by spatial parameters. To further reduce computational cost, rain occurrence probabilities are updated using a simplified criterion involving only identical-index pairs, corresponding to marginal likelihood contributions,

$$J'_{k} = \sum_{n=1}^{N} \sum_{a=(i,i)} w_{i,i} \pi_{n|N}^{(q)}(k) \left( y_{i}^{(n)} \log \lambda_{k,i, y_{i}^{(n-1)}}^{(n)} + (1 - y_{i}^{(n)}) \log(1 - \lambda_{k,i, y_{i}^{(n-1)}}^{(n)}) \right). \quad (25)$$

Empirically, this approximation yields results comparable to those obtained by maximising the full expression, as these terms dominate the likelihood with respect to the rain probability parameters.

Following Bourotte, Allard, and Porcu 2016, the weights in the pairwise expression are defined as $w_{ij} = \mathbf{1}_{d(s_i, s_j) < t_d}$, where $d(s_i, s_j)$ denotes the distance between locations $s_i$ and $s_j$, and $t_d$ is a threshold. The hybrid M-step for the spatial scale parameters involves maximising a sum of $ND(D-1)/2$ bivariate integrals of the bivariate normal CDF. The periodic formulation and use of the complete probability formula reduce the number of integrals that must be precomputed to $4TD(D-1)/2$.

The hybrid EM algorithm performs well in simulation studies. To ensure that the final estimate corresponds to a local maximum of the pairwise likelihood, the hybrid EM can be used as a preliminary step to obtain reliable initial parameter estimates, which are then supplied to the PEM algorithm to efficiently converge to a candidate that is guaranteed to be a local maxima of the pairwise likelihood.

### 2.2.2 Marginal rain intensity inference

**Low rainfall probability** The parameters $p_{k,s}$ in Eq. (6) are estimated by fitting the trigonometric polynomial on the empirical frequencies of low rain calculated over the whole period.

**Penalised likelihood** $\sigma$ and $\kappa$, the EGPD parameters, are not fully independent, as they both capture aspects of the spread of the distribution in the bulk. To mitigate instability a penalisation term on $\kappa$ is incorporated into the inference procedure, following the approach implemented for the constant-parameter case in the ExtendedExtremes.jl Julia package (Gamet and Jalbert 2024). With $L_{EGPD}$ the likelihood of an EGPD sample,

$$p\ell((r_s^{(1)}, \dots, r_s^{(N)}); \kappa_{k,s}, \sigma_{k,s}, \xi_{k,s}) = \sum_{k=1}^{K} \sum_{t=1}^{T} \sum_{\substack{n \in \{1, \dots, N\} \\ t_n = t, z^{(n)} = k, r_s^{(n)} > u}} \left[ \log L_{\text{EGPD}}(r_s^{(n)} - u \mid \kappa_{k,s}^{(t)}, \sigma_{k,s}^{(t)}, \xi_{k,s}^{(t)}) - \frac{1}{N_k} \frac{(\kappa_{k,s}^{(t)} - 1)^2}{0.1} \right], \quad (26)$$

where $N_k$ is the total number of values of $r_s > u$ in state $k$. This penalisation ensures that $\kappa$ is stable and close to 1. Empirical tests showed that this ensures a realistic shape-scale trade-off in the tail and prevents unstable tail fits.

### 2.2.3 Truncated Gaussian inference

To make inference on latent covariance parameters $\theta_{CR}$ tractable, a weighted pairwise likelihood

$$\log L_{R,wp}((r_s^{(n)})_{s \in \{1, \dots, D\}, n \in \{1, \dots, N\}}, \theta_{CR}) = \sum_{s_1, s_2} \sum_{t_1, t_2} w_{s_1, s_2, t_1, t_2} \log L_{R2}(r_{s_1}^{(t_1)}, r_{s_2}^{(t_2)}; \theta_{CR}) \quad (27)$$

is employed, following Obakrim, Benoit, and Allard 2025. The weights are chosen as $w_{s_1, s_2, t_1, t_2} = 1$ if the distance between locations $s_1$ and $s_2$ is less than a threshold $t_d$ and if the temporal distance $|t_2 - t_1|$ is less than a threshold $t_t$.

The computation of each bivariate log-likelihood $\log L_{R2}(r_{s_1}^{(t_1)}, r_{s_2}^{(t_2)}; \theta_{CR})$ depends on whether the
values of $r_{s_1}^{(t_1)}$ and $r_{s_2}^{(t_2)}$ are zero or positive. Each value is computed as

$$ \log L_{R2}(r_{s_1}^{(t_1)}, r_{s_2}^{(t_2)}; \theta_{CR}) = \left\{ \begin{array}{ll} -\frac{1}{2} \left[ \log(1 - \rho^2) + \frac{x_R(s_1, t_1)^2 - 2\rho x_R(s_1, t_1)x_R(s_2, t_2) + x_R(s_2, t_2)^2}{1 - \rho^2} \right], & r_{s_1}^{(t_1)} > 0, r_{s_2}^{(t_2)} > 0, \\ \\ \log \Phi \left( \frac{T(s_2, t_2) - \rho x_R(s_1, t_1)}{\sqrt{1 - \rho^2}} \right), & r_{s_1}^{(t_1)} > 0, r_{s_2}^{(t_2)} = 0, \\ \\ \log \Phi_2(T(s_1, t_1), T(s_2, t_2), \rho), & r_{s_1}^{(t_1)} = 0, r_{s_2}^{(t_2)} = 0, \end{array} \right. \quad (28) $$

where $x_R(s, n) = \Phi_{T(s,n)}^{-1}(F_{z^{(n)},s}^{(n)}(r_s^{(n)}))$ is computed only in the strictly positive rain case and $\rho = C_R((s_1, t_1), (s_2, t_2))$. To simplify the inference and contrary to the other parameters, the covariance function parameters are chosen constant in time and independent on the hidden states. Indeed, estimating the covariance parameters for each hidden state would require dividing the data into K clusters, causing the uncertainty to increase excessively. In particular, only non-zero rainfall observations provide information about the correlation parameters, which constitutes a major limitation for the driest states.

## 2.3 Simulation

The sequences of states and rain occurrences $z, y$ are first sampled according to Algorithm 1. This is achieved by first simulating the Markov chain $z$ and then using the corresponding parameters to generate the rain occurrences $y$. The second step, using the same hidden states parameters, is to

**Algorithm 1** Simulation of the hidden states and rain occurrence $Z, Y$.

**Input:** Rain occurrence parameters $((Q(k, \ell))_{k, \ell \in \{1, \dots, K\}}, (\eta_k)_{k \in \{1, \dots, K\}}, (\lambda_{k, y})_{k \in \{1, \dots, K\}}, y \in \{0, 1\}^D, (\rho_{CY,k})_{k \in \{1, \dots, K\}})$, previous occurrence $y^{(0)}$.

**Output:** Simulated sequence of states $(z^{(n)})_{n \in \{1, \dots, N\}}$ and rain occurrences $(y^{(n)})_{n \in \{1, \dots, N\}}$.
$z^{(1)} \sim \eta$ ▷ Generate sequence of states

for $n = 2$ to $N$ do $z^{(n)} \sim (Q(z^{(n-1)}, \ell))_{\ell \in \{1, \dots, K\}}$

end for

for $n = 1$ to $N$ do $x_Y^{(n)} \sim N(0, \exp(-h/\rho_{CY,z^{(n)}}^{(n)}))$ ▷ Generate Gaussian

$\forall s, y_s^{(n)} = \begin{cases} 1 & \text{if } x_{Y_s}^{(n)} \le \Phi^{-1}(\lambda_{z^{(n)},s, y_s^{(n-1)}}^{(n)}) \\ 0 & \text{else} \end{cases}$ ▷ Truncate to 0 and 1

end for

Return: $(z^{(n)})_{n \in \{1, \dots, N\}}, (y^{(n)})_{n \in \{1, \dots, N\}}$

simulate from the truncated spatio-temporal Gaussian field that is used in the rain intensity $X_R(s, n)$. For this purpose, an approximation of the spatio-temporal Gaussian simulation at time $n$ is used, which depends only on $\ell$ previous days, with $\ell$ chosen large enough to capture the temporal dependence:

$$ X_R(\cdot, n) | X_R(\cdot, n-1), \dots, X_R(\cdot, 1) \approx X_R(\cdot, n) | X_R(\cdot, n-1), \dots, X_R(\cdot, n-\ell). \quad (29) $$

As the vectors $X_R(\cdot, n), \dots, X_R(\cdot, n-\ell)$ are jointly Gaussian, and due to the stationary assumption, $X_R(\cdot, n) | X_R(\cdot, n-1), \dots, X_R(\cdot, n-\ell)$ can be generated as a Gaussian vector with conditional mean $\mu^{(n)}$ and conditional covariance matrix $\Sigma^{(n)}$, computed from the covariance model and the previous $\ell$ days simulations. This is combined with the truncation conditions of Eq.(8) in algorithm 2. The last step is to transform the latent truncated Gaussian field $x_R$ into zero precipitation and non-zero amounts using the marginal parameters according to algorithm 3.

## 2.4 Validation

A set of evaluation indicators is considered to assess model performance. These include indicators closely linked to the model structure, such as correlations and rain occurrence probabilities, as well
**Algorithm 2** Iterative Simulation of the Spatio-Temporal Gaussian Field $X_R | Y, Z$.

**Input:** Covariance function $C_R(h, u)$, temporal lag $\ell$, occurrences $(y^{(n)})_{n \in \{1, \dots, N\}}$, states $(z_{n \in \{1, \dots, N\}}^{(n)})$, dry probabilities $1 - \lambda_{z^{(n)}, s, y_s^{(n-1)}}^{(n)}$

**Output:** Simulated $x_R \in \mathbb{R}^{D \times N}$

Compute truncation bounds $l_n, u_n$ according to Eq. (8).

**Initialisation:** Draw

$$x_R(\cdot, 1) | y^{(1)}, z^{(1)} \sim \text{TMVN}(0, C_R(\cdot, 0), l_1, u_1)$$

**for** $n = 2$ to $\ell$ **do**

Get $\mu^{(n)}, \Sigma^{n}$, conditioned on **all** previous days information.

Sample

$$x_R(\cdot, n) | y^{(n)}, z^{(n)} \sim \text{TMVN}(\mu^{(n)}, \Sigma^{(n)}, l_n, u_n)$$

**end for**

**for** $n = \ell + 1$ to $N$ **do**

Get $\mu^{(n)}, \Sigma^{n}$, conditioned on **the** $\ell$ previous days information.

Sample

$$x_R(\cdot, n) | y^{(n)}, z^{(n)} \sim \text{TMVN}(\mu^{(n)}, \Sigma^{(n)}, l_n, u_n)$$

**end for**

**Return:** $(x_R(s, n))_{s,n}$

**Algorithm 3** Transformation into precipitation amounts $R | X_R, Z$

**Input:** Latent Gaussian field $x_R(s, n)$, dry probabilities $1 - \lambda_{z^{(n)}, s, y_s^{(n-1)}}^{(n)}$, intensity CDFs $F_{z^{(n)}, s}^{(t_n)}$

**Output:** Simulated precipitation field $r \in \mathbb{R}^{D \times N}$

**for** each site $s$ and time step $n$ **do**

if $x_R(s, n) > T(s, n)$ **then** ▷ Equivalent to $y_s^{(n)} = 1$

$$\text{Transform to precipitation intensity } r_s^{(n)} = (F_{z^{(n)}, s}^{(t_n)})^{-1} (\Phi_{T(s,n)}(x_R(s, n)))$$

else

$$r_s^{(n)} = 0$$

end if

**end for**

**Return:** $(r_s^{(n)})_{s,n}$
as more complex, model-agnostic metrics designed to assess large-scale events. These complementary indicators are then used to assess whether the generated rainfall series reproduce the targeted properties.

### 2.4.1 Marginal metrics

**Local dry (wet) episodes** For a location $s$, a dry spell of duration $m$ starting on date $n$ is a sequence fulfilling the conditions

$$Y_s^{(n-1)} = 1, Y_s^{(n)} = 0, \dots, Y_s^{(n+m)} = 0, Y_s^{(n+m+1)} = 1. \quad (30)$$

The empirical distribution of the dry spell duration $m$ computed from the observations is then compared with an envelope of empirical distributions obtained from the simulations.

**Rainfall totals** To assess the ability of the generator to reproduce seasonal rainfall characteristics, the cumulative monthly precipitation is computed for each location $s \in \mathcal{D}$, month $m$ as

$$R_{s,m}^{(\text{month})} = \sum_{n \in \mathcal{M}(m)} R_s^{(n)}, \quad (31)$$

where $\mathcal{M}(m)$ is the set of all days belonging to the month $m$. Compared to daily rainfall, the monthly accumulated precipitation provides a more robust and interpretable measure of model performance, since it integrates both the frequency of wet and dry days and the variability in rain intensity. This aggregated indicator thus captures the combined effect of occurrence and intensity processes, which are often difficult to assess from only daily data. For each location and month, the empirical distribution of $R_{s,m}^{(\text{month})}$ in years is summarised using quantiles for both the observations and the ensemble of simulations.

### 2.4.2 Spatio-temporal metrics definitions

**Spatial continuity** Locations that are close in space tend to have closer values of rain intensity: when one location records no rain, neighbouring locations with non-zero values are unlikely to be experiencing extreme precipitation. This is evaluated by the Continuity Ratio (CR), introduced in D. S. Wilks 1998 and defined for each pair of locations ($s_1, s_2$) as

$$CR(s_1, s_2) = \frac{\mathbb{E}[R_{s_1} | R_{s_1} > 0, R_{s_2} = 0]}{\mathbb{E}[R_{s_1} | R_{s_1} > 0, R_{s_2} > 0]} \quad (32)$$

Locations close in space will tend to have a small $CR(s_1, s_2)$ while far apart locations will have a value close to 1. Reproducing this ratio can be challenging, see Z. Li 2014; Zhou, Y. Meng, and Abbaspour 2019.

**Rain Occurrence Ratio** An objective of the generator is to accurately reproduce large-scale spatial events, such as widespread dry or wet days. To evaluate this property, inspiration is taken from the integrated mean temperature of Miloshevich et al. 2023, and an indicator called the Rain Occurrence Ratio (ROR) is proposed (see Gobet, Métivier, and Parey 2025). It is defined as a function of the observed or simulated rain occurrence $Y$ for each time step $n$ as

$$ROR(n) = \frac{1}{D} \sum_{s=1}^{D} Y_s^{(n)} \quad (33)$$

A low (high) value of $ROR(n)$ denotes a dry (wet) day for many locations at the same time. As a time series, its marginal distribution and autocorrelation function can be studied for the simulations and compared to those obtained for the observations.

To assess the model's ability to reproduce persistent dry and wet conditions, the durations of threshold-based events in the spatially averaged rainfall ratio ($ROR(n)$)$_{n \in \{1, \dots, N\}}$ are analysed. A dry event of length $m$ starting at date $n$ is defined as a consecutive sequence of time steps during which the rainfall ratio remains below a fixed low threshold $q_{\text{dry}}$.
The distribution of such run lengths provides a measure of the generator's ability to produce extended drought episodes, both in time and across space. Similarly, a wet event is defined as a consecutive sequence of exceedances above a high threshold $q_{\text{wet}}$. By comparing the empirical distributions of dry and wet episode lengths between simulations and observations, the ability of the stochastic generator to reproduce large-scale spatio-temporal droughts and wet spells consistent with the observed climatology can be evaluated.

**Quantile Exceedance Ratio** The Quantile Exceedance Ratio (QER) is defined similarly to assess large spatial events with high rainfall. It is defined as a function of the observed or simulated rain value $R$ for each time step $n$ as

$$QER_{\alpha}(n) = \frac{1}{D} \sum_{s=1}^{D} 1_{R(s,n)>q_{\alpha}(s)}. \qquad (34)$$

For a high quantile $q_{\alpha}$, reproducing this indicator well in the simulation means the generator captures effectively the spatial dependence of high rain.

## 2.5 Implementation

Likelihood computations in the E-step of the EM and hybrid algorithm (Eq. (4)) require evaluation of multivariate normal integrals, while the hybrid EM M-step for the pairwise likelihood only involves bivariate integrals. The multivariate integrals are computed using the algorithm of (Genz 1992), which provides accurate quasi-Monte Carlo estimates with controllable error, and the bivariate integrals use an analytic approximation (Tsay and Ke 2023) that is accurate and differentiable. Simulation of the truncated Gaussian field relies on a quasi-Monte Carlo method with minimax tilting of (Z. I. Botev 2017). The implementation is in Julia (Bezanson et al. 2017), see the supplementary material for more details. It is available on Github at https://github.com/caroline-cognot/HMM_Spatial_codepaper/.

# 3 Application to French precipitation data

## 3.1 Data

Most spatial generators apply to a small region with homogeneous climates. Obakrim, Benoit, and Allard 2025 work on the South-East of France, Evin, Favre, and Hingray 2018 on Switzerland, Dawkins et al. 2022 on Eastern England. Ailliot, Thompson, and Thomson 2009 focuses on New Zealand, which only has one climate type (oceanic). The generator proposed in this work is intended to handle a widespread network of stations across the country. It is applied to a set of $D = 37$ French weather stations from the European Climate Assessment & Dataset (ECA&D, Klein Tank et al. 2002). The selected stations have 100% valid data from 1 January 1973 to 31 December 2023, ie, $N = 18627$ days. 6 stations have been selected to represent a diversity of climatic conditions, in order to illustrate the marginal rainfall properties. Station 737 (Lille-Lesquin), located in the North, has a very rainy climate, even in the summer. Station 32 (Bourges), located in the center of France, has a mild climate, and all seasons are equally wet. Station 757 (Nice), located in the South-East, has a warm Mediterranean climate, with a strong contrast in the precipitation that is more important in the winter than in the summer. Location 755 (Embrun), the only mountainous station in the dataset, has a cool climate, with precipitation amounts that are high even in the summer. Station 323 (Strasbourg), located in the north-east, has a continental climate. Station 34 (Bordeaux) is in the south-west with an oceanic climate. Additional information on the dataset can be found in the Supplementary material.

## 3.2 Occurrence model

### 3.2.1 Hyperparameter selection

Two hyperparameters are to be determined to fit the occurrence model described in section 2.1.1: the number of classes $K$ and the degree of the trigonometric polynomials $d$. The Integrated Complete
<figure><img src="image_3.png" /><figcaption>Figure 2: Yearly mean of fitted rain probability parameters $\lambda_{k,s,y}$ for each location $s$ in each of the states $k = 1, 2, 3, 4$. The first row corresponds to the probability of rain given a previous dry day, the second row to the probability of rain given a previous wet day. The selected 6 stations are highlighted by black circles.</figcaption></figure>

Likelihood (ICL) criterion is used for this purpose, as it is suitable for hidden Markov models and, more generally, for models with latent variables (Celeux and Durand 2007). Let $\theta_{K,d}$ be the set of model parameters for $K$ classes and trigonometric degree of $d$, $L_C$ be the complete likelihood used in the EM algorithm. Then

$$ \mathrm{ICL}(K, d) = \log(L_C(\mathbf{y}^{(1:N)}, \hat{z}^{(1:N)}; \theta_{K,d})) - |\theta_{K,d}| \log(N)/2 \quad (35) $$

$\hat{z}^{(1:N)}$ are the most likely values for the hidden variable $z$, obtained using the Viterbi algorithm.

Table 1 shows the values of ICL as a function of $K$ for several values of the degree of seasonality $d$. The criterion favours the choice $K = 4, d = 1$, which is selected as the best model.

### 3.2.2 Rain occurrence parameters

Figure 2 represents the yearly average rain probabilities $\lambda_{k,s,y}^{(t_n)}$ for all the stations according to the four hidden states identified by the HMM and given the occurrence of the previous day. The four states show clear different patterns : the first one displays high rain probabilities in the whole country, in the second one the probabilities are much higher in the North than in the South, in the third one the probabilities are slightly lower in the North than in the South, and in the fourth state, probabilities are low in the whole country. The six stations selected for illustration are highlighted by black circles. Figure 3 shows their rain probability pattern over the year. In most locations and states, the probability of rain conditionally to a rainy previous day (dashed line) is, as expected, significantly higher than the one for a dry previous day (solid line). The yearly evolution of the spatial scale parameter, $\rho_{CY,k}^{(t)}$, is represented in Figure 4(A) together with the spatial scale estimated in the case of a unique state. When four states are identified, the spatial scale parameter is approximately 250 km. When there is no hidden state, the spatial scale averages 400 km, which means the hidden states captured some large-scale spatial correlation while the conditional covariance function captures short-range correlations. Figure 4(B) shows the evolution over the year of the transition parameters $Q^{(t)}(i, j)$. It is worth noticing that states 1 (high rain probabilities everywhere) and 4 (low rain probabilities everywhere) are much more stable than states 2 and 3 (different probabilities in the south and the north of the country).

## 3.3 Intensity model

### 3.3.1 Marginal rain intensity parameters

The estimated marginal rain parameters $p_{k,s}$, $\sigma_{k,s}$, $\xi_{k,s}$, $\kappa_{k,s}$ (see section 2.1.2) are shown in figure 5 for Lille and Nice, along with the parameters estimated in the case of a unique state. The parameters
<figure><img src="image_4.png" /><figcaption>Figure 3: Yearly evolution of fitted rain probability parameters $\lambda_{k,s,y}^{(t_n)}$ for selected locations. The full line represents the probability of rain given a previous dry day, while the dashed line corresponds to the probability of rain given a previous wet day. Each colour denotes a hidden state.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4: (A) Yearly evolution of the spatial scale parameter $\rho_{CY,k}^{(t)}$. Each colour corresponds to a hidden state, consistently with figure 3. The dotted line represents the fitted values for a generator with no hidden state, i.e, $K = 1$. (B) : Yearly evolution of the transition probabilities between states. Each of the four graphs shows the transitions from one state to all others, with colours indicating the arrival state, consistent with Figure 3.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 5: Yearly evolution of fitted rain intensity parameters for two of the selected locations. $p_{k,s}^{(t)}$ is the probability of low rain, $\sigma_{k,s}^{(t)}$, $\kappa_{k,s}^{(t)}$, $\xi_{k,s}^{(t)}$ are the EGPD parameters as defined in 2.1.2. The coloured lines represent the parameters inferred in each of the states determined by the Hidden Markov Model. The black lines correspond to the parameters estimated in the absence of hidden states.</figcaption></figure>

for the other locations are shown in the supplementary material. For all stations, the parameters estimated when considering only a single state are close to those of the wettest state.

The low rain probability $p_{k,s}^{(t_n)}$ is lowest in the wettest state 1, indicating that this state is characterised by higher rainfall than the others, both in terms of occurrence and intensity. The drier states 3 and 4 exhibit a high probability of low rainfall: not only is the probability of rain low, but even when it occurs, the amounts are small. Nice displays a different pattern from the other stations: the low rain probability in state 3 is as low as in state 1, indicating that state 3 represents a specific rainy regime for this Mediterranean location.

The parameter $\sigma_{k,s}^{(t)}$ representing the spread of the rain intensity distribution shows a strong seasonal aspect that depends on the location: it is higher in the summer for locations like Lille, Bourges, and Strasbourg, in the North and Center of France. It is higher in the autumn for Embrun and Nice, which are in the South-East.

$\kappa_{k,s}^{(t)}$ also represents part of the spread, localised at the origin of the distribution. It represents the deformation of the GPD (extreme distribution) in EGPD (full distribution). It is penalised to be close to 1, which makes it vary in the year near this value, except for Nice in state 2 where it is mostly 0.

The shape parameter $\xi_{k,s}^{(t)}$ is close to either one or zero in the driest state (state 4) for Lille, Bourges, Bordeaux, and Strasbourg due to the scarcity of rainy days in this state, which leads the optimisation procedure to produce values close to the boundary. For relatively dry locations such as Nice or Embrun, this effect occurs even for state 3. As a result, it is mostly equal to zero, implying a simpler exponential distribution.

### 3.3.2 Covariance parameters

The estimated spatio-temporal covariance parameters of the truncated Gaussian process in Equation (11), as defined by the covariance function in Equation (12), are reported in Table 2.

The variance parameter is set to 1. The estimated value of the spatial regularity parameter $\nu$ is close to 0.5, which corresponds to an exponential kernel. The estimation of the spatial scale $r$ is very close to the average spatial scale in the precipitation occurrence. The separability parameter $\beta$ is significantly different from zero, indicating clear non-separability. The temporal scale $\tau$ is very low, which means a fast-decaying correlation in time that can be observed on the curves. In particular, after $\ell = 2$ days, the temporal correlation vanishes, which validates the use of the conditional simulation in algorithm 2. It is implemented with $\ell = 10$.

## 3.4 Validation

To validate the generator, 500 rainfall simulations are produced and compared with the observations using the indicators presented in section 2.4. Additional indicators are also presented in the supplementary material.
<figure><img src="image_7.png" /><figcaption>Figure 6: Distribution of temporal length of dry spells for the observations (blue line) and for 500 simulations (orange and grey envelopes) for 6 locations.</figcaption></figure>

### 3.4.1 Distribution properties at each location

Figure 6 shows the distribution of local dry episode lengths as defined in Equation (30) for the 6 chosen locations. The red envelope represents the 25 to 75 % quantiles of the distributions in the simulation. The observations tend to be in the red envelope or to have slightly longer episodes. The gray envelope, covering all simulations, shows that simulations are able to generate longer events than those observed. Similar results hold for simulated wet spells (see Supplementary Material), demonstrating that the generator accurately reproduces local wet and dry events, including those lasting more than one day.

Figure 7 shows cumulative monthly rain characteristics as defined in Equation (31) for several locations for all years. For each location and month, the empirical distribution of $R_{s,m,y}^{(month)}$ across years is summarised using the 0.1, 0.5, and 0.9 quantiles for both the observations and the ensemble of simulations. The chosen locations exhibit quite different seasonal characteristics, with different interannual variability. Strasbourg station shows more cumulative rain in the summer than in the winter, whereas Bordeaux shows the contrary. Lille, Embrun, and Bourges, despite having very different wet/dry characteristics (see figure 3), all exhibit a relatively low variation during the year, while Nice has a very strong seasonal pattern, in particular at the summer-autumn transition in September - October. This shows that the periodic EGPD model has been able to capture well both the seasonal characteristics and the location-specific behaviour of rain intensity. Similar results are found when investigating 5-days cumulative rain, see Supplementary material.

### 3.4.2 Spatial continuity

Figure 8 shows different indicators computed on the rain simulations and compared to the rain observations: lag-0 cross-covariance between pairs of locations, tail index, and continuity ratio as defined in Eq. (32). All indicators show a good fit between simulations and observations. In particular, the continuity ratio is very well reproduced, which validates the joint occurrence/intensity simulation strategy.

### 3.4.3 Dry/wet large-scale events

For the large-scale indicators evaluated on the occurrence model, the generator with four hidden states, spatial correlation and marginal memory, referred to as HMM1_Spa_K4, is compared to several simpler alternatives in order to show the contribution of each component of the generator. First, the generator of Gobet, Métivier, and Parey 2025 under conditional independence, using the same hyperparameters ($K = 4, d = 1$) and marginal memory, is referred to as HMM1_Ind_K4. Comparing its properties with those of the new generator emphasises the role of spatial covariance in the occurrence process, relative to the assumption of conditional independence. Next, the usefulness of the hidden states is assessed.
<figure><img src="image_8.png" /><figcaption>Figure 7: Monthly cumulative rain for selected locations. 0.1 quantile in orange, 0.5 quantile in green, 0.9 quantile in blue. Each simulated envelope includes a central ribbon corresponding to the interquartile range (0.25–0.75 quantiles) and an outer ribbon showing the full range of the simulations (minimum–maximum).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 8: From left to right: lag-0 cross-covariance between pairs of locations, tail index ($q = 0.95$), and continuity ratio as defined in Eq. (32).</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 9: Distribution (first row) and auto-correlation function (second row) of the ROR indicator from Eq. (33). The values in the observations are the blue dots, and each coloured ribbon corresponds to the min-max envelope from simulation of the generators HMM1_Ind_K4, HMM1_Spa_K1, HMM1_Spa_K4.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 10: Distribution of length in days of dry spatial episodes, defined as ROR($n$) < 0.2 for the observations and 4 generators. The observed values are the blue curves, and each coloured ribbon corresponds to the min-max envelope from simulation of the generators HMM0_Spa_K4, HMM1_Ind_K4, HMM1_Spa_K1, HMM1_Spa_K4.</figcaption></figure>

using the generator HMM1_Spa_K1, which corresponds to a spatially correlated distribution without hidden states. This generator is conceptually close to D. S. Wilks 1998, although here a parametric covariance function is used, while the original generator relied on an empirical covariance matrix. This version contains spatial dependence in the emissions and a one-day local memory, but lacks additional temporal structure provided by hidden states. Finally, to examine the contribution of the marginal memory, it is removed in a generator denoted HMM0_Spa_K4.

Figure 9 shows the distribution and auto-correlation function of the ROR indicator from Eq. (33) for the observations (blue dots) and the three generators with marginal memory; the generator HMM0_Spa_K4 is not shown, in order to avoid overcrowding the graph. HMM1_Ind_K4 (gray) clusters the observations into “very wet”, “very dry”, and “intermediate” days (North wet, South dry, and vice versa). However, in the absence of the spatial component, the generator primarily produces very dry, mild, or very wet days across the 37 locations, resulting in the top-plot curves with three modes that do not reflect reality. Remarkably, the correlation function remains close to the observed one (bottom plots), indicating that the temporal structure is still well captured. HMM1_Spa_K1 (green) has a limited temporal dependence and does not separate the complexity of the model into weather types. As shown on the correlation function of the spatial ROR (bottom plots), this is not enough to capture the temporal dependence. The generator also struggles at representing very wet or very dry events, mostly generating days with ROR = 0.5 which are not representative of reality (top plots). Finally, the generator with spatial emission and marginal memory, HMM1_Spa_K4 (red), shows a good fit to the observations, though slightly underestimating the temporal correlation in summer.

Figure 10 shows the distribution of dry episode lengths for the observations and the 4 generators for a threshold that defined a "dry day" set to 20%. As was observed previously, generators that lack spatial or temporal dependence are outperformed by the two spatial generators with four states HMM1_Spa_K4, HMM0_Spa_K4. This means both the temporal dependence in the states and the local memory are necessary to fully capture the large-scale drought events. This is also true when the threshold is decreased to 10% (see Supplementary material).
<figure><img src="image_12.png" /><figcaption>Figure 11: Distribution (first row) and auto-correlation function (second row) of QER$_{0.95}$. The simulations from the generator are represented by its min-max envelope in orange, and the observed values are the blue dots.</figcaption></figure>

### 3.4.4 Joint quantile exceedance

Figure 11 shows the distribution and temporal correlations of the Quantile Exceedance Ratio $QER_\alpha$ (Eq. (34)) for $\alpha = 0.95$, with observed values denoted by blue circles and simulated values by orange ribbons. The case $QER_\alpha = 0$ is the most common (no exceedance), and the probability of exactly one location experiencing heavy rainfall above its quantile $\alpha = 0.95$ (QER close to 0.02) is approximately 0.15 for all seasons (top plots). For exactly 2 locations (QER close to 0.05), this probability is approximately 0.1. This is well represented by the simulations. The generator slightly underestimates the temporal correlation (bottom plots), as it does for other indicators.

## 4 Discussion and conclusion

The stochastic weather generator developed in this work focuses on precipitation and is designed to reproduce spatial and temporal variability over a large, climatically diverse region. The model separates rainfall into occurrence and intensity components, which are linked through a truncated latent Gaussian framework. This decomposition allows each component to be modelled with appropriate statistical structures: a Hidden Markov Model with spatially dependent emissions for the occurrence process, yielding interpretable weather regimes, and an Extended Generalised Pareto model for the intensity, conditional on the truncated Gaussian field. The joint simulation procedure combines both processes, ensuring spatial continuity and realistic temporal dynamics.

Overall, the simulation results demonstrate that the proposed spatio-temporal generator provides a realistic and coherent reproduction of rainfall statistics across temporal and spatial scales. At the local level, the marginal distributions of wet and dry spell lengths and the seasonal variability of cumulative rainfall are well captured. The generator reproduces both the persistence of wet and dry events and the strong regional contrasts between continental, oceanic, and Mediterranean climates. These results confirm the ability of the periodic EGPD parameterisation to account for location-specific seasonality in rain intensity distributions. The auto-correlation analysis of rainfall quantities shows that the temporal dependence is generally well reproduced, although slightly underestimated for some stations, particularly during summer. This suggests that the Gaussian latent process captures short-term persistence well but may benefit from a richer temporal structure or non-stationary covariance to better represent long multi-day events. Pairwise indicators confirm that the joint occurrence and intensity dependencies between stations are realistically simulated. The close alignment between observed and simulated joint probabilities and the accurate reproduction of continuity ratios highlight that the coupling between the occurrence process $Y$ and the intensity process $R$ is consistent and effective. This represents a clear improvement over earlier generators that did not explicitly link rainfall occurrence and amounts. At larger spatial scales, the generator successfully captures the dynamics of widespread wet and dry episodes. Incorporating both the hidden Markov model and the spatial covariance structure enables a realistic representation of large-scale coherence, with the proposed generator clearly outperforming simpler or spatially independent versions for the generation of dry and wet spells. The
use of a thresholded rain occurrence ratio ROR₁₀ further shows that the generator correctly represents the transition between light rain and dry conditions, which is crucial for drought or low-rainfall impact studies. Finally, the joint quantile exceedance analysis demonstrates that the generator reproduces the frequency and clustering of heavy-rainfall events across sites. These large-scale properties make this generator a promising tool for climate impact studies and hydrological applications.

Several strengths of the approach can be emphasized. First, the separation between occurrence and intensity modelling enhances interpretability and allows the use of specialised distributions for each process, while the new conditional simulation procedure couples accurately both processes. Second, the latent Gaussian formulation in both the occurrence and intensity processes ensures scalability, as spatial dependence is controlled through a small number of covariance parameters. Third, an efficient modified version of the Expectation-Maximisation (EM) algorithm, the hybrid EM, has been developed to estimate the HMM parameters with spatial conditional dependence, using appropriate approximations of the *D*-dimensional multivariate normal CDF at both E and hybrid-M steps. All inferences were performed using automatic differentiation, improving both estimation speed and accuracy compared with gradient-free or approximate gradient methods. Finally, the seasonal parameterisation provides smooth transitions across the year, avoiding discontinuities often observed when data are split by month or season.

Some improvements could be worth making. The generator tends to underestimate the probability of completely dry days across all stations, possibly due to the limited spatial range of the dependence structure or the influence of large-scale atmospheric phenomena not captured by the current covariance model. This can maybe be made by enforcing the existence of a "completely dry" state; however, this does not suit the philosophy of this study, where the hidden states were supposed to be completely learned from the data with no prior knowledge. Studying such very dry spells can be done separately; however, the strength of stochastic weather generators, such as the one in this study, is to be able to represent the full range of possibilities. The temporal dependence may also be improved by adding more previous days' memory in the rain probability. Future work will focus on extending this framework to multivariate weather generation. The latent Gaussian structure is naturally compatible with additional variables, such as temperature or wind speed. Incorporating temperature would open the way to hydrological applications, where simulated precipitation and temperature could be combined within hydrological models to assess streamflow and drought risks under various climate scenarios.

<table><caption>Table 1: ICL values for several *K*, *d*, models for the rain occurrence fitted on French station data. The best value for each degree *d* is in italics, and the best overall value is bolded.</caption><thead><tr><th>ICL value</th><th>K = 1</th><th>K = 2</th><th>K = 3</th><th>K = 4</th><th>K = 5</th></tr></thead><tbody><tr><td>d = 0</td><td>-3.12255e5</td><td>-2.95841e5</td><td>-2.94279e5</td><td>-2.92816e5</td><td>-2.929e5</td></tr><tr><td>d = 1</td><td>-3.09932e5</td><td>-2.94661e5</td><td>-2.93450e5</td><td><strong>-2.92643e5</strong></td><td>-2.93197e5</td></tr><tr><td>d = 2</td><td>-3.09927e5</td><td>-2.95501e5</td><td>-2.95018e5</td><td>-2.95164e5</td><td>-2.96406e5</td></tr></tbody></table>

<table><caption>Table 2: Estimated covariance model parameters for the covariance model in Equation (12))</caption><thead><tr><th>Parameter</th><th>β</th><th>r</th><th>ν</th><th>τ</th><th>α</th><th>δ</th></tr></thead><tbody><tr><td>Value</td><td>0.74</td><td>274.55</td><td>0.48</td><td>0.47</td><td>1.00</td><td>0.42</td></tr></tbody></table>

## Acknowledgements

Caroline Cognot was funded by Agence Nationale de la Recherche et de la Technologie (grant no. 2023/0061).
# A Additional information on the statistical model

## A.1 Illustration of the correlated binary-valued distribution

The emission distribution in the HMM model presented in the main article relies on a correlated binary occurrence vector $Y = (Y_1, ..., Y_D)$, which is generated in two steps by

$$X_{Y_1}, ..., X_{Y_D} \sim \mathcal{N}(0, C_Y), \quad \forall s, Y_s = \begin{cases} 1 & \text{if } X_{Y_s} \le \Phi^{-1}(\lambda_s) \\ 0 & \text{else} \end{cases}, \qquad (36)$$

with $\Phi$ the cumulative distribution function (CDF) of the standard normal distribution.

<figure><img src="image_13.png" /><figcaption>Figure 12: Generations of a binary correlated field for constant values of $\lambda = \mathbb{P}(Y_s = 1)$ over space and different values of spatial scale $\rho$ in an exponential covariance function $C_Y(h) = \exp(-h/\rho)$.</figcaption></figure>

Figure 12 illustrates examples of such binary fields for constant threshold $\lambda$ and different values of the spatial scale parameter $\rho$ in the exponential covariance $C_Y(h) = \exp(-h/\rho)$. Black pixels represent $Y_s = 1$ (rain) and white pixels represent $Y_s = 0$ (no rain). For very small $\rho$ (e.g., $\rho = 0.01$), the field resembles independent Bernoulli draws. As $\rho$ increases, larger and smoother spatial patterns emerge.

## A.2 Overview of the model parameters

Figure 13 summarises the different parameters involved in the generator and presented in the model definition section, together with the processes to which they are related. The key and novel point compared to other strategies that model occurrences and intensity independently, is that both processes are related according to the rain probabilities $\lambda_{k,s,y_s}^{(n)}$.

# B Details on the computational choices

## B.1 Integral computation methods

Computations of the likelihood $L_Y$ within the E-step of the EM algorithm require efficient evaluation of multivariate normal integrals, whereas the maximisation of the pairwise likelihood $L_{Y,wp}$ in the M-step only requires bivariate integrals. These bivariate integrals must be computed efficiently enough to allow optimisation and, ideally, automatic differentiation.

For the multivariate case, the algorithm of Genz 1992 is employed. This method uses efficient Quasi-Monte Carlo sampling and variable transformations to achieve high accuracy and provides a parameter to control the approximation error through the number of Monte Carlo samples. More
<figure><img src="image_14.png" /><figcaption>Figure 13: Outline of the model parameters, grouped into occurrence (HMM), intensity marginals and dependence parameters. The parameters are linked to the variables they influence.</figcaption></figure>

recent approaches approximate the target distribution using bivariate or trivariate conditioning (see, for instance, Bhat 2018 and references therein), but they lack an explicit error-control parameter and tend to produce larger errors for highly correlated data, making them unsuitable for the present application. Another alternative could have been to adopt block covariance decomposition methods, such as those proposed by Genton, Keyes, and Turkiyyah 2018, which are designed to efficiently compute the multivariate CDF for very high-dimensional problems ($D \sim 10^4$). However, this approach was not implemented because the current setting involves a moderate spatial dimension with numerous temporal points rather than extremely high-dimensional spatial data.

In the bivariate case, Lin 1995; Genz 2004; Tsay and Ke 2023 have developed efficient analytic approximations of the integral, with the most recent approach offering the best performance: it is accurate both in the central region and in the tails and often matches or outperforms existing methods. The approximation error remains small even for high correlations, and the resulting expression is easily differentiable.

## B.2 Integral simplifications

### B.2.1 Reformulation of integral problems into CDF of multivariate normal

The integrals involved in the likelihoods of the HMM have one finite and one infinite bound in each dimension, depending on the observed value of precipitation occurrence. They can be rewritten as the CDF of a multivariate normal distribution to take advantage of efficient methods.

**Bivariate case** In the bivariate case, any integral of the pdf $f_\rho$ of normal $N\left(\begin{bmatrix} 0 \\ 0 \end{bmatrix}, \begin{bmatrix} 1 & \rho \\ \rho & 1 \end{bmatrix}\right)$ can be written as a function of the bivariate normal cumulative distribution function $\Phi_2$ and the univariate normal cumulative distribution function $\Phi$

$$ \int_{-\infty}^{a} \int_{-\infty}^{b} f_{\rho} = \Phi_{2}(a, b, \rho) \qquad (37) $$

$$ \int_{a}^{\infty} \int_{-\infty}^{b} f_{\rho} = \Phi(a) - \Phi_{2}(a, b, \rho) \qquad (38) $$

**Multivariate case** Let $R(B) = \prod_{i=1}^{D} I_i$ with $I_i = (-\infty, b_i)$ or $(b_i, \infty)$ the integration region in the multivariate case, and $f_\Sigma$ the pdf of multivariate $N(0, \Sigma)$.
Let

$$D = \operatorname{diag}(d_1, \dots, d_D), \quad d_i = \begin{cases} -1, & \text{if the bounds are } (b_i, \infty), \\ 1, & \text{if the bounds are } (-\infty, b_i), \end{cases}$$

and define

$$\tilde{b}_i = d_i b_i, \quad \tilde{\Sigma} = D \Sigma D.$$

Then the mixed-bound integral can be rewritten as

$$\int_{R(B)} f_{\Sigma}(u) du = \int_{-\infty}^{\tilde{b}_1} \dots \int_{-\infty}^{\tilde{b}_D} f_{\tilde{\Sigma}}(v) dv = \Phi_{\tilde{\Sigma}}(\tilde{b}_1, \dots, \tilde{b}_D),$$

### B.2.2 Simplifications in the pairwise maximisation

The sum to maximise with respect to the parameters during the M-step is

$$\tilde{J} = \sum_{n=1}^{N} w_n \sum_{i<j} w_{ij} \log L_{Y2}((y_i^{(n)}, y_j^{(n)}); (\lambda_i^{(n)}, \lambda_j^{(n)}), \rho_{CY}^{(n)}). \quad (39)$$

As the parameters $\lambda, \rho_{CY}$ are periodic, the sum can be simplified as

$$\tilde{J} = \sum_{t=1}^{T} \sum_{n, t_n=t} w_n \sum_{i<j} w_{ij} \log L_{Y2}((y_i^{(n)}, y_j^{(n)}); (\lambda_i^{(t)}, \lambda_j^{(t)}), \rho_{CY}^{(t)}). \quad (40)$$

**No marginal memory** In case where there is no marginal memory; $\lambda_i^{(n)} = \lambda_i^{(t_n)}$. There are only 4 different possible values for the pair $(y_i^{(n)}, y_j^{(n)})$: $P_1 = (0, 0)$, $P_2 = (0, 1)$, $P_3 = (1, 0)$ and $P_4 = (1, 1)$. Then $\tilde{J}$ can be written as

$$\tilde{J} = \sum_{p=1}^{4} \sum_{t=1}^{T} \sum_{n, t_n=t, (y_i^{(n)}, y_j^{(n)})=P_p} w_n \sum_{i<j} w_{ij} \log L_{Y2}(P_p; (\lambda_i(t), \lambda_j(t)), \rho_{CY}^{(t)}). \quad (41)$$

This reduces the computation of $ND(D-1)/2$ integrals to at most $4TD(D-1)/2$. Using the reformulation from Appendix B.2.1, only one of the 4 situations has to be computed, as the others can be easily obtained from marginal probabilities. In the end, only $TD(D-1)/2$ computations are necessary.

**With marginal memory** In the case with marginal memory, $\lambda_i^{(n)} = \lambda_{i, y_i^{(n-1)}}^{(t_n)}$. Similar simplifications allow to group observation pairs by 16 different possible values for $(y_i^{(n)}, y_j^{(n)}, y_i^{(n-1)}, y_j^{(n-1)})$.

This reduces the computation of $ND(D-1)/2$ integrals to at most $16TD(D-1)/2$. Using the reformulation from section B.2.1, only $4TD(D-1)/2$ are needed.

## B.3 Software

All implementation was performed in Julia (Bezanson et al. 2017). For estimation purposes, the automatic differentiation package Forwarddiff.jl (Revels, Lubin, and Papamarkou 2016) was used with the package Optimization.jl (Dixit and Rackauckas 2023) and the LBFGS solver. The autodifferentiation-compatible implementation of the Matérn kernel function of the BesselK.jl package (Geoga et al. 2023) was used. The integrals of the multivariate normal distribution and bivariate normal distribution are computed, respectively, using Julia package MvNormalCDF.jl for the Expectation step (with original code from Genz 1992) and a custom implementation of Tsay and Ke 2023 for the Maximisation step. The Truncated Normal simulation method that is used is the Julia package TruncatedMVN.jl (Z. Botev and Belzile 2015), which uses the recent method of Z. I. Botev 2017 rather than using acceptance-rejection techniques or Gibbs sampling, as implemented in the R package tmvtnorm (Wilhelm and G 2025). Acceptance-rejection techniques are not suitable for high dimensions, and Gibbs sampling also quickly becomes computationally expensive.
<table><caption>Table 3: Locations and characteristics of the 37 stations used in this study.</caption><thead><tr><th>ID</th><th>Name</th><th>Lat</th><th>Lon</th><th>Elevation (m)</th><th>Rain frequency</th><th>Max dry spell duration (days)</th><th>Max daily rainfall</th><th>Mean daily rainfall</th></tr></thead><tbody><tr><td>32</td><td>BOURGES</td><td>+47:03:33</td><td>+002:21:35</td><td>161</td><td>0.47</td><td>36</td><td>79.0</td><td>4.28</td></tr><tr><td>33</td><td>TOULOUSE</td><td>+43:37:16</td><td>+001:22:44</td><td>151</td><td>0.4</td><td>34</td><td>82.7</td><td>4.32</td></tr><tr><td>34</td><td>BORDEAUX</td><td>+44:49:50</td><td>-000:41:29</td><td>47</td><td>0.48</td><td>28</td><td>87.6</td><td>5.33</td></tr><tr><td>36</td><td>PERPIGNAN</td><td>+42:44:14</td><td>+002:52:22</td><td>42</td><td>0.25</td><td>50</td><td>222.0</td><td>5.87</td></tr><tr><td>736</td><td>ABBEVILLE</td><td>+50:08:10</td><td>+001:50:02</td><td>69</td><td>0.53</td><td>31</td><td>65.0</td><td>4.04</td></tr><tr><td>737</td><td>LILLE</td><td>+50:34:12</td><td>+003:05:51</td><td>47</td><td>0.51</td><td>33</td><td>62.8</td><td>3.93</td></tr><tr><td>738</td><td>CAEN</td><td>+49:10:48</td><td>-000:27:22</td><td>67</td><td>0.55</td><td>37</td><td>93.0</td><td>3.66</td></tr><tr><td>740</td><td>ALENCON</td><td>+48:26:44</td><td>+000:06:37</td><td>143</td><td>0.52</td><td>29</td><td>70.0</td><td>3.86</td></tr><tr><td>742</td><td>NANTES</td><td>+47:09:00</td><td>-001:36:32</td><td>26</td><td>0.5</td><td>35</td><td>96.7</td><td>4.44</td></tr><tr><td>745</td><td>DIJON</td><td>+47:16:04</td><td>+005:05:18</td><td>219</td><td>0.46</td><td>35</td><td>92.8</td><td>4.46</td></tr><tr><td>749</td><td>POITIERS</td><td>+46:35:38</td><td>+000:18:52</td><td>123</td><td>0.45</td><td>33</td><td>92.3</td><td>4.16</td></tr><tr><td>750</td><td>CLERMONT</td><td>+45:47:13</td><td>+003:08:58</td><td>331</td><td>0.4</td><td>37</td><td>76.0</td><td>3.97</td></tr><tr><td>786</td><td>MONTELIMAR</td><td>+44:34:52</td><td>+004:43:59</td><td>73</td><td>0.33</td><td>46</td><td>218.4</td><td>7.57</td></tr><tr><td>755</td><td>EMBRUN</td><td>+44:33:56</td><td>+006:30:08</td><td>871</td><td>0.31</td><td>51</td><td>67.2</td><td>6.55</td></tr><tr><td>756</td><td>TARBES</td><td>+43:11:17</td><td>+000:00:00</td><td>360</td><td>0.47</td><td>27</td><td>84.0</td><td>6.32</td></tr><tr><td>757</td><td>NICE</td><td>+43:38:56</td><td>+007:12:32</td><td>2</td><td>0.24</td><td>61</td><td>191.4</td><td>8.64</td></tr><tr><td>758</td><td>BASTIA</td><td>+42:32:26</td><td>+009:29:07</td><td>10</td><td>0.31</td><td>59</td><td>232.4</td><td>7.0</td></tr><tr><td>434</td><td>BREST</td><td>+48:26:39</td><td>-004:24:43</td><td>94</td><td>0.63</td><td>24</td><td>82.6</td><td>5.14</td></tr><tr><td>39</td><td>MARIGNANE</td><td>+43:26:16</td><td>+005:12:58</td><td>9</td><td>0.23</td><td>75</td><td>161.3</td><td>6.26</td></tr><tr><td>322</td><td>RENNES</td><td>+48:04:08</td><td>-001:44:02</td><td>36</td><td>0.5</td><td>29</td><td>70.4</td><td>3.77</td></tr><tr><td>323</td><td>STRASBOURG</td><td>+48:32:58</td><td>+007:38:25</td><td>150</td><td>0.47</td><td>32</td><td>66.3</td><td>3.68</td></tr><tr><td>2184</td><td>ROUEN</td><td>+49:22:59</td><td>+001:10:54</td><td>151</td><td>0.55</td><td>27</td><td>81.3</td><td>4.11</td></tr><tr><td>2190</td><td>TOURS</td><td>+47:26:40</td><td>+000:43:38</td><td>108</td><td>0.46</td><td>31</td><td>62.0</td><td>4.02</td></tr><tr><td>2192</td><td>BASEL</td><td>+47:36:52</td><td>+007:30:36</td><td>263</td><td>0.5</td><td>26</td><td>58.6</td><td>4.14</td></tr><tr><td>2195</td><td>LIMOGES</td><td>+45:51:40</td><td>+001:10:30</td><td>402</td><td>0.49</td><td>32</td><td>77.2</td><td>5.66</td></tr><tr><td>2199</td><td>GOURDON</td><td>+44:44:42</td><td>+001:23:48</td><td>260</td><td>0.46</td><td>33</td><td>102.6</td><td>5.05</td></tr><tr><td>2200</td><td>MILLAU</td><td>+44:07:07</td><td>+003:01:07</td><td>712</td><td>0.39</td><td>41</td><td>115.0</td><td>5.04</td></tr><tr><td>2203</td><td>MONT-DE-MARSAN</td><td>+43:54:35</td><td>-000:30:01</td><td>59</td><td>0.43</td><td>27</td><td>68.8</td><td>5.43</td></tr><tr><td>2205</td><td>ST-GIRONS</td><td>+43:00:19</td><td>+001:06:25</td><td>414</td><td>0.47</td><td>31</td><td>95.4</td><td>5.74</td></tr><tr><td>2207</td><td>MONTEPELLIER</td><td>+43:34:37</td><td>+003:57:47</td><td>2</td><td>0.26</td><td>45</td><td>299.5</td><td>6.65</td></tr><tr><td>2209</td><td>AJACCIO</td><td>+41:55:05</td><td>+008:47:34</td><td>5</td><td>0.28</td><td>89</td><td>147.6</td><td>6.17</td></tr><tr><td>11244</td><td>CHASSIRON</td><td>+46:02:49</td><td>-001:24:42</td><td>11</td><td>0.47</td><td>32</td><td>59.6</td><td>4.01</td></tr><tr><td>11245</td><td>PLOUMANAC'H</td><td>+48:49:33</td><td>-003:28:23</td><td>55</td><td>0.59</td><td>22</td><td>51.4</td><td>3.93</td></tr><tr><td>793</td><td>LA HAGUE</td><td>+49:43:31</td><td>-001:56:23</td><td>6</td><td>0.52</td><td>39</td><td>52.5</td><td>3.96</td></tr><tr><td>11246</td><td>NANCY</td><td>+48:34:52</td><td>+005:57:35</td><td>336</td><td>0.51</td><td>36</td><td>80.0</td><td>4.36</td></tr><tr><td>11247</td><td>BELLE ILE</td><td>+47:17:40</td><td>-003:13:06</td><td>34</td><td>0.49</td><td>32</td><td>69.1</td><td>3.8</td></tr><tr><td>11249</td><td>ORLY</td><td>+48:43:01</td><td>+002:23:04</td><td>89</td><td>0.46</td><td>35</td><td>75.4</td><td>3.68</td></tr></tbody></table>

## C Details on the French station data

Figure 14 (left) shows a map of the different locations, with the stations selected as illustrative examples highlighted in red. Their respective rain characteristics are shown in figure 14 (right). Table 3 gives for each location its ID, name, coordinates, elevation, and summary statistics related to their rain characteristics: rainy day frequency, maximum dry spell duration, maximum daily rainfall, mean daily rainfall. These statistics vary among locations, with locations close to the Mediterranean sea exhibiting longer dry spells but also much higher maximum rain. On the contrary, location in the North-West tend to have more frequent rain, with shorter dry spells, and more moderate maximum daily rainfall.

<figure><img src="image_15.png" /><figcaption>Figure 14: Left: map of the $D = 37$ weather stations used in this study. The six stations highlighted in red are used as illustrative examples in the validation section due to their contrasting climatic conditions. Right: mean daily rain amounts (blue) and daily rain probability (red) for each date of the year for the 6 locations chosen as examples.</figcaption></figure>
# D Details on the estimated parameters

## D.1 Interpretability of the hidden states

### D.1.1 State description

This is an additional description of the states found in the main paper.

The results of the HMM parameter estimation allow the hidden states to be interpreted as weather regimes. The first state corresponds to an overall wet regime, the second to wet conditions in the North and dry conditions in the South, the third to wetter conditions in the South than in the North, and the fourth to dry conditions everywhere.

The states can be ordered according to the rain probability, from state 1 with the highest probability to state 4 with the lowest. This hierarchy is verified for almost all stations, as shown by Figure 3 of the main text. The hierarchy of wettest to driest states is more pronounced for Lille (North, oceanic) than for the southern stations. For Nice in the South-East state 3 is as wet as state 1, with a high difference according to the previous day memory. In Embrun, a South mountainous location, the rainy state 1 and drier state 2 stay similar all year, but states 3 and 4 tend to be more rainy in the summer than in the winter. The opposite behaviour is observed in Lille, located in the north, where these states are relatively wetter in winter than in summer. For Strasbourg (North-East), Bordeaux (South-West), and Bourges (Central France), states 2 and 3 exhibit very similar moderate rain probabilities. Despite this similarity across locations, their seasonal modulation differs, with a more pronounced summer decrease in rain probability at Bordeaux. These two states appear to overlap for these locations, while they show a clear different behaviour (state 2 rainier than state 3) for the other stations.

Transition probabilities provide useful information on the inferred states (see Figure 4(B) in main text). States 1 and 4 are the most persistent, with a high probability of staying in the state. Moreover, it is almost impossible to go directly from state 1 (wettest) to state 4 (driest). States 2 and 3 act as transition states: from state 1 (or 4), there is a much higher transition rate to go to state 2 or 3 than to state 4 (or 1). From state 2, there is a relatively high probability to go to state 1 (and a lower probability, but non-zero, to go into 4). From state 3, it depends on the season: in spring, the most likely transition is into state 3; in autumn, states 2, 3, and 4 have similar probabilities.

The spatial scale also changes between states (see Figure 4(A) in main text). State 4 (the driest) has a lower estimated scale, suggesting this state has shorter-range events. In contrast, states 1 (wettest) and 2 (wet north) have a higher spatial scale, indicating larger-range events.

### D.1.2 Relationship with observed meteorological events

The four states can be illustrated by some remarkable meteorological events in the past. Figure 15 illustrates the most probable sequences of these states, colour-coded consistently with previous figures, for years 1976, 1999, 2003, 2005, and 2023. Major heatwaves are retrieved, as they correspond to dry events. Those occurring in August 2003 and during the summer of 1976 correspond to sustained episodes of the driest state (state 4), as anticipated. Both events were preceded by anomalously dry spring conditions, also captured in the generator by extended sequences of state 4. The year 1999 was "wet" according to the French climatology (see (SDES) 2024). It exhibits a marked dominance of the wettest state (state 1) and comparatively few occurrences of state 4. In contrast, 2005 is identified as a dry year and shows infrequent and short-lived wet episodes, consistent with the limited representation of state 1. In 2023 (Météo-France 2024), three periods show interesting behaviours. There was a precipitation deficit in the winter 2022-2023, with no rain between 21 January and 21 February, which is reflected by the persistence of state 4 over this period. In contrast, the end of the year experienced the longest rain episode since 1988, with a length only seen twice in the historical records, which is characterised by a very long episode of state 1. Perhaps more interestingly, the summer is characterised by a long period of the second state. This state expresses a North/South contrast, with wet north and dry south. Météo-France reported autumn-like conditions in the northern half of the country while record-breaking heat affected the south, consistent with the interpretation of this state.

## D.2 Rain parameters

Figure 16 shows the estimated rain marginal parameters for the remaining 4 out of 6 stations not shown in the main text.
<figure><img src="image_16.png" /><figcaption>Figure 15: Most likely a posteriori sequence of states for several years. The wettest state 1 is in blue, state 2 (rainy north, dry south) in green, state 3 (rainy south) is in purple, and state 4 (driest) is in orange.</figcaption></figure>

# E Additional indicators for the evaluation of the generator

## E.1 Local dry spells

Figure 17 shows the distribution of local wet episode lengths for the 6 chosen locations. The red envelope represents the 25 to 75 % quantiles of the distributions in the simulation. The observations tend to be in the red envelope or to have slightly longer episodes. The gray envelope, covering all simulations, shows that simulations are able to generate longer events than those observed. These results are shared with dry events (see main text), which shows that the generator correctly reproduces local wet/dry events, even for lengths longer than 1 day.

## E.2 5-day cumulative rain

To complement the monthly analysis, the cumulative precipitation over consecutive 5-day periods is also examined, providing insight into short-term rainfall persistence and temporal dependence. For each location s and each day of the year t, the 5-day cumulative rainfall is computed from the daily intensities $R_s^{(n)}$ as

$$R_{s,n}^{(5-day)} = \sum_{j=n-4}^{n} R_{s,j}^{(j)}, \qquad (42)$$

using a moving 5-day window over the full annual cycle. Compared to daily precipitation, this aggregated indicator smooths out day-to-day noise and reflects the occurrence of multi-day rainfall episodes, which are crucial for hydrological applications.
<figure><img src="image_17.png" /><figcaption>Figure 16: Yearly evolution of fitted rain intensity parameters for the 4 selected locations not shown in the main text. $p_{k,s}^{(t)}$ is the probability of low rain, $\sigma_{k,s}^{(t)}$, $\kappa_{k,s}^{(t)}$, $\xi_{k,s}^{(t)}$ are the EGPD parameters.</figcaption></figure>

Figure 18 shows the cumulative 5-day characteristics for each date of the year as defined in the main text. The blue ribbons represent the 90% quantile of that day of the year, green ribbons the median of that day, and orange ribbons the 10% quantile for that day, with dark blue lines representing the observed values. The low quantile is almost always 0 and is not visible on the plot. For locations like Nice or Embrun, the 5-day cumulative rain is mostly close to 0, with the high-quantile days being responsible for most of the monthly precipitation, denoting a low temporal dependence. For the other locations, the median 5-day precipitation is higher. In all cases, the simulation envelopes contain the observations.

## E.3 Temporal dependence in the simulated rain amounts

Figure 19 shows the auto-correlation function of the rain intensity for the same 6 locations for observations and simulations. The observations in blue lie mostly within the gray simulation envelope, although the generator tends to underestimate the temporal correlation, especially for Bordeaux or Nice.

## E.4 Probability of joint rain occurrence

Figure 20 shows joint probabilities of rain occurrence between pairs of stations, for the four seasons. Observed values on the x-axis are compared with simulated values on the y-axis and are found to be aligned on the 1:1 line, both for same-day probabilities and for lag-1 probabilities.

## E.5 Spatial dry spell durations

Figure 21 shows the distribution of dry episode lengths for the observations and the 4 occurrence models, for a dry day threshold set to 10%. The results are similar to the less extreme threshold of 20%, however, a small underestimation of the lengths of dry episodes can be seen in the summer.
<figure><img src="image_18.png" /><figcaption>Figure 17: Distribution of temporal length of wet spells for the observations (blue line) and in 500 simulations (red and grey envelopes) for 6 locations.</figcaption></figure>

## E.6 Relaxed ROR definition

Using an indicator that treats very small rainfall amounts as dry days is more relevant for drought risk assessment than defining a dry day strictly as one with zero precipitation. Therefore, a threshold of $R > 1$ mm can be chosen when defining a wet day, and the ROR$_{1.0}$(n) defined as

$$ \text{ROR}_{1.0}(n) = \frac{1}{D} \sum_{s=1}^{D} 1_{\{R_s^{(n)} > 1.0 \text{ mm}\}} \quad (43) $$

This indicator of dry days can be correctly represented only if the occurrence process Y and the rain intensity process R are properly linked, as both low rainfall amounts and zero rainfall must be jointly captured. Figure 22 shows that the generator successfully reproduces this indicator. The 1.0 mm tolerance increases 3 times the proportion of "dry" days.
<figure><img src="image_19.png" /><figcaption>Figure 18: 5-day cumulative rain for selected locations. 0.1 quantile in orange, 0.5 quantile in green, 0.9 quantile in blue. The lowest quantile is mostly 0, which makes it not visible. Each simulated envelope includes a central ribbon corresponding to the interquartile range (0.25–0.75 quantiles) and an outer ribbon showing the full range of the simulations (minimum–maximum).</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Figure 19: Auto-correlation function of the precipitation intensity for selected locations, for time lags between 0 and 10 days. The values from the observations are in blue, and the simulated envelope is in red.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>Figure 20: Joint rain occurrence probabilities on the same day (columns 1 and 2) and on consecutive days (columns 3 and 4), for the different seasons DJF: Winter, MAM: Spring, JJA: Summer, SON: Autumn, and all pairs of sites.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Figure 21: Distribution of length in days of dry spatial episodes, defined as ROR(n) < 0.1 for the observations and 4 generators. The values in the observations are the blue curves, and each coloured ribbon corresponds to the min-max envelope from simulation of the generators HMM0_Spa_K4, HMM1_Ind_K4, HMM1_Spa_K1, HMM1_Spa_K4.</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>Figure 22: Distribution (first row) and auto-correlation function (second row) of ROR$_{1.0}$ defined as the rain occurrence ratio with a tolerance of 1mm. The simulations from the generator are represented by its min-max envelope in orange, and the observed values are the blue dots.</figcaption></figure>
References

(SDES), Statistical Data and Studies Department (2024). *Observing climate change and its effects*. URL: https://www.statistiques.developpement-durable.gouv.fr/%20edition-numerique/chiffres-cles-du-climat/en/3-observing-%20climate-change-and-its.

Ailliot, Pierre, Denis Allard, et al. (2015). "Stochastic weather generators: an overview of weather type models". en. In: *Journal de la Societe Francaise de Statistique* 156.1, pp. 101-113. ISSN: 2102-6238.

Ailliot, Pierre, Craig Thompson, and Peter Thomson (2009). "Space-time modelling of precipitation by using a hidden Markov model and censored Gaussian distributions". en. In: *Journal of the Royal Statistical Society: Series C (Applied Statistics)* 58.3, pp. 405-426. ISSN: 1467-9876. DOI: 10.1111/j.1467-9876.2008.00654.x. URL: https://onlinelibrary.wiley.com/doi/abs/10.1111/j.1467-%209876.2008.00654.x.

Baum, Leonard E. et al. (1970). "A Maximization Technique Occurring in the Statistical Analysis of Probabilistic Functions of Markov Chains". In: *The Annals of Mathematical Statistics* 41.1, pp. 164-171. ISSN: 0003-4851. URL: https://www.jstor.org/stable/2239727.

Baxevani, Anastassia and Jan Lennartsson (2015). "A spatiotemporal precipitation generator based on a censored latent Gaussian field". en. In: *Water Resources Research* 51.6, pp. 4338-4358. ISSN: 1944-7973. DOI: 10.1002/2014WR016455. URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/%202014WR016455.

Benoit, Lionel, Denis Allard, and Gregoire Mariethoz (2018). "Stochastic Rainfall Modeling at Sub-kilometer Scale". en. In: *Water Resources Research* 54.6, pp. 4108-4130. ISSN: 1944-7973. DOI: 10.1029/2018WR022817. URL: https://onlinelibrary.wiley.com/doi/abs/10.1029/2018WR022817.

Bezanson, Jeff et al. (Jan. 2017). "Julia: A Fresh Approach to Numerical Computing". In: *SIAM Review* 59.1, pp. 65-98. ISSN: 0036-1445. DOI: 10.1137/141000671. URL: https://epubs.siam.org/doi/10.1137/141000671.

Bhat, Chandra R. (Mar. 2018). "New matrix-based methods for the analytic evaluation of the multivariate cumulative normal distribution function". In: *Transportation Research Part B: Methodological* 109, pp. 238-256. ISSN: 0191-2615. DOI: 10.1016/j.trb.2018.01.011. URL: https://www.sciencedirect.com/science/article/pii/%20S0191261517306847.

Botev, Z. I. (Jan. 2017). "The Normal Law Under Linear Restrictions: Simulation and Estimation via Minimax Tilting". en. In: *Journal of the Royal Statistical Society Series B: Statistical Methodology* 79.1, pp. 125-148. ISSN: 1369-7412. DOI: 10.1111/rssb.12162. URL: https://dx.doi.org/10.1111/rssb.12162.

Botev, Zdravko and Leo Belzile (Nov. 2015). *TruncatedNormal: Truncated Multivariate Normal and Student Distributions*. en. Institution: Comprehensive R Archive Network Pages: 2.3. DOI: 10.32614/CRAN/package.TruncatedNormal. URL: https://CRAN.R-project.org/package=TruncatedNormal.

Bourotte, Marc, Denis Allard, and Emilio Porcu (Nov. 2016). "A flexible class of non-separable cross-covariance functions for multivariate space-time data". en. In: *Spatial Statistics*. Spatial Statistics Avignon: Emerging Patterns 18, pp. 125-146. ISSN: 2211-6753. DOI: 10.1016/j.spasta.2016.02.004. URL: https://www.sciencedirect.com/science/article/pii/%20S2211675316000191.

Brisson, N et al. (Jan. 2003). "An overview of the crop model stics". In: *European Journal of Agronomy*. Modelling Cropping Systems: Science, Software and Applications 18.3, pp. 309-332. ISSN: 1161-0301. DOI: 10.1016/S1161-0301(02)00110-7. URL: https://www.sciencedirect.com/science/article/pii/%20S1161030102001107.

Carreau, J., P. Naveau, and E. Sauquet (2009). "A statistical rainfall-runoff mixture model with heavy-tailed components". en. In: *Water Resources Research* 45.10. ISSN: 1944-7973. DOI: 10.1029/2009WR007880. URL: https://onlinelibrary.wiley.com/doi/abs/10.1029/2009WR007880.

Celeux and Durand (2007). *Selecting hidden Markov model state number with cross-validated likelihood — Computational Statistics*. URL: https://link.springer.com/article/10.1007/s00180-007-0097-%201.

Cessie, S. Le and J. C. Van Houwelingen (1994). "Logistic Regression for Correlated Binary Data". In: *Applied Statistics* 43.1, p. 95. ISSN: 00359254. DOI: 10.2307/2986114. URL: https://www.jstor.org/stable/2986114?origin=crossref.

Chandler, Richard E. (Dec. 2020). "Multisite, multivariate weather generation based on generalised linear models". In: *Environmental Modelling & Software* 134, p. 104867. ISSN: 1364-8152. DOI:
10.1016/j.envsoft.2020.104867. URL: https://www.sciencedirect.com/science/article/pii/%20S1364815220309245.

Coles, Stuart (2001). *An Introduction to Statistical Modeling of Extreme Values*. Springer Series in Statistics. London: Springer London. ISBN: 978-1-84996-874-4 978-1-4471-3675-0. DOI: 10.1007/978-1-4471-3675-0. URL: http://link.springer.com/10.1007/978-1-4471-3675-0.

Cox, D. R. and N. Reid (Sept. 2004). "A note on pseudolikelihood constructed from marginal densities". en. In: *Biometrika* 91.3, pp. 729-737. ISSN: 0006-3444, 1464-3510. DOI: 10.1093/biomet/91.3.729. URL: https://academic.oup.com/biomet/article-lookup/doi/10.1093/%20biomet/91.3.729.

Dai, Bin, Shilin Ding, and Grace Wahba (Sept. 2013). "Multivariate Bernoulli distribution". en. In: *Bernoulli* 19.4. ISSN: 1350-7265. DOI: 10.3150/12-BEJSP10. URL: https://projecteuclid.org/journals/bernoulli/volume-19/%20issue-4/Multivariate-Bernoulli-distribution/10.3150/12-%20BEJSP10.full.

Dawkins, Laura C. et al. (Apr. 2022). "The Advanced Meteorology Explorer: a novel stochastic, gridded daily rainfall generator". In: *Journal of Hydrology* 607, p. 127478. ISSN: 0022-1694. DOI: 10.1016/j.jhydrol.2022.127478. URL: https://www.sciencedirect.com/science/article/pii/%20S0022169422000531.

Dempster, A. P., N. M. Laird, and D. B. Rubin (Sept. 1977). "Maximum Likelihood from Incomplete Data Via the EM Algorithm". en. In: *Journal of the Royal Statistical Society Series B: Statistical Methodology* 39.1, pp. 1-22. ISSN: 1369-7412, 1467-9868. DOI: 10.1111/j.2517-6161.1977.tb01600.x. URL: https://academic.oup.com/jrsssb/article/39/1/1/7027539.

Dixit, Vaibhav Kumar and Christopher Rackauckas (Mar. 2023). *Optimization.jl: A Unified Optimization Package*. DOI: 10.5281/ZENODO.7738525. URL: https://zenodo.org/record/7738525.

Durban, M. and C. A. Glasbey (Sept. 2001). "Weather modelling using a multivariate latent Gaussian model". In: *Agricultural and Forest Meteorology* 109.3, pp. 187-201. ISSN: 0168-1923. DOI: 10.1016/S0168-1923(01)00268-4. URL: https://www.sciencedirect.com/science/article/pii/%20S0168192301002684.

Evin, Guillaume, Anne-Catherine Favre, and Benoît Hingray (Jan. 2018). "Stochastic generation of multi-site daily precipitation focusing on extreme events". en. In: *Hydrology and Earth System Sciences* 22.1, p. 655. DOI: 10.5194/hess-22-655-2018. URL: https://hal.inrae.fr/hal-02912352.

Fan, Jianqing et al. (Mar. 2017). "High Dimensional Semiparametric Latent Graphical Model for Mixed Data". en. In: *Journal of the Royal Statistical Society Series B: Statistical Methodology* 79.2, pp. 405-421. ISSN: 1369-7412, 1467-9868. DOI: 10.1111/rssb.12168. URL: https://academic.oup.com/jrsssb/article/79/2/405/7041934.

Furrer, Em and Rw Katz (July 2007). "Generalized linear modeling approach to stochastic weather generators". en. In: *Climate Research* 34, pp. 129-144. ISSN: 0936-577X, 1616-1572. DOI: 10.3354/cr034129. URL: http://www.int-res.com/abstracts/cr/v34/n2/p129-144/.

Gamet, Philémon and Jonathan Jalbert (2022). "A flexible extended generalized Pareto distribution for tail estimation". en. In: *Environmetrics* 33.6, e2744. ISSN: 1099-095X. DOI: 10.1002/env.2744. URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/env.2744.

(2024). *ExtendedExtremes.jl*. URL: https://github.com/JuliaExtremes/ExtendedExtremes.jl.

Gao, Xin and Peter X.-K. Song (2011). "Composite Likelihood Em Algorithm with Applications to Multivariate Hidden Markov Model". In: *Statistica Sinica* 21.1, pp. 165-185. ISSN: 1017-0405. URL: https://www.jstor.org/stable/24309267.

Garçon, R. (Aug. 1996). "Prévision opérationnelle des apports de la Durance à Serre-Ponçon à l'aide du modèle MORDOR. Bilan de l'année 1994-1995". In: *La Houille Blanche* 82.5, pp. 71-76. ISSN: 0018-6368. DOI: 10.1051/lhb/1996056. URL: https://doi.org/10.1051/lhb/1996056.

Genton, Marc G., David E. Keyes, and George Turkiyyah (Apr. 2018). "Hierarchical Decompositions for the Computation of High-Dimensional Multivariate Normal Probabilities". In: *Journal of Computational and Graphical Statistics* 27.2. _eprint: https://doi.org/10.1080/10618600.2017.1375936, pp. 268-277. ISSN: 1061-8600. DOI: 10.1080/10618600.2017.1375936. URL: https://doi.org/10.1080/10618600.2017.1375936.

Genz, Alan (1992). "Numerical Computation of Multivariate Normal Probabilities". In: *Journal of Computational and Graphical Statistics* 1.2, pp. 141-149. ISSN: 1061-8600. DOI: 10.2307/1390838. URL: https://www.jstor.org/stable/1390838.
Genz, Alan (Aug. 2004). "Numerical computation of rectangular bivariate and trivariate normal and t probabilities". en. In: *Statistics and Computing* 14.3, pp. 251-260. ISSN: 0960-3174. DOI: 10.1023/B:STCO.0000035304.20635.31. URL: http://link.springer.com/10.1023/B:%20STCO.0000035304.20635.31.

Geoga, Christopher J. et al. (Feb. 2023). "Fitting Matérn smoothness parameters using automatic differentiation". en. In: *Statistics and Computing* 33.2, p. 48. ISSN: 1573-1375. DOI: 10.1007/s11222-022-10127-w. URL: https://doi.org/10.1007/s11222-022-10127-w.

Gneiting, Tilmann (June 2002). "Nonseparable, Stationary Covariance Functions for Space-Time Data". en. In: *Journal of the American Statistical Association* 97.458, pp. 590-600. ISSN: 0162-1459, 1537-274X. DOI: 10.1198/016214502760047113. URL: http://www.tandfonline.com/doi/abs/10.1198/%20016214502760047113.

Gobet, Emmanuel, David Métivier, and Sylvie Parey (Sept. 2025). "Interpretable seasonal multisite hidden Markov model for stochastic rain generation in France". English. In: *Advances in Statistical Climatology, Meteorology and Oceanography* 11.2, pp. 159-201. ISSN: 2364-3579. DOI: 10.5194/ascmo-11-159-2025. URL: https://ascmo.copernicus.org/articles/11/159/2025/.

Heagerty, Patrick J. and Subhash R. Lele (Sept. 1998). "A Composite Likelihood Approach to Binary Spatial Data". en. In: *Journal of the American Statistical Association* 93.443, pp. 1099-1111. ISSN: 0162-1459, 1537-274X. DOI: 10.1080/01621459.1998.10473771. URL: http://www.tandfonline.com/doi/full/10.1080/%2001621459.1998.10473771.

Hughes, James P., Peter Guttorp, and Stephen P. Charles (Mar. 1999). "A Non-Homogeneous Hidden Markov Model for Precipitation Occurrence". In: *Journal of the Royal Statistical Society Series C: Applied Statistics* 48.1, pp. 15-30. ISSN: 0035-9254. DOI: 10.1111/1467-9876.00136. URL: https://doi.org/10.1111/1467-9876.00136.

Hundecha, Yeshewatesfa, Markus Pahlow, and Andreas Schumann (2009). "Modeling of daily precipitation at multiple locations using a mixture of distributions to characterize the extremes". en. In: *Water Resources Research* 45.12. ISSN: 1944-7973. DOI: 10.1029/2008WR007453. URL: https://onlinelibrary.wiley.com/doi/abs/10.1029/%202008WR007453.

Intergovernmental Panel On Climate Change (IPCC) (July 2023). *Climate Change 2021 – The Physical Science Basis: Working Group I Contribution to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change*. en. 1st ed. Cambridge University Press. ISBN: 978-1-009-15789-6. DOI: 10.1017/9781009157896. URL: https://www.cambridge.org/core/product/identifier/%209781009157896/type/book.

Kirshner, Sergey (2005). "Modeling of multivariate time series using hidden Markov models". PhD thesis. University of California, Irvine.

Kleiber, William, Richard W. Katz, and Balaji Rajagopalan (Jan. 2012). "Daily spatiotemporal precipitation simulation using latent and transformed Gaussian processes: PRECIPITATION SIMULATION USING GAUSSIAN PROCESSES". en. In: *Water Resources Research* 48.1. ISSN: 00431397. DOI: 10.1029/2011WR011105. URL: http://doi.wiley.com/10.1029/2011WR011105.

Klein Tank, A. M. G. et al. (2002). "Daily dataset of 20th-century surface air temperature and precipitation series for the European Climate Assessment". en. In: *International Journal of Climatology* 22.12, pp. 1441-1453. ISSN: 1097-0088. DOI: 10.1002/joc.773. URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/joc.773.

Li, Xin and Vladan Babovic (May 2019). "Multi-site multivariate downscaling of global climate model outputs: an integrated framework combining quantile mapping, stochastic weather generator and Empirical Copula approaches". en. In: *Climate Dynamics* 52.9, pp. 5775-5799. ISSN: 1432-0894. DOI: 10.1007/s00382-018-4480-0. URL: https://doi.org/10.1007/s00382-018-4480-0.

Li, Zhi (Aug. 2014). "A new framework for multi-site weather generator: a two-stage model combining a parametric method with a distribution-free shuffle procedure". en. In: *Climate Dynamics* 43.3, pp. 657-669. ISSN: 1432-0894. DOI: 10.1007/s00382-013-1979-2. URL: https://doi.org/10.1007/s00382-013-1979-2.

Lin, Jinn-Tyan (Apr. 1995). "A Simple Approximation for the Bivariate Normal Integral". en. In: *Probability in the Engineering and Informational Sciences* 9.2, pp. 317-321. ISSN: 1469-8951, 0269-9648. DOI: 10.1017/S0269964800003880. URL: https://www.cambridge.org/core/journals/probability-in-the-%20engineering-and-informational-sciences/article/simple-%20approximation-for-the-bivariate-normal-integral/%209455086EE8DBAB83643B34ADCB9A410D.
Meng, Xiao-Li and Donald B. Rubin (1993). "Maximum Likelihood Estimation via the ECM Algorithm: A General Framework". In: *Biometrika* 80.2, pp. 267-278. ISSN: 0006-3444. DOI: 10.2307/2337198. URL: https://www.jstor.org/stable/2337198.

Météo-France (2024). *Climate report for the year 2023*. URL: https://www.meteofrance.fr/sites/meteofrance.fr/files/files/%20editorial/bilan_2023_web.pdf.

Miloshevich, George et al. (Apr. 2023). "Probabilistic forecasts of extreme heatwaves using convolutional neural networks in a regime of lack of data". In: *Physical Review Fluids* 8.4, p. 040501. DOI: 10.1103/PhysRevFluids.8.040501. URL: https://hal.science/hal-04124796.

Naveau, Philippe et al. (2016). "Modeling jointly low, moderate, and heavy rainfall intensities without a threshold selection". en. In: *Water Resources Research* 52.4, pp. 2753-2769. ISSN: 1944-7973. DOI: 10.1002/2015WR018552. URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/202015WR018552.

Obakrim, Said, Lionel Benoit, and Denis Allard (Sept. 2025). "A multivariate and space-time stochastic weather generator using a latent Gaussian framework". en. In: *Stochastic Environmental Research and Risk Assessment* 39.9, pp. 3677-3701. ISSN: 1436-3259. DOI: 10.1007/s00477-024-02897-8. URL: https://doi.org/10.1007/s00477-024-02897-8.

Pickands, James (Jan. 1975). "Statistical Inference Using Extreme Order Statistics". In: *The Annals of Statistics* 3.1, pp. 119-131. ISSN: 0090-5364, 2168-8966. DOI: 10.1214/aos/1176343003. URL: https://projecteuclid.org/journals/annals-of-statistics/%20volume-3/issue-1/Statistical-Inference-Using-Extreme-Order-%20Statistics/10.1214/aos/1176343003/full.

Ragulina, Galina and Trond Reitan (Apr. 2017). "Generalized extreme value shape parameter and its nature for extreme precipitation using long time series and the Bayesian approach". In: *Hydrological Sciences Journal* 62.6, pp. 863-879. ISSN: 0262-6667. DOI: 10.1080/02626667.2016.1260134. URL: https://doi.org/10.1080/02626667.2016.1260134.

Rajagopalan, Balaji and Upmanu Lall (1999). "A k-nearest-neighbor simulator for daily precipitation and other weather variables". en. In: *Water Resources Research* 35.10, pp. 3089-3101. ISSN: 1944-7973. DOI: 10.1029/1999WR900028. URL: https://onlinelibrary.wiley.com/doi/abs/10.1029/201999WR900028.

Revels, Jarrett, Miles Lubin, and Theodore Papamarkou (July 2016). *Forward-Mode Automatic Differentiation in Julia*. arXiv:1607.07892 [cs]. DOI: 10.48550/arXiv:1607.07892. URL: http://arxiv.org/abs/1607.07892.

Richardson, C. W. (1981). "Stochastic simulation of daily precipitation, temperature, and solar radiation". en. In: *Water Resources Research* 17.1, pp. 182-190. ISSN: 1944-7973. DOI: 10.1029/WR017i001p00182. URL: https://onlinelibrary.wiley.com/doi/abs/10.1029/%20WR017i001p00182.

Rodriguez-Iturbe, Ignacio and Peter S. Eagleson (1987). "Mathematical models of rainstorm events in space and time". en. In: *Water Resources Research* 23.1, pp. 181-190. ISSN: 1944-7973. DOI: 10.1029/WR023i001p00181. URL: https://onlinelibrary.wiley.com/doi/abs/10.1029/WR023i001p00181.

Soubeyroux, Jean-Michel et al. (Nov. 2024). *A quel climat s'adapter en France selon la TRACC ?* Tech. rep. Meteo-France. URL: https://hal.science/hal-04797481.

Sparks, Nathan J. et al. (Mar. 2018). "IMAGE: a multivariate multi-site stochastic weather generator for European weather and climate". en. In: *Stochastic Environmental Research and Risk Assessment* 32.3, pp. 771-784. ISSN: 1436-3259. DOI: 10.1007/s00477-017-1433-9. URL: https://doi.org/10.1007/s00477-017-1433-9.

Stoehr, Julien and Stephane S. Robin (Apr. 2024). *Composite likelihood inference for the Poisson log-normal model*. arXiv:2402.14390 [stat]. DOI: 10.48550/arXiv.2402.14390. URL: http://arxiv.org/abs/2402.14390.

Touron, Augustin (Sept. 2019). "Chapter 5 - Bivariate modeling of temperature and precipitations". English. en. In: *Modélisation multivariée de variables météorologiques*. Université Paris-Saclay. URL: https://theses.hal.science/tel-02319170.

Tsay, Wen-Jen and Peng-Hsuan Ke (Apr. 2023). "A simple approximation for the bivariate normal integral". In: *Communications in Statistics - Simulation and Computation* 52.4, pp. 1462-1475. ISSN: 0361-0918. DOI: 10.1080/03610918.2021.1884718. URL: https://doi.org/10.1080/03610918.2021.1884718.
Vaittinada Ayar, Pradeebane et al. (Feb. 2020). "Space-time simulation of precipitation based on weather pattern sub-sampling and meta-Gaussian model". en. In: *Journal of Hydrology* 581, p. 124451.
ISSN: 0022-1694. DOI: 10.1016/j.jhydrol.2019.124451. URL: https://www.sciencedirect.com/science/article/pii/%20S0022169419311862.

Varin, Cristiano, Gudmund Høst, and Oivind Skare (June 2005). "Pairwise likelihood inference in spatial generalized linear mixed models". In: *Computational Statistics & Data Analysis* 49.4, pp. 1173–1191.
ISSN: 0167-9473. DOI: 10.1016/j.csda.2004.07.021. URL: https://www.sciencedirect.com/science/article/pii/%20S0167947304002403.

Varin, Cristiano, Nancy Reid, and David Firth (2011). "An Overview of Composite Likelihood Methods". In: *Statistica Sinica* 11.1, pp. 5–42.
ISSN: 1017-0405. URL: https://www.jstor.org/stable/24309261.

Verdin, Andrew et al. (2019). "BayGEN: A Bayesian Space-Time Stochastic Weather Generator". en. In: *Water Resources Research* 55.4, pp. 2900–2915.
ISSN: 1944-7973. DOI: 10.1029/2017WR022473. URL: https://onlinelibrary.wiley.com/doi/abs/10.1029/%202017WR022473.

Wilhelm, Stefan and Manjunath B. G (Sept. 2025). *tmvtnorm: Truncated Multivariate Normal and Student t Distribution*. URL: https://cran.r-project.org/web/packages/tmvtnorm/%20index.html.

Wilks, D. S. (Sept. 1998). "Multisite generalization of a daily stochastic precipitation generation model". In: *Journal of Hydrology* 210.1, pp. 178–191.
ISSN: 0022-1694. DOI: 10.1016/S0022-1694(98)00186-3. URL: https://www.sciencedirect.com/science/article/pii/S0022169498001863.

Wilks, Daniel S. (2009). "A gridded multisite weather generator and synchronization to observed weather data". en. In: *Water Resources Research* 45.10. ISSN: 1944-7973. DOI: 10.1029/2009WR007902.
URL: https://onlinelibrary.wiley.com/doi/abs/10.1029/%202009WR007902.

Yiou, P. (Apr. 2014). "AnaWEGE: a weather generator based on analogues of atmospheric circulation". English. In: *Geoscientific Model Development* 7.2, pp. 531–543.
ISSN: 1991-959X. DOI: 10.5194/gmd-7-531-2014. URL: https://gmd.copernicus.org/articles/7/531/2014/.

Yiou, P. et al. (May 2008). "Weather regime dependence of extreme value statistics for summer temperature and precipitation". English. In: *Nonlinear Processes in Geophysics* 15.3, pp. 365–378.
ISSN: 1023-5809. DOI: 10.5194/npg-15-365-2008. URL: https://npg.copernicus.org/articles/15/365/2008/.

Zhou, Lingfeng, Yaobin Meng, and Karim C. Abbaspour (Aug. 2019). "A new framework for multi-site stochastic rainfall generator based on empirical orthogonal function analysis and Hilbert-Huang transform". In: *Journal of Hydrology* 575, pp. 730–742.
ISSN: 0022-1694. DOI: 10.1016/j.jhydrol.2019.05.047. URL: https://www.sciencedirect.com/science/article/pii/S0022169419304925.

Zucchini, Walter and Peter Guttorp (1991). "A Hidden Markov Model for Space-Time Precipitation". en. In: *Water Resources Research* 27.8, pp. 1917–1923.
ISSN: 1944-7973. DOI: 10.1029/91WR01403. URL: https://onlinelibrary.wiley.com/doi/abs/10.1029/91WR01403.