Bayesian trend analysis of extreme wind using observed and hindcast series off the Catalan coast, NW Mediterranean Sea

It has been suggested that climate change might modify the occurrence rate and magnitude of large oceanwave and wind storms. The hypothesised reason is the increase of available energy in the atmosphere–ocean system. Forecasting models are commonly used to assess these effects, given that good-quality data series are often too short. However, forecasting systems are often tuned to reproduce the average behaviour, and there are concerns on their relevance for extremal regimes. We present a methodology of simultaneous analysis of observed and hindcast data with the aim of extracting potential time drifts as well as systematic regime discrepancies between the two data sources. The method is based on the peak-over-threshold (POT) approach and the generalized Pareto distribution (GPD) within a Bayesian estimation framework. In this context, storm events are considered points in time, and modelled as a Poisson process. Storm magnitude over a reference threshold is modelled with a GPD, a flexible model that captures the tail behaviour of the magnitude distribution. All model parameters, i.e. shape and location of the magnitude GPD and the Poisson occurrence rate, are affected by a trend in time. Moreover, a systematic difference between parameters of hindcast and observed series is considered. Finally, the posterior joint distribution of all these trend parameters is studied using a conventional Gibbs sampler. This method is applied to compare hindcast and observed series of average wind speed at a deep buoy location off the Catalan coast (NE Spain, western Mediterranean; buoy data from 2001; REMO wind hindcasting from 1958 on). Appropriate scale and domain of attraction are discussed, and the reliability of trends in time is addressed.


Introduction
Interest on natural hazard prevention, prediction and mitigation has increased along the last decades: strong wind storms, extreme wind gusts, hurricanes and tornados are not an exception.The dangers, that a climate change might induce, add an additional challenge to the statistical analysis of extreme wind data as possible trends in extremal winds might occur, increasing the inherent difficulties of extremal analysis.Performing extremal analysis requires long records spanning decades, even centuries, to characterise the rarest events.This need is exacerbated when the data series are potentially affected by a trend in time.Indeed, whichever model is fitted, the uncertainty of the estimates is large when only data series with few events are available.
Bayesian methods were introduced in the field of natural hazards two decades ago and they have succeeded as a flexible and consistent way of controlling uncertainty (Sánchez-Arcilla et al., 2008a;Coles and Tawn, 1996;Gelman et al., 1995).Simultaneously, the models for analysing extremal events have evolved.The models known as excesses over threshold or peak over threshold (POT) are nowadays in common use for natural hazards and, particularly, in wind hazard modelling (Walshaw, 1994;Palutikof et al., 1999;Della-Marta et al., 2009;Coles, 2001).However, trend analyses require non-stationary models while the standard POT models are based on the assumption of stationarity of the time process.Non-stationary POT models are available, with applications in different frameworks (Beguería et al., 2011;Hundecha et al., 2008;Tramblay et al., 2011).

M. I. Ortego et al.: Bayesian trends of extreme wind
An obvious way of reducing the uncertainty of the estimates of both extremal parameters and their trends consists in analysing data spanning a longer time period.As direct observations cannot be extended beyond the availability of measuring devices at the site of interest, hindcast, forecast or satellite data series may offer an alternative source of information (Winterfeldt et al., 2010;Izaguirre et al., 2010).However, it is known that hindcast wind data seldom fit exactly the directly observed wind (Bolaños et al., 2004).There are several reasons for the misfit: possible mis-calibration of models; the fact that models actually aim at giving a kind of energy average, not quasi-instantaneous (e.g. 10 min average) winds; etc.However, when interest is focused in climatic features as extreme event statistics, hindcast data information is very valuable.Nevertheless, the joint use of directly observed and hindcast data introduces further complication in the models to be used, as a kind of conciliation between the two types of data must be considered.Some authors have assessed this compatibility between types of data in order to assess extremal characteristics.Different approaches have been used, combining information from different buoy networks and hindcast models in different areas: the North Atlantic (Forte et al., 2012;Winterfeldt et al., 2010) and the Mediterranean (Cañellas et al., 2007;Izaguirre et al., 2010), among others.
This contribution aims at the joint analysis of nonoverlapping hindcast wind data and data coming from direct measurements.This analysis evaluates the differences between both types of observations and checks for (linear) trends in the extremal parameters.The procedure is based on a non-stationary POT model which is analysed using Bayesian techniques.The illustration series in this contribution shows no overlap between model hindcast and measured data, and therefore the proposed model does not account for it.An analogous model for overlapping series was presented by Ortego et al. (2012), where it was shown that the results are not significantly altered when accounting for overlap, albeit for significant wave height data.
Section 2 describes the data and the site briefly.Section 3 presents the model and the estimation methods.Finally, results are presented in Sect. 4.

