A Continuous Multisite Multivariate Generator for Daily Temperature Conditioned by Precipitation Occurrence

Abstract

Temperature is one of the most influential weather variables necessary for numerous studies, such as climate change, integrated water resources management, and water scarcity, among others. The temperature and precipitation are relevant in river basins because they may be particularly affected by modifications in the variability, for example, due to climate change. We developed a stochastic model for daily precipitation occurrences and their influence on maximum and minimum temperatures with a straightforward approach. The Markov model has been used to determine everyday occurrences of rainfall. Moreover, we developed a multisite multivariate autoregressive model to represent the short-term memory of daily temperature, called MASCV. The reduction of parameters is an essential factor addressed in this approach. For this reason, the normalization of the temperatures was performed through different nonparametric transformations. The case study is the Jucar River Basin in Spain. The multisite multivariate stochastic model of two states and a lag-one accurately represents both occurrences as well as maximum and minimum temperature. The simulation and generation of occurrences and temperature is considered a continuous multivariate stochastic process. Additionally, time series of multiple correlated climate variables are completed. Therefore, we simplify the complexity and reduce the computational time for the simulation.

1. Introduction

The stochastic modeling approach has been widely used for hydrologic time series analysis, including modeling weather variables [1,2] and flood prediction [3]. Therefore, it is essential to build accurate forecast models in the hydrologic process [4]. The stochastic modeling of temperature conditioned by precipitation is proposed when a day is wet or dry. It is necessary to estimate the temperature for both dry and wet days.

The most common stochastic model is the first-order Markov model with two states. It was introduced by Gabriel and Newman [5] and has since been used and modified by many authors [3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19]. The rainfall occurrence process focuses on representing the dry and wet days. Critical probability ($p_c$) depends on transition probabilities of dry–wet ($p_{01}$) and wet–wet days ($p_{11}$) [8]. Thus, the occurrence is a bivariate function (0,1) that relies on the uniform random process ($u$) and critical probabilities [12].

The temperature displays a near-normal distribution. It is common to consider the normal distribution to minimize the skewness coefficient of observed data. In other cases, root transformation is used [20].

The statistical characteristics of the series change throughout the whole year. The mean, variance, and standard deviation are periodic—these statistical changes are yearly, monthly, daily, and even hourly. This periodicity is commonly modeled by the finite Fourier series [9,16,21,22]. The periodic component is necessary for standardization, which is a process that converts into a series with a mean of zero and a standard deviation of one. Standardization is usually applied to series with no follow-up normal distributions [7,23,24,25]. Moreover, standardization helps remove the seasonal effects [26]. Standardized series are used to calculate the autoregressive parameters of the stochastic model [27].

It is desirable to generate synthetic temperature based on weather data [9]. The process should be capable of representing different statistics from historical data [18]. The basic statistics of the series are the sample mean, standard deviation, skewness coefficient, cross-correlation coefficients, and probability distribution. We must fit a probability distribution function (two or more parameters) or use a nonparametric transformation. Many authors have proposed stochastic weather generators (WG) for different purposes, a few of which are downscaling, prediction, or simulation [28]. Some developed models are Weather GENerator (WGEN) [29], CLImate GENerator (CLIGEN) [30], Agriculture and Agri-Food Canada Weather Generator (AAA-FC) [27], CLIMate GENerator (CLIMGEN) [31], Long Ashton Research Station-Weather Generator (LARS-WG) [21], the weather generator CLIMA [32], the closed skew-normal weather generator (WACS-Gen) [26], École de technologie supérieure Weather Generator (WeaGETS) [22], and the Multi-site Rainfall Simulator (MRS) [18]. We propose the multivariate autoregressive model of climate variables (MASCV).

For precipitation, weather generators analyze the occurrence similarly. For example, CLIMA [32], CLIGEN [30], and WGEN [9] model the occurrence process based on transition probabilities of lag-one and two states for each month (two parameters). LARS-WG [21] is based on an empirical function of wet and dry days for each month (21 parameters). CLIMGEN [31] models the occurrence process through transition probabilities of a second-order Markov model at a monthly scale (four parameters). WeaGETS [22] models the transition probabilities of a third-order Markov model with two states on a biweekly scale (eight parameters). WASC-Gen uses a lag-one Markov model with different states and parameters derived from the Bayesian information criterion [26]. AA-FC-WG models the occurrence process with a second-order Markov model at a monthly scale using two parameters [27]. Other stochastic models have proposed modifications or new methods, such as the Markov renewal model [33], Copula-based models [34], the Hidden Markov model [35,36], the Semi-Markov model [37], generalized chain-dependent process models [38], the autoregressive model [39], and artificial intelligence [40].

