Space-time calibration of wind speed forecasts from regional climate models
Luiz E. S. Gomes Thaís C. O. Fonseca Kelly C. M. Gonçalves Ramiro Ruiz-Cárdenas


Relay-assisted Device-to-Device Networks: Connectivity and Uberization Opportunities



    Quentin Le Gall12, Bartłomiej Błaszczyszyn2, Elie Cali1 and Taoufik En-Najjary1
1Modelling and Statistical Analysis, Orange Labs Networks, Châtillon, France

Email: quentin1.legall@orange.com, elie.cali@orange.com and taoufik.ennajjary@orange.com
2Inria-ENS, Paris, France
Email: quentin.le.gall@ens.fr, bartek.blaszczyszyn@ens.fr
    Received 27 September 2019; accepted 25 May 2020; published 1 July 2020
==============================================================================================================================================================================================================================================================================================================================================




Numerical weather predictions (NWPs) are systematically subject to errors due to the deterministic solutions used by numerical models to simulate the atmosphere. Statistical postprocessing techniques are widely used nowadays for NWP calibration. However, time-varying bias is usually not accommodated by such models. The calibration performance is also sensitive to the temporal window used for training. This paper proposes space-time models that extend the main statistical postprocessing approaches to calibrate NWP model outputs. Trans-Gaussian random fields are considered to account for meteorological variables with asymmetric behavior. Data augmentation is used to account for censoring of the response variable. The benefits of the proposed extensions are illustrated through the calibration of hourly 10-meter height wind speed forecasts in Southeastern Brazil coming from the Eta model.







§ INTRODUCTION



Numerical weather predictions (NWPs) are often based on mathematical models that make deterministic predictions from current atmospheric conditions. Such models rely on complex systems of differential equations that do not have an analytical solution. Therefore, numerical integration is used to simulate the physical, dynamic, and thermodynamic processes of the atmosphere depending on its current state. From this initial solution, it is possible to solve the system for any future time of interest <cit.>.

These numerical systems are solved in a discrete grid which means that the model assumes uniform predictions for any point in the region belonging to a given cell of the grid. In particular, forecasts at each grid cell are based on average data in its region (e.g., average elevation and predominant vegetation). Under these assumptions, the representativeness of predictions in cells with, for example, complex orography, dense vegetation, or presence of water bodies, is inadequate due to differences between the actual characteristics of the surface and the homogenization made by the model. Therefore, outputs from NWP models may not be representative at specific locations <cit.>, thus producing systematic errors.

Meteorological phenomena are often not well described by a single numerical prediction. As an alternative, ensembles of forecasts, i.e., groups of outputs coming from multiple runs of either, the same or several models under varying initial conditions and model physics, are considered to allow quantification of uncertainty. The ensembles can be interpreted as a Monte Carlo experiment aiming to produce a range of future states of the atmosphere from different initial conditions <cit.>. Moreover, the simulation of scenarios can indicate extreme events that would not be identified by just one run of the numerical model <cit.>. However, according to <cit.>, ensemble forecasts are often under-dispersive, i.e., the ensemble spread is too narrow to account for all the uncertainty. As a solution, calibration is often used to correct these error patterns.

Statistical postprocessing techniques have been widely used in recent years to minimize the above limitations of numerical models and to enhance both, the reliability and statistical consistency of numerical weather predictions. Many methods for the statistical calibration of ensemble forecasts are now available. Recent reviews on this subject can be found in <cit.> and <cit.>. The pioneer postprocessing technique in this context was the application of model output statistics <cit.> which considers a multiple linear regression relating the responses (e.g., observed wind) to the set of ensemble members, assuming constant variance. Several extensions of this precursor method have been proposed, such as updatable MOS <cit.> which varies the size of the training period; MOC <cit.> which directly models the prediction error; and the generalized MOS <cit.> which considers generalized linear models instead of Gaussian distributions for the response vector. An extension of the MOS technique known as ensemble MOS <cit.> allows for a spread-skill relationship between the dispersion of ensemble members and the response variance <cit.>. In the context of spatial calibration, the geostatistical output perturbation <cit.> allows spatial dependence, resulting in calibrated weather fields for a fixed temporal horizon. Moreover, the spatial EMOS method <cit.> combines the EMOS and GOP techniques.

Estimation of parameters in statistical postprocessing usually occurs in a subset of the data, called a training period, defined through a moving time window that accounts for the effect of both past observations and numerical predictions. If the training period is reasonably long, it becomes easier to estimate the uncertainty of predictions <cit.>. However, longer training periods may introduce distortions in the calibration due to seasonal effects. This challenging tradeoff suggests that the window size must be tailored for the specific application to achieve good results. To illustrate, <cit.> analyzed the effects of the training window size on the uncertainty estimation of parameters in the calibration of temperature and sea level pressure predictions coming from the same numerical model. The authors reported gains with the use of a time window as large as 25 days for both case studies and highlighted the need for an automatic way to choose the length of the training period.

The above shortcomings of current statistical postprocessing approaches can be addressed by the incorporation of statistical techniques such as spatially structured calibration models. This approach may be useful in correcting NWP errors in regions where numerical models have smoothed important characteristics of the terrain, while the Kalman filter <cit.> and Bayesian dynamic models <cit.> are natural alternatives to account for seasonality and temporal dynamics of bias in postprocessing models.

The proposed models should also be able to accommodate the spatial characteristics of the terrain, as well as fairly skew and censor the behavior of the response, which can be essential for prediction and uncertainty quantification of certain variables, such as wind speed. In this context, the transformed Gaussian dynamic model has been useful in applied settings, such as precipitation modeling <cit.>.  <cit.> and <cit.> also followed the latent Gaussian approach of <cit.> to predict weather variables. In particular, <cit.> proposed a spatiotemporal model for rainfall data based on the convolution of covariance functions. Note that in this case, the resulting covariance structure is usually not obtained in closed form. <cit.> proposed a stochastic partial differential equation-based model and applied it to postprocessing of precipitation forecasts. The model is based on the spectral domain. For further details about the spectral approach to spatiotemporal modeling, see <cit.>.

In general, the proposed model is a known spatial dynamic regression model, which we have extended to include Box-Cox transformation to adapt to the wind speed characteristics such as asymmetry and censoring. Developing asymmetric models for spatiotemporal applications is non-trivial, as discussed in recent papers. <cit.> mentioned that the parameters in asymmetric spatial processes are not well identified, even if the number of locations is large. <cit.> indicated that for skew-normal spatial models when the asymmetry is high, the spatial correlation between two locations approaches 1 regardless of the distance between these two locations. In this context, the proposed model is an interesting alternative that accounts for non-stationarity in time, asymmetry, censoring, and spatiotemporal correlation. The Gaussian model and conditional Gaussian extensions are usual choices for the spatiotemporal setup due to the validity of the process and analytical results available for fast predictions. In the context of spatial generalized linear models <cit.>, the model would be hierarchical, and the spatiotemporal random effects would be assumed in the second level of this hierarchy. We did not follow this approach, as it would be computationally even more demanding than the conditional Gaussian alternative with Box-Cox transformation.