Wind data
The statistical properties of a hindcast model and a series of observations are compared.A series of wind speed has been constructed, joining actual data measured at a deep buoy location (10 min average) (Bolaños et al., 2009) and hindcast data provided by a REMO model within the HIPOCAS project (hourly average at 10 m height) (Sotillo et al., 2005;Guedes Soares et al., 2002).Wind speeds at the Tarragona deep buoy (XIOM network, longitude 40.68 • N, latitude 1.47 • E) are intermittently available between the years 2004 and 2012, totalling 7 years and 158 days with observations.The series of wind speed at the nearest HIPOCAS grid point (longitude 40.50 • N, latitude 1.50 • E) is available for the period 3 January 1958 to 31 December 2001.Both locations are shown in Fig. 1.Admitting that the statistical parameters of the model may evolve over time, it is possible to assess evidence of differences between series (h for hindcast REMO and b for buoy) as well as evidence of linear trends for the parameters.
The criteria to define a wind event have been chosen to meet three hypotheses which are used in the proposed methodology.The first and second ones assume that wind events are conditionally independent both in time occurrence and in size, i.e. once the parameters of the model are fixed the occurrence and size are independent from event to event.These two assumptions are difficult to test using limited time series, but obvious autocorrelations should be avoided.The third assumption is that events are punctual in time.This is not true as wind events can be characterised by a duration over a threshold and a definition of punctual event is then required.Events have been defined in the following manner: an event starts when the recorded average wind speed in the reference time series is greater than 15 m s −1 .It finishes when the recorded average wind speed is less than 15 m s −1 and remains at this level for at least 3 days.Therefore, the minimum gap between events is 3 days and we ensure that events are approximately independent.It is debatable whether 3 days separation between events is enough for approximate independence from event to event.For instance, the so-called twin storms observed on the east coast of the Iberian Peninsula (Jiménez et al., 1997) are candidates to violate the assumed independence.However, recent studies on these phenomena are not conclusive in this respect (Sánchez-Arcilla et al., 2008b;Wojtanowicz, 2010).The magnitude of the event is defined as the maximum average wind speed recorded during the event, and the corresponding time of occurrence as the instant of recording of this value.The h 0 = 15 m s −1 threshold ensures that the excesses have approximately a generalized Pareto distribution (GPD) distribution.As wind speed measurements may have a ratio (relative) scale, the studied magnitude is log wind speed (Egozcue et al., 2006).The data set is displayed in Fig. 2. Gaps in the buoy segment represent a loss of information but do not represent an inconvenience for the methodology.However, we are not considering the distribution of events within a year and gaps covering only one season may cause small increments/decrements of the number of events or a distortion of their magnitude.

The peak over threshold model
Models of excesses over a threshold, often named peak over threshold (POT; Embrechts et al., 1997) extensively used in hazard analysis of natural phenomena (Davison and Smith, 1990;Egozcue and Ramis, 2001;Coles et al., 2003;Egozcue et al., 2006;Leadbetter, 1991).Independently of the specific estimation method, the POT framework consists of a system of modelling assumptions.The standard assumptions for stationary process (Embrechts et al., 1997) can be summarised as follows: -events are identified as points occurring in time; -the elapsed times between consecutive events are random, identically distributed and mutually independent; -each event has an associated random magnitude; -excesses of magnitudes over a given threshold are identically distributed and mutually independent; -elapsed times and magnitude excesses are mutually independent.
In this study of wind events we adopt a slightly modified POT model in order to cope with non-stationarity.The assumption of identically distributed elapsed times and magnitude excesses is changed: both random variables have distributions in a family which parameters can change with time and with the data source, buoy or hindcast.Also the statements of independence are modified to conditional independence -i.e. for fixed parameters of the distribution of elapsed times and magnitude excesses, random observations are independent.
In the stationary case the elapsed times are usually assumed to be exponentially distributed, thus the number of excesses over the selected threshold of magnitude is an homogeneous Poisson process.Also, a common approach is to assume that excesses over the threshold have a generalized Pareto distribution (GPD), as proposed in Davison and Smith (1990) or in Embrechts et al. (1997).In the present approach, the occurrence of events with magnitude over the threshold is assumed to be an inhomogeneous Poisson process.Excesses over the threshold are modelled by a GPD, although considering that its parameters might depend on time and data source.

