Julie Bessac

Valérie Pierre Ailliot

V Monbet

Ms J Bessac

Gaussian linear state-space model for wind fields in the North-East Atlantic

Keywords: Multisite wind generators, Space-time model, State-space model, EM algorithm, Identifiability

Gaussian linear state-space model for wind fields in the

North-East Atlantic

Introduction

Many natural phenomena and human activities depend on wind conditions. However meteorological data are often available over periods of time that are not long enough to estimate reliably probabilities of complex events. In order to overcome this insufficiency, stochastic weather generators have been developed. Those stochastic weather generators are statistical models that simulate sequences of meteorological variables with statistical properties similar to the ones of the observations. They have been adopted in impact studies as a computationally inexpensive tool that generates quickly as many synthetic time series of unlimited length as desired, see for instance [START_REF] Srikanthan | Stochastic generation of annual, monthly and daily climate data: A review[END_REF] and references therein. Stochastic weather generators can be adapted to in-filling tools that simulate missing data [START_REF] Yang | Spatialtemporal rainfall simulation using generalized linear models[END_REF] or to downscaling global climate models, see for instance [START_REF] Maraun | Precipitation downscaling under climate change: Recent developments to bridge the gap between dynamical models and the end user[END_REF] and references therein. Wind generators have in particular been used to assess various quantities related to wind power production [START_REF] Brown | Time series models to simulate and forecast wind speed and wind power[END_REF][START_REF] Castino | Stochastic modelling of wind velocities time series[END_REF][START_REF] Hofmann | Nowicob-a tool for reducing the maintenance costs of offshore wind farms[END_REF], drift of objects in the ocean (Ailliot et al., 2006a) or coastal erosion [START_REF] Skidmore | Stochastic wind simulation for erosion modeling[END_REF].

A review of stochastic models for wind time series can be found in [START_REF] Monbet | Survey of stochastic models for wind and sea state time series[END_REF]. Most of the existing models are designed for wind time series at a single location. The most classical approach consists in using the Box-Jenkins methodology, where an ARIMA model is fitted after achieving stationarity and applying a marginal transformation to obtain Gaussian like margins. Nonlinear models have also been proposed and, in particular, weather type models with a discrete latent variable, see [START_REF] Ailliot | Markov-switching autoregressive models for wind time series[END_REF] and references therein.

Generalizations to space-time models have been explored recently. Multisite wind models have to deal with the temporal and spatial dependence and it is known that these two components are generally not separable when air masses are moving in a prevailing direction [START_REF] Gneiting | Nonseparable, stationary covariance functions for spacetime data[END_REF]. Black-box models such as artificial neural networks may be fitted but they lead to noninterpretable models [START_REF] Lei | A review on the forecasting of wind speed and generated power[END_REF]. A first alternative is based on Gaussian fields [START_REF] Gneiting | Nonseparable, stationary covariance functions for spacetime data[END_REF][START_REF] Rychlik | A spatial-temporal model for wind speeds variability[END_REF] where non-separable parametric covariance functions can be considered to take into account the mean displacement of the air masses. Another approach consists in using vector AutoRegressive-Moving-Average models [START_REF] Haslett | Space-time modelling with long-memory dependence: Assessing ireland's wind power resource[END_REF][START_REF] De Luna | Predictive spatio-temporal models for spatially sparse environmental data[END_REF] where the wind dynamic is described by the autoregressive matrices. Motions can be introduced using covariates or latent variables. For example, in [START_REF] Ailliot | An autoregressive model with time-varying coefficients for wind fields[END_REF]) the autoregressive coefficients depend on a latent process that describes the motion of the air masses. In [START_REF] Šaltytė Benth | Spatial-temporal model for wind speed in Lithuania[END_REF], a latent field describes the spatial structure of the autoregressive parameters at each station. Following similar ideas, various authors developed models that aim at embedding physical insights into a probabilistic model. The Bayesian framework is very convenient to deal with such coupling [START_REF] Wikle | Spatiotemporal hierarchical bayesian modeling tropical ocean surface winds[END_REF]. For instance, in [START_REF] Milliff | Ocean ensemble forecasting. part i: Ensemble mediterranean winds from a bayesian hierarchical model[END_REF], classical partial differential equations for the wind at the sea surface are perturbated by adding a white noise and the parameters are estimated following a Bayesian inference method.

In the present paper, a structural model that aims at simulating wind speed at several locations is investigated. The main idea consists in introducing a latent variable that aims at describing regional wind conditions and the observed local wind is modeled as a function of the regional wind at different lags in order to reproduce the mean displacement of the air masses. The framework is kept simple with a linear Gaussian model used to describe both the dynamics of the latent process and the link between the latent and the observed process.

It leads to an interpretable model with efficient numerical procedures available for parameter estimation and simulation. Despite its simplicity, the model leads to non-separable and anisotropic covariance functions. No physical equations were embedded because their resolution is generally computationally too expensive for a stochastic generator but the suggested model involves quantities that have a physical meaning in the proposed context. It could be used as a surrogate of the atmospheric model (emulator) for data assimilation or data fusion.

The data considered in this paper are presented in Section 2. The model is described in Section 3. Parameter estimation and fitting procedures are also discussed in this section. The model is validated in Section 4 and it is shown in particular that the fitted model is able to reproduce the anisotropy and non-separability of the data. However, the model includes a large number of parameters and various reduced models are introduced in Section 5. Conclusions are given in Section 6. Parameter identifiability and non full-symmetry are proven in the supplementary materials.

The wind dataset

In situ data are neither available on a long time period nor on a large area offshore Brittany in France. Reanalysis data, which are obtained by combining observations with numerical weather prediction models, provides a relevant alternative for meteorological or climatological studies. In this paper we consider wind speed at 10 meters above sea level extracted from the ERA Interim Full dataset produced by the European Center of Medium-range Weather Forecast (ECMWF). It can be freely downloaded and used for scientific purposes at the URL http://data.ecmwf.int/data/. This dataset is available on a regular space-time grid with a temporal resolution of 6 hours and a spatial resolution of 0.75 • . However the methodology introduced in this paper could easily be adapted to handle datasets with a more complicated space-time sampling such as the one obtained when considering networks of meteorological stations.