Artificial intelligence (AI) has been used to predict hydrological data, including rainfall, rainfall runoff, and temperature [4]. Various machine learning methods have been applied for rainfall and temperature prediction, such as the support vector machine (SVM) classifier [41], ANN [42], long short-term memory (LSTM), statistical and multiple linear regression (MLR) [43,44], and classification and regression trees [45].

Stochastic weather generators are used to represent different variables such as precipitation and temperature. Long temperature series are commonly required, but historical data are sometimes short in length. Therefore, it is necessary to apply stochastic weather generators for hydrological design [12]. The integrated water resources management (IWRM) depends on water availability [46]. For the hydric balance and various decisions at a river basin scale, the rainfall is fundamental [10], primarily to generate long synthetic series.

The issue with synthetic series is preserving the low-frequency variability in temperature and rainfall [47]. Daily temperature series have a strong correlation [48]. To preserve low-frequency variability, different authors have proposed the fast Fourier transform (FFT) algorithm [49], using crop variables [50], annual and monthly autoregressive models [51], k-nearest neighbor bootstrap [52], and the copulas approach [53].

The objective of this paper was to develop an adequate stochastic model for daily maximum and minimum temperature with a straightforward approach. We propose a multisite multivariate stochastic model for precipitation occurrences and temperature capable of reproducing spatial and temporal dependence. We corrected the daily stochastic model using an annual multisite multivariate autoregressive model. For precipitation, we applied a multivariate multisite lag-one Markov model with two states (wet and dry) and few parameters. A normal distribution was used to define the precipitation occurrence process. For the maximum and minimum temperature, we implemented the first-order multisite multivariate autoregressive model (MAR(p)) for daily and annual temperature. Our approach simplifies the temperature simulation (continuous and nonparametric), which implies a considerable advantage and versatility concerning other stochastic generators (parametric distribution function for each month or biweekly).

The normalized mean and variance periodicity was modeled as a continuous daily scale through the Fourier series. The model validation was performed by synthetic series of precipitation, which accurately represent the main statistics of observed data for different climates. MASCV was programmed in MATLAB and is capable of simulating with various parameters for both occurrences and temperature. Moreover, we proposed a multisite multivariate stochastic model in which different wet thresholds and nonparametric transformations can be applied. A relevant advantage is that we modeled continuously for the whole series. The method developed in this paper generates constant day-to-day rainfall occurrences and temperature.

2. Materials and Methods

The purpose of generating long series of temperatures is to evaluate the effects of hydrologic changes [4,9] and analyze different scenarios, such as environmental, agricultural [28], or climate change. We focused on developing long multivariate synthetic series of maximum and minimum temperature. The lag-one Markov model has been applied for plenty of weather modelers due to its accurate representation of dependence [12]. Said model is based on terms of occurrence given the previous day. We modeled the rainfall occurrence and the temperature separately (Figure 1), so we defined the model according to Equation (1):

$$Y_t=T_t{X}_t$$

(1)

where $T_t$ is the temperature model, $X_t$ is the precipitation occurrence model, and $Y_t$ is the whole stochastic process.

2.1. Multivariate Precipitation Occurrence (Dry–Wet)

First, we must determine whether a day is wet or dry according to the day prior. We used a bivariate function ($X_t$), and if the precipitation is greater than a given threshold, the t day is wet ($X_t=1$); otherwise, the t day is dry $\left(X_t=0\right)$. The first-order Markovian approach only depends on the previous day, whether it was wet or dry. The high-order Markov model has been notably studied [54,55] and is recommended for vast persistence [10]. The two- and three-order Markov models significantly improve the fit [22,56,57]. In other cases, the results are nearly equivalent to the first-order Markov model [54]. One disadvantage of the high-order Markov model is the increase in the number of parameters [28]. The threshold for precipitation occurrence could change according to different data. This depends on the minimum value of precipitation amount that the station in the study can report. Common threshold values are 0 [9,21,29,30], 0. [22,26,58], 0.2 [25], 0.254, and 0.3 mm [18]. We used the wet thresholds (0.001, 0.01, 0.1, and 0.25) and identified the ones that provided the best results. A disadvantage in the case of stochastic modeling is the limited data from the historical period.

For the occurrence process, we used the Wilks approach [12]. The conditional probabilities for a one-order Markov model are a dry day followed by a wet day, $p_{01}$, a wet day followed by a wet day, $p_{11}$, a wet day followed by a dry day, $p_{10}$, and a dry day followed by a dry day, $p_{00}$. The conditional probabilities have a complementary relationship between them.