Occurrence of wind events
Wind events were defined in Sect. 2 and their wind speed is assumed to be larger than the threshold 15 m s −1 .The main assumption is that their occurrence in time corresponds to an inhomogeneous Poisson process (Grandell, 1997;Coles, 2001).In a homogeneous Poisson process the parameter is the Poisson rate, interpreted as the mean number of events per year.When this mean is not constant, the Poisson rate is replaced by the Poisson intensity λ(t).If T is the random time from an origin to the next event, and (t, t + dt) denotes a short enough time interval, then λ(t) is for t > 0, where F T is the cumulative distribution function (CDF) of T .Therefore, λ(t)dt is the probability of the event occurring in (t, t + dt], conditional to the event not occurring before t.As a consequence of Eq. ( 1), the CDF and probability density function (PDF) of T , for t > 0, are respectively.The Poisson intensity is modelled as depending on time in two ways.First, we assume a linear trend on time of the Poisson rate, in order to discuss evidence in favour or against the presence of a trend in time.On the other hand, wind events come from two different sources, REMO hindcasts and observations on a buoy.The Poisson intensity corresponding to the two sources might differ.This possibility is modelled as a systematic difference, which can be viewed as a step function on time or as an indicator function of the data source.Suppose that wind events detected by REMO took place at times t 0 , t 1 , . . ., t n , and those observed by the buoy at times τ 0 , τ 1 , . . ., τ m , with t n < τ 0 .The available data are the elapsed times between consecutive events, i.e. t i −t i−1 , i = 1, 2, . . ., n (REMO events) and τ i − τ i−1 , i = 1, 2, . . ., m (buoy events).Then, the whole time interval covered by observations is (t 0 , τ m ).The Poisson intensity can be written as where λ h is the Poisson intensity at the first event at t 0 which corresponds to a REMO observation and δ λ is the increment of Poisson intensity due to the fact that the observation comes from the buoy.The symbol I b is an indicator with value equal to 1 for event times in which events are recorded by the buoy, and 0 otherwise.The parameter α λ is the total increase of the Poisson intensity from t 0 to τ m due to the linear trend of the Poisson intensity irrespective to the type of observation.
In order to proceed to a Bayesian estimation of the parameters the likelihood function of (λ h , α λ , δ λ ) given the elapsed times is required.Since the elapsed times are assumed independent, given the parameters, the required likelihood is where, for simplicity, {t i }, {τ j } represents all available data on elapsed times between consecutive events.When interruptions of observations are present, as is the case in the buoy data set, some time differences τ j − τ j −1 do not correspond to elapsed times between consecutive events; these time intervals are ignored in the likelihood (see further details in Ortego et al., 2012), i.e. the likelihood (Eq.2) is valid when there are no interruptions in the observations.After grouping integrals within the exponentials, the log-likelihood The log-likelihood in Eq. ( 3) is essentially that proposed in Ortego et al. (2012) for ocean waves, although a bit simpler due to the fact that in the present data set there is no overlapping in time between the hindcast and buoy observations.Once the model for occurrences in time is specified, its capabilities and limitations can be discussed.The main ability of the model is to capture larger (smaller) concentrations of wind events along the observation period.As the change of Poisson intensity is considered linear, what is controlled is only the long-term change of the mean number of events per year.Annual or few-year periodic changes are not represented by the model.Moreover, occurrence of events is assumed to be independent of the size of the events in the inhomogeneous Poisson process.As a consequence, possible changes in the Poisson intensity λ(t) along the observation time do not also tell us about possible changes in the size of events.For instance, an increasing linear trend of λ(t) would imply a larger mean number of events per year, but, in the model, this is not related to the possibility of changes of the size of the events.

