Probabilistic Gaussian Copula Regression Model for Multisite and Multivariable Downscaling

By Ouarda, T B M J

Proquest LLC

ABSTRACT

Atmosphere-ocean general circulation models (AOGCMs) are useful to simulate large-scale climate evolutions. However, AOGCM data resolution is too coarse for regional and local climate studies. Downscaling techniques have been developed to refine AOGCM data and provide information at more relevant scales. Among a wide range of available approaches, regression-based methods are commonly used for downscaling AOGCM data. When several variables are considered at multiple sites, regression models are employed to reproduce the observed climate characteristics at small scale, such as the variability and the relationship between sites and variables. This study introduces a probabilistic Gaussian copula regression (PGCR) model for simultaneously downscaling multiple variables at several sites. The proposed PGCR model relies on a probabilistic framework to specify the marginal distribution for each downscaled variable at a given day through AOGCM predictors, and handles multivariate dependence between sites and variables using a Gaussian copula. The proposed model is applied for the downscaling of AOGCM data to daily precipitation and minimum and maximum temperatures in the southern part of Quebec, Canada. Reanalysis products are used in this study to assess the potential of the proposed method. Results of the study indicate the superiority of the proposed model over classical regression-based methods and a multivariate multisite statistical downscaling model.

(ProQuest: ... denotes formulae omitted.)

1. Introduction

Atmosphere-ocean general circulation models (AOGCMs) are commonly used to simulate large-scale climate evolution. Information provided by these models is widely used to produce future climate projections. AOGCM data are generally produced on regular grids with a low horizontal resolution around 2.58 longitude and latitude (approximately 250-300 km). However, this resolution is coarse for regional and local climate studies. Downscaling techniques have been developed to refine AOGCM data and provide information at more relevant scales. These techniques can be classified into two main categories: dynamic methods and statistical methods (Wilby and Wigley 1997; Herrera et al. 2006). Dynamic methods use regional climate models (RCMs), which have the same basic principles as AOGCMs, with a high resolution between 25 and 50 km. RCMs only cover a limited portion of the globe. Dynamic methods require large computational capabilities and substantial human resources. Statistical methods, on the other hand, consider statistical relationships between large-scale variables (predictors) and small-scale variables (predictands). The main advantages of these methods are their sim- plicity and their low computational costs. Thus they represent a good alternative to dynamic methods in the case of limited resources.

The performance of a statistical downscaling model depends on its ability to reproduce the observed statis- tical characteristics of local climate (Wilby 1998; Wilby et al. 2002; Gachon et al. 2005; Hessami et al. 2008). Statistical downscaling methods are required to take into account the climatic characteristics of predictands, such as their observed variability, in order to provide reliable meteorological information at the local scale (Wilby and Wigley 1997). The proper reproduction of the variability in downscaling applications is a very im- portant issue, since a poor representation of the vari- ability could lead to a poor representation of extreme events. In addition, if data are required at multiple sta- tions, then a spatial or multisite model is employed to better represent the observed correlation between sites and predictands.

Once a downscaling model has been developed, downscaled data can be entered as input into an envi- ronmental model-for example, a hydrological model for streamflow in a watershed that requires climate in- formation at a finer scale (Cannon 2008). Precipitation and temperature are commonly considered as predictands in a downscaling problem. In hydrology, streamflows depend strongly on the spatial distribution of precipi- tation in a watershed, and on the interactions between temperature and precipitation, which determines whether precipitation falls as rain or snow (Lindstroeurom et al. 1997). Therefore, maintaining realistic relationships between sites and variables in downscaled results is particularly important for a number of applications such as hydro- logical modeling.

Statistical downscaling techniques can be grouped into three main approaches: stochastic weather generators (Wilks and Wilby 1999), weather typing (Conway et al. 1996), and regression methods (Wilby et al. 2002). Re- gression based methods are commonly used for down- scaling AOGCM data. They can find a direct relationship between large-scale predictors and local predictands. A variety of regression methods have been used in the literature for downscaling purposes, such as multiple linear regression (MLR) (Wilby 1998; Hellstroeurom et al. 2001; Huth 2004; Hessami et al. 2008; Jeong et al. 2012a), principal component analysis (PCA) with MLR (Huth 2004), canonical correlation analysis (CCA) combined with MLR (Huth 2002; Palutikof et al. 2002), artificial neural networks (ANNs) (Schoof and Pryor 2001; Cannon 2008, 2011), and singular value decomposition (SVD) (Widmann et al. 2003). Regression-based methods usually perform well for downscaling purposes, but their major drawback is that they generally provide only the mean or the central part of the predictands (Cawley et al. 2007). Therefore the variance of the modeled mean will typically be smaller than the variance of the observed series. The reason is in part that regression models cannot represent the influence of small-scale phenomena. On the other hand, stochastic weather generators use a random num- ber conditioned upon large-scale model output state and can be applied at a single site as well as at multiple sites (Wilks1998,1999;Qianetal.2002;Palutikofetal.2002). However, the main limitation of weather generators is the difficulty in adjusting the parameters in a physically re- alistic and consistent manner under future climate states (Wilby et al. 1998).