This paper proposes a unifying approach to statistical postprocessing by accounting for both spatial and temporal dependencies overcoming the need for large time windows in the definition of training sets. Our proposal is inspired by current postprocessing approaches, which are extended based on both, spatial dynamic linear models and data augmentation techniques. An optimized Markov chain Monte Carlo (MCMC) scheme based on the robust adaptive Metropolis algorithm <cit.> is used to perform efficient statistical inference and prediction.

The primary motivation of this study is the calibration of 10-meter height wind speed forecasts since that meteorological variable presents some particular features that might not be properly addressed by the usual postprocessing techniques. According to <cit.>, some of these features include: intermittent atmospheric regimes with a predominance of a certain wind direction in given regions; spatial and temporal correlation; non-Gaussianity; non-stationarity; conditional heteroscedasticity, i.e., the variance of the wind speed changes frequently in time; seasonal annual and diurnal components due to the effects of the sun and seasons; and possible trends.

The remainder of this article is organized as follows. Section <ref> presents the wind speed data which motivates our proposed modeling approach. Section <ref> briefly describes the main postprocessing approaches used to calibrate numerical weather forecasts; Section <ref> presents the proposed extensions to enhance current postprocessing techniques and describes the inference procedure. The performance of the proposed extensions for applications using real data and simulated experiments is evaluated and compared with the current approaches in Sections <ref> and <ref>, respectively. The article concludes with a discussion in Section <ref>.




§ WIND SPEED DATA FOR THE STATE OF MINAS GERAIS



The main motivation for the proposed models presented in this paper is the calibration of wind speed forecasts generated by the Eta regional climate model for the state of Minas Gerais, Brazil. The state, located in the Southeast region of Brazil, is formed by abundant plateaus interspersed with rugged relief, ranging from 100 to 2800 meters above sea level, which provides to the state exceptional water resources. The predominant vegetation, known as Cerrado(savanna), is characterized by large variations in the landscape between the rainy and dry seasons. This results in a seasonal influence of the surface roughness on the displacement of the winds, which are stronger during winter and spring. The climate in Minas Gerais varies from hot semiarid to humid mesothermal. Rainfall distribution is uneven, with the northern region having long periods of drought and the highest temperatures, while the southern region (high elevation areas) concentrates the highest total annual rainfall <cit.>.

Two sources of hourly 10-meter height wind speed data for the region of interest, covering the two years from 1 October 2015 to 30 September 2017, were used in this application: ground measurements and numerical forecasts. The measurement data are recorded by an irregular network of weather stations operated by the Brazilian National Institute of Meteorology (INMET), available at < <http://www.inmet.gov.br> >. A set of 68 stations from this network, spread over Minas Gerais and its surroundings were selected with the requirement of at least 70% of data availability. In these data, distances in space range from 7.9 to 1145 km, and elevation of stations varies from 60 to 1359 m. The numerical accuracy of the data is to one decimal place.

Numerical forecasts are generated by the Eta model <cit.>. It is a regional climate model run at the Center for Weather Forecasting and Climate Studies (CPTEC). Its extent covers most of South America and Central America. The model results are disclosed twice a day at 0000 UTC and 1200 UTC. We used hourly Eta model outputs from the run began at 1200 UTC with a horizontal resolution of 15 km, 50 vertical layers, and a lead time of up to 264 hours (11 days). This lead time allows obtaining up to 10 different wind speed forecasts for each hour, coming from the daily runs of the model. These groups of forecasts composed our ensembles in this application. This was because just one version of the model is issued at each running time. Hence, ensembles of forecasts for several numerical models were not available. Figure <ref> presents the distribution of selected stations in the region of interest as well as the discrete grid of the numerical model.



Gridded forecasts were bilinearly interpolated to the locations of the 68 stations, to form the dataset with hourly 10-meter height wind speed forecasts for calibration at these locations. According to <cit.>, more complex interpolation methods can be used but these are unlikely to produce considerable gains if the grid is very fine.

To illustrate the local characteristics of wind speed, Figure <ref> presents some histograms with the wind speed distribution during the seasons at three selected stations: Viçosa, Muriaé, and Pampulha (Belo Horizonte). The histograms clearly show an asymmetric distribution with a frequent point mass at zero. 
In general, lower wind speeds are recorded during the fall and summer and this pattern extends to the vast majority of available weather station readings in the dataset. Particularly, the wind speed recorded during the spring has the highest threshold and the smallest point mass at zero. This wind regime is even more evident at the Pampulha (Belo Horizonte) station, as illustrated in Figure <ref>(c). These specific aspects of wind speed distribution motivate the space-time postprocessing models proposed in this paper, as well as the need to consider the left censoring and data transformation.










§ STATISTICAL POSTPROCESSING MODELS



We now describe the main postprocessing models used in the probabilistic calibration of forecasts generated by NWP models. Our proposed methods based on dynamic spatial models are presented as an automatic bias correction alternative for ensemble forecasting in longer training periods.



 §.§ Static calibration models


The model output statistics <cit.> defines a relation between
an ensemble with m forecasts F_1,...,F_m, individually distinguishable, for an univariate quantity of interest Y through a multiple linear regression model given by:

    Y = θ_0 + θ_1F_1 + ... + θ_mF_m + ε,

where the error term ε assures that E(ε)=0 and  Var(ε) = σ^2. 

An extension of MOS called ensemble MOS <cit.>, also known as non-homogeneous regression, allows the response variance to depend on the dispersion of ensemble members. This improvement is called the spread-skill relationship. It is based on the premise of a positive relationship between ensemble spread and forecast absolute error.
The EMOS method assumes that ε is such that:

    E(ε) = 0,    Var(ε) = σ^2_* = β_0 + β_1S^2,

where S^2 is the sample variance of ensemble members and = (β_0,β_1)' are non-negative coefficients. <cit.> assumed Gaussianity and estimated model parameters by minimization of continuous ranked probability score <cit.>. Other applications of this method in the context of wind speed, strong winds and wind direction can be found in  <cit.>, <cit.>, and <cit.>, respectively.

A spatial extension of MOS called geostatistical output perturbation <cit.> is the pioneer statistical postprocessing model that considers the correlation between the measurements of a meteorological variable at distinct locations. This technique produces calibrated forecasts for entire weather fields. Let { Y(), ∈ S ⊂ℝ^2 } be a random weather field and Y = (y(_1),...,y(_n))', an observed sample at n locations. Considering a set of m ensemble members for these locations, represented by F_1 = (F_1(_1),...,F_1(_n))', … , F_m = (F_m(_1),...,F_m(_n))', and assuming = ( ε(_1),...,ε(_n))', a vector of observations from a Gaussian process {ε(), ∈ S } is obtained, such that:

    E() = 0_n,    Cov(ε(_i),ε(_j)) = Σ_i,j = σ^2C(_i,_j),    i,j = 1,…,n,