We focus on 18 gridded locations between latitudes 48 • N and 49.5 • N and longitudes 6.25 • W and 9 • W (see Figure 1). The dataset consists of 33 years of wind data from 1979 to 2011 and we focus on the month of January. Further, the statistical inference is based on the assumption that the 33 months of January are 33 independent realizations of a common stationary stochastic process. This assumption is usual for meteorological processes but it does not take into account low frequency variations such as the North Atlantic Oscillation (NAO).

In the studied area prevailing air masses moves are generally eastward.

It induces non-separability and non full-symmetry properties of the space- time covariance function of the wind speed as for the dataset of wind speed in Ireland considered in [START_REF] Haslett | Space-time modelling with long-memory dependence: Assessing ireland's wind power resource[END_REF][START_REF] Gneiting | Nonseparable, stationary covariance functions for spacetime data[END_REF]. The lagged by 1 cross-correlations shown in Figure 2 highlight this phenomenon.

Indeed, the asymmetry with respect to the difference of longitude shows that the correlation between y t (p) and y t+1 (p ) is higher when location p is more westerly with respect to p than when p is easterly with respect to p ; in average western locations see the meteorological events before the eastern locations.

This asymmetry is less pronounced in latitude but reveals flows from north to south. Furthermore, the correlations reveal some anisotropy as dependences in latitude and longitude differ (see Figure 2).

Wind speed distribution is known to be skewed. It is often modeled as a Weibull distribution [START_REF] Brown | Time series models to simulate and forecast wind speed and wind power[END_REF] but other distributions such as the skew normal distribution have also been considered [START_REF] Flecher | A stochastic daily weather generator for skewed data[END_REF].

A classical method to handle such asymmetry in time series analysis consists in applying a Box-Cox transformation in order to get a time series with approximately Gaussian marginal distribution. This method has been extensively used for analyzing wind time series at a single location, see for example [START_REF] Brown | Time series models to simulate and forecast wind speed and wind power[END_REF]. In [START_REF] Rychlik | A spatial-temporal model for wind speeds variability[END_REF]) a different power transformation λ i is used at each location. More precisely, let us denote

     y λ i ,i,t = y λ i i,t -1 λ i if λ i > 0 y λ i ,i,t = log(y i,t ) if λ i = 0,
with y i,t the wind speed at time t and location i. Following [START_REF] Hinkley | On quick choice of power transformation[END_REF], λ i can be estimated by searching the roots of the asymmetry measure

S(λ i ) = mean(y λ i ,i,t ) -median(y λ i ,i,t ) var(y λ i ,i,t ) . ( 1
)
The resulting estimates are shown in Figure 1 with values ranging from about 1 (Gaussian distribution) in the north-west to 0.7 closer in the south-east.

Despite this spatial variability, we have chosen to use the same power transformation at all sites in order to preserve the spatial structure of the wind fields as done in [START_REF] Haslett | Space-time modelling with long-memory dependence: Assessing ireland's wind power resource[END_REF]. The value λ = 0.85 is used in the sequel. It is the average value of the λi shown on Figure 1. The simulation results given in Section 4 (see Figure 5) indicate that this simple transformation permits to reproduce the marginal distributions of the wind data considered in this study.

A linear Gaussian state-space model for wind speed

State-space models first appeared in engineering and have then been extensively used in many domains. State-space representations bring a very flexible framework for modeling time series [START_REF] Durbin | Time series analysis by state space methods[END_REF][START_REF] Brockwell | Time series: theory and methods[END_REF] and space-time processes [START_REF] Wikle | A general science-based framework for dynamical spatio-temporal models[END_REF]. The model introduced in this section is a linear Gaussian state-space model. One of the main advantages of this class of models is that estimation, forecasting and smoothing can be processed through general and efficient procedures.

Model

The observed wind fields are generally smooth, which leads to a high correlation between the different sites. Although the smoothness observed here is inherent to reanalysis data that are known to be smoother than observations [START_REF] Milliff | Ocean ensemble forecasting. part i: Ensemble mediterranean winds from a bayesian hierarchical model[END_REF], it is coherent with the considered spatio-temporal scale.

This regularity suggests to explain an important part of the multisite wind by using a common scalar process (the 'regional wind condition'). This scalar process, denoted by {X t } in the sequel, can not be observed directly and is thus introduced as a latent (or 'hidden') process. In order to model the prevail-ing motion of the air masses we propose to let the wind conditions at western locations depend more on the leading one-lag X t+1 and X t signals than on the lagged signal X t-1 with the reverse phenomenon at eastern locations. More precisely, the Gaussian state-space model, which is considered in this paper, is defined as

(M)      X t+1 = ρX t + σ t+1 , Y t = α 1 X t+1 + α 0 X t + α -1 X t-1 + Γ 1/2 η t for t ≥ 0, Y t ∈ R K
is the observed process, its K coordinates correspond to the meancorrected transformed wind speed at the K = 18 locations. { t } and {η t } are independent Gaussian white noises with zero-means and identity covariance matrices. α 1 , α 0 and α -1 are K-dimensional vectors that link the lagged values of the regional process {X t } to local wind conditions. The covariance matrix Γ ∈ R K×K models the spatial structure of small-scale fluctuations. In finance and economics or when high dimensional data are considered this covariance matrix is often assumed to be diagonal [START_REF] Wikle | A general science-based framework for dynamical spatio-temporal models[END_REF].

Here it would imply that the local wind conditions are conditionally independent given the regional conditions, which is a very strong assumption. As a first step, we have chosen to work with a full non-parametric covariance matrix but reduced parametric models are explored in Subsection 5.1. In the sequel, we denote Λ = (α 1 |α 0 |α -1 ) ∈ R K×3 and θ = (ρ, σ, Λ, Γ) the unknown parameters.

The temporal dynamics of the observed process is mainly contained in the latent process {X t } and explained by the coefficient ρ. The model thus imposes the same long-term temporal dynamics at each location. Under the assumption |ρ| < 1, the AR(1) process {X t } is stationary and so is the process

{Y t }.

Second-order structure and identifiability

Identifiability is required to get sensible and reliable parameter estimates. The introduction of a latent process {X t } is a source of non-identifiability since the unknown parameters need to be identified uniquely from the distribution of the observed {Y t } and Gaussian linear state-space models are known to be nonidentifiable without additional constraints [START_REF] Hannan | The statistical theory of linear systems[END_REF][START_REF] Ljung | System identifiability[END_REF][START_REF] Bai | Identification and estimation of dynamic factor models[END_REF][START_REF] Bork | Macro factors, Monetary policy analysis and affine term structure models[END_REF]. Identifiability of linear Gaussian state-space models was initially investigated in control theory and has been largely explored during the last decades. However we could not find any result that applies directly to the model considered in this paper.