To correctly estimate the temporal variance of down- scaled data series, three main approaches have been proposed in the literature: inflation (Huth 1999), ran- domization (von Storch 1999; Clark et al. 2004), and ex- pansion (Burger and Chen 2005). Inflation is usually carried out by multiplying the downscaled data by a con- stant factor, randomization consists in adding a random noise, and in the expansion approach the predicted vari- ances are constrained to match the variance of the ob- served data by adding a constraint to the regression cost function. A problem with the inflation approach is that the spatial correlations between sites can be mis- represented. However, randomization can also be applied in a multisite downscaling framework. In a randomiza- tion procedure, regression and stochastic weather generator approaches can be combined in a single hy- brid model. Thus, the resulting stochastic hybrid model can overcome weaknesses of both approaches. For in- stance, Jeong et al. (2012b) employed a randomization procedure to reproduce the cross-site correlation of pre- cipitation occurrence and amount among the observation sites using the multivariate normal distribution. As well, Jeong et al. (2012c) proposed a multivariate multisite statistical downscaling model (MMSDM) for simulta- neous downscaling of climate variables including daily maximum and minimum temperatures for multiple ob- servation sites. The MMSDM employs multivariate multiple linear regression (MMLR) to simulate deter- ministic series from large-scale reanalysis data and adds spatially correlated random series to the deterministic series of the MMLR to complement the underestimated variance and to reproduce the spatial correlation of var- iables from multiple sites and an at-site temporal corre- lation between variables. In the same way, Khalili et al. (2013) proposed a hybrid approach by combining a linear regression component with a stochastic component based on a spatial moving average process to reproduce the observed spatial dependence between extremes of tem- peratures at different sites. However, Burger and Chen (2005) indicated that randomization in a hybrid approach is based on a static noise model and thus it failed to represent local changes in atmospheric variability in a climate change simulation, which is well explained by expended downscaling. The expanded downscaling can be applied to multiple predictands by constraining the covariance matrix of the predicted series to be equal to the observed covariance matrix. For example, Cannon (2009) introduced a multivariate ridge regression with negative ridge parameters for improving the covariance structure of multivariate linear downscaling models. This procedure is conceptually similar to expanded down- scaling since both force the covariance structure of the predictions to match that of the observations. Thereby, the observed variability is reproduced using a single de- terministic regression component, and thus, the repro- duced variability is not static and may change in a climate change simulation. On the other hand, von Storch (1999) suggested that the inflation and expansion approaches are inappropriate techniques, because the implicit assump- tion that all local variability can be traced back to the large-scale using a deterministic regression model is im- proper and is not the case in reality. For this reason, von Storch (1999) suggests using randomization approaches that are more realistic than expanded approaches.

As an alternative to the three existing techniques for specifying predictand variance, the variability of pre- dictand can also be captured by modeling the whole distribution. In this regard, probabilistic approaches have provided significant contributions in downscaling applications (e.g., Bates et al. 1998; Hughes et al. 1999; Bellone et al. 2000; Vrac and Naveau 2007; Cannon 2008; Fasbender and Ouarda 2010). They allow re- producing the whole distribution by modeling the effect of the AOGCM predictors on the parameters of the predictand distributions. For example, Williams (1998) employs an ANN to model parameters of a mixed Bernoulli-gamma distribution for precipitation at a sin- gle site. For this purpose, the problem that arises is how to extend probabilistic approaches in multisite down- scaling tasks. In this context, Cannon (2008) applied the principles of expanded downscaling in a probabilistic modeling framework to allow for a realistic represen- tation of spatial relationships between precipitation at multiple sites. On the other hand, as indicated pre- viously, expanded approaches are not appropriate in a climatic downscaling problem.

Given the disadvantages of these three existing tech- niques for specifying predictand variances and covari- ance structure, we propose in this paper a copula-based approach as an alternative solution to extend the prob- abilistic modeling framework in multisite downscaling tasks. The methodology proposed here presents another advantage in that multiple climatic variables are down- scaled at multiple sites simultaneously and consistently to produce realistic relationships between both sites and variables.

Multivariate dependence structures can be modeled using classical distributions such as the multivariate normal distribution. However, the multivariate normal approach cannot adequately reproduce the dependence structure of hydrometeorological data when existing asymmetries are significant, such as for precipitation var- iables. To address this limitation, the use of copula-based approaches can be beneficial since copula functions are more flexible and may be better suited to the data (e.g., Schoeurolzel and Friederichs 2008). Copulas have recently become very popular, especially in fields such as econo- metrics, finance, risk management, and insurance. In re- cent years, the application of copulas has also made significant contributions in the field of hydrometeorology. Schoeurolzel and Friederichs (2008) provide a brief overview of copulas for applications in meteorology and climate research. Models based on copulas have been introduced for multivariate hydrological frequency analysis (Chebana and Ouarda 2007; El Adlouni and Ouarda 2008), risk assessment, geostatistical interpolation, and multivariate extreme values (e.g., De Michele and Salvadori 2003; B^ardossy 2006; Renard and Lang 2007; Kazianka and Pilz 2010). In addition, copulas allow one to describe the de- pendence structure independently from the marginal dis- tributions, and thus use different marginal distributions at the same time without any transformations (Sklar 1959; Dupuis 2007).

The goal of the present study is to develop and test a probabilistic Gaussian copula regression (PGCR) model for multisite and multivariable downscaling. The model can reproduce observed spatial relationships between sites and variables and specify at each site and for each day, the conditional distributions of each vari- ables. To this end, PGCR uses a probabilistic frame- work to address the limitations of regression-based approaches, namely (i) the poor representation of ex- treme events, (ii) the poor representation of observed variability, and (iii) the assumption of normality of data. The PGCR model specifies a marginal distribution for each predictand through AOGCM predictors as well as a Gaussian copula to handle multivariate dependence between margins. The proposed model can be consid- ered as a hybrid approach combining a probabilistic regression-based downscaling model with a stochastic weather generator component. The main advantage of this model compared to conventional hybrid approaches is that the temporal variability may change in future climate simulations.

The present paper is structured as follows: After a presentation of the MMLR model and the MMSDM as a classical regression-based method in statistical down- scaling, the proposed PGCR model is presented. The PGCR model is then applied to the case of daily pre- cipitations and maximum and minimum temperatures in the southern part of the province of Quebec, Canada. Reanalysis data are used in order to assess the potential of the proposed method. After the calibration of PGCR model, an independent dataset is used to assess the downscaling quality. Results are compared with those obtained with MMLR and MMSDM using statistical criteria and climatic indices that describe the frequency, intensity, and duration of the variables of interest. Fi- nally, a discussion and conclusions are provided.

2. Methodology

The MMLR model, the MMSDM, and the PGCR model are presented respectively in sections 2a, 2b,and 2c. The probabilistic framework for the PGCR model is presented with a description of the different selected marginal distributions for each predictand. Then, a sim- ulation procedure is presented using Gaussian copulas to produce the dependence structure between several pre- dictands at multiple sites.

a. Multivariate multiple linear regression

Statistical downscaling from multiple AOGCM pre- dictors to numerous meteorological observation sites is one of the situations where predictions of several de- pendent variables are required from a set of independent variables. Multivariate regression approaches have been used in many scientific areas in order to analyze rela- tionships between multiple independent variables and multiple dependent variables (Jeong et al. 2012b).