The normal critical probability $(p_c)$ depends on the previous day $\left(X_{t-1}\right)$; if it is dry, $p_c=p_{01}$, otherwise $p_c=p_{11}$. The boundary for the precipitation occurrence process is determined for a random spatially correlated normal matrix $n_t={\varnothing }^{-1}{u}_t \left[0,1\right]$. Only one day is computed if the random variable is equal to or less than a critical probability.

$$if{n}_t\le p_c\therefore X_t=1;n_t>p_c\therefore X_t=0$$

(2)

The transition probability vectors vary on a monthly and daily scale. Therefore, the occurrence process was simulated through different parameters. The lag-one Markov model with two states enables the interaction between wet and dry days. Moreover, it allows the residence time and periodicity in both the first (dry) and second (wet) states.

For the daily occurrence process, one of the disadvantages is the limited number of data, primarily for dry seasons and arid river basins. These data limit the quality of results, and it is common to model the transition probabilities monthly or biweekly [9,22,26,27,30,31].

Our approach was to simulate the transition probabilities on a daily scale, similar to Woolhiser and Pegram [59]. It is imperative to preserve most of the characteristics for days, months, and years. The number of parameters increases from a monthly to a daily scale; therefore, it is essential to determine the most efficient model. In this study, we focused on the relevance of evaluating different analyses, which depend on the number of parameters, ${\mathrm{p}}_{01\mathsf{\tau }}={\mathrm{P}}_{01}{|}_{\mathsf{\tau }}$ and ${\mathrm{p}}_{01\mathsf{\tau }}={\mathrm{P}}_{11}{|}_{\mathsf{\tau }}$.

Due to the distinct temporal scales of analysis, the number of parameters will increase according to the variability of the occurrence process. The transition probabilities follow a uniform distribution applicable in the case of univariate modeling. However, in the case of multivariate modeling, a transformation to the normal distribution is more efficient [15]. A transformation was performed from uniform to normal, $p_{n11\mathsf{\tau }}={\varnothing }^{-1}\left(p_{11\mathsf{\tau }}\right) \mathrm{and} p_{\mathrm{n}01\mathsf{\tau }}={\varnothing }^{-1} \left(p_{01\mathsf{\tau }}\right)$. The spatial correlation between probability vectors with lag-k is shown in Equation (3):

$$r_k{}^{\left(i,j\right)}\left(p\right)={(p_{\mathrm{n}11}{}^{\left(i,j\right)}-p_{\mathrm{n}01}{}^{\left(i,j\right)})}^{k}i\ne j$$

(3)

where $r_k{}^{\left(i,j\right)}$ is the spatial correlation with lag-k, $p_{\mathrm{n}11}{}^{\left(i,j\right)}$ is the normalized probability vector of two consecutive days of rainfall occurrence, and $p_{\mathrm{n}01}{}^{\left(i,j\right)}$ is the normalized probability vector of a dry day followed by a wet day. Based on these vectors, a spatially correlated matrix was proposed in Equation (4):

$$M_k=\left[\begin{array}{ccc}{r}_k{}^{\left(1,1\right)}\left(p\right)& \cdots & r_k{}^{\left(1,n\right)}\left(p\right)\\ ⋮& \ddots & ⋮\\ r_k{}^{\left(n,1\right)}\left(p\right)& \cdots & r_k{}^{\left(n,n\right)}\left(p\right)\end{array}\right]$$

(4)

where $M_k$ is the cross-correlation matrix. The Cholesky factorization was used to determine the spatial dependence of the series [15,19,60]. For a positive definite matrix, the Cholesky factorization is a lower triangular matrix, $D=\left[M\right]{\left[M\right]}^{\prime }$. Multiplying this lower triangular matrix by a random normal matrix results in a random spatially correlated normal matrix, n, in Equation (5):

$$n={\left[M\right]}^{\prime }\left[N\right]=\left[\begin{array}{ccc}{M^{\prime }}^{\left(1,1\right)}& \cdots & 0\\ ⋮& \ddots & ⋮\\ {M^{\prime }}^{\left(n,1\right)}& \cdots & {M^{\prime }}^{\left(n,n\right)}\end{array}\right]x\left[\begin{array}{ccc}{N}^{\left(1,1\right)}& \cdots & N^{\left(1,m\right)}\\ ⋮& \ddots & ⋮\\ N^{\left(n,1\right)}& \cdots & N^{\left(n,m\right)}\end{array}\right]$$

(5)

