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

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.


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 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 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.

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): where T t is the temperature model, X t is the precipitation occurrence model, and Y t is the whole stochastic process.

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 (X t = 0). 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 ( ) depends on the previous day ( ); if it is dry, = , otherwise = . The boundary for the precipitation occurrence process is determined for a random spatially correlated normal matrix = ∅ [0,1]. Only one day is computed if the random variable is equal to or less than a critical probability. The normal critical probability (p c ) depends on the previous day (X t−1 ); 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 = ∅ −1 u t [0, 1]. Only one day is computed if the random variable is equal to or less than a critical probability.

of 22
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, p 01τ = P 01 τ and p 01τ = P 11 τ .
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τ = ∅ −1 (p 11τ ) and p n01τ = ∅ −1 (p 01τ ). The spatial correlation between probability vectors with lag-k is shown in Equation (3): where r k (i,j) is the spatial correlation with lag-k, p n11 (i,j) is the normalized probability vector of two consecutive days of rainfall occurrence, and p n01 (i,j) 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): 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 = [M] [M] . Multiplying this lower triangular matrix by a random normal matrix results in a random spatially correlated normal matrix, n, in Equation (5): where n is the random spatially correlated normal matrix, M (n,n) is the lower triangular matrix for n series, and N (n,m) 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.

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]. 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.
where u is the normalized temperature mean for wet (µ 1τ ) and dry days (µ 0τ ) and standard deviation for wet (s 1τ ) and dry days (s 0τ ). 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 (µ 1τ , µ 0τ , s 1τ , s 0τ , p 01τ , p 11τ ). 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 (µ 1τ , µ 0τ ) and standard deviation (s 1τ , s 0τ ). We determined the standardized series (z τ ) at a daily scale, canceling the mean and standard deviation of the normalized series according to Equation (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): 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):

Evidence for the Goodness of Fit
The residual series must satisfy the residual normality (ε ∼ = 0, s ε ∼ = 1), which is neither correlated (r k (ε) ∼ = 0) nor has a biased judgment (g ε ∼ = 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].

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 ≤ p c , it reaches the wet state. The stochastic temperature process starts by generating a random number with a normal distribution (ε) and then obtains the coupled standardized series (z t ) and corrected using the annual low-frequency multivariate stochastic model. Finally, the inverse normalization (y τ −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.

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 2 . 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 3 ) and Contreras (852 hm 3 ). 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 ( Figures A1-A4). the residual series must be within the confidence limits [63]. Another evaluated aspect in the stochastic process is the Akaike information criterion (AIC) [65].

Generation of Multivariate Synthetic Series
For the generation of synthetic series, the stochastic model is divided into two states, the precipitation occurrence ( ) and the maximum temperature and temperature range ( ). The precipitation occurrence will appear if the transition vector ( ) is a random normally distributed number greater than the random normal distribution ( ); in other words, if ≤ , it reaches the wet state. The stochastic temperature process starts by generating a random number with a normal distribution ( ) and then obtains the coupled standardized series ( ) and corrected using the annual low-frequency multivariate stochastic model. Finally, the inverse normalization ( ) 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.

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 2 . 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 3 ) and Contreras (852 hm 3 ). 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 ( Figures A1-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].  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].

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.

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τ were in a range of 0.05 (minimum) to 0.45 (maximum). On the other hand, the transition probabilities p 11τ were between 0.3 and 0.95 (according to Figure 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.

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 were in a range of 0.05 (minimum) to 0.45 (maximum). On the other hand, the transition probabilities 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 and can be obtained ( Figure 3). On the other hand, the confidence limits for the vectors and were determined from the approximation to the t-distribution with 95% confidence. For the Fourier probability , 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 , which had no considerable fluctuations between October and May, decreases in June and 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τ and p 01τ can be obtained (Figure 3). On the other hand, the confidence limits for the vectors p 11τ and p 01τ were determined from the approximation to the t-distribution with 95% confidence. For the Fourier probability p 01τ , 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τ , 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.

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).

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 ( , ), 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.   Fourier series was applied to the mean and standard deviation (µ τ , s τ ), 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).   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. 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). 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 ±0.01. The skewness coefficient and 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 ±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). On the other hand, the series was also considered stationary since it complied with ∅ 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 ±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.

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 ±5 days from the 1:1 line (Figure 8a). The monthly rainfall occurrences varied by ±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) 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. 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)    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.
(a) (b) Figure 9. Scatter plots for observed mean versus generated temperature (mean 66 years observed and 1000 simulated series) for two models: (a) maximum temperature for each calendar day (blue) and (b) temperature range for each calendar day (red).
(a) (b) Figure 10. Scatter plots for observed standard deviations versus generated (mean 66 years observed and 1000 simulated series) for two models: (a) maximum temperature for each calendar day (blue) and (b) temperature range for each calendar day (red).
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. 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. 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  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 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). 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.
corrected by the annual model, we can simulate low frequency. The interannual variability represents the autocorrelation and cross-correlations of observed data ( Table 4). 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.

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

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.

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.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
The spatial and temporal variation of the characteristics of meteorological elements in the Jucar River Basin are presented in Figures A1-A4, showing the spatial distribution of annual elements based on coordinates and 16 meteorological stations in Jucar River Basin. The grid interpolation uses inverse distance weight interpolation (IDW). Figure A5 shows the daily lag correlation for residual series.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
The spatial and temporal variation of the characteristics of meteorological elements in the Jucar River Basin are presented in Figures A1-A4, showing the spatial distribution of annual elements based on coordinates and 16 meteorological stations in Jucar River Basin. The grid interpolation uses inverse distance weight interpolation (IDW). Figure A5 shows the daily lag correlation for residual series.