We consider s meteorological stations where, for each site ( j 5 1, ..., s), we consider three predictands: max- imum daily temperature Tmaxj, minimum daily tem- perature Tminj, and daily precipitation Precj. The precipitation data Precj is of particular interest since the responses are the product of an occurrence process, which decides whether or not there is any rainfall on a particular day, and an amount process, which governs the amount of rainfall, given that some rainfall is ob- served. Therefore, Precj can be decomposed into two separate variables: one for precipitation occurrence Pocj and one for wet-day precipitation amount Pamj. To define wet days, the occurrence was limited to events with a precipitation amount larger than or equal to 1mmday21 to avoid problems associated to trace mea- surements and low daily values. A dry day is defined as a day having less than 1 mm of precipitation (Jeong et al. 2012b).

The precipitation amount vector Pamj for a site j is not normally distributed. Indeed, the gamma distribution has often been found to provide a good fit to rainfall amounts in many studies (e.g., Stephenson et al. 1999; Yang et al. 2005). Thus, appropriate transformation should be performed before developing a regression- based precipitation amount model. Yang et al. (2005) proposed the Anscombe transformation for transforming the precipitation amount to a normal distribution. If the vector Pamj is gamma distributed, the distribution of Rij5 Pami1j/3 on a day i at a site j,wheretheR is the Anscombe residuals, is normal (e.g., Yang et al. 2005; Jeong et al. 2012b).

We denote A the matrix grouping all the predictand vectors Tmax, Tmin, Poc, and R, of dimension n 3 m where m 5 4s:

... (1)

We define X as the n 3 l dimensional matrix that con- tains the multiple predictor variables. One can estimate the parameter matrix B of dimensions l 3 m, which can define the linear relationship between the two matrices X and A . Therefore, the MMLR can be expressed as

... (2)

where E is the residual matrix of dimensions n 3 m. The parameter matrix B can be estimated using the ordinary least squares (OLS) method, which is given by

... (3)

Under the assumption that the errors are normally dis- tributed, B^ by the OLS is the maximum likelihood es- timator. Then the deterministic series of predictands can be obtained using the following MMLR equations and atmospheric predictors:

... (4)

Note that the wet day was determined when the de- terministic series of daily probability of precipitation occurrence by the MMLR occurrence model was larger than 0.5. The parameter matrix B^ is affected by multi- collinearity, which produces large standard errors of estimated parameters (Jeong et al. 2012b). A number of methods have been employed in order to limit the in- fluence of multicollinearity, such as ridge regression, principal component analysis (Fasbender and Ouarda 2010), canonical correlation regression (Huth 2004), stepwise regression, and lasso regression (Hammami et al. 2012). In this study, principal component analysis was employed to deal with the multicollinearity problem.

b. Multivariate multisite statistical downscaling model

The MMLR predicts only deterministic components explainable by linear regression and the independent atmospheric variables X for the different multiple sites. As mentioned previously, this deterministic compo- nent underestimates the temporal variability of each pre- dictand and cannot adequately reproduce the correlation between sites and variables. To this end, a classical sto- chastic randomization procedure is commonly employed by adding correlated random series to the deterministic component. The resulting model is a MMSDM (Jeong et al. 2012c). The MMSDM employs MMLR to simu- late deterministic series from large-scale reanalysis data and adds spatially correlated random series to the deterministic series of the MMLR to complement the underestimated variance and to reproduce the corre- lation between sites and variables.

The residual (or error) matrix E [n 3 m] of the MMLR model is described as

... (5)

For the MMSDM, correlated random noise among the predictands at multiple sites was generated from multi- variate normal distribution and added to the deter- ministic series of the MMLR.

The cross-correlated error matrix H [n 3 m] is generated from a multivariate normal distribution having zero error mean and an error covariance matrix S 5 SCS equal to that of the residual matrix E [H ; Nm (0, S)], where S is a diagonal matrix of standard deviations and C is a corre- lation matrix of the residual matrix E. Generated residuals were then added to the downscaled predictand as

... (6)

The residual vectors for precipitation amount at each site may be not normally distributed and may be skewed. To overcome this problem, a probability distribution mapping technique was adapted and the generated precipitation amount was adjusted using the gamma distribution.

c. Probabilistic Gaussian copula regression

1) PROBABILISTIC REGRESSION

In most applications, regression models are performed to describe a mapping that approximates the conditional mean of the prediction and data. This mapping is ap- propriate if the data are generated from a deterministic function that is corrupted by a normally distributed noise process with constant variance (Cannon 2008). When the noise process has nonconstant variance or is nonnormal, it is more appropriate to use a model that fully describes the conditional density of the predictand in a probabilistic framework. Thus, the distribution of each predictand at the observed sites must be represented by an appro- priate probability density function (PDF), and then we employ a regression model with outputs for each pa- rameter in the assumed PDF of noise process. In this paper, the normal distribution is chosen for the tem- perature variables. According to Dorling et al. (2003), for a normally distributed noise process with non- constant variance, the conditional density regression would have two outputs: one for the conditional mean and one for the conditional variance. For Tmax and Tmin, the model is described by

... (7)

... (8)

... (9)

and

... (10)

where E(.) is expectation and x(t) is the value of pre- dictors at the day t, and coefficients amax j, aminj, bmaxj, and bmin j are estimated separately. Then the conditional normal PDF of Tmaxj for a day t is given by

... (11)

And the conditional normal PDF of Tminj for a day t is given as

... (12)

Note that for Pocj, a standard problem is that the dry-wet dichotomy leads to a Bernoulli process. In this case we use a logistic regression given by

... (13)

where pj (t) is the probability of precipitation occur- rence at a site j on a day t and cj isthecoefficientofthe logistic model. Thus the conditional distribution of Pocj is given by

... (14)

The Pamj is modeled through a conditional gamma dis- tribution with shape parameters aj(t) and scale parame- ter bj (t) given by (Cannon 2008)

... (15)

where dj and ej are the coefficients of the model. Thus, the conditional gamma PDF for Pamj on a day t is given by

... (16)

where G(.) is the gamma function.

Finally, let us define the random matrix Y of dimen- sion n 3 m grouping each predictand at each site with

... (17)