where $n$ is the random spatially correlated normal matrix, $M^{{}^{\prime }\left(n,n\right)}$ is the lower triangular matrix for n series, and $N^{\left(n,m\right)}$ is the random normal matrix for n series and m days. We used this matrix to generate multivariate precipitation occurrences. In addition, they were used to determine synthetic wet and dry days.

2.2. Multivariate Maximum and Minimum Temperature

The maximum and minimum temperatures were modeled using the original maximum temperature data and the difference between the maximum and minimum temperature, Trange in Equation (6). The target of modeling the temperature range is to avoid negative values of the synthetic series [61].

$$Trange=Tmax-Tmin$$

(6)

The temperature process is commonly modeled considering normal distribution. However, according to our experience, in some cases, the maximum temperature and temperature range follow a normal distribution. We focused on the nonparametric transformations. These are practical solutions to model the temperature [25,56,62].

Seasonal variability is one of the most significant characteristics of the stochastic temperature process. The daily temperature shows recurrent changes within the year. Parameters such as the mean and standard deviation are periodic components [1]. The periodicity was analyzed through the Fourier series [59]. The number of parameters is reduced when these are periodic or seasonal. The Fourier series [23,59] is present in Equation (7). Several harmonics were used to represent 90% of the explicative variance of the observed data.

$$v=\overline{u}+{\displaystyle \sum }_{j=1}^h\left[A_j\cos \left(\frac{2\pi \tau }{w}\right)+B_j\sin \left(\frac{2\pi \tau }{w}\right)\right]$$

(7)

where $\overline{u}$ is the normalized temperature mean for wet (${\mu }_{1\tau }$) and dry days (${\mu }_{0\tau }$) and standard deviation for wet ($s_{1\tau }$) and dry days ($s_{0\tau }$). A and B are the Fourier coefficient vectors, $j$ is the harmonic, and $h$ is the total number of harmonics, equal to (w − 1)/2 for even numbers and equal to w/2 for uneven numbers. For example, in a daily simulation, we have 365 days, and the maximum number of harmonics is 182. Important harmonics were selected according to the accumulated period-gram, defined as the ratio of mean standard deviation (MSD) of the harmonics to the observed series. We accepted 90% of the original data representation to select the significant harmonics applied at the mean, standard deviation, and transition probabilities $({\mu }_{1\tau }, {\mu }_{0\tau }, s_{1\tau }, s_{0\tau }, p_{01\tau }, p_{11\tau }$).

Once the periodic component was modeled, we standardized both temperatures (maximum and range), allowing the analysis of temporal dependence. In addition, a standardized series served to generate residual series. Standardization removes the periodicity of the series based on the mean (${\mu }_{1\tau },{\mu }_{0\tau }$) and standard deviation (${\mathrm{s}}_{1\tau }, {\mathrm{s}}_{0\tau }$). We determined the standardized series ($z_{\tau }$) at a daily scale, canceling the mean and standard deviation of the normalized series according to Equation (8):

$$z_{\tau }=\frac{y_{\tau }-{\mu }_{0\tau }}{s_{0\tau }}; y_{\tau }=0; z_{\tau }=\frac{y_{\tau }-{\mu }_{1\tau }}{s_{1\tau }} y_{\tau }>0;$$

(8)

A multivariate autoregressive model, MAR(1), with constant parameters was applied, in which the best fit needs to represent the conditions of temporal and spatial dependence. The first order of the multivariate autoregressive model was determined based on the Cholesky factorization. In the same way, temporal dependence was conditioned by precipitation occurrence in Equation (9):

$${\left[\varphi \right]}_1=M_0{M}_1; {\left[\varphi \right]}_0{\left[\varphi \right]}_0{}^T=M_0-{\left[\varphi \right]}_1{M}_1{}^T; {\left[\varphi \right]}_0{\left[\varphi \right]}_0{}^T=D$$

(9)

where ${\left[\varphi \right]}_1$ is the lag 1 autoregressive coefficient matrix, ${\left[\varphi \right]}_0$ is the lag 0 autoregressive coefficient matrix, ${\left[\varphi \right]}_0{}^T$ is the transposed matrix, $M_0$ and $M_1$ are the cross-correlation matrices, and $D$ is the positive definite matrix. Finally, we determined the white noise based on the multivariate autoregressive coefficient matrix using Equation (10):

$${\epsilon }_{\tau }={\left[\varphi \right]}_0{}^{-1}\left({\left\{z\right\}}_{\tau }-{\left[\varphi \right]}_1{\left\{z\right\}}_{\tau -1}\right)$$

(10)

2.3. Evidence for the Goodness of Fit