{Y t } is a zero-mean stationary Gaussian process that is thus characterized by its second-order structure given below

cov θ (Y t , Y t ) = σ 2 1 -ρ 2 α 1 (α 1 + ρα 0 + ρ 2 α -1 ) t + α 0 (ρα 1 + α 0 + ρα -1 ) t + α -1 (ρ 2 α 1 + ρα 0 + α -1 ) t + Γ, (2)
cov θ (Y t , Y t+1 ) = σ 2 1 -ρ 2 α 1 (ρα 1 + α 0 + ρα -1 ) t + α 0 (ρ 2 α 1 + ρα 0 + α -1 ) t + ρα -1 (ρ 2 α 1 + ρα 0 + α -1 ) t , ( 3
) cov θ (Y t , Y t+k ) = σ 2 1 -ρ 2 ρ k-2 (α 1 + ρα 0 + ρ 2 α -1 )(ρ 2 α 1 + ρα 0 + α -1 ) t , (4)
for all k ≥ 2.

The study of this space-time covariance function leads to the following Proposition, which is proven in the supplementary materials.

Proposition 1 Assume that (M) holds. Assume further that σ 2 1-ρ 2 = 1 and that the vectors α 1 , α 0 and α -1 are linearly independent. Then the parameters can be identified from the distribution of the process {Y t }.

These identifiability constraints are interpretable and were always satisfied when fitting the model to the data. The first one implies that X t has a unit variance, the variance of the wind at the different locations being explained by the scaling matrix Λ. The second one implies that Y t actually depends on the three lagged values X t-1 , X t and X t+1 and not only on one or two lagged values.

We will see in Section 4 that the proposed model enables to reproduce various complex properties of the observed space-time covariance. Under constraints of Proposition 1, the covariance defined by (2-4) is neither fullsymmetric nor separable (see the supplementary materials). Other non-symmetric space-time covariance models have been proposed in the literature.

Some of them have been fitted to the Irish wind dataset, see for instance [START_REF] Gneiting | Nonseparable, stationary covariance functions for spacetime data[END_REF]. They generally rely on strong assumptions such as spatial stationarity and isotropy, which are not realistic for our dataset. A noticeable exception is the model proposed in (de Luna and Genton, 2005) which is based on the specification of a vector autoregressive process and captures a part of the anisotropy that is observed on the Irish dataset.

Parameter estimation

Two methods of estimation have been implemented and compared. The first one is a method of moments based on the second-order structure of the process {Y t } given by (2-4). It consists in numerically minimizing the following objective function

θ → cov(Y t , Y t ) -cov θ (Y t , Y t ) 2 2 + cov(Y t , Y t+1 ) -cov θ (Y t , Y t+1 ) 2 2 (5) + cov(Y t , Y t+2 ) -cov θ (Y t , Y t+2 ) 2 2 + cov(Y t , Y t+3 ) -cov θ (Y t , Y t+3 ) 2 2 ,
where cov denotes the empirical covariance function and . 2 stands for the matrix Frobenius norm. This method, denoted by GMM for Generalized Method of Moments in the sequel, is usual in geostatistics [START_REF] Cressie | Statistics for spatial data[END_REF]. We have chosen to consider only the first four lags of the autocovariance function when building the objective function ( 5). It corresponds to the minimal number of terms needed to identify the parameters (see the supplementary materials).

Simulation results indicate that including more lags does not lead to more accurate estimates.

The second method performs Maximum Likelihood (ML) estimation using the Expectation-Maximization (EM) algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the em algorithm[END_REF]. The EM algorithm aims at maximizing the incomplete log-likelihood function

θ → E(log(p(X 1 , ..., X T , Y 1 , ..., Y T ; θ))|Y T 1 = y T 1 )
by performing recursively two steps (E-step and M-step). For linear Gaussian state-space models efficient numerical procedures exist for both steps. In the E-step, the Kalman recursions lead to an exact computation of the various conditional expectations involved and in the M-step analytical expressions of the maximizers of the intermediate function are available. More details about the Kalman recursions and EM-algorithm can be found in the supplementary materials.

Both methods are sensitive to the initial parameter values, which need to be chosen carefully. We used the following procedure that involves the properties of the second-order structure of {Y t }:

-ρ = cov(Y t , Y t+3 ) i,j cov(Y t , Y t+2 ) i,j
for all i, j ∈ {1, ..., K} is initialized as the empirical

mean of cov(Y t , Y t+3 ) i,j cov(Y t , Y t+2 ) i,j .
-Λ is estimated by minimizing

θ Λ → cov(Y t , Y t+1 ) -cov θ (Y t , Y t+1 ) 2 2 + cov(Y t , Y t+2 ) -cov θ (Y t , Y t+2 ) 2 2
as a function of Λ with ρ being fixed to the value obtained in the previous step. Note that this function does not depend on Γ according to (3) and ( 4).

-Γ is determined by minimizing

θ Γ → cov(Y t , Y t ) -cov θ (Y t , Y t ) 2 2
as a function of Γ with ρ and Λ being fixed to the value obtained in the previous steps.

These rough estimates are used as initial conditions of the numerical optimization of the function (5) to compute the GMM estimates, which in turn are used to initialize the EM algorithm. An extra step could be added to refine the output of the EM algorithm with a numerical optimization of the likelihood function, which is known to be more efficient close to local maxima [START_REF] Durbin | Time series analysis by state space methods[END_REF]). However we did not find any improvements in practice with such a procedure.

Properties of the estimates

Under suitable conditions, GMM [START_REF] Newey | Large sample estimation and hypothesis testing[END_REF] and ML [START_REF] Newey | Large sample estimation and hypothesis testing[END_REF][START_REF] Shumway | Time series analysis and its applications[END_REF][START_REF] Hannan | The statistical theory of linear systems[END_REF][START_REF] Caines | Linear stochastic systems[END_REF] estimates are consistent and asymptotically Gaussian. In order to assess the performances of the estimates for the practical application considered in this paper, we perform a simulation study. N = 100 independent sets of the size of the studied data are simulated for the parameters set estimated by ML on the wind data.