where 0_n is a zero n-vetor and C(·,·) is a valid spatial correlation function <cit.> for any pair of locations in S. This method is a multivariate extension of the model described in (<ref>) and therefore also does not allocate the ensemble information in its structure.

Finally, considering simultaneously the ensemble information and the spatial structure of meteorological variables, the spatial EMOS <cit.> is presented combining the EMOS and GOP models. It consists of the GOP model formulation but assumes that:

    Cov(ε(_i),ε(_j)) = Σ^*_i,j = D_i,iC(s_i,s_j)D_j,j,    i,j = 1,...,n,

where D = diag(√(β_0 + β_1 S^2_1),...,√(β_0 + β_1 S^2_n)) is an n× n diagonal matrix with S^2_i representing the sample variance of ensemble members for location i. The models presented here, in their canonical form, are not appropriate for calibration of asymmetric variables (e.g., wind speed and rainfall).



 §.§ Proposed dynamic calibration models



The postprocessing models proposed here improve previous approaches in several ways, accounting for both, temporal dynamics and asymmetric distributions with mass at zero through a spatiotemporal modeling approach with censoring. These novel methods are presented next. In particular, the truncated Gaussian model which has been used in precipitation prediction <cit.> is presented in the context of forecast calibration from NWP models.



  §.§.§ Spatiotemporal trans-Gaussian model