The residual series must satisfy the residual normality ($\epsilon \cong 0, s_{\epsilon }\cong 1$), which is neither correlated $(r_k\left(\epsilon \right)\cong 0)$ nor has a biased judgment $(g_{\epsilon }\cong 0)$. The probability of the residual series is verified by the confidence interval limits as well as for the mean, skewness coefficient, standard deviation, and correlations to corroborate the residual series normality. It must satisfy the mean and correlation of the entire series within 95% of the confidence limits [63]. In the case of the standard deviation of the chi-squared test, it should comply with a 95% confidence for a normal distribution [64]. The skewness coefficient of the residual series must be within the confidence limits [63]. Another evaluated aspect in the stochastic process is the Akaike information criterion (AIC) [65].

2.4. Generation of Multivariate Synthetic Series

For the generation of synthetic series, the stochastic model is divided into two states, the precipitation occurrence $(X_t)$ and the maximum temperature and temperature range $(T_t).$ The precipitation occurrence will appear if the transition vector $(p_c)$ is a random normally distributed number greater than the random normal distribution $(n_t)$; in other words, if $n_t \le p_c$, it reaches the wet state. The stochastic temperature process starts by generating a random number with a normal distribution ($\epsilon$) and then obtains the coupled standardized series ($z_t$) and corrected using the annual low-frequency multivariate stochastic model. Finally, the inverse normalization ($y_{\tau }{}^{-1}$) was used.

The obtainment of synthetic series through the stochastic process allows for validating the developed model. The multivariate series were generated with the same characteristics as the observed data. The statistical parameters of both series were assessed and should not be significantly different; consequently, the developed model can be validated. The main metrics were determined. Mean absolute error (MAE), root mean square error (RMSE), and percent error estimate (PE) were defined.

2.5. Study Area

The Jucar River Basin is part of the Jucar Basin Agency (JBA), located in the eastern portion of the Iberian Peninsula, Spain. The basin covers an area of approximately 22,291 km². Information regarding the zone was obtained from the official website of the confederation (www.chj.es). The most relevant surface runoff is the Jucar River, which captures the surface runoff of all sub-basins [66]. The most significant reservoirs are Alarcon (1088 hm³) and Contreras (852 hm³). The river rises from the Tragacete (1600 ms.n.m) and subsequently arrives at reservoirs Toba, Alarcon, Molinar, and Tous. The study area’s limit ends where the Mediterranean Sea is reached (Figure 2). Rainfall in the Jucar River Basin has decreased since 1980 [67,68]. Temporal and spatial variation characteristics of meteorological elements in the Jucar River Basin are presented in Appendix A (Figure A1, Figure A2, Figure A3 and Figure A4).

The Jucar River Basin is divided into five sub-basins: Alarcon, Contreras, Molinar, Tous, and Huerto Mulet. The precipitation data for the study area were obtained from the Spain 02 database [69], a regular grid (20 × 20 km). The historical data have information from 1950 to 2015. The observed data were interpolated using the inverse distance weighting (IDW) method to generate a rainfall and temperature series for each sub-basin. The IDW method is one of the most common interpolation techniques [46,70,71,72].

3. Results

According to the information of the Jucar River Basin, the average annual temperature was between 17.5 °C (Contreras) and 21 °C for the years 1950–2015. In the present study, we defined four wet day thresholds (0.001, 0.01, 0.1, and 0.25 mm). The range of rainy days from October to May was between 7.3 and 12.17, and for June to September, the average precipitation occurrences were between 2.4 and 8.72 days. The months with few precipitation occurrences were between June and September, an important factor because there were little data for the stochastic modeling process. In the case of the months from October to May, the information was essential for the stochastic modeling to perform with better confidence.

3.1. Multivariate Occurrence Synthetic Series

The occurrence process was developed through a Markov model of two states. The transition probability vectors were identified, and the daily noise level can be observed. The transition probabilities for all sub-basins $p_{01\mathsf{\tau }}$ were in a range of 0.05 (minimum) to 0.45 (maximum). On the other hand, the transition probabilities $p_{11\mathsf{\tau }}$ were between 0.3 and 0.95 (according to Figure 3).

The transition vectors were evaluated with four parameters for the Fourier series. The objective of analyzing the transition probability vectors was to select the best representation of the wet–dry event. Simulations were carried out from two harmonics in the Fourier series to reach approximately 90% of the explicative variance. However, the process of rainfall occurrence can be represented by only a few parameters [9,22,29,73]; therefore, using few harmonics is acceptable. The frequency of the precipitation for four parameters represents a smoothed probability of occurrences, and with these few parameters, a good approximation of $p_{11\tau }$ and $p_{01\tau }$ can be obtained (Figure 3). On the other hand, the confidence limits for the vectors $p_{11\tau }$ and $p_{01\tau }$ were determined from the approximation to the t-distribution with 95% confidence. For the Fourier probability $p_{01\tau }$, four main patterns were observed, predominantly the Fourier probability increase in March, April, and May, after a decrease in June and July, an increase again in August and September, and finally, from October to February, the smoothed Fourier probabilities were nearly constant. A similar pattern was present for the Fourier probability $p_{11\tau }$, which had no considerable fluctuations between October and May, decreases in June and July, and increases once more in August and September. These Fourier series are sufficient for proper stochastic performance of rainfall occurrence.