Results

In order to validate the proposed model we check its physical realism and its ability to generate artificial wind conditions with statistical properties similar to the ones of the dataset. We compare the GMM and ML estimates through this validation in order to investigate their robustness in a practical context.

Interpretability

The loading matrix Λ links the latent process to the observed wind conditions.

The values of α 1 and α -1 shown on Figure 3 reveal the site-dependent relations with the latent process. Western locations depend more on X t+1 than on X t-1 and the reverse is true for eastern locations. This was expected since western locations are the first locations affected when meteorological events enter in the studied region.

Since large-scale variability is supposed to be contained in the latent process, Γ should contain only small-scale variations. This is confirmed when comparing the spatial sill and range of Γ with the ones of the original covariance function of the data (see Figure 4). The shape of Γ has a block structure that is induced by the geometry of the domain and the numbering of the sites (see Figure 1). The level sets of the blocks, except the top right corner (and by symmetry bottom left corner), look like saddle point level sets: the model better explains the wind observed at the central locations of the domain than at the locations that are close to the boundary. The top right corner has elliptical level sets. These geometrical differences raise problems when trying to develop simple parametric models for Γ (see Section 5.1).

Realism of simulated sequences

In order to further validate the model, we have checked its ability to simulate realistic wind conditions. For that, artificial time series are simulated with the fitted models and their statistics are compared with the ones of the original data. According to the quantile-quantile plots shown on Figure 5, the model is able to reproduce the general shape of the marginal distribution of the process at the central station 9 except for very low wind speed. Similar results were obtained at other locations.

Figure 7 shows that the cross-correlations at lags 0 and 1 are well reproduced by the fitted models with a slightly better fit for the GMM estimates.

This was not unexpected since the GMM is designed to make the first lags of the empirical autocovariance function coincide with the one of the fitted model. Figure 6 shows however that the fit is better for lags greater than one day with the ML estimates, which take into account longer term dynamics and leads to a higher value of ρ (0.76 for ML against 0.70 for GMM). The better fit of the ML estimates is also coherent with Table 1. Note also that the models reproduce the time shift between locations 13 and 18 that is induced by the prevailing westerly flow (see Figure 6).

Some improvements of the model

In this section we explore reduced models for Γ and Λ with the aim of reducing the number of parameters involved in the model.

Parameterization of Γ

The spatial structure of the estimated Γ shown on Figure 4 suggests to model the covariance between locations i and j in {1, ..., K} as a function of the distance d i,j between these locations. In the sequel, we consider two different models, one with Gaussian correlation function Γ i,j = σ i σ j (exp(-λ 1 d 2 i,j ) + λ 2 δ i,j ) for i, j ∈ {1, ..., K}, and the other with wave correlation function

Γ i,j = σ i σ j sin(λ 1 d i,j ) λ 1 d i,j
+ λ 2 δ i,j for i, j ∈ {1, ..., K}, where (σ 1 , ..., σ K , λ 1 , λ 2 ) are positive parameters and δ i,j denotes the Kronecker delta. λ 1 and λ 2 are respectively the range and nugget parameters, and σ 2 i (1 + λ 2 ) represents the variance of the field at location i. These models are well defined covariance functions [START_REF] Cressie | Statistics for spatial data[END_REF][START_REF] Abrahamsen | A review of Gaussian random fields and correlation functions[END_REF] and are denoted respectively (M Γ∼Gauss ) and (M Γ∼Sinus ) hereafter.

The difference in dependence on latitude and longitude of Γ suggests the use of an anisotropic distance [START_REF] Refice | On the use of anisotropic covariance models in estimating atmospheric dinsar contributions[END_REF][START_REF] Haskard | An anisotropic matérn spatial covariance model: Reml estimation and properties[END_REF][START_REF] Šaltytė Benth | Spatial-temporal model for wind speed in Lithuania[END_REF] d i,j = ∆Lat(i, j) 2 + θ 1 ∆Long(i, j) 2 + θ 2 ∆Lat(i, j)∆Long(i, j)

where ∆Lat(i, j) and ∆Long(i, j) denote respectively the difference in latitude and longitude between locations i and j expressed in kilometers. The

constraint θ 1 > θ 2 2
4 is imposed to ensure the positive-definiteness of the distance.

These covariance structures have first been fitted by least square estimation to the estimated Γ shown on Figure 4. The fit is globally good for the wave covariance whereas the Gaussian shape can not cope with the negative correlations observed between western and eastern locations. However the covariance between the northern and southern locations are poorly reproduced. As mentioned in Section 4.1 these blocks have a particular elliptical shape that is difficult to reproduce with parametric models. Estimated anisotropy coefficients for the sinus and the Gaussian structures are respectively ( θ1 , θ2 ) = (0.2, 0.04) and ( θ1 , θ2 ) = (0.23, 0.005). For both models θ 1 is lower than one and θ 2 is close to zeros and thus the spatial range is maximum in the west-east direction.

In a second step, the parameters have been re-estimated using the GMM and ML methods. A numerical optimization needs to be performed in the M-step of the EM algorithm to update the values of (σ 1 , ..., σ K , λ 1 , λ 2 ). Note that the function to minimize can be expressed in a compact way (see supplementary materials) that leads to an efficient numerical procedure. The Similar results were obtained for the marginal distributions and the temporal correlation functions but the description of the spatial structure was deteriorated when using a (M Γ ) model instead of (M). This miss-specification is also confirmed by the Bayes Information Criterion (BIC) values given in Table 2 where BIC = -2 log L + N p log(N obs ) with L the likelihood of the model, N p the number of parameters and N obs the number of observations. The reduced models (M Γ ) are clearly outperformed by the full model (M). Other parametric models such as the Matérn one have been tried without more success and it seems difficult to find a simple reduced model that can reproduce all the complexity of the covariance matrix Γ of the observation error.

Parameterization of Λ

The structure of α 1 , α 0 and α -1 reveals a quadratic dependence in longitude and the dependence in latitude suggests the use of an intercept depending on latitude (see Figure 8). The following parameterization is then proposed The parameterization is easily handled in the GMM procedure whereas a numerical optimization is again needed to update Λ in the M-step of the ML procedure. Moreover a joint optimization on Λ and Γ should be done since both of them are involved in the same part of the log-likelihood function. In main innovation, with respect to the other space-time models that have been proposed for meteorological variables, is the introduction of a continuous latent process describing regional conditions. The proposed models are interpretable and can reproduce the marginal distribution of the wind speed and important properties of the space-time covariance structure such as the asymmetries induced by prevailing motions of the air masses.

