## 2 Materials and Methods
### GAMLSS Approach  
The GAMLSS proposed by [PERSON] and [PERSON] [30] were chosen to describe the relationship between drought factors and explanatory covariates. These are widely used nonstationary models for hydro-meteorological elements [28; 31; 35; 39]. GAMLSS are univariate distributional regression models, wherein all the parameters of the assumed distribution for the response can be modelled as additive functions of the explanatory variables through linear or nonlinear functions.  
\begin{table}
\begin{tabular}{c c} \hline \hline
**Threshold Value** & **Description** \\ \hline SPEI \(>-0.50\) & Normal \\ \(-1.00<\) SPEI \(\leq-0.50\) & Mild drought \\ \(-1.50<\) SPEI \(\leq-1.00\) & Moderate drought \\ \(-2.00<\) SPEI \(\leq-1.50\) & Severe drought \\ SPEI \(\leq-2.00\) & Extreme drought \\ \hline \hline \end{tabular}
\end{table}
Table 1: Threshold values of the SPEI index classification.  
The principle of the GAMLSS model is as follows.  
In the GAMLSS, the times series \(y_{i}\), where \(i=1,2,\ldots,n\) follows the probability density function \(\mathrm{f}\big{(}y_{i}|\theta^{i}\big{)}\) conditional on \(\theta^{i}=(\theta_{1i},\theta_{2i},\ldots,\theta_{ai})\), where \(\theta^{i}\) is the vector of the o para meter of the distribution. Let \(y^{T}=(y_{1},y_{1},\ldots,y_{n})\) be the \(n\)-length vector of the response variables in the observations. For \(\mathrm{k}\ =1,2,\ldots,p,\ J_{k}=0\), a fully parametric model for \(\mathrm{g}_{k}(\cdot)\) is a known monotonic link function relating \(\theta_{k}\) to the explanatory variable and random effects through an additive model given by:  
\[\mathrm{g}_{k}(\theta_{k})=\eta_{k}=X_{k}\beta_{k}+\sum_{j=1}^{k}h_{jk}\Big{(} x_{jk}\Big{)} \tag{14}\]  
where \(\theta_{k}\) and \(\eta_{k}\) are vectors of length \(n\), \(X_{k}\) is a matrix (\(n\times m\)) of the explanatory variables, and \(\beta_{k}\) is a parameter vector of length \(m\). The function \(h_{jk}(\cdot)\) is an unknown function of the explanatory variable \(x_{jk}\) for \(j=1,2,\ldots,k\) and \(k=1,2,\ldots,p\). Equation (14) represents a semiparametric model, wherein the first two population parameters \(\theta_{1}\) and \(\theta_{2}\) are usually characterized as location and scale parameters, denoted by \(\upmu\) and \(\upsigma\), respectively, whereas the remaining parameters, if any, are characterized as shape parameters.  
Two basic algorithms were used for maximizing the penalized likelihood: the Conjugate Gradient (CG) algorithm [51] and the Rigby and Stasinopoulos (RS) algorithm [52; 53]. The RS algorithm was easy to implement and has been successfully used for fitting all the distributions, including the logistic distribution. The Akaike information criterion (AIC) [54] was used to evaluate the goodness of fit and avoid model overfitting.  
\[\mathrm{AIC}\ =2\mathrm{l}\ -2\mathrm{ln}(\hat{L}) \tag{15}\]  
where \(\mathrm{l}\) is the number of estimated parameters in the model, and \(\hat{L}\) is the maximum value of the likelihood function of the model. All computations in this study were performed with the GAMLSSS R package available from the Comprehensive R Archive Network (CRAN) of the R library.