Event magnitude model
The wind magnitude associated with each event can be selected in different ways.Traditionally, wind speed magnitude is taken as velocity in m s −1 without any additional consideration.However, the natural scale and domain of wind speed should be taken into account.Large wind speeds do rather behave in a ratio scale (also known as relative scale), as shown by the thresholds chosen for most cyclone classification systems and the Beaufort scale levels of 7 or more.In fact, the absolute scale is near to be meaningless: a difference of 1 m s −1 on a reference wind speed of 2 m s −1 represents a factor 3/2, whereas the same difference on a reference wind speed of 20 m s −1 is considered almost irrelevant.A simple way of considering these issues of scale is to take logarithms on wind speeds, and thus consider that the magnitude associated with events is the logarithm of the maximum average wind speed.Accordingly, the magnitude of a wind event is herein taken as the natural logarithm of the measured wind speed in m s −1 .
As we consider that events occur following a nonhomogeneous Poisson process, magnitudes are assumed to be conditionally independent from event to event.In order to model magnitude excesses over a threshold of h 0 = log(15 m s −1 ), we assume that they follow a generalized Pareto distribution (GPD).The GPD provides a suitable asymptotic model for excesses over a high enough threshold (Pickands III, 1975;Davison and Smith, 1990;Dupuis, 1998).Furthermore, we also consider that excesses over h 0 must be GPD in the Weibull domain of attraction, i.e. the support of the excesses has an unknown finite upper limit.This assumption is based on the fact that, on Planet Earth, wind speed is physically limited and, accordingly, the existence of such upper bound is granted, even though the limit itself is not known.This assumption has been successfully applied in the extremal analysis of other weather magnitudes, such as rainfall and ocean-wave heights (Pawlowsky-Glahn et al., 2005;Egozcue et al., 2005Egozcue et al., , 2006;;Sánchez-Arcilla et al., 2008a).
For a magnitude X, denote the excess by Y = X − h 0 conditional on X > h 0 .One of the standard parameterisations of the GPD is (Embrechts et al., 1997) where β > 0 is a scale parameter, ξ is a shape parameter, and y sup is the upper limit of the distribution support.When ξ ≥ 0, the GPD has an unlimited tail, i.e. y sup = +∞, and belongs to the Fréchet domain of attraction of maxima (ξ > 0); for ξ = 0, the GPD is in the Gumbel domain of attraction.
Under the assumption that the GPD has an upper bound, the shape parameter must be ξ < 0 (Weibull domain of attraction).In order to take into account the restriction ξ < 0, other parameterisations of GPD in the Weibull domain of attraction have been proposed (Ortego et al., 2010(Ortego et al., , 2012)); we adopt here the parameterisation in Ortego et al. (2012).We use instead of the classical parameters ξ , β, given that the GPD distribution is in the Weibull domain of attraction (ξ < 0).The old parameters β and ξ < 0 can be obtained from the new ones as −ξ = exp(ν) and β = exp(ν)/ exp(µ).Introducing these parameters in Eq. ( 4), the CDF and the corresponding PDF are for 0 ≤ y < y sup , respectively.We have introduced two assumptions -namely that the magnitude to be treated is the log-wind speed and that the excesses over h 0 must have a limited distribution in the Weibull domain of attraction.The compatibility of both assumptions can be checked on the available data. Figure 3 shows the likelihood contours of raw and log-transformed wind speed for both data series.
In both data sources, the likelihood corresponding to the raw data cover a substantial region with ξ > 0 (Fréchet domain of attraction), whereas for log-transformed data the coverage is considerably shifted towards values of ξ < 0 (Weibull domain of attraction).Table 1 shows the probabilities of the Weibull domain of attraction for raw and logtransformed wind speed at the location for both data series.
Buoy data are likely to correspond to a GPD with ξ < 0 both in the raw and log scales, but this is not the case for hindcast data.In a raw scale, the Fréchet domain of attraction is more likely for this second data series.These results suggest that after taking logarithms both series are likely to correspond to the Weibull domain of attraction, a further argument in favour of considering a ratio scale for large wind speeds.Figures 4 and 5 show the Q-Q plots of the HIPOCAS/buoy data using the GPD with median values of ξ , β as reference.This approach is quite restrictive because the GPD parameters are assumed to be uncertain and possibly changing in time.The figures show that GPD is a suitable model for both data series.
The upper bound parameter µ = log y sup is so large that it only depends on universal physical laws and geometric aspects of the Earth.We consider that human-scale climatic q q q q qq q qq q q q q qq q q q qq q qq q q q q qq q q q q q qq qq q q q q q q qq q q q q q qq qq q q q q qqqq q q q qqq q q q qq q q q qq q qq qqq qqq q q q qq q q qq qq q qq q qq q q q q q q q q q q q q 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.0 0.5 1.0 1.5