Let { Y_t(), ∈ S ⊂ℝ^2, t=1,...,T } be a random spatial field in discrete time t. The observed response vector at n locations Y_t = (Y_t(_1),...,Y_t(_n))' is composed by censored random variables, Y_t(_i) ≥ c, i=1,…,n, t=1,…,T. Assume that Y_t follows a truncated Gaussian distribution as follows:


    Y_t()={[    BC^-1( X_t(), λ) , BC^-1( X_t(), λ) > c,;                    c, BC^-1( X_t(), λ) ≤ c, ].

where c is a known constant, λ is the unknown transformation parameter, X_t() is a Gaussian process and BC(·,λ) represents the family of Box-Cox transformations <cit.> defined as:

    BC(y,λ)={[    ( y^λ-1 ) /λ ,  λ≠ 0      y > 0,;            log y, λ = 0      y > 0.;                   ].


Therefore, X_t() is a latent Gaussian field which allows for asymmetry in the resulting process of interest Y_t(). That is, after transforming the possibly asymmetric process Y_t(), the resulting field will follow a Gaussian process. The trans-Gaussian model considered adds flexibility to the Gaussian model usually assumed for random spatiotemporal fields. This strategy was considered due to its simple specification as opposed to other alternatives to develop non-Gaussian processes which may result in invalid processes <cit.>. The proposed trans-Gaussian field accounts for asymmetry (e.g., the log transformation) which is a known characteristic of wind speed data modeling.

Furthermore, censoring is admitted through a known constant c. In the context of wind speed calibration, c represents the minimum allowed value for observed wind speed (usually 0 m/s). However, depending on the application context, this value can be greater than zero. Figure <ref> illustrates the practical effect of this modeling setup on simulated wind speed values at a single location. Note that values below c from transformed and truncated Gaussian latent process induce a mass at c in the process of interest Y_t(). Application of the complete model with Box-Cox transformation and censoring is facilitated by the use of data augmentation techniques <cit.>. It also naturally allows for missing data. 

Spatiotemporal trans-Gaussian models with different transformation functions are widely used, especially in the context of precipitation forecasting. Specific transformation functions are defined for each application. In general, the power transform function <cit.> and a quadratic form of the power transform function <cit.> are the most used due to their simplicity. However, other more complex functions can be found, such as the power-exponential transformation <cit.>. For the wind speed prediction context, we selected the Box-Cox transformation, a power transform function extension, often used in temporal and spatial applications as discussed, respectively, in <cit.> and <cit.>.



In the usual setup for statistical calibration of numerical models time windows are defined to account for observed and predicted temporal variation, the so-called training period. Although larger training windows result in smaller uncertainty, they introduce distortions due to seasonal effects. In general, the seasonal patterns of meteorological variables are well defined (e.g., solar forcing – 24 hours and seasons – 3 months). Thus, the Kalman filter <cit.> allows for both longer temporal training periods and the inclusion of temporal dynamics in the bias parameters. In the following, we consider the Bayesian dynamic approach which considers the Kalman filter for forward filtering and backward posterior smoothing for a fully Bayesian inference procedure as proposed by <cit.>. This approach allows for time-varying regression coefficients while accounting for spatial dependence in the region of interest, overcoming the over-smoothing of numerical forecasts.



  §.§.§ Dynamic geostatistical output perturbation


The proposed model adds a temporal dynamic to the GOP formulation by allowing the coefficients to vary over time. The covariance modeling is stochastic over time with Beta-Gamma evolution depending on discount factors. Assuming the structure stated in (<ref>) for the transformation of the original process Y_t(s), the dynamics are introduced by defining the observation and state equations as follows:


    X_t= F'_t_t + _t,      _t ∼ N(0_n,Σ_t),   Σ_t = φ_t^-1H, 
    _t= G_t_t-1 + _t,      _t ∼ T_n_t-1(0_p,W_t),
    φ_t = γ_t φ_t-1/δ_V,      γ_t ∼Beta(κ_t,κ̅_t),

where _t = ( ε_t(_1),...,ε_t(_n))' follows a zero mean multivariate normal distribution with covariance matrix Σ_t = φ_t^-1H such that φ_t=1/σ^2_t is the precision and H is the correlation matrix with elements H_i,j = C(_i,_j),   i,j = 1,…,n. C(·,·;ϕ) is a valid correlation function depending on an unknown parameter ϕ. For instance, a simple correlation function is C(_i,_j;ϕ) = exp(-ϕ_i - _j), with ϕ > 0 representing the spatial
decay parameter and _i - _j, the Euclidean distance between locations _i and _j, i,j=1,...,n. Other possible choices are the Matérn class <cit.> and the Cauchy class <cit.>. X_t = ( X_t(_1),...,X_t(_n) )', F'_t is a matrix with dimension n × r (r ≥ m) composed by covariates (e.g., predicted ensembles, latitude, longitude, and elevation) and _t represents the state-space vector of variables with dimension r.

For the purely temporal components in (<ref>) and (<ref>), G_t is an evolution matrix with dimension r, _t are mutually independent and follow a zero mean multivariate Student-t distribution with n_t-1 degrees of freedom and unknown scale matrix W_t which can be estimated using discounting factors. The degrees of freedom n_t-1 and the shape parameters κ_t and κ̅_t are defined through a Beta-Gamma stochastic evolution <cit.>. The parameter δ_V ∈ (0,1] operates as a discount factor, so the larger the discount is, smaller will be the random shock in the observational covariance. When δ_V = 1, the variance is static over time, i.e., σ^2_t = σ^2, ∀ t. The initial information at time t = 0 assumes _0|D_0 ∼ T_n_0(m_0,C_0) and φ_0|D_0 ∼ G(n_0/2,d_0/2).



  §.§.§ Spatiotemporal ensemble model output statistics


Analogously to the model in (<ref>), this proposed model combines the Spatial EMOS with DLMs. Also assuming the structure stated in (<ref>), the spatiotemporal Gaussian model for X_t(s) is given by:


    X_t= F'_t_t + _t,      _t ∼ N(0_n,Σ^*_t), 
    _t= G_t_t-1 + _t,      _t ∼ N(0_p,W_t),

where _t = ( ε_t(_1),…,ε_t(_n))' follows a zero mean multivariate normal distribution with covariance matrix Σ_t^* with elements Σ_t_i,j^* =  D_t_i,iH_i,jD_t_j,j,   i,j = 1,…,n, also with H_i,j = exp(-ϕ_i - _j) and 
    D_t = diag(√(β_0 + β_1 S^2_1,t),...,√(β_0 + β_1 S^2_n,t)),
 an n-dimensional diagonal matrix, such that S^2_i,t is the sample variance of the ensemble forecast for location i at time t. Different from the model in (<ref>), _t is normally distributed with unknown covariance matrix W_t which can also be estimated using discounting factors. The initial information at time t = 0 assumes _0|D_0 ∼ N(m_0,C_0).



 §.§ Inference procedure


Let y=(y_1,...,y_T)' be the collection of T observed time series at n spatial locations over S ⊂ℝ^2 and, Θ=(_0:T,σ^2_0:T, ϕ, λ)'
and ^*=(_0:T,,ϕ,λ)' be the parameter vector in (<ref>) and (<ref>), respectively,
 such that _0:T=(_0,...,_T) and σ^2_0:T=(σ^2_0,...,σ^2_T).

The inference procedure is performed under the Bayesian paradigm and the model specification is complete after assigning a prior distribution to the parameter vectors Θ and Θ^*. An advantage of following the Bayesian paradigm is that the inference procedure is performed under a single framework with uncertainty about parameter estimation being naturally accounted for. Moreover, uncertainty about spatial interpolations and temporal predictions is naturally described through the credible intervals (CI) of the respective posterior predictive distributions. Assuming that the components of Θ and Θ^* are independent a priori, from Bayes' theorem we obtain the following posterior distribution of Θ:

    p(Θ | y)    ∝ p(y|Θ) × p(Θ) 
       ∝∏_t=1^T |Σ_t|^-1/2exp{ -1/2∑_t=1^T (x_t - F_t'_t)'Σ_t^-1(x_t - F_t'_t)}
       ×exp{ -1/2∑_t=1^T (_t - G_t_t-1)'W_t^-1(_t - G_t_t-1)}
       ×∏_{ (i,t): y_t(s_i)>c }{ y_t(s_i)}^λ-1× p(Θ).


The posterior distribution of Θ^* is obtained analogously. The prior distributions assigned to, respectively, Θ and Θ^*, are described in Section <ref>.

The kernel of this distribution does not result in a known closed-form expression. Markov chain Monte Carlo (MCMC) methods are used to obtain samples from the posterior distribution of interest. In particular, the Gibbs sampling algorithm is used when there are missing or censored observations, i.e., respectively, Y_t(s) is missing or Y_t(s) ≤ c; the forward filtering backward sampling algorithm <cit.> for _0:T and σ^2_0:T; and the robust adaptive Metropolis algorithm <cit.> for the remaining parameters.




§ APPLICATION TO WIND SPEED FORECASTS IN THE STATE OF MINAS GERAIS



Following the motivation presented in Sections <ref> and <ref>, the proposed postprocessing extensions are used in this section to calibrate hourly wind speed forecasts (24 hours ahead) for the state of Minas Gerais generated by the Eta regional climate model. In particular, we consider a subset of 59 weather stations and a random selection of 20 days per season from the temporal range of the available dataset. Thus, we obtain the overall and per season results in aggregate predictions for 20 days (480 hours) and 80 days (1920 hours), respectively. The results are compared with those obtained through current postprocessing approaches. In this application, we assess the overall performance of the postprocessing approaches in terms of minimizing systematic errors of wind speed forecasts in the study region.



 §.§ Model settings





For each fitted model, the same mean structure is defined depending on the mean of the ensemble of wind speed forecasts and a set of auxiliary variables defined at each location: roughness length (see calculation details in Appendix <ref>), latitude, longitude, and elevation. Thus, matrix F'_t in (<ref>) has row components 
    (1,f̅(s_i),roughness(s_i),elevation(s_i),latitude(s_i),longitude(s_i), 1,0 )'.
 The ensemble mean f̅ is used to avoid possible multicollinearity problems. This procedure does not entail a great loss of information, since there is a high linear correlation between the members, as shown in Figure <ref>, implying reasonable stability of the numerical forecasts from the Eta model. According to <cit.>, the ensemble mean can capture a possible point anomaly (e.g., cold front), so its use is also indicated by theoretical aspects. The seasonal pattern of solar forcing is represented by a Fourier harmonic through evolution matrix G_t, which is nested in dynamic linear models. In particular, this matrix is given by G_t = diag( G_1,  G_2), where G_1 is an identity matrix of order 6 and 

    G_2 = (
    [ +cos(2 π/24)  sin(2 π/24); -sin(2 π/24)  cos(2 π/24) ]).


To simplify the inference procedure established in Section <ref>, the covariance matrix of the evolution equation W_t in (<ref>) and (<ref>) is estimated using discounting factors <cit.>. We set the discount factors δ_T for trend referring to the intercept, the numerical prediction and the geographic location components, and δ_S for seasonality component. The discount factor δ_V is exclusive to the DGOP model and works as an artificial instrument to insert random shocks in the observational covariance along time, as shown in (<ref>). This analysis considers the exponential correlation function to model spatial dependence. This structure was chosen due to its simplicity. The number of spatial locations in the dataset is not very large (n=68). Several papers have indicated that there would not be enough information coming from data to estimate smoothness parameters in this context, such as <cit.>, which indicates that the parameters in the Matérn model cannot all be estimated consistently. Note, for instance, that if the shortest distance in space is too large, smoothness parameters will be poorly estimated even for large datasets. Given these difficulties for this analysis, we fixed the smoothness and considered the exponential correlation.

Moreover, we assigned reasonably vague priors to model parameters. Specifically, for exclusive parameters of DGOP: _0 ∼ T_1(0,I_8) where I_8 is an identity matrix of order 8 and φ∼ G(1, 1/10); for exclusive parameters of STEMOS: _0 ∼ N(0_8,I_8) and ∼ NT_(0,∞)(0_2,10I_2); and for common parameters: ϕ∼ G( 2, max(d)/6) based on the assumption that the practical range (when the spatial correlation is equal to 0.05) is reached at half of the maximum distance (max(d)/2) between weather station locations and λ∼ N(1,10), such that the expected value is 1, representing the null effect of the Box-Cox transformation. The robustness of these prior distributions in simulated experiments is displayed in Appendix <ref> (Figure <ref>).



To evaluate the improvements in the calibration of wind speed forecasts for the state of Minas Gerais, six models were fitted to the dataset described in Section <ref>. The models were obtained from distinct settings of the proposed models, which are described in Table <ref>. A sensitivity analysis, described in Appendix <ref> (Table <ref>), was performed to choose the best combination of discount factors (δ̃_T, δ̃_S, δ̃_V)' for (δ_T, δ_S, δ_V)'. This was used as an alternative to directly estimate this parameter which avoids over-smoothing the observed data and aims to achieve good predictive power.

All setups considered left censoring with Box-Cox transformation as described in (<ref>). Establishing certain parameters of the proposed models on specific values, we obtain simpler models as particular cases. For instance, if the discount factors are set to 1, the corresponding dynamic parameters become static over time. Setting the transformation parameter λ equal to 1, the Box-Cox transformation keeps the response in its original scale. Finally, when a large value is assigned to the decay parameter of the exponential correlation function ϕ, the correlation matrix H in (<ref>) and (<ref>) tends to an identity matrix resulting in spatial independence between observed locations. We used a moving window training dataset of 10 days (240 hours) as a training period. 



 §.§ Computational details



The MCMC algorithm was implemented in the  programming language, version 3.4.1 <cit.>. To speed up this simulation process, some particular functions were implemented in  language using the library  <cit.> through the package  <cit.>. The analyses were carried out using a laptop computer with a 2.70 GHz Intel Core i7-7500 processor, 32 GB RAM, running with a Microsoft Windows 7 Professional operating system.

For each fitted model over time, convergence tests were performed to check agreement and convergence of two parallel chains starting from different values. The convergence diagnostic was given through the dependence factor <cit.> and the R̂ statistics <cit.>. With chains convergence ensured by these tests, we ran operationally a single chain with 12,500 iterations, discarding the first 500 and sampling at every 5^th step. 






 §.§ Results



Model comparison is performed through the out-of-sample predictive performance using the root-mean-square error (RMSE), the mean absolute error (MAE), the index of agreement <cit.>, and the interval score <cit.>. The goodness of fit also is determined through the well-known Bayesian criteria: the deviance information criterion <cit.> and the logarithm of the pseudo marginal likelihood <cit.>. Each of these criteria is described in more detail in Appendix <ref>. The computational efficiency is also used as a secondary model comparison criterion. 



Table <ref> reports the average of the MAE, RMSE, d, IS, DIC, LPML, and computational time (in minutes) for 24-hour wind speed forecasts at 10-meter height obtained under the Eta model and the six fitted models. All the fitted models present better results than the ensemble mean for all the criteria considered in common. In particular, models considering temporal dynamics performed better than the others according to all predictive performance criteria. Considering the MAE criterion, they were better than ensemble mean and the best static model (MOS), respectively, by 43% and 13%. Similar results were obtained in terms of the RMSE. Concerning the index of agreement d, while the ensemble mean performed better than the static models, the dynamic models outperformed the ensemble mean by 20% and the static models by 32%. Thus, considering the prediction intervals of probabilistic forecasts, the static models obtained the highest values of IS. Moreover, the IS criterion allowed identifying clearer differences between the dynamic models than the other criteria considered. The models that considered the purely spatial component performed significantly worse than the rest in terms of the DIC and LPML. In contrast, the addition of the temporal component benefited both predictive performance and goodness of fit. Remarkably, the DMOS presented better results for all criteria.

Moreover, the boxplots of the average MAE, RMSE and d of 24-hour deterministic forecasts for wind speed at 10-meter height during the seasons are provided as additional figures in Appendix <ref> (Figures <ref> – <ref>). It shows that the forecast evaluation criteria considered, along with wind speed, have different patterns across seasons. Regarding computational time, the DGOP-based models were more efficient on average, taking about 2.4% to 8.9% less time than STEMOS-based models. Its better performance is related to the sequential estimation of dynamic parameters.

Among the model comparison criteria considered, the RMSE has a particular feature that can be used to discriminate different error sources in the context of calibrating wind speed forecasts. In general, error sources are related to either local conditions or general properties of the NWP models. According to <cit.>, this criterion can be decomposed into two additive parts: the amplitude and the phase errors. Figure <ref> shows illustrative panels for each of these additive components recorded at Viçosa weather station. In some cases, these components may be not explicitly well-defined generating clear visualization. Generally, both occur but with different intensities. The amplitude errors are related to forecasts with correct temporal evolution but with a systematic difference from the actual measure. An observation of this error type is shown in Figure <ref>(a). 
 
Note the similarity of the trajectory between both predictions and actual values, except for the level. On the other hand, the phase errors are related to forecasts with correct amplitude but with a mismatch of the temporal evolution of the actual measure. Figure <ref>(b) shows a case for which the numerical forecast does not respect the duration of the intraday seasonal cycle. This can be considered a phase error. Specific corrections are indicated for each component of the error. Statistical postprocessing methods apply linear corrections and can minimize amplitude errors. However, this process is ineffective for phase errors, which are associated with cross-correlation between the forecast and actual time series. Thus, this error type is invariant under linear transformations. Panels of Figure <ref> display the decomposition components of RMSE for deterministic forecasts from probabilistic models during the seasons. Figure <ref>(a) shows the portion of RMSE associated with phase errors. For all seasons, the forecast errors of dynamic models present lower values for this component than the static models. Among the dynamic models, the full-featured models associate their RMSE values slightly less to phase errors than the simplest model DMOS. In this way, since the decomposition elements of RMSE are complementary, the models DGOP and STEMOS associate greater parts of the RMSE value to the amplitude errors compared to the other models, as shown in Figure <ref>(b).


A larger portion of RMSE associated with amplitude errors still allows substantial improvements through linear corrections. Therefore, the dynamic models have greater potential for better calibrated probabilistic predictions among all fitted models. These benefits can be obtained by including more covariate information about the local pattern of the weather phenomenon or more detailed geographical features of relief (e.g., climatology, vegetation type, and presence of water bodies).

Furthermore, our dynamic proposal described in (<ref>) assumes the dynamic variance in time. Therefore, this feature induces fat tails of residuals distribution. In that case, the integrated process is no longer Gaussian but Student-t distributed. Also, Q-Q plots of residuals are presented in Appendix <ref> (Figure <ref>) as an indication of the good fit of the assumed sampling distribution. The STEMOS presents better results in terms of tail modeling, indicating that this model is able to better capture large observations compared to the other fitted models.

An overall calibration assessment of forecasts was carried out through the verification rank histogram <cit.> which evaluates the forecast reliability. An ideal calibration is obtained when its bins are uniformly distributed. Figure <ref> shows the verification rank histograms for forecasts from all fitted models considering the aggregate result for all periods and weather stations in Minas Gerais. The skewed aspect of the histogram for forecasts obtained by the ensemble mean implies that the measurements at weather stations are lower than the predictions, i.e., the ensemble mean commonly overestimates the actual values. All tested calibration models significantly improve the numerical forecasts. Specifically, the forecasts from statistical postprocessing models present similar histograms. These histograms have aspects closer to a uniform distribution than the histogram for uncalibrated forecasts. This result indicates greater reliability. Indeed, the calibration process removes the systematic error resulting in more accurate forecasts as shown in Table <ref>. The forecasts from dynamic models have slightly less underdispersed distribution. Among these, the DGOP and the STEMOS present the best calibration overall. Nevertheless, the observed ∪-shape can represent a lack of variability. The irregularities of boundary bins (bins 1 and 10) suggest that the calibrated forecasts are not ideal at the outermost values. According to wind speed measurements at all weather stations in Minas Gerais, we interpret the outermost values 0.2 and 3.8, respectively, as 0.1 and 0.9-quantiles.




The models with spatial component (see Table <ref>) allow straightforward spatial interpolation through kriging <cit.>, assuming that the latent process X_t() in (<ref>) is a realization of the Gaussian random fields. Therefore, potential improvements in site forecasts may be obtained from the proposed statistical postprocessing processes over the entire weather field. Figure <ref> shows the 24-hour weather field forecast of wind speed at 10-meter height in Minas Gerais starting on 21 July 2016, 1200 UTC. Figure <ref>(a) shows the numerical prediction for the entire weather field obtained using the bilinearly interpolated mean of ensemble forecasts. Figure <ref>(b) shows the deterministic DGOP weather field forecast defined as the median of the predictive distribution. Note that the calibrated forecasts are smoothly spatially distributed in contrast to the numerical prediction. As shown in Figure <ref>, the ensemble members from the Eta model, specifically the ensemble mean, generally overestimate the actual wind speed at 10-meter height. Thus, bias-corrected smooth distribution over the weather field has greater reliability. Lastly, Figure <ref>(c) shows the margin of error of the 90% prediction interval defined as half the width of the credible interval. This plot is useful to indicate regions where there is more uncertainty of predictions. Three regions on the map, closer to the center, are spotted as having larger variability.



Regarding model parameter fitting and interpretation, Figure <ref> shows the posterior distribution behavior of selected parameters from DGOP. Figure <ref>(a) shows the posterior distribution of the Box-Cox transformation parameter λ over the daily range. The estimated value of this parameter ranges between 0.4 and 0.7. The log transformation (λ=0), commonly indicated for wind speed data normalization <cit.>, is not supported by this range. Although the square root transformation (λ=0.5), also indicated, is supported, this parameter performs a cyclical movement related to the seasons. For the summer and spring, which have the highest wind speeds, higher values of λ are suggested. The opposite occurs for seasons with low wind profile (fall and winter). As expected, the behavior of wind speed data is less skewed in fields where the records have a higher level <cit.>. Therefore, values closer to 1 are suggested for the Box-Cox transformation parameter, implying a weak or null (λ=1) effect for this transformation. For this application, these conditions are not supported by the estimated range at any time. These results indicate that the data transformation should not be restricted to a single and static transformation when the temporal windows are wide. If the model is going to be fitted for a whole year, then a λ_t process might be needed. In our case, it is sufficient to assume a constant λ over the training period of 10 days (240 hours). Figure <ref>(b) shows the posterior distribution of the spatial correlation decay parameter ϕ over the daily range. The practical range represents the distance where the spatial correlation between the weather stations decays until almost zero. For example, this distance is estimated to range between 20 and 90 km. In practice, the weather stations within an approximate 2 to 6 grid cell radius are correlated with each other since the Eta model grid has a resolution of 15 km (see Figure <ref>). This parameter also performs a cyclic movement. These values are evidence that the decay parameter can vary over seasons. If a wide temporal window is considered then ϕ_t can be estimated. Note that this approach would be computationally demanding since the precision matrices would be inverted at each time point. Since our aim is the practical use of the proposed models in the context of wind power prediction, we did not pursue this avenue further. Moreover, the parameter ϕ is well estimated only if the spatial correlation decreases to (approximately) zero in the observed spatial domain. If this is not the case for cross-sectional data, then temporal observations may aid in the adequate estimation of this parameter. Thus, fixing ϕ to be the same for each station in our analysis resulted in a consistent estimation of this parameter and satisfactory interpolations.
 

 



§ APPLICATION TO SIMULATED EXPERIMENTS



Following the results for hourly calibration of wind speed forecasts in the state of Minas Gerais presented in Section <ref>, in this section we apply the proposed postprocessing extensions to evaluate the improvements that can be produced under different conditions in simulated experiments. The simulated weather fields have similar characteristics to the real data application: only non-negative values are allowed with left censoring, and skew behavior. The extension of the simulated region is also considered, preserving approximately the same distance between the measurement locations.

We designate the simulated experiments as scenarios. The three scenarios are implemented as follows. The first scenario is a simplistic representation of a weather field generated by the MOS model. This experiment aims to explore the use of full-featured models when an elementary field is available. Although these methods are not parsimonious in this context, we explore how the complexity of the models can impact the calibration process and its results in practice. The remaining scenarios are compound, involving temporal and spatial statistical components, and are generated by DGOP and STEMOS models, respectively. The difference between these experiments is the same as that between the generating models, i.e., the latter accounts for the spread-skill relationship. Furthermore, we investigate the loss of calibration performance when a simpler method is applied in situations where more robust methods are required. In this way, we intend to validate and encourage the use of the proposed postprocessing extensions in applications ranging from simple to challenging. The experiments use the same model settings already considered including the mean structure of the models, the prior distributions, and the training period of 10 days (240 hours). The calibration performance is evaluated as before (see Table <ref>). 

The average of the MAE, RMSE, d, IS, DIC, and LPML for 24-hour simulated wind speed forecasts are shown in Table <ref>. All experiments were generated from a fixed simulated ensemble. The difference between the results obtained by the ensemble mean is related to different structures of the generating models in each scenario. We expected that these results to be similar but not necessarily identical. Scenario 1 is the simplest one regarding reduction of systematic errors according to the comparison criteria. Note that a comprehensible field facilitates the calibration performed by the models. The MOS model, generator of this scenario, obtained the best criteria regarding the adjustment to the data (DIC and LPML), as expected. However, its direct extensions, GOP, DMOS, and DGOP, resulted in similar values. These results confirm the good ability of the proposed postprocessing extensions to adapt parameters to limit values when necessary. The ability to nest precursor models proved to be functional. Thus, the operational use of our extensions is encouraged independently of the weather field characteristics. The gains of scenarios 2 and 3 are smaller compared to the previous scenario due to the greater complexity involved in this simulation. Given these conditions, the performance improved with greater number of featuresin the model. In both scenarios, each proposed extension, DGOP and STEMOS, obtained the best results even when it was not the generating model. For example, these gains ranged from 48% in scenario 3 to 54%, in scenario 2 concerning the ensemble mean for the MAE criterion. The implementation of the MOS model in scenario 3 resulted in worse performance than the ensemble mean in this same criterion, defeating the purpose of postprocessing models. From these experiments, the importance of choosing the model to apply in a given weather field is evident. Applying a simple model in complex scenarios can result in a poor gain or even loss in the calibration process, as seen in scenario 3, while the opposite does not occur, as seen in scenario 1.









§ DISCUSSION



This work presents two new statistical calibration models for weather fields forecasts obtained from regional climate models. The proposed models generalize the well-established spatial postprocessing techniques GOP and SEMOS by combining them with Bayesian dynamic models. Thus, the proposed methods take into account the spatial and temporal correlations, allowing the spatiotemporal calibration of forecasts obtained from NWP models. 

As discussed in Section <ref>, the homogenization of topographic relief achieved by the NWP models results in outputs with potential spatially correlated errors. Thus, the calibration through geostatistical models preserves the inherent spatial correlation structure of the weather field. Additionally, the intrinsic seasonal patterns of meteorological phenomena, discussed in Section <ref>, motivate the use of postprocessing techniques that include temporal correlation in the calibration process as the Bayesian dynamic models do. In contrast with the well-established calibration methods that result in a single forecast for a fixed horizon, the proposed methods allow a single sequential calibration over time, more suitable to hourly data, producing a bias-corrected forecast window output. Moreover, the dynamic structure avoids the issue of empirically optimizing the length of training period due to the intrinsic ability to weight past observations over time <cit.>. As <cit.> suggested, the proposed methods provide an automatic way of choosing the length of the training period. This is balanced such that, parameter estimation is not affected starting from a minimum length threshold, thus avoiding distortions of the calibration process. Also, other advantages of our methods are the natural manipulation of missing data and the versatility of models.
Through the data augmentation technique, the missing data are interpreted as latent data and parameter estimation is performed by a fully Bayesian approach, without a large computational cost.

Our proposed model assumes the Box-Cox transformation parameter λ to be fixed over time and, for our application, this is reasonable for each season. However, if the model is going to be used for the whole year, a moving λ_t might provide better predictions, since asymmetry can vary with the seasons. This issue is a subject of future research.

In particular, we investigate the calibration process of wind forecasts at 10-meter height, which are well known to be locally predicted with systematic errors. As discussed in Section <ref>, these errors are associated with the complex landscape of Minas Gerais and its strong influence on wind speed behavior at low heights. The study region can be considered a low wind speed area. Figure <ref> shows a common pattern recorded at the 68 weather stations: a large number of zero observations, rare wind speeds above 6 m/s, skewed distribution, and different wind speed behavior from other seasons during spring. These characteristics are also commonly observed in many precipitation modeling and forecasting applications. In line with the application performed by <cit.>, the proposals introduce data transformation within the dynamic model, resulting in a flexible sampling distribution for the errors, which can potentially be asymmetric. Hence, our approach leads to narrower predictions in comparison with simpler models without temporal dependence. In this application, the space-time models overall result in slightly better calibration performance, as shown in Table <ref> of Section <ref>. The decomposition of the RMSE, explored also in Figure <ref>, shows that the proposed models associate greater portions of their RMSE with amplitude errors that can still be corrected through linear transformations. This association implies that potential improvements can still be made by including additional information using auxiliary variables. Among the probabilistic models, DGOP and STEMOS produced more reliable forecasts, as shown in Figure <ref>. The persistence of ∪-shape of verification rank histograms indicates a lack of variability. However, the specific irregularities in boundary bins can also be an indicator of conditional biases. <cit.> suggested a local exploratory analysis of model fit and forecast variability. Given the large variability of Minas Gerais relief, the simplistic geostatistical models (see Table <ref>) may also underestimate local characteristics. The exponential correlation function was used due to its good functionality and applicability in the calibration context, as seen in <cit.>, <cit.> and <cit.>. More general correlation functions were not considered due to the difficulty of estimating the smoothness parameters, as discussed in <cit.>.

The inference procedure was carried out through a fully Bayesian approach. The required prior distributions assigned for model specification were weakly informative. This approach has attractive features such as naturally taking into account the uncertainty about parameter estimation, but it is still expensive computationally as shown in Table <ref>. Therefore, it is only recommended for medium-term and long-term forecasts. 

Finally, useful structural improvements can be achieved, for example, by weighting each available ensemble member, as proposed by <cit.>. We work with a particular case in which the same weight is assigned to each member, i.e., the ensemble (arithmetic) mean. <cit.> also proposed the inclusion of local climatological information through auxiliary variables. The proposed models can also be applied to space-time calibration of forecasts from grid-based NWP of other asymmetric censored meteorological variables in which their regime influences seasonal patterns (e.g., precipitation). Additional research on this subject is currently underway.
 



§ ACKNOWLEDGEMENTS

We thank Dr. Chou Sin Chan (CPTEC/INPE) for providing the Eta model outputs used in this study and for her advice during the research. The first author also thanks the financial support from the partnership between Fundação de Amparo à Pesquisa do Estado de Minas Gerais (FAPEMIG) and Companhia Energética de Minas Gerais S.A. (CEMIG) under Project APQ-03813-12.


tocsectionAppendices



§ CALCULATION OF THE ROUGHNESS LENGTH
 

According to <cit.>, the vegetation present on a surface influences the aerodynamic roughness characteristics encountered by the mean wind flow over that surface, affecting both the mean wind speed and direction predicted by numerical models and various other atmospheric parameters. Therefore, the surface roughness length, z_0, defined as the height at which the wind speed equals zero, has an important role in the modeling of atmospheric processes.

The aerodynamic roughness length parameter, z_0, used as a covariate in the application, was estimated for each calibration site from key atmospheric variables, following some principles of the Monin-Obukhov similarity theory <cit.>. Particularly, assuming a logarithmic wind profile, the averaged wind speed u_i at height z_i (10 meters), given by:


    u_i= u_*/k[ ln( z_i/z_0) -Ψ(ζ_i) ]

was used to derive the roughness length z_0 as 


    z_0= z_iexp(-u_i k/u_* - Ψ(ζ_i)),

where k is the von Karman constant, u_* is the friction velocity, Ψ(ζ_i) is the stability correction function of the wind profile and ζ_i=z_i/L is the dimensionless stability parameter given by the height above ground, z_i, normalized by the Obukhov length, L.

Measurements of air temperature, air pressure, sensible heat flux, momentum flux, and wind stress can be used to derive the parameters u_*, L and Ψ(ζ_i). In this application hourly air temperature and air pressure data were obtained from the available meteorological stations, while hourly reanalysis data from the CFSV2 model <cit.> were used as the source of heat and momentum fluxes, after interpolation from the CFSV2 regular grid to the calibration sites.

Hourly values of the roughness length, z_0, were first obtained from (<ref>) for the two years of available data. Then, median values of z_0 by month and by hour within each month (288 values in total for each calibration site) were calculated, considering only those z_0 values that were estimated during neutral conditions of atmospheric stability (such condition is achieved when |L|>500).





§ ROBUSTNESS OF PRIOR DISTRIBUTIONS



In this section, we verify the robustness of prior distributions assigned to the parameters of the proposed postprocessing models applied in Section <ref>. Based on these prior distributions in simulated experiments, the posterior and prior distributions associated with the parameter vectors from DGOP and STEMOS are presented in Figure <ref>. For clarity of exposition, only the results of a subset of common parameter _0 (vector with dimension r) are exhibited. Specifically, the intercept parameter at moment t=0 represented by θ_0,0 is shown. 

Even when using this specific set of prior distributions in a distinct application from that previously reported, Figure <ref> shows that the posterior distributions differ significantly from the prior distributions for all parameters. This is evidence that the gain of information is mostly provided by the data, ensuring that non-informative properties are preserved in general applications.







§ SENSITIVITY ANALYSIS OF DISCOUNT FACTORS



This sensitivity analysis was performed to choose the best combination of discount factors (δ̃_T, δ̃_S, δ̃_V)' for (δ_T, δ_S, δ_V)'. In particular, we also consider a subset of 59 weather stations and a random selection of 5 days per season from the temporal range of the available dataset. Following the same model comparison criteria described in Section <ref>, Table <ref> reports the sensitivity analysis of discount factors through average MAE, RMSE, d, IS, DIC, and LPML values for 24-hour wind speed forecasts at 10-meter height. The setups of the dynamic models differ only in the discount factor values. The results indicate greater sensitivity for δ_T, which significantly reduces the predictive performance when assuming lower values. In contrast, the variation of values for δ_S and δ_V only slightly alters the predictive performance. Thus, we determine (δ̃_T, δ̃_S, δ̃_V)' = (0.99, 0.95, 0.99)'.







§ MODEL COMPARISON CRITERIA
 

In this section, we briefly describe the model comparison criteria used to compare the prediction of the fitted models in Section <ref>. The first three criteria (RMSE, MAE and index of agreement) are appropriate to compare numerical predictions from the Eta model, which provides only deterministic estimates, with the proposed postprocessing models. The probabilistic forecasts are evaluated through IS, which takes into account the amplitude and coverage of the prediction intervals in a parsimonious way.



 §.§ Mean absolute error and root-mean-square error

Standard measures of goodness of fit were also entertained in this
study for comparison purposes. The root-mean-square error
(RMSE) and the mean absolute error (MAE) are given by:

    =1/nT∑_i=1^n∑_t=1^T(y_t(s_i)-ŷ_t(s_i))^2
         =1/nT∑_i=1^n∑_t=1^T|y_t(s_i)-ŷ_t(s_i)|,

where ŷ_t(s_i) is obtained through a Monte Carlo estimate of the posterior mean of the predictive distribution, that is, E[y_t( s_i) | y], across N draws. Smaller values of RMSE and MAE indicate better model fit.



 §.§ Index of agreement
 
<cit.> introduced a standard measure for assessing the quality of forecasts. The index of agreement (d) ranges between 0 (absence of agreement) and 1 (perfect agreement), and is given by:

    d = 1 - ∑_i=1^n ∑_t=1^T(y_t(s_i) - ŷ_t(s_i))^2/∑_i=1^n ∑_t=1^T(|ŷ_t(s_i) - y̅|+|y_t(s_i) - y̅|)^2,

where y̅ = 1/n∑_i=1^n∑_t=1^T y_t(s_i).



 §.§ Interval score

The interval score <cit.> is a scoring rule for interval predictions considering the symmetric prediction interval with level (1-α)×100%. The score is rewarded by accurate intervals and penalized when there is no coverage of the forecast. If actual values are contained in the prediction interval, this measure is reduced to the range amplitude. The average IS is given by:

    IS =      1/nT∑_i=1^n∑_i=1^T (û_t(s_i)-l̂_t(s_i)) 
        +2/α(l̂_t(s_i)-y_t(s_i))1{ y_t(s_i)<l̂_t(s_i) }
        +2/α(y_t(s_i)-û_t(s_i))1{ y_t(s_i)>û_t(s_i) }

where l̂_t(s_i) and û_t(s_i) are, respectively, the lower bound obtained by the α/2 quantile, and the upper bound, obtained by the 1-α/2 quantile based on the predictive distribution. The indicator function is represented by 1.

Smaller IS values indicate more efficient probabilistic forecasts.



 §.§ DIC


Particularly useful in Bayesian model selection problems, the deviance information criterion <cit.> is a hierarchical modeling generalization of the Akaike information criterion <cit.>. This criterion is negatively oriented implying that models with smaller DIC should be preferred to models with larger DIC. The DIC is given by:

    DIC = -2 ∫log{p(y|Θ)}p(Θ|y)dΘ + P_D,

where y is a vector of observed values and Θ is the parameter vector. Thus, p(y|Θ) represents the likelihood function and p(Θ | y), the posterior distribution. Defined as a Bayesian measure of model complexity, P_D is given by:

    P_D = 2 log{p(y|Θ̃)} -2 ∫log{p(y|Θ)}p(Θ|y)dΘ,

with Θ̃ denoting the Bayes estimator. 

A Monte Carlo approximation of DIC is given by:

    DIC = -2 1/M∑_m = 1^Mlog{p(y|Θ^(m))} + P̂_D,

where 
    P̂_D = 2 log{p(y|Θ̃)} -2 1/M∑_m = 1^Mlog{p(y|Θ^(m))},

with Θ^(m) denoting the m-th MCMC sample of Θ from posterior distribution p(Θ|y), m = 1, ..., M.



 §.§ LPML


A component of the Bayes factor, the logarithm of the pseudo marginal likelihood <cit.> is given by:

    LPML = 1/n∑_i=1^nlog{CPO_i},


where CPO_i represents the conditional predictive ordinate (CPO) for location i and is given by:

    CPO_i =       p(y(s_i)| y(s_-i)) 
    
    =      (∫1/p(y(s_i)|Θ)p(Θ|y)dΘ)^-1,

with y(s_i)=(y_1(s_i), ..., y_T(s_i))' and y(s_-i)=(y(s_1),...,y(s_i-1), y(s_i+1), ..., y(s_n))', i = 1,...,n. Note that the CPO_i is based on the leave-one-out-cross-validation process and estimates the probability of y(s_i) given the observation of y(s_-i). In particular, a Monte Carlo approximation of CPO_i is given by:

    CPO_i = (1/M∑_m = 1^M1/p(y(s_i)|Θ^(m)))^-1,

with Θ^(m) denoting the m-th MCMC sample of Θ from posterior distribution p(Θ|y), m = 1, ..., M.

Finally, the Monte Carlo approximation of LPML is given by: 


    LPML = 1/n∑_i=1^nlog{CPO_i}.


The preferred model maximizes this criterion.





§ SUPPLEMENTARY RESULTS














jasa