Λ = 1 | Long | Long 2       β Lat 1 β Lat 4 β Lat 7 β 2 β 5 β 8 β 3 β 6 β 9       351
An important advantage of Gaussian linear state-space models is that efficient and easy to implement procedures of estimation are available. Two estimation procedures, one based on a method of moments (GMM) and the other on the likelihood function (ML) have been compared. GMM yields to better results when looking at the short-term space-time structure but ML is better in reproducing the long-term dynamics.

According to the BIC values given in Table 2, the ranking of the model coincides with the complexity of the model and the quality of the model is systematically worsened when the number of parameters is reduced. Note that higher-order autoregressive models have been considered for modeling the dynamics of the hidden state but they led to very slight improvements and are not further discussed here (the model with autoregressive models of order 2, denoted (M 2 ), is given in Table 2). In order to check the relevance of the BIC criterion, we have performed a cross-validation study (see supplementary materials) that confirmed the ranking of the models given by BIC. Similar results were obtained on the Irish wind dataset considered in [START_REF] Haslett | Space-time modelling with long-memory dependence: Assessing ireland's wind power resource[END_REF][START_REF] Gneiting | Nonseparable, stationary covariance functions for spacetime data[END_REF], which has a different space-time resolution with daily data and stations on an irregular spatial grid. This highlights the difficulty to find parsimonious and realistic models for describing the spacetime evolution of wind.

Figure 1 :

 1 Figure 1: Left panel: mean wind speed at the 18 numbered points under study in the North-East Atlantic. Right panel: estimated values of the power in the Box-Cox method at the 18 locations.

Figure 2 :

 2 Figure 2: Lagged-one cross-correlations against differences of latitude (left) and longitude (right).

Figure 3 :

 3 Figure 3: ML estimates of α 1 (left panel) α 0 (middle panel) and α -1 (right panel).

Figure 4 :

 4 Figure 4: Empirical covariance matrix of the wind data (left) and ML estimates of Γ (right).

Figure 5 :

 5 Figure 5: Quantile-Quantile plot at location 9 for the model (M) and the parameters estimated by GMM (left) and ML (right). The dashed lines correspond to 90% prediction intervals computed by simulation.

Figure 6 :Figure 7 :

 67 Figure6: Observed (full lines) and simulated (dashed lines) cross-correlations between locations 13 and 18 (upper row) and auto-correlation at location 9 (lower row) for the model (M) with parameters estimated by GMM (left) and ML (right). 90% prediction intervals are computed from 100 independent simulated samples of the size of the original data.

Figure 8 :

 8 Figure 8: Estimated α 1 (top), α 0 (middle) and α -1 (bottom) against longitude at latitude 48 • N (left) and against latitude at longitude 6.75 • W (right). Solid line: ML estimation of Λ for model (M) , dashed line: parametric structure fitted by least square.

Table 1 :

 1 Bias, standard deviation and RMSE of parameters estimates. For the multidimensional parameters, minimal and maximal values are given in brackets.

		Bias		Sd		RMSE
	Parameters	GMM	ML	GMM	ML	GMM	ML
	ρ	0.036	0.004	0.022	0.017	0.042	0.017
	α 1	[-0.11;-0.009] [-0.069;-0.019] [0.065;0.108] [0.071;0.097] [0.067;0.149] [0.068;0.127]
	α 0	[-0.047;-0.234] [0.054;0.144]	[0.11;0.182] [0.11;0.144] [0.125;0.292] [0.127;0.228]
	α -1	[-0.080;0.022] [-0.035;0.012] [0.078;0.114] [0.062;0.104] [0.086;0.139] [0.079;0.117]
	Γ	[-0.199;0.007] [-0.108;0.013] [0.058;0.367] [0.029;0.368] [0.053;0.199] [0.053;0.115]

Table 2 :

 2 Table of log-likelihoods and BIC indexes for the different models. models have been validated in the same way as the model (M) (see Section 4).

	Model	Parameters Log-likelihood BIC
	(M 2 )	209	-24849	52040
	(M)	208	-24954	52238
	(M Λ )	186	-25399	52895
	(M Γ∼Gauss )	78	-29110	59082
	(M Γ∼Sinus )	78	-35615	72094

Acknowledgments

We thank Peter Thomson, Statistics Research Associates, for stimulating discussions and comments on the modeling and on identifiability problems. We thank the reviewers for their helpful comments.

order to avoid a numerical optimization in a high-dimensional space, separate optimizations in Λ and in Γ have been performed leading to a so-called Generalized EM algorithm (see the supplementary materials for more details). The reduced (M Λ ) and the full (M) models give similar results for the marginal distribution and the autocorrelation function. (M Λ ) leads also to an accurate description of the spatial structure of the data (see Figure 9). Again, lagged-one correlations are better reproduced by GMM parameters than by ML parameters. The model (M Λ ) is slightly inferior to the full model (M) in terms of BIC according to Table 2. Nevertheless, it clearly outperforms the models (M Γ ). It seems easier to find an appropriate reduced model for the loading matrix Λ than for the covariance matrix of the observation error Γ.

General discussion

Several multisite models, all based on Gaussian linear state-space models, are proposed to generate synthetic multivariate time series of wind speed. The

Supplementary materials

This file contains the proof of Proposition 1, a description of the Expectation-Maximization (EM) algorithms used to fit the models introduced in the paper together with a description of the Kalman recursions involved in the EMalgorithm (supp estimation.pdf).


Áttila L Rodrigues

Christophe Chatelain

Tânia Tomé

Mário J De Oliveira

Critical behavior in lattice models with two symmetric absorbing states

Keywords: numbers: 05.70.Ln, 05.50.+q, 64.60.Ht

We analyze nonequilibrium lattice models with up-down symmetry and two absorbing states by mean-field approximations and numerical simulations in two and three dimensions. The phase diagram displays three phases: paramagnetic, ferromagnetic and absorbing. The transition line between the first two phases belongs to the Ising universality class and between the last two, to the direct percolation universality class. The two lines meet at the point describing the voter model and the size ℓ of the ferromagnetic phase vanishes with the distance ε to the voter point as ℓ ∼ ε, with possible logarithm corrections in two dimensions.