reference GPD (param given)
observed quantile reference quantile q q q q qq q qq q q q q qq q q q qq q qq q q q q qq q q q q q qq qq q q q q q q qq q q q q q qq qq q q q q qqqq q q q qqq q q q qq q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Figure 4. Q-Q plot of wind REMO data using the GPD with median values of ξ and β as reference.This approach is quite restrictive because the GPD parameters are assumed to be uncertain and possibly changing in time.GPD is a suitable model for this data series.Green lines are 0.05-significance contours for Kolmogorov-Smirnov test of goodness of fit.
changes cannot change it, and it is therefore constant over time (Ortego et al., 2012).We also assume that µ may be different for hindcast data and data observed at the buoy.The proposed model for µ is then where the difference δ µ represents systematic regime differences between buoy observations and the hindcast data.
On the other hand, the GPD parameter ν might be affected by a (linear) trend in time, apart from the possible systematic regime differences between the buoy and the REMO series.Accordingly, the model proposed is where δ ν is the difference in ν between the series h (hindcast REMO) and b (buoy).The parameter α ν is the total drift of ν from t 0 to t N , from the start of the hindcast series to the end of buoy observations.The excess magnitude Y has been recorded for each event, for both series h and b.The data set of pairs of occurrence times and excesses is {(t i , y i ), i = 1, . . ., N}.It is not necessary to consider the distinction of occurrence times of series h and b.Times are generically denoted as t i and the notation in Sect.3.2 is no longer necessary: N = n + m and the τ i are a subset of the t i in the present section.The likelihood of the q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.05 0.10 0.15 0.20 0.25 0.0 0.1 0.2 0.3

reference GPD (param given)
observed quantile reference quantile q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Figure 5. Q-Q plot of buoy data using the GPD with median values of ξ and β as reference.This approach is quite restrictive because the GPD parameters are assumed to be uncertain and possibly changing in time.GPD is a suitable model for this data series.Green lines are 0.05-significance contours for Kolmogorov-Smirnov test of goodness of fit.
model parameters, given the sample, is The corresponding log-likelihood, for 0 < y i < y sup , is which will be used in the Bayesian estimation of the parameters.

Bayesian estimation
Bayesian methods (e.g.Gelman et al., 1995) are very useful in contexts where data are scarce, such as hazard assessment problems.Bayesian methods allow us to assess the uncertainty of the estimates of parameters, usually large in these situations.Parameters are assumed to be random variables, and uncertainty is described through their distribution.
The likelihood of the parameters given the data, φ, θ, factorises into two terms depending only on φ and θ respectively: L(θ|D) × L(φ|D), where D denotes the sample of occurrence times and excess magnitudes for the two series h and b.
The posterior distribution of the parameters is proportional to the product of the prior and the likelihood: The likelihood of the parameters is concentrated in a finite range (mainly because of the assumption of a GPD model with an upper bound).An improper uniform prior has been assumed for both sets of parameters, according to this finiterange concentration feature (Box and Tiao, 1973, p. 21).The prior density for the occurrence parameters, π 0 (φ) has been assumed to be (improper) uniformly distributed on λ h > 0, The prior density for the magnitude parameters, π 0 (θ), has been assumed to be uniformly distributed on µ 0 , δ µ , ν 0 , δ ν , α ν .The whole shape of the posterior density gives an assessment of the uncertainty in the estimation procedure.The posterior π(φ, θ) may also be used to derive a point estimate, e.g. the most likely value (posterior mode) or the expected value.
The posterior in Eq. ( 11) has a quite complex form when considering all parameters simultaneously.A first simplification comes from the factorisation of the posterior density into terms containing the occurrence parameters φ and the excesses parameters θ.Therefore, the estimation of φ and θ can be carried out separately.For both terms in the posterior density, fixing all but one of the involved parameters, the conditional log-posterior density becomes a tractable univariate log-density, which can be satisfactorily sampled using a Gibbs sampler (Robert and Casella, 2000).