3.2. Stochastic Multisite Multivariate Temperature Series

In the case of the maximum temperature and temperature range, the skewness coefficient of the historical series was near normal. For this reason, we did not consider normalization. Normality was assessed based on the skewness coefficient test for the 95% confidence limit. The daily temperature skewness coefficient was near the normal distribution, according to the confidence limits (Figure 4).

Fourier series was applied to the mean and standard deviation (${\mu }_{\tau }, s_{\tau }$), and the objective was to reduce the number of parameters (Figure 5). The series were carried out with four parameters, which provide a good fit for the model for the Jucar River Basin. In Figure 5, we can observe the mean, standard deviation, and the Fourier series for the stochastic models for wet and dry days with a 95% confidence interval. The results from fitting Fourier series for the model in Figure 5a,b reflect the smoothed mean for wet and dry days. In case of the standard deviation presented in Figure 6a,b, the fitted curve includes the original data noise.

We calculated the standardization based on the normal series, mean, and standard deviation of the fitted series. The objective of the standardization is to remove the series periodicities and to obtain a mean of zero and variance of one. On the other hand, the standardized series were determined by the multivariate autoregressive model. For the Jucar River Basin, we defined the autoregressive parameters and residual series for different wet thresholds.

We selected the wet threshold of 0.001 for determining the temperature multisite multivariate stochastic generator. First, we evaluated the normality of the residual series. Mean, deviation, skewness coefficient, and lag-one autocorrelation were calculated (

Table 1. AIC for different wet thresholds and stochastic occurrence model.

Sub-Basin	Wet Day Threshold (mm)
	0.001 *	0.01	0.10	0.25
Alarcon	−590.2	−515.3	−362.2	−310.9
Contreras	−681.5	−551.2	−459.5	−421.3
Molinar	−562.3	−420.7	−261.2	−215.1
Tous	−587.4	−463.6	−380.5	−340.0
Huerto Mulet	−610.5	−554.7	−467.3	−427.3

* Best performance.

).

The autocorrelation for this series was also determined within the 95% confidence limits for both models’ maximum temperature and temperature range. In addition, we applied different tests to confirm that the residual series can be considered a normal distribution with a mean of zero, variance of one, and skewness coefficient of zero (Figure 7). The residual series of the two developed stochastic models were very similar.

The residual series of the MAR(1) for the maximum temperature had a mean near zero, the variance was around 0.85, the skewness coefficient was between −0.245 and 0.026, and the lag-one autocorrelation was around $\pm$0.01. The skewness coefficient and autocorrelation are within the confidence limits, the average is within 99% of the confidence limits, and therefore we assumed the normality of the residual series (

Table 2. Normality analysis for residual series for M1 and M2 (wet threshold 0.001).

Model	Statistical/Sub-Basin	Alarcon	Contreras	Molinar	Tous	Huerto Mulet
1 *	Mean	−6.7 × 10−5	−2.0 × 10−4	7.5 × 10−5	−2.6 × 10−4	5.6 × 10−5
	Deviation	0.842	0.821	0.834	0.812	0.850
	Skewness coefficient	−0.185	−0.242	−0.245	−0.089	0.026
	Lag-one autocorrelation	0.005	0.003	0.009	−0.041	−0.042
	AIC	−8369	−9508	−8814	−10,152	−7958
2 **	Mean	−3.09 × 10−4	−5.98 × 10−4	−6.83 × 10−5	−2.70 × 10−4	−2.65 × 10−4
	Deviation	0.937	0.920	0.942	0.900	0.942
	Skewness coefficient	−0.009	−0.038	−0.044	0.112	−0.028
	Lag-one autocorrelation	0.036	0.043	−0.030	−0.082	−0.039
	AIC	−6471	−8233	−5957	−10,294	−5896

* Maximum temperature. ** Temperature range.

).

On the other hand, the series was also considered stationary since it complied with ${\varnothing }_1<1$ for all sub-basins. The AIC value was −10,152 with the examined parameters.