Introduction

Studies of phase transitions and critical phenomena in equilibrium and nonequilibrium systems have shown that the critical behavior can be organized into universality classes, which are identified by a small number of features among them the symmetry. The most widely studied universality class, due to its experimental implications, is the Ising universality class which includes equilibrium as well as nonequilibrium systems [1,2]. The main feature of the Ising universality class is the up-down symmetry. Another class, that includes only nonequilibrium systems, is the direct percolation (DP) universality class [START_REF] Marro | Nonequilibrium Phase Transitions in Lattice models[END_REF][START_REF] Henkel | Non-Equilibrium Phase Transitions[END_REF][START_REF] Hinrichsen | [END_REF]. The main feature that distinguishes this class from the others is the presence of a single absorbing state. That is, if a system displays a continuous phase transition to a single absorbing state with no extra symmetry or conserving laws then its critical behavior belongs to the DP universality class [6,7].

Our purpose here is to analyze models that possess features pertaining to the two universality classes. Here, we focus on models belonging to the Ising and DP universality classes, that is, systems with up-down symmetry and an absorbing state [8,9,10,11,12]. Due to the up-down symmetry, these systems, in fact, have two symmetric single absorbing states. Lattice models of this type have been studied numerically [8] revealing the existence of two transition points when a parameter is varied: a symmetry breaking point belonging to the Ising universality class and another belonging to the DP universality class. Another numerical study in two dimensions by the method of Langevin equation [9], with two parameters, has shown that the two critical points evolve into two transition lines that meet at a certain point, given rise to a third transition line. This study confirm the Ising and DP character of the two lines and the existence of a point on the third line corresponding to the voter model [2], which is not the meeting point. It has been claimed that the symmetry breaking transition is not of the Ising type [11], but Monte Carlo simulations on a two-dimensional interacting monomers model [12] confirmed the Ising transition. In the present study, we wish to clear up this point and analyze the role played by the voter model. This is achieved by studying lattice models with up-down symmetry and two absorbing states in two and three dimensions by numerical simulations and mean field approach. Our results confirmed the existence of the Ising and DP lines and revealed that two lines meet at the voter point and that this is the only point at the third line which is critical.

The phase diagram in the two-dimensional space of parameters displays three phases: paramagnetic (P), ferromagnetic (F) and absorbing (A), as shown in figure 1a. The PF line that separates the P and F phases is a critical line belonging to the Ising universality class. The FA line that separates the A phase and the F phase, understood as an active phase, is a critical line belonging to the DP universality class. The two lines, PF and FA, meet at the point corresponding to the voter model. As to the PA line, that separates the P and A phases, it is a discontinuous transition line that ends at the voter point. Along the PA line the jump in the order parameter is a constant and equals the jump observed in the voter model. For an appropriate description of the critical properties around the voter point, it is convenient to introduce two parameters, ε and r, that define a reference frame centered with origin at the voter point. The ε-axis and r-axis are parallel and perpendicular to the direction defined by the discontinuous transition line, as shown in figure 1a.

The order parameter m of the F phase, shown in figure 1b is assumed to obey the following scaling relation, around the voter point,

m = ψ ± (εr -1/φ ), ( 1
)
where φ is the crossover exponent and ψ + and ψ -are valid for ε > 0 and ε < 0, respectively. The universal function ψ -(x) is a step function. The universal function ψ + (x) vanishes according to ψ + (x) ∼ (a -x) β I , where β I is the Ising order parameter exponent and a is a positive constant so that, near the Ising critical line, the order Ising parameter m behaves as

m = A(ε I -ε) β I , (2)
where ε I = ar 1/φ and the amplitude A diverges as A ∼ r -β I /φ . The universal function ψ + (x) approaches its maximum value 1 according to 1 -ψ + ∼ (x + b) β DP where β DP is the DP order parameter exponent and b is a constant so that, near the DP critical line, the DP order parameter ρ = 1 -m behaves as

ρ = B(ε + ε DP ) β DP , (3)
where ε DP = br φ and the amplitude B diverges as B ∼ r -β DP /φ . The Ising and DP critical lines occur when m → 0 and ρ → 0, respectively, and are given by ε = ar 1/φ and ε = -br 1/φ so that the size of the F phase behave near the voter point as

ℓ ∼ ε φ . ( 4
) Table 1. Probability of transition σ i → -σ i . σ i -6 -4 -2 0 +2 +4 +6 +1 p 0 p 1 p 2 1/2 1 -p 2 1 -p 1 1 -p 0 -1 1 -p 0 1 -p 1 1 -p 2 1/2 p 2 p 1 p 0

Model

We consider stochastic variables that take two values, σ i = ±1, called spin variables.

They are located on the sites of a regular lattice which is either a two-dimensional triangular lattice or a three-dimensional cubic lattice. Notice that both lattices have the same coordination number, which is six. The system evolves in time according to a continuous time Markov process with a single-site change. The single-site transition rate w i (σ) from σ = (σ 1 , . . . , σ i , . . . , σ N ) to σ i = (σ 1 , . . . , -σ i , . . . , σ N ) is set up in such way as to hold the properties: a) it has the up-down symmetry, that is, it is invariant under the transformation σ → -σ; b) it depends on the neighboring spin variables only through the sum of these spin variables. A transition rate that fulfils these two properties is given by

w i (σ) = 1 2 {1 -σ i f (s i )}, s i = δ σ i+δ , (5)
where the sum in δ extends over the six nearest neighbor sites, and f (s) is an odd function of s.

The most general form of f (s) has three parameters. Let us denote by σ j , j = 1, 2, 3, 4, 5, 6 the spins of the six nearest neighbor sites of a central site j = 0, so that

s 0 = σ 1 + σ 2 + σ 3 + σ 4 + σ 5 + σ 6 . (6)
Taking into account that σ j = ±1, then s 0 may take only the values 0, ±2, ±4, ±6 from which follows that the most general form of an odd function f (s) can be written as

f (s) = a 2 s + b 8 s 3 + c 32 s 5 , (7)
where a, b and c are parameters. Here, it is convenient to define new parameters

p 0 = [1 + f (6)]/2, p 1 = [1 + f (4)]/2 and p 2 = [1 + f (2)]/2, that is, p 0 = 1 2 (1 + 3a + 27b + 243c), (8)
p 1 = 1 2 (1 + 2a + 8b + 32c), (9)
p 2 = 1 2 (1 + a + b + c). (10)
In terms of these parameters, the transition probabilities are as shown in table 1.