It is important to mention that the matrix of predictand Y is different from the matrix A of Eq. (1), since the latter contains the transformed variables for the pre- cipitation amount. The conditional PDF ftk [yk (t) j x(t)] at the time t of the kth element of Y, where k 5 1, ..., m, is then given by Eq. (9) if yk is a Tmax, Eq. (10) if it is a Tmin, Eq. (12) if it is a Poc, and Eq. (14) if it is a Pam. Then all coefficients amaxj, aminj, bmaxj, bminj, cj, dj, and ej for all sites are set following the method of maximum likelihood by minimizing the negativelogpredictiveden- sity (NLPD) cost function (Haylock et al. 2006; Cawley et al. 2007; Cannon 2008):

... (18)

This is carried out via the simplex search method (Lagarias et al. 1998). This is a direct search method that does not use numerical or analytical gradients.

2) CONDITIONAL SIMULATION USING GAUSSIAN COPULA

Once the proposed probabilistic regression model has been trained, it can be used to estimate the PDF of each predictand for a given day when we have the AOGCM predictors. Then, it is possible to create synthetic pre- dicted series of each predictand by sampling in the obtained PDF on each day. In these steps, it is important to maintain realistic relationships between sites and var- iables. Indeed, consistency of predictions between sites and variables is very important, particularly in hydro- logical modeling.

Reproducing the relationships of multiple random variables when each variable is normally distributed is possible by using the multivariate normal distribution. However, this is not the case in a number of studies such as the present one (Pam and Poc). Indeed, Pam vari- ables are not usually normally distributed. In this regard, Pitt et al. (2006) proposed a Gaussian copula framework that can describe the dependence part of the model but allows the margins to be normal or not, discontinuous or continuous.

A copula is a multivariate distribution whose marginals are uniformly distributed on the interval [0, 1]. The multivariate function C is called a copula if it is a contin- uous distribution function and each marginal is a uniform distribution function on [0, 1]; that is, C[0, 1]q /[0, 1] with

... (19)

in which each Ui ;Un(0, 1) and u5(u1, . . . , uq). If C is a Gaussian copula, then

... (20)

where F is the standard normal cumulative distribution function and Fq (w; C) is the cumulative distribution function for a multivariate normal vector w having zero mean and covariance matrix C.FollowingPitt et al. (2006), we use latent variables to transform the marginal distributions of each predictand to a standard normal distribution. The dependence structure between predic- tands is reproduced by assuming a multivariate Gaussian distribution for the latent variables z(t) ' Nm (0, C), using the following equations:

... (21)

Note that Poc is a discrete variable for which the cumu- lative distribution function Ftk is discontinuous. Thus, in ordertomapzk (t) onto the full range of the normal distribution, the cumulative probabilities Ftk [yk (t)] for Poc are randomly drawn from a uniform distribution on [0, 1 2 p(t)] for dry days and [1 2 p(t), 1] for wet days. Finally, the goal in this step of calibration is the estima- tion of the copula parameter, which is the correlation matrix C of the latent variables z(t). Then, for a new day t0 it is possible to generate y 5 [y1, ..., y m ]thatcan be randomly generated by the Gaussian copula C(z; C). Figure 1 illustrates how the probabilistic regression model and the Gaussian copula are combined to produce one simulation at a day t0 . The value of the synthetic predictand time series y^k(t0)isgivenbyFt2k 1[yk], except for precipitation occurrences. For precipitation occur- rences, y^k (t0 ) 5 1ifyk is greater than 1 2 pk (t0 ), and y^k(t0)50ifyk is lessthan 12pk(t0).

For the application of a copula model, Vogl et al. (2012) and Laux et al. (2011) mentioned that it is an indispensable prerequisite that the marginals are in- dependent and identically distributed (iid). If this is not the case, an appropriate transformation has to be ap- plied to the data to generate iid variates. Thereby, for the PGCR model, all marginal distributions obtained from the probabilistic regression model for a given day areassumedtobeiid.

3. Data and study area

The study area is located in Quebec (Canada), in latitudes between 458 and 608N and longitudes between 658 and 758W (see Fig. 2). Observed daily Tmax, Tmin, and Prec are selected as predictands at four stations: Cedars, Drummondville, Seven Islands, and Bagotville. These stations are located around the St. Lawrence River and the St. Lawrence Gulf (see Fig. 2). The correspond- ing series are provided by Environment Canada weather stations and cover the period between 1 January 1961 and 31 December 2000.

The reanalysis products from the National Centers for Environmental Prediction (NCEP) and the National Center for Atmospheric Research (NCAR) are used to evaluate the potential of the downscaling method. All the NCEP-NCAR data are averaged on a daily basis from 6-h data for on the original regular grid of 2.58 latitude by 2.58 longitude. The obtained predictors are linearly interpolated on the Gaussian grid of 3.75833.758 from the Canadian Centre for Climate Modelling and Analysis (CCCma) Coupled General Circulation Model, version 3 (CGCM3),. These predictors (except the wind direction) are then normalized to the reference period 1961-90. Six grid points covering the study area are se- lected (see Fig. 2) and, for each grid point, 25 NCEP predictors are provided (see Table 1). For each day, 150 variables are available for the downscaling process. The total dataset is divided into two independent sets: a cali- bration period between 1961 and 1990 and a validation period between 1991 and 2000. A PCA is first performed in order to reduce the number of AOGCM predictors and to deal with the multicollinearity problem. The l first components that preserve more than 97% of the variance of the original NCEP predictors are then used in the proposed study as predictor variables.

4. Results

The PGCR model was trained for the calibration period, using Tmax, Tmin, and Prec data from the four stations and the 40 predictors obtained by the PCA. All the coefficients amaxj, aminj, bmaxj, bminj, cj, dj,andej for each site were set following the maximum likelihood estimator by minimizing the NLPD cost function for each predictand. Then, once the parameters of the PDF have been estimated for each day t and for each pre- dictand, all obtained conditional marginal distributions are used to transform the observed data to have a new dataset for the calibration period in the open interval (0, 1). These new datasets are then used to obtain the latent variables z(t) and to fit the parameter of the Gaussian copula. Thereafter, 100 realizations are gen- erated of the precipitation and maximum and minimum temperature series for the validation period (1991- 2000), as shown in Fig. 1.