Similar results were obtained for the stochastic model for temperature range: the residual series had a mean of around −0.0005, variance of at least 0.92, skewness coefficient between −0.044 and 0.112, and lag-one autocorrelation of about $\pm$0.08 (Figure A5). For this stochastic model, we assumed normality and stationarity of the residual series as well as for stochastic Model 1 (Table 2). According to the AIC, the stochastic maximum temperature model (M1) was similar to the stochastic temperature range (M2) for all sub-basins.

3.3. Generation of Multivariate Synthetic Temperature Series

For the process of generating synthetic series, 1000 series were created considering the same length as the sample (66 years). The statistical sum, mean, standard deviation, and skewness coefficient were determined for both the synthetic and historical series. The occurrence depends on the correlated multivariate precipitation probabilities, critical probability, and normal distribution. Multivariate synthetic series of rainfall occurrences were generated for the five sub-basins using the stochastic process. The same occurrence series was applied for both Model 1 (M1) and Model 2 (M2) to avoid the bias from the MAR(1) model.

For the historical series, the sum of the number of rainy days in 66 years was calculated. In the same way, the monthly occurrences of the synthetic and the historical series were determined. Several statistical tests were performed to validate that the generated series are not substantially different from the historical period. The tests applied to the precipitation occurrence were the k-s test to verify that the results come from the same distribution, the t-test for equality of means, and the Wilcoxon test for equality of medians. The tests were applied considering a 95% reliability, and it was concluded that there is insufficient evidence to refute that the generated precipitation occurrences and historical series, both monthly and daily, are significant.

The scatter plots of daily mean occurrences (66 years) of historical versus simulated data can be seen in Figure 8 for the five sub-basins. The daily rainfall occurrences varied $\pm$5 days from the 1:1 line (Figure 8a). The monthly rainfall occurrences varied by $\pm$0.2 days. Monthly occurrences provide better results than daily occurrences due to the number of parameters applied (four in total). We used the same simulated rainfall occurrence (four parameters) to generate the maximum temperature and temperature range.

The statistical mean, standard deviation, and skewness coefficient of daily data were computed. The MAR(1) for mean daily temperatures had a deviation of ±1 °C, which is more accurate than the stochastic Model 2 for the daily average temperature range (±1.5 °C). Model 1 achieved better results on both a daily and a monthly scale (Figure 9). The observed and generated series were not significantly different according to the k-s test, which indicates that they originate from the same distribution, in addition to sharing the same average according to the t-test and the same median according to the Wilcoxon test. These tests were applied to the daily average temperatures and the monthly averages. The best approximations were provided for the monthly data with higher reliability. Moreover, the RMSE, MAE, and PE (

Table 3. Performance analysis for Model 1 and Model 2.

Parameter	Model 1 (M1)					Model 2 (M2)
	1	2	3	4	5	1	2	3	4	5
RMSE (°C/day)	1.881	1.107	1.392	1.156	0.767	1.786	1.019	1.290	1.078	0.780
MAE (°C/day)	1.503	0.820	1.092	0.896	0.572	1.455	0.773	1.021	0.852	0.582
PE (%)	0.043	0.027	0.017	0.031	0.021	0.035	0.049	0.040	0.025	0.012

) display the adequate performance of Model 1 (M1) and Model 2 (M2).

The f-test was used for the standard deviation (Figure 10), which indicates the equality of deviations for both models and showed minor differences. Both the maximum temperature and temperature range performed well in the daily standard deviation. It is worth mentioning that the deviations were overestimated concerning the observed data. For a monthly scale, both models effectively reproduced the standard deviation.

Regarding the skewness coefficient, MAR(1) was underestimated for the observed data. The same dispersion was present for the two models (M1 and M2) and for the daily and monthly averages (Figure 11).

Even though the normalization’s function can be adjusted on average for the confidence limits, this offers a variation for the observed skewness. The skewness of observed daily data was between −1.2 and 1.2, and the simulated was between −0.5 and 0.5. Therefore, the multivariate stochastic model underestimated the skewness coefficient by less than −0.5 and more than 0.5. On a monthly scale, the skewness of the simulated precipitation distribution was underestimated similarly to the daily scale.

For monthly temperature, the multisite multivariate model preserved the main statistics. Figure 12 presents all months for the 5 sub-basins and 66 years, in which the temporal and spatial dependence was adequately performed. For maximum temperature, the values were between 5.5 and 35 °C, with variability of ±2 °C (Figure 12a). In the case of temperature range, the values were between 5 and 25 °C with the same variability. For the monthly mean of all sub-basins, the error was only ±0.1 °C.