Two well known models are special cases. When p 0 = 1, p 1 = 5/6 and p 2 = 2/3, (a = 2/3, b = c = 0) we recover the so called voter model [2,9] whose transition rate is given by

w voter i (σ) = 1 2 {1 - 1 3 σ i s i }, s i = δ σ i+δ . (11)
The Glauber-Ising model, whose transition rate is given by

w GI i (σ) = 1 2 {1 -σ i tanh(Ks i )}, s i = δ σ i+δ , (12)
is recovered when the parameter are connected by

p 0 = 1 2 (1 + tanh 6K), (13)
p 1 = 1 2 (1 + tanh 4K), (14)
p 2 = 1 2 (1 + tanh 2K). ( 15
)
In this case, detailed balance is obeyed and the stationary distribution is the Gibbs distribution associated to the Ising Hamiltonian. In other cases there is no detailed balance as can be seen by the following sequence of transitions. Consider two nearest neighbor sites initially in states (++). All the other nearest neighbors of this pair of spins are positive except one of them which is negative. Next, they change their states according to the cyclic sequence (++) → (+-) → (--) → (-+) → (++). Using the rules of table 1 we find that the ratio between the transition probabilities of this path and of its reverse equals (1 -p 1 ) 2 p 2 p 0 /(1 -p 0 )(1 -p 2 )p 1 p 1 , which in general is distinct from the unit. Thus, there is no microscopic reversibility and no detailed balance.

To meet the two features of the Ising and DP universality class, that is, up-down symmetry and two equivalent single absorbing states, the transition rate should lead to an absorbing state. This is accomplished by imposing the restriction p 0 = 1 because in so doing the transition rate will vanish whenever all spins are up or all spins are down, which are the two states identified as the equivalent absorbing states. In this case, the transition rates do not obey detailed balance, which amounts to say that the two features mentioned above are only possible in systems that in the stationary states are out of equilibrium. When p 0 = 1, the phase diagram is restricted to the plane (p 1 , p 2 ), and the voter model corresponds to the point (5/6, 2/3) of this diagram. As we shall see, the phase diagram displays a paramagnetic phase, a ferromagnetic phase and an absorbing state, as explained before.

Mean-field approximation

A qualitative phase diagram is obtained by the use of a mean-field approximation. To this end, we start by writing the function f (s) in the form where the first summation runs over i from 1 to 6, and has 6 terms; the second summation runs over i, j, k with the restriction 1 ≤ i < j < k ≤ 6 and has 20 terms; and the third summation runs over i, j, k, ℓ, n with the restriction 1 ≤ i < j < k < ℓ < n ≤ 6 and has 6 terms. The parameters A, B, and C are related to p 0 , p 1 , and p 2 by

f (s) = A 6 i σ i + B 20 i,j,k σ i σ j σ k + C 6 i,j,k,ℓ,n σ i σ j σ k σ ℓ σ n (16)
p 0 = 1 2 (1 + A + B + C) (17)
p 1 = 1 2 (1 + 2 3 A - 2 3 C) (18)
p 2 = 1 2 (1 + 1 3 A - 1 5 B + 1 3 C) (19)
The time evolution of the magnetization σ 0 is given by

d dt σ 0 = -2 σ 0 w 0 (σ) = -σ 0 + f (s 0 ) (20) or d dt σ 0 = -σ 0 + A σ i + B σ i σ j σ k + C σ i σ j σ k σ ℓ σ n (21)
where translation invariance has been used. Next, we use the simplest mean-field approximation, for which σ i σ j σ k = m 3 and σ i σ j σ k σ ℓ σ n = m 5 , where

σ i = m, to write dm dt = -(1 -A)m + Bm 3 + Cm 5 . (22)
Here, we are interested in the case p 0 = 1 which is equivalent to 1

-A = B + C so that dm dt = -(B + C)m + Bm 3 + Cm 5 , (23)
which is identical to the deterministic part of the Langevin equation considered in reference [9]. In the stationary state, the possible solutions are as follows. a) m = 0, corresponding to the paramagnetic (P) phase; it is stable as long as B + C > 0, which is equivalent to 12p 1 + 15p 2 < 20. b) m = ±1, corresponding the absorbing (A) phase; it is stable as long as B + 2C < 0, which is equivalent to p 1 > 5/6. c) m = 0 and m = ±1, corresponding to the ferromagnetic (F) phase. The phase diagram in the plane (p 1 , p 2 ) is shown in figure 2. The P-F line is given by B + C = 0, or 12p 1 + 15p 2 = 20, and the F-A line is given by B + 2C = 0, or p 1 = 5/6. The point corresponding to the voter model is located at the point where the two lines meet, that is, at B = C = 0, or p 1 = 5/6 and p 2 = 2/3. The P-F and F-A transitions are continuous whereas the P-A transition is discontinuous. The discontinuous P-A transition line shown in figure 2 was obtained as follows. We start by writing equation ( 23) in the form dm/dt = -∂f /∂m where f = (B + C) m 2 /2 -Bm 4 /4 -Cm 6 /6 + K. Interpreting f (m) as a free energy then the discontinuous transition occurs when f (0) = f (±1), which gives 3B + 4C = 0 or 5p 2 + 8p 1 = 10.

Inside the F phase, the magnetization m is given by

m = 15p 2 + 12p 1 -20 15p 2 -12p 1 1/2 . ( 24
)
Near the critical line PF, p 1 = p FP = (20 -15p 2 )/12, we may write

m = A(p 1 -p FP ) 1/2 , A = 5 2 p 2 - 2 3 -1/2
(25) so that the amplitude A diverges with an exponent 1/2 as one approaches the voter point. Analogously, near the FA transition, p 1 = p FA = 5/6, the order parameter ρ behaves as

ρ = B(p FA -p 1 ), B = 5 4 p 2 - 2 3
-1

(26) so that the amplitude B diverges with an exponent 1 as one approaches the voter point.

It is worth to notice that the size ℓ = p FA -p PA of the ferromagnetic phase, given by

ℓ = 5 4 p 2 - 2 3 ( 27
)
vanishes linearly as one approaches the voter point.

Simulations