Figure 3 illustrates an example of the obtained result using the PGCR model at Cedars station during 1991 for both Tmax and Tmin. Results indicate that the estimated series are close to the true observed series for both Tmax and Tmin. Moreover, the majority of the observations are within the 95% confidence intervals based on the estimated standard deviations. This indicates that the proposed model adequately depicts the natural process and its fluctuations. Figure 4 illustrates PGCR results for precipitations at Cedars station during 1991. Figures 4a and 4b show respectively the estimated series of the shape and the scale for conditional gamma distribution. The estimated series of the probability of precipitation occur- rences is shown in Fig. 4c, and the synthetic precipitation series and the observed series are shown in Fig. 4d.Wecan see that the PGCR model provides interesting results for both Poc and Pam.

a. Univariate results

For assessing the downscaling quality, data between 1991 and 2000 are used. Two approaches are considered to validate the PGCR model. The first approach is based on a direct comparison between the estimated and ob- served values using statistical criteria, while the second approach is based on calculating climate indices. In the two validation approaches, the PGCR model results are compared to those obtained using the MMLR and the MMSDM models.

In the first validation approach, three statistical criteria are used for model validation. These criteria are given by

... (22)

... (23)

... (24)

where n denotes the number of observations, yobst refers to the observed value, yestt is the estimated value, t de- notes the day, and s is the standard deviation. The first criterion is the mean error (ME), which is a measure of accuracy. The second criterion is the root-mean-square error (RMSE), which is given by an inverse measure of the accuracy and must be minimized, and the last cri- terion D measures the difference between observed and modeled variances. This criterion evaluates the performance of the model in reproducing the observed variability.

The PGCR and MMSDM models give probabilistic predictions. Point forecasts can be made by estimating the conditional mean for each day. For each predictand, values of the RMSE, ME, and D for all the three models PGCR, MMLR, and MMSDM are given in Table 2. The RMSE and the ME for PGCR and MMSDM were cal- culated using the conditional mean for each day. How- ever, the differences between observed and modeled variances are obtained using the mean D values of 100 realizations. From Table 2 it can be seen that all the three models give similar results in terms of RMSE and ME for maximum and minimum temperatures. However, for downscaled precipitations, PGCR shows the best per- formance, since it has lower RMSE and close to zero ME compared to both the MMSDM and MMLR. Table 3 indicates also that MMSDM performs better than MMLR in terms of both RMSE and ME for downscaled precipitations. This result is because the MMLR model is in reality biased for precipitation variables. The main reason for this bias is that zero precipitation amounts were included to calibrate MMLR amount model. Fur- thermore, the Anscombe residuals R from the observed precipitation amount may not be normally distributed. For this reason, the MMSDM model employs a proba- bility mapping technique to correct this bias. Moreover, from Table 2 it can also be noted that, in terms of D, the PGCR model reproduces better the temporal variability compared to MMLR and MMSDM. In the second vali- dation approach, we consider a set of several climate indices that have been proposed for northern climates for maximum and minimum temperatures. The defini- tions of these climate indices are presented in Table 3. These indices are chosen to evaluate the performance of downscaling models and reflect temperature charac- teristics including the frequency, intensity, and duration of temperature extremes (Wilby 1998; Wilby et al. 2002; Gachon et al. 2005; Hessami et al. 2008). Likewise, for downscaled precipitation several climate indices defined in Table 4 are considered to assess the downscaling quality of precipitation. PGCR and MMLR are then compared by computing the RMSE for each of the cli- matic indices for both precipitation and temperature. For the PGCR model the RMSE of these indices are calcu- lated using the mean RMSE values of 100 realizations.

The results for temperature indices are presented in Table 5. These results indicate that PGCR performs better than both MMLR and MMSDM. MMLR gives better results only for the FSL at Bagotville, the GSL at Drummondville, and the 90th percentile for minimum temperature at Seven Islands. Similarly, Table 6 sum- marizes the RMSE of downscaled precipitations cli- matic indices for both PGCR and MMLR models and for the four weather stations during the validation period (1991-2000). For precipitation amounts five indices are considered: the mean precipitation of wet days (MPWD), the 90th percentile of daily precipitation (P90), the maximum 1-day precipitation (PX1D), the maximum 3-day precipitation (PX3D), and the maximum 5 day precipitations (PX5D). Additionally, three other in- dices are considered for precipitation occurrences: the maximum number of consecutive wet days (WRUN), the maximum number of consecutive dry days (DRUN), and the number of wet days (NWD). In terms of RMSE, re- sults indicate that PGCR gives better results than MMLR and MMSDM for all precipitation indices except for WRUN at Cedars, Pmax90 at Seven Islands, and DRUN at Bagotville, for which the MMSDM gives the best re- sults. This result shows the role of the gamma distribution in the PGCR model to replicate adequately the charac- teristic of precipitation amounts. Furthermore, the use of the logistic regression allows the PGCR model to better reproduce observed monthly characteristics of precipi- tation occurrences based on the indices WRUN, NWD, and DRUN.

b. Interstation and intervariable results

To evaluate the ability of the conditional Gaussian copula in the PGCR model to replicate the observed cross-site correlations for each predictand, the scatter- plots of observed and modeled cross-site correlations of each predictand for PGCR, MMLR, and MMSDM are plotted (Fig. 5). The correlation values of the PGCR and MMSDM model were obtained using the mean of the correlation values calculated from 100 realizations. For all predictands, MMLR overestimates the cross-site cor- relations as shown in Fig. 5, and PGCR and MMSDM reproduce well these cross-site correlations. On the other hand, PGCR outperformed the MMSDM for precipi- tation occurrences. Similarly, Fig. 6 shows the scatterplots of observed and modeled cross-predictand correlations for PGCR, MMLR, and MMSDM during the validation period. This figure shows that both PGCR and MMSDM are able to reproduce more adequately these cross- predictand correlations. This is a great achievement of the two multivariable models, in comparison with the univariate MMLR model. Basically, both PGCR and MMSDM precisely simulate these cross-predictand correlations and there is no clear difference between them except when precipitation amount is present. Indeed, MMSDM has difficulty reproducing the cross- correlation predictand. Indeed, these cross-correlation values can be affected by the probability mapping step that is used to correct the bias and to reproduce the adequate distribution of precipitation. However, when evaluating the PGCR model there is no need to rely on transformation steps or on bias correction pro- cedures and the mapping in the conditional distribu- tion is automatic using its probabilistic regression component.