Results and discussion
The wind data set described in Sect. 2 has been modelled using the POT-GPD framework defined in Sects.3.2 and 3.3.A sample of the posterior density in Eq. ( 11) has been obtained using a Gibbs sampling algorithm, with three chains, 10 000 draws and a thinning ratio of 1 : 10.A burn-in of 50 % has been applied.The convergence of the joint chain has been assessed using the Gelman criterion (Gelman et al., 1995).
The log wind speed magnitude (excesses over log 15 m s −1 ) has been modelled using a GPD with the proposed (ν, µ) parameterisation.Figures 4 and 5 show that the hypothesis of a GPD (Weibull domain of attraction) is not rejected at 0.05 significance for both series.The presence of time trends and differences between REMO and deep buoy data can be assessed using Figs.7 and 8.
Figure 7 shows the joint posterior PDF of δ ν and α ν (lower-left panel).This joint PDF is characterised by its large dispersion thus pointing out the need of larger records to reliably estimate possible time trends in the shape of the GPD and the differences between the two data series.The marginal and conditional mode of α ν (lower-right panel) differ substantially from 0. But a criterion based on Bayesian discrepancy p value has been used (Gelman et al., 1996) and the value α ν = 0 is placed on the central part of the posterior marginal, thus making the possible positive trend in ν nonsignificant.We can conclude that there is some evidence of positive trend in the shape parameter but it should considered doubtful as there is no strong evidence against no trend during the observed time interval.Similarly, the change in ν from hindcast data to buoy data has the mode placed at a positive value near to 1 (upper panel), pointing out differences between the GPD shape for the two series.However, the value δ ν = 0 is fairly centred in the posterior marginal distribution of δ ν , thus meaning that the change in ν is not significant and should be considered carefully.These results are consistent with the decrease in the size of the extreme events during the buoy period observed in Fig. 2. The trend in ν (α ν ) is positive but not significant, and the difference between the REMO and buoy periods (δ ν ) is negative but non-significant.Altogether this leads to smaller values of ν for the buoy period than for hindcast.In the classical parameterisation this would mean a smaller ξ parameter for the buoy period than for the REMO.The upper limit of the GPD distribution in the classical parameterisation is −β/ξ .If the β parameters were the same for both periods, this would lead to a lower upper limit of the GPD distribution for buoy period than for REMO.
In the proposed parameterisation, δ µ accounts for the difference between the upper limit of the distribution of both periods.Figure 8 (lower-left panel) shows the posterior joint PDF of total drift in the parameter µ, α µ and the difference in µ corresponding to the two series.The value δ µ = 0, although not centred with respect to the marginal PDF (upper panel), is in a 90 % credible interval.The hypothesis of no difference between REMO and buoy data is plausible.The mode of the posterior marginal of δ µ is negative.Therefore, there is only a weak evidence against the upper limit of the wind speed being the same for hindcast and buoy observations.With regard to µ 0 , the posterior median estimate for µ 0 is about 0.53 (lower-right panel) which corresponds to about y sup = 82 m s −1 .The µ 0 density conditional to the mode of δ µ , µ 0 ) (red line, lower-right panel) approximately corresponds to µ 0 = 0.82, i.e. an upper limit of wind speed of 145 m s −1 .Figure 6 shows a kernel estimate of the predictive posterior density of the upper limit of the distribution for hindcast and buoy observations.
Figure 9 shows the estimated posterior density for the intensity difference λ b − λ h and the linear trend parameter α λ which is the total drift of λ along the observation time.The marginal histograms for these parameters are shown in the upper and lower-right panels.The equality of initial Poisson intensities for REMO and buoy series is assessed graphically by means of the line λ b − λ h = 0.The line lies in the lower tail of the posterior PDF, leading to a small Bayesian p value when testing for the equality of both values, i.e. the difference of the reference Poisson intensities is significant.The deep buoy series predicts about two events per year more than those hindcast by the REMO model.Regarding the trend, α λ , its histogram (lower-left panel) is approximately centred at 0. This provides a Bayesian p value near 0.5 when testing no linear trend in the Poisson intensity.Therefore, there is a non-significant trend in the intensity of the Poisson process.
It is still a matter of debate to what extent the frequency and intensity of windstorms may change as a consequence of the hypothetical climate change in the future.The results obtained for λ(t) are non-contradictory with other author's works, mainly devoted to investigating the changes in extreme winds, with methods based on global or regional climate models (e.g.Bolaños et al., 2004;Rockel and Woth, 2007).The slight and non-significant positive time trend observed for λ(t), corresponding to an increase of events, is in agreement to the hypothesis of climate change considered in IPCC reports (IPCC, 2007).However, the consistency of the locally obtained results with IPCC projections should be considered carefully, since these IPCC projections are global and are not necessarily coincident with local ones such as the ones considered in the paper.Actually, any assessment of single series cannot be expected to deliver any sensible result about global/regional climate change effects.The simultaneous analysis of many of these series of different nature and data quality may be needed for that.
We consider that the hindcast model may have a stronger inertia than the buoy measurements; actually, model data (being hourly averages) must exhibit less variability, i.e. have more "inertia" (in a figurative sense), than buoy data (10 min averages).Under steady, non-extremal stormy conditions, hindcast winds would have more energy than true winds, leading to an overestimation of winds.However, after the analysis of this data set, no significant change of either upper limit or shape of wind speed excesses distribution has been detected.The lack of overlap between buoy and hindcast series is not a major problem for this methodology -as shown by Ortego et al. (2012), the Bayesian estimation procedure makes it unnecessary to have overlapping series, though of course the obtained estimates show less uncertainty if the series are larger, and are better if they overlap.