The results obtained from numerical simulations on the cubic and triangular lattice allow us to conclude that the topology of the phase diagram is the same as the one obtained from the mean-field approximation. The numerical simulations on these lattices were carried out as follows. At each time step a site of the lattice is chosen at random and updated according to the transition probabilities in table 1, restricted to the case p 0 = 1. The averages of several quantities were obtained in the stationary state. Defining

m = |s| , (28) the susceptibility,
χ = N{ s 2 -|s| 2 }, (29)
the Binder cumulant

U = 1 - s 4 3 s 2 2 . ( 30
)
Examples of these quantities, for the cubic lattice with N = L 3 sites, calculated along the line p 2 = 1, are shown in figures 3 and 4. At the P-F transition point, the Binder cumulant U is a constant, a result that allows to located the P-F transition line. The critical behavior of m and χ allows to say that, with respect to the ferromagneticparamagnetic transition, the model belongs to the universality class of the Ising model.

To analyze the ferromagnetic-absorbing transition, where the ferromagnetic phase is identified as the active phase, we define the quantity ρ = 1 -m, which is the order parameter that characterizes the transition from active to absorbing. In figure 5, we show ρ versus p 1 along the line p 2 = 1. The location of the transition is obtained by the extrapolation ρ → 0, what allows us to draw the F-A line. The critical behavior shows that with respect to this transition, the model belongs to direct percolation universality class.

Using the data of the P-F and F-A transitions obtained for several values of p 2 we draw the phase diagram of figure 6a. As expected the two lines, P-F and F-A, meet at the voter point located at p 1 = 5/6 and p 2 = 2/3. In figure 6b we show the size ∆p 1 = p DP -p I of the F phase as a function of the distance from the voter point ∆p 2 = p 2 -2/3. As can be seen in this figure, the data points has a finite slope at the origin which amounts do say that the scaling relation

∆p 1 ∼ (∆p 2 ) φ , (31)
is fulfiled with φ = 1.0. We have also determined the nature of the phase transition between the P and A phases. As we cross the transition line the order parameter jumps from zero, which is the value of m within the paramagnetic phase, to the maximum value m = 1, which is the value of m inside the absorbing phase. The susceptibility χ increases as one approaches the transition line from the paramagnetic phase and reaches a finite value χ 0 (∆p 2 ) which depends on the distance ∆p 2 = p 2 -2/3 from the voter point. As one approaches the voter point χ 0 increases without limits and diverges with an exponent γ = 1. Notice that, according to our numerical results for the susceptibility obtained along p 2 = 2/3, the susceptibility χ diverges as one approaches the voter point from the paramagnetic phase with the same exponent γ = 1.

We have also studied the model defined on a triangular lattice by using the same methods. As can be seen in figure 7, the phase diagram also displays the three phases found in the cubic lattice. However, the ferromagnetic phase is very narrow, becoming very difficult to extract the crossover exponent but the two lines, PF and FA, seems to meet tangentially in consistency with an exponent φ = 1 with logarithm corrections as we will argue below.

Discussion

We have studied the phase diagram of a system belonging to the Ising and DP universality classes. The phase diagram in the two parameters has an Ising critical line separating the paramagnetic and the ferromagnetic phase, a DP critical line separating the ferromagnetic phase and the absorbing phase, and a discontinuous phase transition line separating the paramagnetic and the absorbing phase. The three lines meet at the voter point in such as way that the size of the ferromagnetic phase vanishes as one approaches the voter point with a crossover exponent φ. The crossover exponent should be understood as the ratio between the two exponents related to the two relevant scaling fields ε and r around the voter point. Above the critical dimension d c of the voter model, these two exponents do not depend on dimension so that the same should happen to the crossover exponent. Taking into account that the upper critical dimension of the voter critical point is d c = 2, we expect the crossover exponent φ to have the same value for d ≥ 2 with possible logarithm corrections when d = 2. The results coming from mean-field theory gives φ = 1 which allows to conclude that φ = 1 for d ≥ 2, a result confirmed by our calculations in d = 3. As to the two-dimensional case, the size of the ferromagnetic phase is too narrow for a precise numerical calculation but is consistent with the result φ = 1 with logarithm corrections.

The results for the susceptibility around the Ising and DP lines are in accordance with these two universality classes. Concerning the discontinuous transition from the paramagnetic to the absorbing transition, we found that the susceptibility is finite at the transition line but diverges as one approaches the voter point.

It would be very interesting if the picture that we have drawn concerning the critical behavior of systems belonging to the Ising and DP universality classes could be confirmed by a renormalization group calculation by means of a Langevin equation formulation.

Figure 1 .

 1 Figure 1. (a) Phase diagram in the field parameters ε and r. The full circle, located at the origin, corresponds to the voter model (V). The phases are: paramagnetic (P), ferromagnetic (F) and absorbing (A). The solid lines are continuous phase transitions belonging to the Ising and DP universality classes and the dashed line is a discontinuous phase transition. (b) Magnetization m as a function of ε and r.

Figure 2 .

 2 Figure 2. Phase diagram in the plane (p 1 , p 2 ) according to mean-field approximation. The full circle, located at (p 1 , p 2 ) = (5/6, 2/3), corresponds to the voter model. The phases are: paramagnetic (P), ferromagnetic (F) and absorbing (A). The solid lines are continuous phase transitions and the dashed line is a discontinuous phase transition.

Figure 3 .

 3 Figure 3. (a) Magnetization m versus p 1 along p 2 = 1. The transition occurs at p I = 0.5930 as can be inferred from the plot (b) of m 1/β versus p 1 , where β = 0.326.

Figure 4 .

 4 Figure 4. (a) Susceptibility χ and (b) Binder cumulant U versus p 1 along p 2 = 1. The transition occurs at p I = 0.5933.

Figure 5 .Figure 6 .

 56 Figure 5. (a) Order parameter ρ of the F-A transition as a function of p 1 for p 2 = 1, The critical point occurs at p DP = 0.7670 as can inferred from the plot (b) where β = 0.81

Figure 7 .

 7 Figure 7. Phase diagram in the plane (p 1 , p 2 ) according to simulations on a triangular lattice. The full circle, located at (p 1 , p 2 ) = (5/6, 2/3), corresponds to the voter model. The phases are: paramagnetic (P), ferromagnetic and absorbing (A). The inset shows an enlargement around showing the ferromagnetic phase between the P-F and F-A transition lines.