For precipitation, joint probabilities of the events that two sites are both dry or both wet on a given day are displayed in Fig. 7. PGCR and MMSDM adequately simulate these joint probabilities and there is no clear difference between PGCR and MMSDM simulations, both having almost better results compared to the MMLR model. This is a great finding of PGCR and MMSDM, in comparison with the joint probabilities from the single-site MMLR model.

Finally, to evaluate the consistency of local weather variables, Fig. 8 compares differences of mean temper- atures on wet and dry days (mean temperature on wet days minus the corresponding value on dry days) in the synthetic datasets with the observed ones. The values in the plots are for each site and each month (48 data points in total). It appears that the self-consistency between local precipitation and temperatures is reproduced very well in synthetic datasets, from either PGCR or MMLR.

5. Discussion and conclusions

A PGCR model is proposed in this paper for the downscaling of AOGCM predictors to multiple predic- tands at multiple sites simultaneously and for preserving relationships between sites and variables. This model relies on a probabilistic framework in order to describe the conditional density of each predictand for a given day. The PGCR model uses a Gaussian distribution for max- imum and minimum temperature, a Bernoulli distribution for precipitation occurrences, and a gamma distribution for precipitation amounts. In the probabilistic framework, PGCR adopts a regression model with outputs for each parameter in the specified probability density function. To maintain realistic relationships between sites and vari- ables, the PGCR model uses a Gaussian copula that de- scribes dependences between all predictands.

The developed model was then applied to generate daily maximum and minimum temperatures and preci- pitations of four observation sites located in the south- ern part of the province of Quebec (Canada). NCEP reanalysis data were used as predictors in order to assess the potential of the method, although the final objec- tive is to use AOGCM predictors. Application results of the PGCR model were compared with those of the regression-based MMLR model and an MMSDM as a classical multisite and multivariable model.

Results show that all the three models give similar results in terms of RMSE and ME for maximum and minimum temperatures. On the other hand, PGCR per- forms better for downscaled precipitations in terms of ME and RMSE. In addition, the comparison based on temperature and precipitation indices shows that the PGCR model is more able to reproduce extremes and observed variability on a seasonal and interannual ba- sis for both temperature and precipitation. In terms of reproducing spatial and intervariable properties, both PGCR and MMSDM models provide interesting re- sults without significant differences.

Reproduction of the temporal variability in both precipitation and temperature fields is among the most important achievements for the proposed PGCR model. The MMLR model showed difficulty in reproducing the observed variability. Indeed, regression models gener- ally reproduce the mean of the process conditionally to the selected independent variables. As a consequence, the variability of the regression is always smaller than the initial variability. In this regard, von Storch (1999) men- tioned that predictors generated from synoptic-scale fields cannot represent all variability at the subgrid scale. Hence, in this way, the PGCR model presents the ad- vantage of modeling the entire conditional distribution, and thereby gives the mean of the downscaled pre- dictands as well as their variability. The reproduction of the variability is not a typical characteristic of PGCR; indeed, stochastic and hybrid models, such as the MMSDM model, have this characteristic. Nevertheless, these hybrid approaches are based on static noises that are not dependent on the predictors. Thus, hybrid ap- proaches may not represent local change in atmospheric variability in a climate change simulation. Thereby, one advantage of the proposed PGCR model compared to hybrid approaches is that the temporal variability could change in the future, which may explain the fact that PGCR gives better results compared to the MMSDM model, as shown in Table 2. This finding shows the role of the probabilistic regression component model of PGCR that allows predicting not only the conditional mean but also the whole conditional distribution including vari- ability. Therefore, the real gain when using the proposed PGCR model consists not only in the dependence mod- eling through Gaussian copula but also including the advantages of probabilistic regression in a multisite multivariable framework. This is attributable to the characteristic of the marginal effect elimination through the copula. This attractive characteristic helps to model and understand the dependence structure effectively, as it has no relationship with the marginal behavior. However, it should be mentioned that it is more difficult to re- produce correlations when the precipitation amount is considered, since precipitation amount results are af- fected by the precipitation occurrence results. For this reason it would be more judicious to model the occur- rence and amount of precipitation simultaneously by employing only one PDF that describes the discrete- continuous behavior of precipitation, for example by us- ing a Bernoulli-gamma or Poisson-gamma PDF.

An important fact that has not been considered in this work is that the time structure of the downscaled pre- cipitations must be poorly simulated, such as lag-1 cor- relation for downscaled precipitations. The number of AOGCM grid points in this study can also be increased, which would improve the precision of both MMLR and PGCR. In that context, the use of regional-scale predictors from regional climate models (RCMs) instead of coarse-scale AOGCMs will be strongly beneficial for improving downscaled results. Finally, the results of the PGCR model can be improved by developing the PGCR model for each month. This would allow taking into ac- count the seasonal variability of each predictand.

REFERENCES

Bardossy, A., 2006: Copula-based geostatistical models for groundwater quality parameters. Water Resour. Res., 42, W11416, doi:10.1029/2005WR004754.

Bates, B. C., S. P. Charles, and J. P. Hughes, 1998: Stochastic down- scaling of numerical climate model simulations. Environ. Modell. Software, 13, 325-331, doi:10.1016/S1364-8152(98)00037-1.

Bellone, E., J. P. Hughes, and P. Guttorp, 2000: A hidden Markov model for downscaling synoptic atmospheric patterns to pre- cipitationamounts. Climate Res.,15,1-12,doi:10.3354/cr015001.

Burger, G., and Y. Chen, 2005: Regression-based downscaling of spatial variability for hydrologic applications. J. Hydrol., 311, 299-317, doi:10.1016/j.jhydrol.2005.01.025.

Cannon, A. J., 2008: Probabilistic multisite precipitation down- scaling by an expanded Bernoulli-gamma density network. J. Hydrometeor., 9, 1284-1300, doi:10.1175/2008JHM960.1.

_____, 2009: Negative ridge regression parameters for improving the covariance structure of multivariate linear downscaling models. Int. J. Climatol., 29, 761-769, doi:10.1002/joc.1737.

_____, 2011: Quantile regression neural networks: Implementation in R and application to precipitation downscaling. Comput. Geosci., 37, 1277-1284, doi:10.1016/j.cageo.2010.07.005.