Conclusions
A wind speed time series (REMO) has been analysed, together with a wind speed data set registered in a deep buoy in front of the Tarragona coast.A non-stationary Poisson/GPD model accounting for linear time trends and differences between the hindcast and buoy series has been assessed.The wind speed was log-transformed to deal with its ratio scale.The parameterisation of GPD of excesses over 15 m s −1 has been adopted to restrict distributions to be have finite tail, i.e. to be within the Weibull domain of attraction of maxima.The model was fitted using a Bayesian procedure.
The results confirm that the joint analysis of hindcast and directly observed wind speeds is useful to enlarge existing records used in extremal analysis.No significant time trends have been detected in occurrence rates of events and shape parameter of GPD.Importantly, there was no evidence in favour of the existence of differences in shape and upper limit of the GPD for excesses between the two sources of information, thus supporting the idea of using hindcast data for extremal analysis.Nonetheless, there were significant differences in the rate of occurrence of wind events recorded by hindcast and directly observable events, the latter being substantially higher, at about two events a year.
Although the total time of observation has been substantially increased by incorporating hindcast data, the uncertainty of the estimates is too large to attain conclusive results.This is the case of the time trend on the shape of GPD, represented by the parameter α ν , with a marginal distribution suggesting a positive trend, but without a clear statistical significance.

Figure 3 .Figure 3 .
Figure3.Joint likelihood contours for the (ξ, β) parameters of the classical parametrisation of the GPD distribution for the raw data (left) and for the log-transformed data (right).Orange contours for buoy; brown contours for REMO.

Figure 6 .
Figure 6.Kernel estimate of the predictive posterior density of the upper limit of the GPD distribution for hindcast and buoy observations.Orange contours for buoy; brown contours for REMO.The curves are reproduced in a reduced scale in order to show the mode of the velocity upper limit.

Figure 7 .Figure 8 .
Figure7.Lower-left panel shows the joint posterior density for (δ ν , α ν ), the difference between the two data series and the total drift along the observed time interval for the shape parameter ν (lowerleft panel).The upper and lower-right panels show the marginal histograms for the two parameters and their PDF conditional on the posterior joint mode.

Figure 9 .
Figure 9. Lower-left panel shows the joint posterior density for (λ b − λ h ), the difference on Poisson intensity for both signals, and α λ ), the total drift of the Poisson intensity (lower-left panel).The upper and lower-right panels show the marginal histograms for the two parameters and their PDF conditional on the posterior joint mode.
, have been

Table 1 .
Probability of GPD belonging to the Weibull attraction domain, ξ < 0, for the buoy and REMO series, for raw (in m s −1 ) and log-transformed data.