The multisite multivariate stochastic autoregressive Model 1 was selected to compare temperature years with regard to the observed data. The stochastic model can represent the temporal tendency of the results, providing an adequate indication of the yearly temperatures (Figure 13). Due to the design of the multisite multivariate stochastic model corrected by the annual model, we can simulate low frequency. The interannual variability represents the autocorrelation and cross-correlations of observed data (

Table 4. Annual cross-correlation matrix for all sub-basins.

Maximum Temperature Cross-Correlation
Sub-Basin	Alarcon	Contreras	Molinar	Tous	Huerto M	* Alarcon	* Contreras	* Molinar	* Tous	* Huerto M
Alarcon	1.000					1.000
Contreras	0.833	1.000				0.834	1.000
Molinar	0.760	0.797	1.000			0.765	0.819	1.000
Tous	0.239	0.492	0.598	1.000		0.243	0.489	0.610	1.000
Huerto M	0.223	0.371	0.390	0.802	1.000	0.230	0.377	0.398	0.800	1.000
Temperature Range Cross-Correlation
Sub-Basin	Alarcon	Contreras	Molinar	Tous	Huerto M	* Alarcon	* Contreras	* Molinar	* Tous	* Huerto M
Alarcon	1.000					1.000
Contreras	0.719	1.000				0.718	1.000
Molinar	0.891	0.592	1.000			0.895	0.590	1.000
Tous	0.563	0.494	0.754	1.000		0.565	0.499	0.759	1.000
Huerto M	0.504	0.331	0.710	0.918	1.000	0.502	0.333	0.708	0.917	1.000

* Simulated cross-correlation.

). The results indicate that variability was well-reproduced for the stochastic process. Moreover, the maximum and minimum values were performed adequately. For annual temperature, the stochastic model can produce the variability of both temperatures. The variability expressed for the model (sim) was greater than the observed data. Accordingly, the MAR(1) can define the maximum and minimum temperature values.

4. Discussion

The multisite multivariate autoregressive stochastic model (MASCV) was developed using MATLAB and was verified for different stations within the same basin with similar results. A Markov model of two states for the multivariate precipitation occurrence and the conditioned multivariate stochastic model for temperature can represent spatial and temporal parameters to study on-site conditions on daily, monthly, and annual scales. The Markov model of two states with few parameters has been able to depict the precipitation occurrence process. On the other hand, the temperature was accomplished considering normal distribution. Therefore, a reduction of parameters was achieved in generating the temperature. This represents a critical simplification in obtaining the daily temperature, which according to the performed parameterization, provides acceptable results for the Jucar River Basin. Moreover, it simplifies the complexity and reduces computational time.

The stochastic multivariate autoregressive model with few parameters adequately reproduced the daily and monthly temperatures. The tests showed insignificant differences between the observed and generated temperatures.

The primary objective of this stochastic model is to determine the monthly runoff and incorporate it into the integrated water resources management. It is noteworthy that this validation is performed on a monthly scale. Therefore, MASCV must be able to represent the statistics of the precipitation occurrence and temperatures in different timescales. The developed stochastic MAR(1) adequately reproduced the main statistics.

5. Conclusions

The multivariate autoregressive model of climate variables (MASCV) is a daily stochastic weather generator, programmed in MATLAB, with several user advantages, i.e., wet day selection, number of harmonics for each case, normalization type, and synthetic series to generate and automatize graphs. Moreover, MASCV can generate extreme temperatures over extended periods of time. This is unique in stochastic weather generators because only a few can reproduce extreme events. Furthermore, MASVC can be used for bias correction for climate change studies, in this case, perturbing parameters according to climate models. Finally, the results of MASCV can be incorporated for environment analysis.

MASCV presents the completion of a multivariate model for precipitation occurrences, i.e., a Markov model of two states and the dependence of temperature with rainfall occurrence process. This multisite multivariate stochastic model is meant to become a beneficial tool in a simplified manner, which may allow the incorporation of different climatic and hydrological variables.

A Markov model of first order can reflect the time dependence of the precipitation occurrence, preserving comparative statistical data with the historical series. Moreover, the spatial and temporal structure using a stochastic multisite multivariate model and coupling daily and annual temperature correction reproduced adequately in the different timescales.

Models M1 and M2 were suitably performed for different temperatures. For wet and dry days, the multisite multivariate stochastic model can adjust to real dependence with precipitation occurrence and spatial and temporal dependence of daily, monthly, and annual temperatures.

This approach greatly simplifies the process of simulating precipitation, which implies a considerable advantage and versatility over other stochastic generators. The reduction of parameters is an important factor addressed in this approach for determining the temperature and considering continuous modeling for days, months, and years.