Cawley,G. C.,G.J. Janacek,M.R. Haylock, andS.R.Dorling,2007: Predictive uncertainty in environmental modelling. Neural Networks, 20, 537-549, doi:10.1016/j.neunet.2007.04.024.

Chebana, F., and T. B. M. J. Ouarda, 2007: Multivariate L-moment homogeneity test. Water Resour. Res., 43, W08406, doi:10.1029/ 2006WR005639.

Clark, M., S. Gangopadhyay, L. Hay, B. Rajagopalan, and R. Wilby, 2004: The Schaake shuffle: A method for recon- structing space-time variability in forecasted precipitation and temperature fields. J. Hydrometeor., 5, 243-262, doi:10.1175/ 1525-7541(2004)005,0243:TSSAMF.2.0.CO;2.

Conway, D., R. Wilby, and P. Jones, 1996: Precipitation and air flow indices over the British Isles. Climate Res., 7, 169-183, doi:10.3354/cr007169.

De Michele, C., and G. Salvadori, 2003: A generalized Pareto intensity-duration model of storm rainfall exploiting 2-copulas. J. Geophys. Res., 108, 4067, doi:10.1029/2002JD002534.

Dorling, S. R., R. J. Foxall, D. P. Mandic, and G. C. Cawley, 2003: Maximum likelihood cost functions for neuralnetwork models of air quality data. Atmos. Environ., 37, 3435-3443, doi:10.1016/ S1352-2310(03)00323-6.

Dupuis, D. J., 2007: Using copulas in hydrology: Benefits, cau- tions, and issues. J. Hydrol. Eng., 12, 381-393, doi:10.1061/ (ASCE)1084-0699(2007)12:4(381).

El Adlouni, S., and T. B. M. J. Ouarda, 2008: E^ tude de la loi conjointe de^bit-niveau par les copules: Cas de la rivi^ere Ch^ateauguay (Study of the joint law flow-level by copulas: Case of the Chateauguay River). Can. J. Civil Eng., 35, 1128- 1137, doi:10.1139/L08-054.

Fasbender, D., and T. B. M. J. Ouarda, 2010: Spatial Bayesian model for statistical downscaling of AOGCM to minimum and maximum daily temperatures. J. Climate, 23, 5222-5242, doi:10.1175/2010JCLI3415.1.

Gachon, P., and Coauthors, 2005: A first evaluation of the strength and weaknesses of statistical downscaling methods for sim- ulating extremes over various regions of eastern Canada. Sub-component, Climate Change Action Fund (CCAF), Environment Canada Final Rep., 209 pp. [Available online at http://www.cccsn.ec.gc.ca/?page5r-and-p.]

Hammami, D., T. S. Lee, T. B. M. J. Ouarda, and J. Le, 2012: Predictor selection for downscaling GCM data with LASSO. J. Geophys. Res., 117, D17116, doi:10.1029/2012JD017864.

Haylock, M. R., G. C. Cawley, C. Harpham, R. L. Wilby, and C. M. Goodess, 2006: Downscaling heavy precipitation over the United Kingdom: A comparison of dynamical and statistical methods and their future scenarios. Int. J. Climatol., 26, 1397- 1415, doi:10.1002/joc.1318.

Hellstroeurom, C., D. Chen, C. Achberger, and J. Reuroaiseuroanen, 2001: Comparison of climate change scenarios for Sweden based on statistical and dynamical downscaling of monthly precipita- tion. Climate Res., 19, 45-55, doi:10.3354/cr019045.

Herrera, E., T. B. M. J. Ouarda and B. Bob^ee 2006: Me^thodes de d^esagre^gation applique^es aux Mod^eles du Climat Global Atmosph^ere-Oc^ean (MCGAO) [Downscaling methods applied to atmosphere-ocean general circulation models (AOGCM)]. J. Water. Sci., 19, 297-312, doi:10.7202/014417ar.

Hessami, M., P. Gachon, T. B. M. J. Ouarda, and A. St-Hilaire, 2008: Automated regression-based statistical downscaling tool. Environ. Modell. Software, 23, 813-834, doi:10.1016/ j.envsoft.2007.10.004.

Hughes, J. P., P. Guttorp, and S. P. Charles, 1999: A non- homogeneous hidden Markov model for precipitation occur- rence.J.Roy.Stat.Soc.,48C,15-30,doi:10.1111/1467-9876.00136.

Huth, R., 1999: Statistical downscaling in central Europe: Evalua- tion of methods and potential predictors. Climate Res., 13, 91- 101, doi:10.3354/cr013091.

_____, 2002: Statistical downscaling of daily temperature in central Europe. J. Climate, 15, 1731-1742, doi:10.1175/ 1520-0442(2002)015,1731:SDODTI.2.0.CO;2.

_____, 2004: Sensitivity of local daily temperature change estimates to the selection of downscaling models and predictors. J. Cli- mate, 17, 640-652, doi:10.1175/1520-0442(2004)017,0640: SOLDTC.2.0.CO;2.

Jeong, D. I., A. St-Hilaire, T. B. M. J. Ouarda, and P. Gachon, 2012a: CGCM3 predictors used for daily temperature and precipitation downscaling in southern Qu^ebec, Canada. Theor. Appl. Climatol., 107, 389-406, doi:10.1007/s00704-011-0490-0.

_____, _____, _____, and _____, 2012b: Multisite statistical downscaling model for daily precipitation combined by multivariate mul- tiple linear regression and stochastic weather generator. Cli- matic Change, 114, 567-591, doi:10.1007/s10584-012-0451-3.

_____, _____, _____, and _____, 2012c: A multivariate multi-site sta- tistical downscaling model for daily maximum and minimum temperatures. Climate Res., 54, 129-148, doi:10.3354/cr01106.

Kazianka, H., and J. Pilz, 2010: Copula-based geostatistical mod- eling of continuous and discrete data including covariates. Stochastic Environ. Res. Risk Assess., 24, 661-673, doi:10.1007/ s00477-009-0353-8.

Khalili, M., V. T. Van Nguyen, and P. Gachon, 2013: A statistical approach to multisite multivariate downscaling of daily extreme temperature series. Int. J. Climatol., 33, 15-32, doi:10.1002/ joc.3402.

Lagarias, J. C., J. A. Reeds, M. H. Wright, and P. E. Wright, 1998: Convergence properties of the Nelder-Mead simplex method in low dimensions. SIAM J. Optim., 9, 112-147, doi:10.1137/ S1052623496303470.

Laux, P., S. Vogl, W. Qiu, H. Knoche, and H. Kunstmann, 2011: Copula-based statistical refinement of precipitation in RCM simulations over complex terrain. Hydrol. Earth Syst. Sci., 15, 2401-2419, doi:10.5194/hess-15-2401-2011.

Lindstroeurom, G., B. Johansson, M. Persson, M. Gardelin, and S. Bergstroeurom, 1997: Development and test of the distributed HBV-96 hydrological model. J. Hydrol., 201, 272-288, doi:10.1016/S0022-1694(97)00041-3.

Palutikof, J. P., C. M. Goodess, S. J. Watkins, and T. Holt, 2002: Generating rainfall and temperature scenarios at multiple sites: Examples from the Mediterranean. J. Climate, 15, 3529-3548, doi:10.1175/1520-0442(2002)015,3529:GRATSA.2.0.CO;2.

Pitt, M., D. Chan, and R. Kohn, 2006: Efficient Bayesian inference for Gaussian copula regression models. Biometrika, 93, 537- 554, doi:10.1093/biomet/93.3.537.

Qian, B., J. Corte-Real, and H. Xu, 2002: Multisite stochastic weather models for impact studies. Int. J. Climatol., 22, 1377- 1397, doi:10.1002/joc.808.

Renard, B., and M. Lang, 2007: Use of a Gaussian copula for multivariate extreme value analysis: Some case studies in hydrology. Adv. Water Resour., 30, 897-912, doi:10.1016/ j.advwatres.2006.08.001.

Schoeurolzel, C., and P. Friederichs, 2008: Multivariate non-normally distributed random variables in climate research-Introduction to the copula approach. Nonlinear Processes Geophys., 15, 761- 772, doi:10.5194/npg-15-761-2008.

Schoof, J. T., and S. C. Pryor, 2001: Downscaling temperature and precipitation: A comparison of regression-based methods and artificial neural networks. Int. J. Climatol., 21, 773-790, doi:10.1002/joc.655.

Sklar, A., 1959: Fonctions de r^epartition ^a n dimensions et leurs marges. Publ. Inst. Stat. Univ. Paris, 8, 229-231.

Stephenson, D. B., K. Rupa Kumar, F. J. Doblas-Reyes, J. F. Royer, F. Chauvin, and S. Pezzulli, 1999: Extreme daily rain- fall events and their impact on ensemble forecasts of the In- dian monsoon. Mon. Wea. Rev., 127, 1954-1966, doi:10.1175/ 1520-0493(1999)127,1954:EDREAT.2.0.CO;2.

Vogl, S., P. Laux, W. Qiu, G. Mao, and H. Kunstmann, 2012: Copula-based assimilation of radar and gauge information to derive bias-corrected precipitation fields. Hydrol. Earth Syst. Sci., 16, 2311-2328, doi:10.5194/hess-16-2311-2012.

von Storch, H., 1999: On the use of ''inflation'' in statistical downscaling. J. Climate, 12, 3505-3506, doi:10.1175/ 1520-0442(1999)012,3505:OTUOII.2.0.CO;2.

Vrac, M., and P. Naveau, 2007: Stochastic downscaling of pre- cipitation: From dry events to heavy rainfalls. Water Resour. Res., 43, W07402, doi:10.1029/2006WR005308.

Widmann, M., C. S. Bretherton, and E. P. Salathe^ Jr., 2003: Sta- tistical precipitation downscaling over the northwestern United States using numerically simulated precipitation as a predictor. J. Climate, 16, 799-816, doi:10.1175/1520-0442(2003)016,0799: SPDOTN.2.0.CO;2.

Wilby, R. L., 1998: Statistical downscaling of daily precipitation using daily airflow and seasonal teleconnection indices. Cli- mate Res., 10, 163-178.

_____, and T. M. L. Wigley, 1997: Downscaling general circulation model output: A review of methods and limitations. Prog. Phys. Geogr., 21, 530-548, doi:10.1177/030913339702100403.

_____, H. Hassan, and K. Hanaki, 1998: Statistical downscaling of hydrometeorological variables using general circulation model output. J.Hydrol.,205,1-19,doi:10.1016/S0022-1694(97)00130-3.

_____, C. W. Dawson, and E. M. Barrow, 2002: SDSM-A decision support tool for the assessment of regional climate change impacts. Environ. Modell. Software, 17, 145-157, doi:10.1016/ S1364-8152(01)00060-3.

Wilks, D. S., 1998: Multisite generalization of a daily stochastic precipitation generation model. J. Hydrol., 210, 178-191, doi:10.1016/S0022-1694(98)00186-3.

_____, 1999: Multisite downscaling of daily precipitation with a sto- chastic weather generator. Climate Res., 11, 125-136, doi:10.3354/ cr011125.

_____, and R. L. Wilby, 1999: The weather generation game: A re- view of stochastic weather models. Prog. Phys. Geogr., 23, 329-357, doi:10.1177/030913339902300302.

Williams, P. M., 1998: Modelling seasonality and trends in daily rainfall data. Proc. Conf. on Advances in Neural Information Processing Systems, Boulder, CO, Neural Information Pro- cessing Systems Foundation, 985-991.

Yang,C.,R.E.Chandler,V.S.Isham,andH.S.Wheater,2005: Spatial-temporal rainfall simulation using generalized linear models. Water Resour. Res., 41, W11415, doi:10.1029/ 2004WR003739.

M. A. BEN ALAYA AND F. CHEBANA

Centre Eau Terre Environnement, Institut National de la Recherche Scientifique, Quebec, Quebec, Canada

T. B. M. J. OUARDA

Institute Center for Water and Environment (iWater), Masdar Institute of Science and Technology, Abu Dhabi, United Arab

Emirates, and Centre Eau Terre Environnement, Institut National de la Recherche Scientifique, Quebec, Quebec, Canada

(Manuscript received 31 May 2013, in final form 9 December 2013)

Corresponding author address: M. A. Ben Alaya, INRS-ETE, 490 rue de la Couronne, Quebec QC G1K 9A9, Canada.

E-mail: [email protected]

DOI: 10.1175/JCLI-D-13-00333.1

Copyright:	(c) 2014 American Meteorological Society
Wordcount:	8326