Sylvia Frühwirth-Schnatter

Helga Wagner

Stochastic model specification search for Gaussian and partial non-Gaussian state space models

Keywords: auxiliary mixture sampling, Bayesian econometrics, noncentered parameterization, Markov chain Monte Carlo, variable selection

Model specification for state space models is a difficult task as one has to decide which components to include in the model and to specify whether these components are fixed or timevarying. To this aim a new model space MCMC method is developed in this paper. It is based on extending the Bayesian variable selection approach which is usually applied to variable selection in regression models to state space models. For non-Gaussian state space models stochastic model search MCMC makes use of auxiliary mixture sampling. We focus on structural time series models including seasonal components, trend or intervention. The method is applied to various well-known time series.

Introduction

State space models are widely used in time series analysis to deal with processes which gradually change over time. Model specification, however, is a challenge for these models as one has to specify which components to include and to decide whether they are fixed or time-varying. For state space models, like for many other complex models, this often leads to testing problems which are non-regular from the view-point of classical statistics.

Thus, a classical approach toward model selection which is based on hypothesis testing such as a likelihood ratio test or information criteria such as AIC or BIC cannot be easily applied, because it relies on asymptotic arguments based on regularity conditions that are violated in this context.

Consider, for example, modeling a time series y = (y 1 , . . . , y T ) through the dynamic linear trend model, defined for t = 1, . . . , T as:

y t = µ t + ε t , ε t ∼ N 0, σ 2 ε , (1)
where µ t follows a random walk with a random drift starting from unknown initial values µ 0 and a 0 :

µ t = µ t-1 + a t-1 + ω 1t , ω 1t ∼ N (0, θ 1 ) , (2)
a t = a t-1 + ω 2t , ω 2t ∼ N (0, θ 2 ) . (3)
A typical specification problem arising for this model is to decide if the drift a t is timevarying rather than constant. However, testing θ 2 = 0 versus θ 2 > 0 results in a nonregular testing problem, because the null hypothesis lies on the boundary of the parameter space.

A similar specification problem is deciding which components are present in this times series model. For instance, is it necessary to include a dynamic drift term a t or should a t be removed because the level µ t follows a simple random walk? This is another non-regular problem, because again the null hypothesis can be rephrased as testing a 0 = θ 2 = 0.

The Bayesian approach is, in principle, able to deal with such non-regular testing problems. Suppose that K different models M 1 , . . . , M K are considered to be candidates for having generated the time series y. In a Bayesian setting each of these models is assigned a prior probability p(M k ) and the goal is to derive the posterior model probability p(M k |y) for each model M k , k = 1, . . . , K.

ACCEPTED MANUSCRIPT

disappears:

µ t = µ t-1 + ω 1t , ω 1t ∼ N (0, θ 1 ) . (4)
Another interesting special case is the linear trend model, where

y t = µ 0 + ta 0 + ε t , ε t ∼ N 0, σ 2 ε . (5)
The practical implementation of this approach is innovative in two respects. First, we employ a new prior for the process variances of the state space model by assuming that the square root of each process variance follows a normal distribution centered at 0. It is well-known that variable selection is, in general, sensitive to the choice of the prior, see e.g. [START_REF] Fernández | Benchmark priors for Bayesian model averaging[END_REF]. We show both for simulated as well as for real data that this prior is less influential on posterior inference when the true process variance is close to 0 than the usually applied inverted Gamma prior. This is in line with [START_REF] Gelman | Prior distributions for variance parameters in hierarchical models[END_REF] who came to similar conclusions for the related random-effects model.

Second, we derive an MCMC method for Gaussian as well as partially Gaussian state space models that performs stochastic model specification search by sampling the indicators simultaneously with the models parameters. The sampler is based on a noncentered parameterization of the state space model which generalizes previous work in this area such as [START_REF] Pitt | Analytic convergence rates and parameterization issues for the Gibbs sampler applied to state space models[END_REF] and [START_REF] Frühwirth-Schnatter | Efficient Bayesian parameter estimation for state space models based on reparameterizations[END_REF]. In combination with the normal prior on the square root of each process variance this leads to Gibbs sampler that is easily implemented. This is in contrast to [START_REF] Shively | Variable selection and function estimation in additive nonparametric regression using a data-based prior[END_REF] who consider a data-based prior which is motivated by the BIC criterion and leads to a sampling scheme where numerical integration has to be performed for each sweep of the MCMC scheme in order to sample indicators and parameters jointly.

To implement this approach for non-Gaussian state space modeling of times series of count data or binary, categorical, or multinomial data we make use of auxiliary mixture sampling [START_REF] Frühwirth-Schnatter | Auxiliary mixture sampling for parameterdriven models of time series of counts with applications to state space modelling[END_REF][START_REF] Frühwirth-Schnatter | Auxiliary mixture sampling with applications to logistic models[END_REF][START_REF] Frühwirth-Schnatter | Improved auxiliary mixture sampling for hierarchical models of non-Gaussian data[END_REF] which is a simple MCMC method for estimating a broad class of discrete-valued models.

Throughout the paper we focus on structural time series models including seasonal components, trend and an intervention effect and apply the method to various well-known time series.

Our method is based on a noncentered parameterization of the dynamic linear trend model which is discussed in the next subsection.

A Noncentered Parameterization

Define two independent random walk processes μt and ãt with standard normal independent increments as well as an integrated process Ãt : μt = μt-1 + ω1t , ω1t ∼ N (0, 1) , (6) ãt = ãt-1 + ω2t , ω2t ∼ N (0, 1) ,

Ãt = Ãt-1 + ãt-1 , (7)
which all are assumed to start at zero: μ0 = ã0 = Ã0 = 0. Combine the state equations ( 6) to ( 7) with following observation equation:

y t = µ 0 + ta 0 + θ 1 μt + θ 2 Ãt + ε t , ε t ∼ N 0, σ 2 ε , (8)
where µ 0 and a 0 are equal to the initial values for the level and the drift component and θ 1 and θ 2 are equal to the variances in the dynamic linear trend model defined in (1) to

(3). The resulting state space model is a noncentered parameterization of the dynamic linear trend model. To verify this define

a t = a 0 + θ 2 ãt , µ t = µ 0 + ta 0 + θ 1 μt + θ 2 Ãt . Then a t -a t-1 = θ 2 (ã t -ãt-1 ) = θ 2 ω2t = ω 2t , ω 2t ∼ N (0, θ 2 ) , µ t -µ t-1 = θ 1 (μ t -μt-1 ) + a 0 + θ 2 ( Ãt -Ãt-1 ) = θ 1 ω1t + a 0 + θ 2 ãt-1 = ω 1t + a t-1 , ω 1t ∼ N (0, θ 1 ) ,
which corresponds to the state equations ( 2) and (3).

ACCEPTED MANUSCRIPT

The noncentered parameterization of the dynamic linear trend model has a representation as a state space model with a state vector of dimension 3:

x t = Fx t-1 + w t , w t ∼ N (0, Q) , (9)
y t = Hx t + z f t α + ε t , ε t ∼ N 0, σ 2 ε , (10)
where x 0 = 0 3×1 and

x t =     μt ãt Ãt     , F =     1 0 0 0 1 0 0 1 1     , Q =     1 0 0 0 1 0 0 0 0     , H = √ θ 1 0 √ θ 2 , z f t = 1 t , α = µ 0 a 0 ′ .
This state space form could be used to perform Kalman filtering and to compute the

integrated likelihood p(y|ϑ) for ϑ = ( √ θ 1 , √ θ 2 , σ 2 ε , µ 0 , a 0 ).
The noncentered parameterization of the dynamic linear trend model, however, is not identified, because in the observation equation ( 8), the sign of √ θ 1 and the sequence {μ t } T 1 may be changed by multiplying all elements with -1 without changing the distribution of y 1 , . . . , y T . If we define a state vector x ⋆ t = (-μ t , ãt , Ãt ) ′ and a parameter ϑ ⋆ = (-√ θ 1 , √ θ 2 , σ 2 ε , µ 0 , a 0 ), then ϑ ⋆ and ϑ, although being different, define the same integrated likelihood:

p(y|ϑ) = p(y|x 1 , . . . , x T , θ 1 , θ 2 , σ 2 ε , µ 0 , a 0 )p(x 1 , . . . , x T )d(x 1 , . . . , x T ) = p(y|x ⋆ 1 , . . . , x ⋆ T , -θ 1 , θ 2 , σ 2 ε , µ 0 , a 0 )p(x ⋆ 1 , . . . , x ⋆ T )d(x ⋆ 1 , . . . , x ⋆ T ) = p(y|ϑ ⋆ ).
Similarly, the sign of √ θ 2 and the sequences {ã t } T 1 and { Ãt } T 1 may be changed without changing the distribution of y 1 , . . . , y T and ϑ

⋆ = ( √ θ 1 , - √ θ 2 , σ 2 ε , µ 0 , a 0 ) and ϑ define the same integrated likelihood, p(y|ϑ) = p(y|ϑ ⋆ ).
As a consequence, the likelihood function p(y|ϑ) is symmetric around 0 in the direction of √ θ 1 and √ θ 2 and therefore multimodal. If the data are generated by a dynamic linear trend model with true parameters (θ tr 1 , θ tr 2 , ξ tr ), where ξ tr = (σ 2,tr ε , µ tr 0 , a tr 0 ), then with increasing number of observations T , the modes of the likelihood function will be close to ( θ tr 1 , θ tr 2 , ξ tr ), (-θ tr 1 , θ tr 2 , ξ tr ), ( θ tr 1 ,θ tr 2 , ξ tr ), and (-θ tr 1 ,θ tr 2 , ξ tr ). If the true variances θ tr 1 and θ tr 2 are positive, then the likelihood function concentrates around

ACCEPTED MANUSCRIPT

four modes. If one the true variances is equal to 0 while the other is positive, two of those modes collapse and the likelihood is bimodal with increasing T . If both variances θ tr 1 and θ tr 2 are equal to zero, then the likelihood function becomes unimodal as T increases.

For illustration, Figure 1 shows contour and surface plots of the (scaled) likelihood function p(y| √ θ 1 , √ θ 2 , σ 2,tr ε , µ tr 0 , a tr 0 ) for a time series of length T = 1000 simulated from a dynamic linear trend model with µ tr 0 = 0.3, a tr 0 = -0.1 and σ 2,tr ε = 1 and four different combinations of θ tr 1 and θ tr 2 . There are clearly four modes, if both process variances are positive, two modes, if one of the variances is restricted to zero and a single mode, if both variances are restricted to 0.

Thus by considering the non-centered parameterization and allowing for nonidentifiability we gain important information about the hypothesis whether the variances of the state space model are zero.

The Parsimonious Dynamic Linear Trend Model

The noncentered parameterization of the dynamic linear trend model is very useful for model selection both for the components and the dynamics. The observation equation (8) of the noncentered parameterization represents the level of the time series y t as a superposition of the components at time t = 0 and the random processes μt and Ãt . Note that neither μt nor Ãt degenerate to a static component. A static component is obtained by setting the appropriate variance equal to 0. For instance, if the variance θ 1 is equal to 0, then √ θ 1 = 0 and μt is not used to explain y t . Similarly, if the variance θ 2 is equal to 0, then √ θ 2 = 0 and Ãt is not used to explain y t . This suggests to consider the choice of the variances θ 1 and θ 2 as a variable selection problem in regression model (8).

To this aim we introduce two binary indicators γ 1 and γ 2 , where √ θ i , and consequently

θ i , is equal to 0, if γ i = 0. If γ i = 1, then √ θ i is an unconstrained unknown parameter
which is estimated from the data under a suitable prior. Evidently, the indicators γ 1 and γ 2 decide if a certain component of the state vector is fixed or changes over time. If both γ 1 = 0 and γ 2 = 0, then the model reduces to a regression model with a linear trend,

given by (5).

To include or delete the trend, an additional indicator δ is introduced which decides, if the initial slope a 0 is equal to 0 or not. If δ = 0, then a 0 is equal to 0; otherwise, if δ = 1,

ACCEPTED MANUSCRIPT

then a 0 is an unknown parameter which is estimated from the data under a suitable prior.

This leads to the following parsimonious dynamic linear trend model:

μt = μt-1 + ω1t , ω1t ∼ N (0, 1) , ( 11
) ãt = ãt-1 + ω2t , ω2t ∼ N (0, 1) , (12)
Ãt = ãt-1 + Ãt-1 , (13)
y t = µ 0 + δta 0 + γ 1 θ 1 μt + γ 2 θ 2 Ãt + ε t , ε t ∼ N 0, σ 2 ε . ( 14
)
For a direct comparison with the usual dynamic linear trend model it is useful to rewrite the parsimonious model in the centered parameterization. Define

a t = δa 0 + θ 2 ãt , (15)
µ t = µ 0 + δta 0 + γ 1 θ 1 μt + γ 2 θ 2 Ãt . (16)
Then ( 11) to ( 14) may be rewritten as:

µ t = µ t-1 + δa 0 + γ 2 (a t-1 -δa 0 ) + γ 1 ω 1t , ω 1t ∼ N (0, θ 1 ) , (17)
a t = a t-1 + ω 2t , ω 2t ∼ N (0, θ 2 ) , (18)
y t = µ t + ε t , ε t ∼ N 0, σ 2 ε . (19)
Evidently, (δ, γ 1 , γ 2 ) = (1, 1, 1) corresponds to the unrestricted dynamic linear trend model

(2). The combination (δ, γ 1 , γ 2 ) = (0, 1, 0) leads to the local level model (4) which is also known as exponential smoothing, the combination (δ, γ 1 , γ 2 ) = (1, 0, 0) leads to a regression model with a deterministic linear trend, given by ( 5) and (δ, γ 1 , γ 2 ) = (0, 0, 0)

leads to i.i.d. normal data, y t ∼ N (µ 0 , σ 2 ε ).
The indicators δ, γ 1 and γ 2 have to be introduced carefully into the centered parametrization. Consider the following alternative choice which appears more natural than ( 17) and ( 18), but leads to nonindentifiability:

µ t = µ t-1 + δa t-1 + γ 1 ω 1t , ω 1t ∼ N (0, θ 1 ) , a t = a t-1 + γ 2 ω 2t , ω 2t ∼ N (0, θ 2 ) .
After recursive substitution we get following representation of the model as a normal linear mixed model:

y t = µ 0 + δta 0 + γ 1 t j=1 ω 1j + δγ 2 t-1 j=1 (t -j)ω 2j + ε t , A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
with fixed effects (µ 0 , a 0 ) and random effects (ω 1j , ω 2j ), j = 1 . . . , t. Only 6 models among the 8 possible combinations of the indicators (δ, γ 1 , γ 2 ) are identifiable, because γ 2 is not identified, if δ = 0. In contrast to that, model ( 17) to ( 19) has the representation

y t = µ 0 + δta 0 + γ 1 t j=1 ω 1j + γ 2 t-1 j=1 (t -j)ω 2j + ε t .
Evidently, all 8 combinations of the indicators (δ, γ 1 , γ 2 ) are identifiable.

The noncentered parameterization of the parsimonious dynamic linear trend model

given by ( 11) to ( 14) has the following representation as a state space model:

x t = Fx t-1 + w t , w t ∼ N (0, Q) , ( 20
)
y t = H(γ 1 , γ 2 )x t + z f t (δ)α + ε t , ε t ∼ N 0, σ 2 ε ,
where x t , F, Q and α are the same as in (9), while H and z f t depend on the model indicators:

H(γ 1 , γ 2 ) = γ 1 √ θ 1 0 γ 2 √ θ 2 , z f t (δ) = 1 δt .

Prior Distributions

To perform Bayesian estimation one has to choose a prior distribution p(δ, γ 1 , γ 2 ) for all possible combinations of indicators. Subsequently, we assume a uniform distribution over all 8 combinations of the indicators. A more flexible distribution is discussed in Section 5.

As common for dynamic linear trend models, we assume that apriori µ 0 and a 0 are independently normally distributed, µ 0 ∼ N (y 1 , P 0,11 σ 2 ε ) and a 0 ∼ N (0, P 0,22 σ 2 ε ). Furthermore we assume an inverted Gamma prior G -1 (c 0 , C 0 ) for the observation variance

σ 2 ε .
In contrast to previous work, we do not use the usual inverted Gamma priors θ 1 ∼ G -1 (d 0,1 , D 0,1 ) and θ 2 ∼ G -1 (d 0,2 , D 0,2 ). We employ a new prior for the process variances of the state space model by assuming that the square root of each process variance follows a normal distribution centered at 0, i.e. √ θ 1 ∼ N (0, B 0,1 σ 2 ε ) and √ θ 2 ∼ N (0, B 0,2 σ 2 ε ). It be should noted that the two priors are equivalent only under the limiting case of

A C C E P T E D M A N U S C R I P T

ACCEPTED MANUSCRIPT

following improper priors: an inverted Gamma prior where d 0,i = -0.5 and D 0,i = 0, i.e. p(θ i ) ∝ √ θ i , and a normal prior where B -1 0,i = 0, i.e. p(

√ θ i ) ∝ constant.
It is well-known, that the hyperparameters of the inverted Gamma prior θ i ∼ G -1 (d 0,i , D 0,i ) strongly influence the posterior density of θ i , if the true value of θ i is close to 0. In contrast to that, the normal prior appears to be less influential and more suitable under model specification uncertainty than the inverted Gamma prior.

Consider, for example, a local level model,

µ t = µ t-1 + ω 1t , ω 1t ∼ N (0, θ 1 ) , y t = µ t + ε t , ε t ∼ N 0, σ 2 ε , (21)
where θ 1 is unknown and σ 2 ε is assumed to be known.

To compare the inverted Gamma prior to the normal prior we consider the posterior density of the parameter ± √ θ 1 which is obtained from θ 1 by multiplying the square root of θ 1 with a random sign. We added the ± sign to emphasize that the sign of this parameter is not identified. The posterior of ± √ θ 1 allows to explore the hypothesis that θ 1 = 0.

Due to the symmetry of the likelihood discussed in Subsection 2.1, the posterior density of ± √ θ 1 is symmetric around zero as long as the prior is also symmetric around 0. If the unknown variance θ tr 1 is significantly different from zero, then the posterior density of ± √ θ 1 is likely to be bimodal with the modes being close to ± θ tr 1 . Otherwise, if θ tr 1 is close to or equal to zero, then the posterior density of ± √ θ 1 is likely to be centered around zero.

For illustration, we consider posterior inference for T = 100 observations simulated from the local level model ( 21) with σ 2 ε = 0.01 and two different values for θ 1 , namely θ tr 1 = 0.01 and θ tr 1 = 0. For any type of prior distribution, the posterior distribution of ± √ θ 1 is derived using numerical integration. For both values of θ tr 1 , the posterior of ± √ θ 1 is plotted in Figure 2 under the inverted Gamma prior θ 1 ∼ G -1 (0.5, D 0 ) as well as under the normal prior ± √ θ 1 ∼ N (0, B 0 σ 2 ε ) for various scale parameters D 0 and B 0 .

Whereas the posterior is fairly robust to the choice of the hyperparameter B 0 in the normal prior, it turns out to be rather sensitive to the hyperparameter D 0 of the inverted Gamma prior. Both posteriors are roughly the same for θ tr 1 = 0.01 and clearly indicate that θ tr 1 > 0. A remarkable difference, however, occurs if θ tr 1 = 0. Under the normal prior, the posterior of ± √ θ 1 is centered at 0 strongly supporting the hypothesis that θ tr 1 = 0. The inverted Gamma density, however, shrinks the posterior of ± √ θ 1 away from 0, falsely indicating that θ tr 1 > 0.

MCMC Estimation

An MCMC approach is implemented to sample jointly the indicators (δ, γ) = (δ, γ 1 , γ 2 ), the unrestricted elements of the parameter β = (µ 0 , a 0 , √ θ 1 , √ θ 2 ), the observation variance σ 2 ε , and the latent state process x = (x 1 , . . . , x T ), where x t is the state vector defined in (9).

When sampling the indicators (δ, γ) we marginalize over the parameters for which variable selection is carried out, as suggested by [START_REF] Geweke | Variable selection and model comparison in regression[END_REF] and [START_REF] Smith | Nonparametric regression using Bayesian variable selection[END_REF], see also [START_REF] George | Approaches for Bayesian variable selection[END_REF]. To make this feasible, we use the noncentered parameterization of the dynamic linear trend model. Conditional on the state process x = (x 1 , . . . , x T ), the observation equation ( 14) defines a standard regression model

y t = z δ,γ t β δ,γ + ε t , ε t ∼ N 0, σ 2 ε . ( 22
)
If all indicators take the value one, then β δ,γ = β and z δ,γ t = z t , where z t = (1, t, μt , Ãt ).

Otherwise the restricted parameter β δ,γ and the corresponding predictors z δ,γ t contain only those elements of β and z t , respectively, for which the corresponding indicator is equal to 1. Under the conjugate prior

β δ,γ ∼ N a δ,γ 0 , A δ,γ 0 σ 2 ε , σ 2 ε ∼ G -1 (c 0 , C 0 ) , (23)
the posterior p(δ, γ|x, y) is obtained from Bayes' theorem:

p(δ, γ|x, y) ∝ p(y|δ, γ, x)p(δ, γ), (24)
where p(y|δ, γ, x) is equal to the marginal likelihood of the regression model ( 22):

p(y|δ, γ, x) = 1 (2π) T /2 |A δ,γ T | 1/2 |A δ,γ 0 | 1/2 Γ(c T )C c 0 0 Γ(c 0 )(C δ,γ T ) c T . ( 25
)
Here A δ,γ T , c T and C δ,γ T denote the posterior moments of β δ,γ and σ 2 ε given below in ( 26) to (28). It should be noted that such a closed form expression for p(y|δ, γ, x) is not available if any of the indicators γ 1 and γ 2 is equal to 1 and an inverted Gamma prior is chosen for θ 1 and θ 2 . The MCMC scheme reads:

(a) Sample the indicators (δ, γ) = (δ, γ 1 , γ 2 ), the initial values µ 0 and a 0 , all variance parameters √ θ 1 and √ θ 2 and the observation variance σ 2 ε jointly in one block:

(a1) Sample the indicators from p(δ, γ|x, y) given in ( 24);

(a2) sample σ 2 ε from G -1 c T , C δ,γ
T , and, conditional on σ 2 ε , sample µ 0 , a 0 (if unrestricted), and all unrestricted variance parameters √ θ 1 and

√ θ 2 jointly from the normal posterior N a δ,γ T , A δ,γ T σ 2 ε where A δ,γ T = (Z δ,γ ) ′ Z δ,γ + (A δ,γ 0 ) -1 -1 , ( 26
)
a δ,γ T = A δ,γ T (Z δ,γ ) ′ y + (A δ,γ 0 ) -1 a δ,γ 0 , c T = c 0 + T /2, (27)
C δ,γ T = C 0 + 1 2 y ′ y + (a δ,γ 0 ) ′ (A δ,γ 0 ) -1 a δ,γ 0 -(a δ,γ T ) ′ (A δ,γ T ) -1 a δ,γ T , (28)
and Z δ,γ is the regressor matrix with rows equal to z δ,γ t ;

(a3) set all restricted initial values and all restricted variances equal to 0.

(b) Sample x = (x 1 , . . . , x T ) from the state space form (20).

(c) Perform a random sign switch for √ θ 1 and {μ t } T 1 . Thus with probability 0.5 the draws of these parameters remain unchanged, while they are substituted by -√ θ 1

and {-μ t } T 1 with the same probability. Perform another random sign switch for

√ θ 2 , {ã t } T 1 and { Ãt } T 1 .
A few comments are in order. The dimension of the normal distribution appearing in step (a2) depends on the number of unrestricted components and is equal to 1 + δ + γ 1 + γ 2 .

In step (b), forward-filtering-backward-sampling (FFBS, Frühwirth-Schnatter (1994);

Carter and Kohn (1994); De Jong and [START_REF] De | The simulation smoother for time series models[END_REF]) is used to sample x = (x 1 , . . . , x T ).

To speed up sampling, a reduced state space form is used if γ 1 or γ 2 is 0. If, for instance, γ 1 = 0, then the observation equation is independent of {μ t } T 1 . FFBS is applied to the reduced state vector x t = (ã t , Ãt ) ′ , while μ1 , . . . , μT is sampled from (11). A similar method

ACCEPTED MANUSCRIPT

applies, if γ 2 = 0, with reduced state vector x t = μt . If both indicators γ 1 and γ 2 are equal to 0, then no FFBS is needed, as sampling of x from the prior is straightforward.

Sampling of the state process in step (b) is based on the noncentered parameterization. The unknown components a t and µ t in the centered parameterization are easily reconstructed from the MCMC draws using ( 15) and ( 16).

We found it useful to start from an unrestricted model and to run the first say 1000 draws of burn-in without variable selection. This allows to generate sensible starting values for the state process and the parameters of the unrestricted model before variable selection actually sets in.

3 Extension to the Basic Structural Model

The Parsimonious Basic Structural Model

In the basic structural model, a seasonal component is added to the dynamic linear trend model discussed in Section 2, see e.g. [START_REF] Harvey | Forecasting, Structural Time Series Models and the Kalman Filter[END_REF]:

s t = -s t-1 -• • • -s t-S+1 + ω 3t , ω 3t ∼ N (0, θ 3 ) , ( 29
)
y t = µ t + s t + ε t , ε t ∼ N 0, σ 2 ε , (30)
where µ t is the same as in ( 2) and (3) and S is the number of seasons. The initial seasonal pattern is given by s 0 = (s -S+1 , . . . , s 0 ) with s -S+1 + . . . + s 0 = 0. In addition to the model specification problems discussed in Section 2, a decision has to be made if a seasonal pattern is present and if this pattern is fixed or dynamic. To this aim, two additional binary stochastic indicators δ 3 and γ 3 are introduced. δ 3 decides, if the initial seasonal pattern is equal to 0, whereas γ 3 controls if it changes over time. As before, the indicators are introduced into the noncentered version of the model.

Combine the following stochastic difference equation:

st = -s t-1 -• • • -st-S+1 + ω3t , ω3t ∼ N (0, 1) , (31)
where s-S+1 = . . . = s0 = 0 with the state equations ( 6) to ( 7) and following observation equation:

y t = µ t + δ 3 s 0,q(t) + γ 3 θ 3 st + ε t , ε t ∼ N 0, σ 2 ε , (32)

ACCEPTED MANUSCRIPT

where θ 3 is equal to the variance of the error term in (29), µ t is the same as in ( 16) and s 0,q(t) with q(t) = 1 + (t -1) mod S is the seasonal component corresponding to time t.

The resulting state space model is a noncentered parameterization of the basic structural model.

If γ 3 = 0, we define θ 3 = 0 and the resulting seasonal pattern is fixed. If δ 3 = 0, we set the initial seasonal pattern to zero, s 0 = 0. If both indicators are equal to 0, then no seasonal pattern is present in the time series and the model reduces to the dynamic linear trend model studied in Section 2.

The non-centered parameterization (32) could be written as

y t = µ 0 + δta 0 + δ 3 s 0,q(t) + γ 1 t j=1 ω 1j + γ 2 t-1 j=1 (t -j)ω 2j + γ 3 t j=1 ω 3j + ε t ,
with fixed effects µ 0 , a 0 and s 0,q(t) and random effects ω 1j , ω 2j and ω 3j . Evidently, all 2 5 = 32 combinations of indicators are identifiable.

As before, the noncentered model is not identified, as the sign of √ θ 3 and the sequence {s t } T 1 may be changed without changing the likelihood function. As a consequence, the likelihood function p(y|ϑ) where ϑ = (

√ θ 1 , √ θ 2 , √ θ 3 , µ 0 , a 0 , s 0 , σ 2 ε ) is symmetric around 0 in the direction of √ θ i , i = 1, 2, 3.
With an increasing number of observations T , the modes of the likelihood function will be close to all combinations of (± θ tr 1 , ± θ tr 2 , ± θ tr 3 , ξ tr ), where ξ tr = (µ tr 0 , a tr 0 , s tr 0 , σ 2,tr ε ). Thus with an increasing number of observations, the likelihood function has eight modes as long as in the data generating process the true variances θ tr 1 , θ tr 2 and θ tr 3 are positive. If one of the true variances is equal to 0 while the others are positive, half of those modes are identical leaving four modes. If two of the true variances are equal to 0 while the other is positive, only two modes are different leaving a bimodal likelihood with an increasing number of observations T . If all variances are equal to zero, then the likelihood function will be unimodal with an increasing number of observations T .

It is easy to verify that in the centered parameterization the parsimonious model is equivalent to combining (2) and (3) with state equation ( 29) and following observation equation:

y t = µ t + δ 3 s 0,q(t) + γ 3 (s t -δ 3 s 0,q(t) ) + ε t , ε t ∼ N 0, σ 2 ε . (33) θ 1 , √ θ 2 , √ θ 3
), the observation variance σ 2 ε , and the latent state process x = (x 1 , . . . , x T ), where x t is following state vector:

x t = μt ãt Ãt st . . . st-S+2 ′ .
The MCMC sampling scheme introduced in Subsection 2.4 is easily modified to deal with a basic structural model. Conditional on the state process x = (x 1 , . . . , x T ), the observation equation ( 32) of the non-centered parameterization of the basic structural model is a standard regression model as in ( 22) with appropriate regressors z δ,γ t . Under the same conditionally conjugate prior for β δ,γ and σ 2 ε as in ( 23), the marginal likelihood p(y|δ, γ, x) and all posterior moments are then computed exactly as in Subsection 2.4. This leads to following MCMC scheme: As in Subsection 2.4, FFBS is applied to a reduced state vector, if any of the indicators γ i = 0 is equal to 0, while the remaining components are sampled from the prior.

Prior Specification

To run the MCMC schemes, prior distributions have to be defined. As before, we assume a uniform prior distribution over all possible indicators δ and γ.

For the observation variance σ 2 ε we choose a hierarchical prior where σ 2 ε ∼ G -1 (c 0 , C 0 ) and C 0 ∼ G (g 0 , G 0 ) with c 0 = 2.5, g 0 = 5 and G 0 = g 0 /(0.75Var(y)(c 0 -1)). For this hierarchical prior it is necessary to add an additional sampling step were C 0 is sampled

ACCEPTED MANUSCRIPT conditional on σ 2 ε from the conditional Gamma posterior C 0 |σ 2 ε ∼ G (g 0 + c 0 , G 0 + 1/σ 2 ε ) at each sweep of the sampler.
Prior (23) assumes normality not only for the initial values µ 0 , a 0 , and s 0 , but also for all remaining parameters. For the same reasons as in Subsection 2.3, we do not use inverted Gamma priors for the variances θ 1 , . . . , θ 3 as usual in the basic structural model, but assume that the parameters ± √ θ 1 , ± √ θ 2 and ± √ θ 3 follow a normal prior.

In our case studies, we found the following prior choices useful for variable selection.

First, we use a partially proper prior which combines the improper prior p(µ 0 ) ∝ 1 for µ 0 with a proper prior N 0, B δ,γ 0 σ 2 ε on the remaining unrestricted elements of β δ,γ , where B δ,γ 0 = B 0 I. This prior corresponds to choosing a δ,γ 0 = 0 and

A δ,γ 0 -1 = 0 (B δ,γ 0 ) -1 . ( 34
)
Under this prior, the sampling scheme described above has to be changed slightly, because the marginal likelihood p(y|δ, γ, x) and the posterior parameter c T read:

p(y|δ, γ, x) = 1 (2π) (T -1)/2 |A δ,γ T | 1/2 |B δ,γ 0 | 1/2 Γ(c T )C c 0 0 Γ(c 0 )(C δ,γ T ) c T , c T = c 0 + (T -1)/2.
Another prior commonly used in model selection is the fractional prior (O'Hagan, 1995).

In the present context, this is a conditional fractional prior for regression model ( 22) which depends on the state vector x and is defined as

p(β δ,γ |σ 2 ε ) ∝ p(y|β δ,γ , σ 2 ε ) b = 1 2πσ 2 ε T b/2 exp - b 2σ 2 ε (y -Z δ,γ β δ,γ ) ′ (y -Z δ,γ β δ,γ ) .
The fractional prior can be interpreted as posterior of a non-informative prior and a fraction b of the data y. It reads

β δ,γ |σ 2 ε ∼ N a δ,γ T , A δ,γ T σ 2 ε /b ,
where a δ,γ T and A δ,γ T are the posterior moments under a non-informative prior:

A δ,γ T = (Z δ,γ ) ′ Z δ,γ -1 , a δ,γ T = A δ,γ T (Z δ,γ ) ′ y. ( 35
)

ACCEPTED MANUSCRIPT

In the MCMC sampling scheme all posterior moments as well as the marginal likelihood p(y|δ, γ, x) have to be modified according to:

c T = c 0 + (1 -b) 2 T, C δ,γ T = C 0 + (1 -b) 2 (y ′ y -(a δ,γ T ) ′ (A δ,γ T ) -1 a δ,γ T ), p(y|δ, γ, x) = b q/2 Γ(c T )C c 0 0 (2π) T (1-b)/2 Γ(c 0 )(C δ,γ T ) c T ,
where q is the dimension of β δ,γ , while a δ,γ T and A δ,γ T are the same as in (35).

UK coal consumption data

We reconsider the series of UK coal consumption, analyzed in [START_REF] Harvey | Forecasting, Structural Time Series Models and the Kalman Filter[END_REF][START_REF] Frühwirth-Schnatter | Data augmentation and dynamic linear models[END_REF] and Frühwirth-Schnatter (1995), among others. Data are quarterly from 1/1960 to 4/1982, see Figure 3, panel (a). We model the series on the log scale by a basic structural model.

All subsequent implementation are carried out using Matlab (Version 7.2.0) on a notebook with a 2.0 GHz processor.

Bayesian Inference for the Unrestricted Basic Structural Model

For illustration, we start with inference for the unrestricted basic structural model without variable selection and compare the priors θ i ∼ G -1 (-0.5, 10 -7 ) with the priors ± √ θ i ∼ N (0, 1). The remaining priors are µ 0 ∼ N (0, 100σ 2 ε ) and σ 2 ε ∼ G -1 (0, 0).

Estimated state components are plotted for the inverted Gamma prior in Figure 3.

The posterior densities of the transformed process variances ± √ θ i , i = 1, . . . , 3 are plotted in Figure 4. Evidently, the posterior density of any parameter ± √ θ i has to be symmetric around zero. If the unknown variance θ i is systematically different from zero, then the posterior density of ± √ θ i is likely to be bimodal; otherwise, if θ i is close to zero, the posterior density of ± √ θ i will be centered around zero. This should allow to explore the hypothesis that θ i = 0.

For the normal prior ± √ θ i ∼ N (0, 1), the posterior densities of ± √ θ 1 and ± √ θ 3 are unimodal and centered at 0, while the posterior of ± √ θ 2 is bimodal. This indicates

ACCEPTED MANUSCRIPT

that θ 1 and θ 3 are equal to 0, while θ 2 > 0. This finding is confirmed by stochastic model selection search in Subsection 3.4.2. Under the inverted Gamma prior, all posterior densities are bimodal and ± √ θ i is bounded away from 0, providing spurious evidence for an unrestricted model.

The MCMC draws underlying Figure 4 were obtaining by running MCMC sampling for 40,000 iterations after a burn-in of 10,000. For the normal prior we use the new MCMC scheme presented in Subsection 3.2 which is based on the non-centered parameterization of the basic structural model. For the inverted Gamma prior we use the usual two-step Gibbs sampler which is based on the centered parameterization, see e.g. [START_REF] Durbin | Time Series Analysis by State Space Methods[END_REF]. MCMC draws for ± √ θ i are obtained by multiplying the square root of the MCMC draws θ (m) i with a random sign.

An interesting difference between the centered and the non-centered parameterization lies in the mixing properties of the corresponding MCMC draws. If some variances are equal to or close to 0, the corresponding MCMC draws mix badly under the centered parameterizaton, while mixing is perfect under the noncentered parameterizaton, see Figure 5.

Stochastic Model Specification Search

Stochastic model specification search is carried out using the prior described in Subsection 3.3. We compare partially proper priors with different prior variances B 0 with fractional priors with different fractions b. For each prior, MCMC sampling was carried out for M = 100, 000 draws after a burn-in of 20, 000 draws. The first 1000 draws of the burn-in were drawn from the unrestricted model, model selection began after these first 1000 draws. Depending on the prior, running MCMC takes between 43 and 52 minutes, see also Table 2.

Results of the variable selection procedure are summarized in Table 1 and2. The most frequently visited model in Table 1 is robust against the prior choice, only the frequency with which this model is selected varies. The same model results for all priors, if in Table 2 an indicator is estimated to be 1, if the corresponding posterior inclusion probability is greater than or equal to 0.5.

As expected from panel (a) and (d) in Figure 3, a seasonal pattern is present in the selected model (δ 3 = 1), but it is fixed and does not change over time (γ 3 = 0). The drift a t is stochastic (γ 2 = 1), but the initial value a 0 is selected to be 0 (δ = 0). This is plausible from panel (c) in Figure 3, where the pointwise confidence band covers a t = 0 at t = 0, but does not contain the restricted line where a t = 0 for all t. Finally, no additional noise ω 1t is added in (2), since γ 1 = 0.

Comparison to Marginal Likelihood Computation

For this case study, it is possible to compare variable selection to the more traditional approach of computing marginal likelihoods p(y|M) of all 32 models M = (δ, δ 3 , γ 1 , γ 2 , γ 3 ), obtained by all possible combinations of indicators. Because we are dealing with a linear

Gaussian state space model it is possible to integrate out all states and the initial values by running a Kalman filter and we are left with a low dimensional integration over the parameter (θ M , σ 2 ε ) where θ M contains all unconstrained variances θ 1 , θ 2 and θ 3 present in model M.

Each marginal likelihood p(y|M) is computed by importance sampling with 20000 draws from the following importance density. For each model, MCMC for the noncentered parameterization with fixed indicators (δ, γ) was run for 40, 000 iterations after a burn-in of 10, 000 draws. The density of a t ν (m, S)-distribution with ν = 10 was fitted to the MCMC draws of (log(θ M ), log(σ 2 ε )) by matching moments, i.e. m is equal to the mean of the transformed MCMC draws, while S is equal to the covariance matrix times (ν -2)/ν.

Computation time for a single model lies between 15 and 16 CPU minutes. Computation for all 32 models took 494 CPU minutes.

The resulting estimators, together with their standard errors, are reported in Table 3 and indicate that importance sampling yields very precise estimates of the marginal likelihood.

Marginal likelihood computations are based on following priors. We choose the same hierarchical prior on the observation variance σ 2 ε as in Subsection 3.4.2 and marginalize over the random hyperparameter C 0 to facilitate computation of the marginal likelihood:

p(σ 2 ε |c 0 , g 0 , G 0 ) = G g 0 0 Γ(c 0 + g 0 ) Γ(c 0 )Γ(g 0 )(σ 2 ε ) c 0 +1 (G 0 + 1/σ 2 ε ) -(c 0 +g 0 ) .

ACCEPTED MANUSCRIPT

Since marginal likelihood computation is not possible for improper priors, we modify the partially proper used in Subsection 3.4.2 for variable selection slightly by assuming µ 0 ∼ N (0, 10000σ 2 ε ), whereas we use the same N (0, B 0 σ 2 ε ) prior with B 0 = 1 and B 0 = 100, respectively, for the initial value of all other components and the process standard deviations. We did not consider the fractional prior, because it is defined as a conditional prior given the all states, which makes it impossible to integrate out the high-dimensional state vector in the fashion described above.

The same prior is used for the stochastic model specification search algorithm which was run for 100,000 iterations after a burn-in of 20000. Computation time for running variable selection is about 52 CPU minutes which is considerably faster than computing the marginal likelihood for all models.

Table 3 reports the number h(M) of times the model M was visited during the 100,000 MCMC iterations. We find that the relative frequency h(M)/100, 000 is close the model probability p(M|y) which is estimated from the marginal likelihood p(y|M) under the uniform prior p(M) = 1/32 underlying variable selection.

In Table 1 we found that changing the prior had no effect on the best model, however, the ranking of the remaining models changed. From Table 3 we see that this effect is also present for the marginal likelihood. This sensitivity of the marginal likelihood to prior choice was even more pronounced when we considered the conventional inverted Gamma prior for the process variances θ i .

Model Selection for Non-Gaussian State Space Models

The investigations in [START_REF] Shively | Variable selection and function estimation in additive nonparametric regression using a data-based prior[END_REF] show that variable selection in state space models is also feasible for binary data. In this section we show how the variable selection approach developed in Section 2 and Section 3 for Gaussian state space model may be extended to nonnormal state space models using auxiliary mixture sampling [START_REF] Frühwirth-Schnatter | Auxiliary mixture sampling for parameterdriven models of time series of counts with applications to state space modelling[END_REF][START_REF] Frühwirth-Schnatter | Auxiliary mixture sampling with applications to logistic models[END_REF][START_REF] Frühwirth-Schnatter | Improved auxiliary mixture sampling for hierarchical models of non-Gaussian data[END_REF]. This allows variable selection for state space modelling of times series of small counts based on the Poisson distribution and of binary as well as categorical time

ACCEPTED MANUSCRIPT

series based on the logit transform. We provide an illustrative application to two time series of small counts.

A Basic structural model for Count Data including Intervention

For count data the basic structural model reads [START_REF] Harvey | The effects of seat belt legislation on British road casualties: A case study in structural time series modelling[END_REF]:

y t ∼ P (e t λ t ) , log λ t = µ t + s t , ( 36
)
µ t = µ t-1 + a t-1 + ω 1t , ω 1t ∼ N (0, θ 1 ) (37) a t = a t-1 + ω 2t , ω 2t ∼ N (0, θ 2 ) , (38)
s t = -s t-1 -• • • -s t-S+1 + ω 3t , ω 3t ∼ N (0, θ 3 ) . ( 39
)
To account for the intervention at t = t int , equation ( 38) is modified in the following way:

µ t = µ t-1 + a t-1 + ∆ + ω 1t .

Stochastic model specification search

Indicators δ, δ 3 , γ 1 , γ 2 and γ 3 are introduced as in Section 3 to select the structural components, and an additional indicator δ 4 is introduced for the intervention effect. In the centered parameterization, equations ( 36) and ( 37) are modified in the following way:

log λ t = µ t + δ 3 s 0,q(t) + γ 3 (s t -δ 3 s 0,q(t) ), ( 40
)
µ t = µ t-1 + δa 0 + γ 2 (a t-1 -δa 0 ) + δ 4 I {t=t int } ∆ + γ 1 ω 1t , (41)
while ( 38) and ( 39) are unaffected. For MCMC estimation, the noncentered version of this model is required which reads:

log λ t = µ 0 + δta 0 + δ 3 s 0,q(t) + δ 4 I {t≥t int } ∆ + γ 1 θ 1 μt + γ 2 θ 2 Ãt + γ 3 θ 3 st ,
where μt and Ãt are defined as in ( 11) to ( 13), while st is defined as in (31).

ACCEPTED MANUSCRIPT

(d2) Sample the component indicator r tj for j = 1, 1 + min(y t , 1) from the following discrete distribution where k = 1, . . . , R(ν tj )

Pr(r tj = k|τ tj , λ t ) ∝ w k (ν tj ) f N (-log τ ij -log λ t -log e t ; m k (ν tj ), s 2 k (ν tj )).
Note that in step (a) marginalizing over the variables and components which are subject to model selection would not be possible for non-Gaussian state space models without the use of auxiliary mixture sampling or another augmentation scheme that leads to a conditionally Gaussian model. Such data augmentation schemes which enable variable selection in non-Gaussian models have been applied earlier by [START_REF] Holmes | Bayesian auxiliary variable models for binary and multinomial regression[END_REF] for binary and multinomial regression model and by Tüchler ( 2008) for binary and multinomial regression models with random effects.

The partially proper normal prior and the fractional prior considered in Subsection 3.3 are easily adjusted for non-Gaussian state space models. A partially proper normal prior combines p(µ 0 ) ∝ 1 with a proper prior N 0, B δ,γ 0 on the remaining unrestricted elements of β δ,γ where B δ,γ 0 = B 0 I corresponds to a δ,γ 0 = 0 and A δ,γ 0 being the same as in (34). The marginal likelihood for this prior reads

p(y|δ, γ, x, u) = |Σ| -1/2 |A δ,γ T | 1/2 (2π) (T -1)/2 |B δ,γ 0 | 1/2 • exp - 1 2 ỹ′ Σ -1 ỹ -(a δ,γ T ) ′ (A δ,γ T ) -1 a δ,γ T .
For a fractional prior, derived as in Subsection 3.3, the marginal likelihood is given as

p(y|δ, γ, x, u) = b q/2 |Σ| -1 (2π) T (1-b)/2 • exp - (1 -b) 2 (ỹ ′ Σ -1 ỹ -(a δ,γ T ) ′ (A δ,γ T ) -1 a δ,γ T ) , (48)
where

(A δ,γ T ) -1 = (Z δ,γ ) ′ Σ -1 Z δ,γ and a δ,γ T = A δ,γ T (Z δ,γ ) ′ Σ -1 ỹ.

Road Safety Data

We analyze a time series consisting of monthly counts of killed or injured pedestrians, aged 6-10, from 1987-2005 in Linz, which is the third largest town in Austria. 2 The observations are a series of small counts not exceeding 5, see Figure 6. A new law intended to increase

ACCEPTED MANUSCRIPT

road safety came into force in Austria on October 1, 1994, since when pedestrians who want to use a pedestrian crossing have to be allowed to cross. Of interest is the effect of this law on the (monthly) risk of being killed or seriously injured in a road accident as a child living in Linz.

The basic structural model with intervention effect for Poisson counts defined in Subsection 4.1 is fitted to the number y t of children killed or seriously injured in time period t, y t ∼ P (e t λ t ), where e t is the number of children living in Linz. Model specification search is carried out to identify an appropriate model.

Bayesian Inference for the unrestricted basic structural model with intervention

For illustration, we start with Bayesian inference without variable selection for the unrestricted basic structural model with intervention. For these data, we were not able to estimate the model under the completely centered parameterization as MCMC did not converge. For this reason, we compare a parameterization where only the season is noncentered [START_REF] Frühwirth-Schnatter | Auxiliary mixture sampling for parameterdriven models of time series of counts with applications to state space modelling[END_REF] under the priors θ i ∼ G -1 (0.1, 0.001) , i = 1, 2, ± √ θ 3 ∼ N (0, 1) with a fully noncentered model with priors ± √ θ i ∼ N (0, 1) , i = 1, 2, 3. For both parameterizations we assume that µ 0 ∼ N (log(y 1 /e 1 ), 1) = N (-9.0084, 1)

and that the unknown initial values of the other components and the intervention effect follow a standard normal prior distribution.

Figure 8 shows histograms of the MCMC draws for ± √ θ i , i = 1, . . . , 3 for both priors.

For the normal prior the posterior of all parameters ± √ θ i , i = 1, . . . , 3 is clearly centered at 0, suggesting that the state space model is overfitting and the data may be explained by a simply Poisson regression model. This finding is confirmed by stochastic model specification search in Subsection 4.2.2. As in Subsection 3.4.1 the inverted Gamma density is very influential and shrinks the posterior densities of ± √ θ 1 and ± √ θ 2 away from 0, spuriously suggesting that θ 1 > 0 and θ 2 > 0.

Figure 9 shows the various components of the unconstrained model like the smoothed level µ t with pointwise 95% credibility intervals. The estimated monthly risk λ t for a child to be seriously injured or killed seems to decrease at the time of intervention. The drift a t is not significantly different from 0 over the whole observation period. The seasonal

ACCEPTED MANUSCRIPT

component has significantly lower values than the annual average in the holiday months July and August and higher values in June and October.

We used 20,000 iterations after a burn-in of 5,000 for each parameterization. As in Subsection 3.4.1, we observe much better mixing behavior of the MCMC sampler under the non-centered parametrization, see Figure 7.

Stochastic Model Specification Search

Stochastic model specification search is carried out using the prior described in Subsection 3.3. We compare partially proper priors with different prior variances B 0 with fractional priors with different fractions b. For each prior, MCMC sampling was carried out for M = 100, 000 draws after a burn-in of 20, 000 draws. The first 1000 draws of the burn-in were drawn from the unrestricted model, model selection began after these first 1000 draws. Depending on the prior, running MCMC takes between 125 and 165 minutes, see also Table 5.

Results of the variable selection procedure are summarized in Table 4 and5 and γ 2 = 0 imply that a t = a 0 = 0 for the whole observation period. The initial seasonal pattern is significant (δ 3 = 1), but does not change over time (γ 3 = 0). The level of the model is constant before and after intervention, because γ 1 = 0. Most importantly, the intervention effect is significant, because δ 4 = 1 is selected. Interestingly, the selected model is no longer a state space model (γ 1 = γ 2 = γ 3 = 0), but a simple Poisson regression model with monthly seasonal dummies and an intervention effect. This finding is confirmed by the marginal likelihoods computed in [START_REF] Frühwirth-Schnatter | Marginal likelihoods for non-Gaussian models using auxiliary mixture sampling[END_REF].

In Table 6 and Figure 10, we compare posterior inference for the intervention effect for the unconstrained basic structural model and the model obtained by variable selection.

We observe here an impressive gain of statistical efficiency for this parameter of interest.

For the unconstrained basic structural model, making the level dynamic before and after

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

the intervention causes quite a loss of information, leading to an intervention effect that is not significant.

The seasonal pattern disappears, if B 0 is rather large in the partially proper prior and if b is rather small in the fractional prior, see Table 4 and5

Purse snatching in Hyde Park, Chicago

For further illustration, we reanalyze a time series of cases of purse snatching y t in the Hyde park neighborhood in Chicago [START_REF] Harvey | Forecasting, Structural Time Series Models and the Kalman Filter[END_REF] reported for the period from January 1968 to September 1973. We consider a simplified version of the model introduced in Subsection 4.1, where no seasonal and no intervention effect is present and the exposures e t are equal to 1.

First, Gibbs sampling was run without variable selection for 15,000 iterations after a burn-in of 10,000. We selected the normal prior N (0, 10) both for µ 0 and a 0 and compare the inverted Gamma prior θ i ∼ G -1 (-0.5, 0.0001) with the normal prior ± √ θ i ∼ N (0, 1), for i = 1, 2. Figure 11 shows histograms of ± √ θ i under both priors. The posterior for ± √ θ 1 is roughly the same under both priors and clearly indicates that θ 1 > 0. Again, the inverted Gamma density is too influential for θ 2 and shrinks the draws away from 0, while for the normal prior the posterior is clearly centered at 0, suggesting that θ 2 = 0.

Second, stochastic model specification search is carried out using the prior described in Subsection 3.3. We compare partially proper priors with different prior variances B 0 with fractional priors with different fractions b. For each prior, MCMC sampling was carried out for M = 100, 000 draws after a burn-in of 20, 000 draws. The first 1000 draws of the burn-in were drawn from the unrestricted model, model selection began after these first 1000 draws. Depending on the prior, running MCMC takes between 25 and 28 minutes, see also Table 8.

Results of the variable selection procedure are presented in

ACCEPTED MANUSCRIPT

drift disappears because δ = 0 and γ 2 = 0 imply that a t ≡ a 0 = 0 for all t. This finding confirms model selection by the marginal likelihoods in [START_REF] Frühwirth-Schnatter | Marginal likelihoods for non-Gaussian models using auxiliary mixture sampling[END_REF].

Concluding remarks

The model space MCMC approach discussed in this paper could be easily adapted to other state space models. An important extension we are investigating currently is searching for fixed and time-varying coefficients in a regression model.

It is possible to extend our approach to several other non-Gaussian time series models.

The method works for any partial Gaussian state space model in the sense of [START_REF] Shephard | Partial non-Gaussian state space[END_REF], i.e. for any state space model that is conditionally Gaussian given a set of auxiliary latent variables. Thus it is straightforward to extend the approach to robust state space modeling based on the Student-t distribution as discussed e.g. in [START_REF] Carlin | A Monte Carlo approach to nonnormal and nonlinear state-space modeling[END_REF] or to modelling binary time series based on probit-link state space models as considered e.g. in [START_REF] Czado | State space mixed models for longitudinal observations with binary and binomial responses[END_REF].

Auxiliary mixture sampling as discussed in [START_REF] Frühwirth-Schnatter | Auxiliary mixture sampling with applications to logistic models[END_REF] allows to consider state space modelling of binary and categorical time series based on the logit link. [START_REF] Frühwirth-Schnatter | Improved auxiliary mixture sampling for hierarchical models of non-Gaussian data[END_REF] shows how to run auxiliary mixture sampling for data from the negative binomial distribution. In combination with the idea of the present paper, this allows variable selection for time series of overdispered count data. Furthermore, Frühwirth-Schnatter et al. ( 2009) discuss efficient auxiliary mixture sampling for data from the binomial and the multinomial distribution where the dimension of the auxiliary variable u t is limited to 2m, with m + 1 being the number of alternatives, even if the number y t of observed counts is increasing. This makes it feasible to extend variable selection to state space modelling of binomial and multinomial data.

A couple of modifications of our approach are worth being considered. First, the uniform prior over all models may be substituted by a more flexible prior which is obtained by assuming that the prior occurrence of δ i = 1 and γ i = 1 is different:

Pr(δ i = 1|α δ ) = α δ , Pr(γ i = 1|α γ ) = α γ . A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
(µ 0 ) ∝ 1, B 0 = 1 0 1 0 1 0 44825 1 1 1 0 0 16323 0 1 0 1 1 11803 p(µ 0 ) ∝ 1, B 0 = 10 0 1 0 1 0 64018 1 1 1 0 0 14267 0 1 1 0 0 10655 p(µ 0 ) ∝ 1, B 0 =
( √ θ 1 , √ θ 2 )/ max(l( √ θ 1 , √ θ 2 )), where l( √ θ 1 , √ θ 2 ) = p(y| √ θ 1 , √ θ 2 , σ 2
,tr ε , µ tr 0 , a tr 0 ) for simulated data with (θ tr 1 , θ tr 2 ) = (0.15 2 , 0.02 2 ) (first row), (θ tr 1 , θ tr 2 ) = (0.15 2 , 0) ( second row), (θ tr 1 , θ tr 2 ) = (0, 0.02 2 ) (third row), and (θ tr 1 , θ tr 2 ) = (0, 0) (last row) 38 1987 1989 1991 1993 1995 1997 1999 2001 2003 1987 1989 1991 1993 1995 1997 1999 2001 2003 1987 1989 1991 1993 1995 1997 1999 2001 2003 2005 -10 -9.5 -9 -8.5 -8 -7.5 (c) 1987 1989 1991 1993 1995 1997 1999 2001 2003 N (0, 1) prior for ± √ θ 1 and ± √ θ 2 ; bottom: G -1 (-0.5, 0.0001)-prior for θ 1 and θ 2

ACCEPTED MANUSCRIPT

  (a) Sample the indicators (δ, γ), the observation variance σ 2 ε and the initial values µ 0 , a 0 , and s 0 and all variance parameters √ Sample x = (x 1 , . . . , x T ) from the state space form corresponding to (32).(c) Perform two random sign switches as in step (c) in Subsection 2.4. Perform a third random sign switch for √ θ 3 and {s t } T 1 .

  . The most frequently visited model is fairly robust against the choice of the prior, only the partially proper prior with the largest B 0 and the fractional prior with the smallest b lead to a more parsimonious model. No trend is present in the selected model, because δ = 0

  . Increasing B 0 and decreasing b forces more parsimonious models. Not surprisingly the seasonal pattern disappears first in a more parsimonious model, because in the selected model only a few seasonal dummies are different from 0, see also Figure 9.

  Figure 1:Contour and surface plots of the (scaled) profile likelihoodl( √ θ 1 , √ θ 2 )/ max(l( √ θ 1 , √ θ 2 )), where l( √ θ 1 , √ θ 2 ) = p(y| √ θ 1 , √ θ 2 , σ 2,trε , µ tr 0 , a tr 0 ) for simulated data with (θ tr 1 , θ tr 2 ) = (0.15 2 , 0.02 2 ) (first row), (θ tr 1 , θ tr 2 ) = (0.15 2 , 0) ( second row), (θ tr 1 , θ tr 2 ) = (0, 0.02 2 ) (third row), and (θ tr 1 , θ tr 2 ) = (0, 0) (last row) 38

Figure 2 :Figure 3 :

 23 Figure 2: Prior and posterior densities for ± √ θ 1 under the prior distributions θ 1 ∼ G -1 (0.5, D 0 ) (left hand side) and ± √ θ 1 ∼ N (0, B 0 σ 2 ε ) (right hand side). Top: priors, middle: posterior, if true value of θ 1 is equal to θ tr 1 = 0.01; bottom: posterior, if true value of θ 1 is equal to θ tr 1 = 0

Figure 5 :

 5 Figure 5: UK coal consumption; MCMC draws for ± √ θ 1 (left), ± √ θ 2 (middle) and ± √ θ 3 (right) under different parameterizations; top: noncentered parameterization witha N (0, 1)-prior for ± √ θ i ; bottom: centered parameterization with a G -1 (-0.5, 10 -7 )prior for θ i

Figure 6 :

 6 Figure 6: Road safety data; (a) counts of killed or injured children, (b) number of children exposed

Figure 7 :Figure 8 :

 78 Figure 7: Road safety data; top: MCMC draws for ± √ θ 1 (left), ± √ θ 2 (middle) and ± √ θ 3 (right) for the noncentered parameterization; bottom: MCMC draws for θ 1 (left), θ 2 (middle) and ± √ θ 3 (right) under the partially noncentered parameterization

Figure 9 :Figure 10 :Figure 11 :

 91011 Figure 9: Road safety data; (a) posterior mean of the risk λ t ; posterior means and pointwise 95% credible regions of (b) the level µ t , (c) the drift a t and (d) the seasonal component s t in year 2005 under the centered parameterization

  Table 7 and 8. Model selection is extremely robust to the prior choice and clearly picks a local level model. The

Table 1

 1

	: UK coal consumption; the three most frequently visited models (among 100000
	MCMC iterations) for various prior distributions
	prior	δ δ 3 γ 1 γ 2 γ 3 frequency
	p

Table 2 :

 2 UK coal consumption; posterior inclusion probability for each indicator under

	various priors
	Prior	δ	δ 3	γ 1	γ 2	γ 3	time (min)
	p(µ 0 ) ∝ 1, B 0 = 1 0.2933 1.0000 0.3712 0.7500 0.2223 p(µ 0 ) ∝ 1, B 0 = 10 0.1716 1.0000 0.3136 0.7410 0.0328 p(µ 0 ) ∝ 1, B 0 = 100 0.0697 1.0000 0.2344 0.7836 0.0101 b = 10 -3 0.1827 1.0000 0.2154 0.9210 0.0677	52.2 51.6 52.0 43.8
	b = 10 -4	0.0630 1.0000 0.0736 0.9728 0.0216	43.6
	b = 10 -5	0.0219 1.0000 0.0600 0.9576 0.0069	43.2

Altenbergerstraße 69, A-4040 Linz, Austria; Phone: ++43 732

8295; e-mail: sylvia.fruehwirthschnatter@jku.at

A shorter version of this time series ranging from 1987-2002 was analyzed in Frühwirth-Schnatter and Wagner (2006).

ACCEPTED MANUSCRIPT

MCMC Estimation

MCMC estimation is implemented using auxiliary mixture sampling for count data. For each t, the distribution of y t |λ t is regarded as the distribution of the number of jumps of an unobserved Poisson process with intensity e t λ t , having occurred in the time interval [0,1]. [START_REF] Frühwirth-Schnatter | Auxiliary mixture sampling for parameterdriven models of time series of counts with applications to state space modelling[END_REF] create such a Poisson process for each observation y t and introduce the (y t + 1) interarrival times of this Poisson process as latent variables, yielding a total of T + T t=1 y t latent variables.

Frühwirth-Schnatter et al. ( 2009) suggest a more efficient method based on introducing for each observation y t at most two latent variables, namely the interarrival time τ t1 between the y t th jump and the next one and, if y t > 0, the arrival time τ t2 of the y t th jump. Since τ t1 ∼ E (e t λ t ) and τ t2 ∼ G (y t , e t λ t ) we have for j = 1, 1 + min(y t , 1):

where ε tj =log ξ tj and ξ tj ∼ G (ν tj , 1) with integer shape parameter equal to ν tj = max(1, (j -1)y t ). The distribution of the negative log Gamma distribution is approximated in [START_REF] Frühwirth-Schnatter | Improved auxiliary mixture sampling for hierarchical models of non-Gaussian data[END_REF] for each integer value ν by a mixture of normal distributions with component indicator r tj :

The number of components R(ν) depends on ν, as do the weights w r (ν), the means m r (ν) and the variances s 2 r (ν), see Frühwirth-Schnatter et al. (2009, Appendix A) for more details.

Introducing the auxiliary variables u = (u 1 , . . . , u T ), where u t = (τ tj , r tj , j = 1, 1 + min(y t , 1)), leads to a conditionally Gaussian state space model:

In (43) we are dealing with a state space model that is conditionally Gaussian with the state vector x t being the same as in Subsection 3.2. The MCMC scheme introduced in Subsection 3.2 for Gaussian state space models needs only a few modifications. First, an

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

additional step has to be added to draw the auxiliary variables u. Second, conditional on the state vector, we are dealing with a regression model with heteroscedastic normal errors with known error variance:

where ỹ denotes the collection of the auxiliary variables (-log

and Σ is a diagonal matrix with elements s 2 r tj (ν tj ). Under the normal prior

, the marginal likelihood in this regression model defines p(y|δ, γ, x, u):

The MCMC scheme reads:

(a1) Sample δ and γ from p(δ, γ|x, u, y) ∝ p(y|δ, γ, x, u)p(δ, γ) conditional on the state process x and the auxiliary variables u using the marginal likelihood (45) obtained from regression model ( 43).

(a2) Sample all unrestricted elements of the initial values of x 0 and all unrestricted variance parameters √ θ i jointly from the multivariate normal distribution N a δ,γ T , A δ,γ T conditional on x and u using the moments ( 46) and ( 47 (d1) Sample ξ t ∼ E (e t λ t ). If y t = 0, set τ t1 = 1 + ξ t . If y t > 0, sample τ t2 from a B (y t , 1)-distribution and set τ t1 = 1τ t2 + ξ t .

ACCEPTED MANUSCRIPT

In this prior, α δ and α γ may be chosen as fixed values, if prior information on the occurrence probabilities is available. If this is not the case, a hyperprior may be put on α δ and α γ as in [START_REF] Smith | Parsimonious covariance matrix estimation for longitudinal data[END_REF] and [START_REF] Frühwirth-Schnatter | Bayesian parsimonious covariance estimation for hierarchical linear mixed models[END_REF]. If both hyperparameters α δ and α γ are iid Uniform on [0,1], then p(δ, γ) = B(1

where B(•, •) is the Beta function. This prior leads to a uniform distribution over model size and outperforms the uniform prior over all models in variable selection for large regression models, see Ley and Steel (2007). In our applications, where model size is small, posterior inference under both priors is virtually the same.

Second, sampling the indicators could be modified. In our MCMC schemes, the indicators (δ, γ) are sampled jointly from the discrete posterior p(δ, γ|x, y) by evaluating the right hand side of ( 25) for all combinations of indicators at each sweep of the sampler. This multi-move sampling is rather time-consuming and may be substituted by single-move Gibbs sampling, i.e. sampling recursively from p(δ j |δ -j , γ, x, y) and p(γ j |γ -j , δ, x, y) as in [START_REF] George | Variable selection via Gibbs sampling[END_REF].

An open issue of our approach is the influence the prior on the initial values and the process variances exercises on final model selection. We demonstrated that the normal prior put on the signed square root of the process variances is far less influential than the usual inverted Gamma for the process variances themselves. The sensitivity analysis carried out for all of our case studies revealed a surprising robustness of the finally selected model against variation in the normal prior. A concise statement which prior scale leads to model consistency in the sense of [START_REF] Casella | Consistency of Bayesian procedures for variable selection[END_REF], however, is far beyond the scope of the present paper.

ACCEPTED MANUSCRIPT

Tüchler, R. ( 2008). Bayesian variable selection for logistic models using auxiliary mixture sampling. Journal of Computational and Graphical Statistics 17, 76-94. Zellner, A. (1971). An Introduction to Bayesian Inference in Econometrics. New York:

Wiley.

ACCEPTED MANUSCRIPT

ACCEPTED MANUSCRIPT


MD, PhD Theodoros Dimitroulas

MD Georgios Giannakoulas

MD Hariklia Dimitroula

MD Tilemahos Sfetsios

MD Despina Parcharidou

MD Haralambos Karvounis

MD Loukas Settas

Significance of serum uric acid in pulmonary arterial hypertension due to systemic sclerosis: a pilot study

Keywords: uric acid, systemic sclerosis, pulmonary arterial hypertension, six minute walking test

Systemic sclerosis (SSc) is a connective tissue disease, which may lead to elevated pulmonary arterial pressure due to pulmonary arterial hypertension and/or left ventricular diastolic dysfunction. Uric acid (UA) has been shown to be elevated in patients with pulmonary hypertension (PH) and heart failure. We aimed to investigate the potent relationship between serum UA and pulmonary pressure as well as functional capacity in patients with SSc. We studied 66 patients (mean age 57.7±12.1 years, 63 women), presenting with SSc. Systolic pulmonary artery pressure assessed by echocardiography, lung function tests, six-minute walk test (6MWT) and serum UA levels were recorded in all patients. In 24 (36%) patients the diagnosis of PH was established by echocardiography (systolic pulmonary artery pressure ≥40mmHg). PH patients had higher UA serum levels compared to patients without PH (5.1±2.1 mg/dl vs 4.2±0.9 mg/dl, p = 0.04). Amongst PH patients, UA values were inversely correlated with the SMWT distance (r = -0.51, p = 0.01). Serum UA values increased in proportion to the functional capacity PH in patients with scleroderma. Further investigations in prospective studies will unfold in detail the pathophysiological significance of UA in SSc patients with PH, and determine its role as a prognostic marker in the assessment and monitoring of the disease.

Introduction

Systemic sclerosis (SSc) is a connective tissue disease characterized by excessive collagen deposition in the skin and the internal organs and vascular involvement with generalized disturbance of the microcirculation [START_REF] Teh | Systemic sclerosis in Sarawak: a profile of patients treated in the Sarawak General Hospital[END_REF][START_REF] Toya | The many faces of scleroderma sine scleroderma: a literature review focusing on cardiopulmonary complications[END_REF][START_REF] Rosato | Regional diastolic function by tissue Doppler echocardiography in systemic sclerosis: correlation with clinical variables[END_REF]. Systemic sclerosis related pulmonary arterial hypertension (SScPAH) is a serious and devastating complication of the disease affecting 8-15% of SSc patients [START_REF] Hachulla | Early detection of pulmonary arterial hypertension in systemic sclerosis: a French nationwide prospective multicenter study[END_REF][START_REF] Mukerjee | Prevalence and outcome in systemic sclerosis associated pulmonary arterial hypertension: application of a registry approach[END_REF].

Pulmonary hypertension (PH) in these patients may also be caused by left ventricular diastolic dysfunction due to myocardial involvement [START_REF] Hachulla | Itinér AIR-Sclérodermie Study Group. The three-year incidence of pulmonary arterial hypertension associated with systemic sclerosis in a multicenter nationwide longitudinal study in France[END_REF]. The presence of elevated pressure in the pulmonary circulation conveys bad prognosis in this population irrespective of the etiology of PH [START_REF] Steen | Changes in causes of death in systemic sclerosis, 1972-2002[END_REF].

Uric acid (UA) is the final oxidation product of purine metabolism. Serum UA levels have been shown to be elevated in conditions of impaired oxidative metabolism such as chronic heart failure [START_REF] Leyva | Serum uric acid as an index of impaired oxidative metabolism in chronic heart failure[END_REF], and obstructive pulmonary disease [START_REF] Braghiroli | Continuous-wave Doppler echocardiographic assessment of severity of calcific aortic stenosis: a simultaneous Doppler-catheter correlative study in 100 adult patients[END_REF]. Elevated UA values have been described in various rheumatic diseases [10] including lupus nephritis [11] and gout [12]. Several studies have demonstrated that increased UA serum concentrations in patients with PH of various etiologies correlate with the severity and the prognosis of the disease [13-15].

There is very little information available on the clinical significance of UA in SSc patients with PH. Therefore, we hypothesised that 1) uric acid levels are elevated in patients with SSc and PH and 2) elevated pulmonary pressure are associated to functional capacity in this population.

Patients and Methods

Study population

This study enrolled 66 patients with SSc as defined by the revised ACR [16] criteria who visited our center for follow-up care between February 2007 and February 2008. Patients were selected after exclusion of hypertrophic, dilated or restrictive cardiomyopathy, left ventricular ejection fraction (LVEF) < 55%, any left and right sided valvular disease more than moderate degree severity, symptomatic lower limb vascular disease, severe interstitial lung fibrosis defined as functional volume capacity (FVC) < 60% of predicted, medical history of hypertension and renal failure determined by serum creatinine levels > 1.2 mg/dl. Smoker patients were also excluded, since smoking has been shown to affect pulmonary function in SSc patients [17]. None of the patients received treatment with diuretics which may affect UA serum levels. All patients were evaluated for cardiac and pulmonary involvement.

Exercise capacity was assessed by six minute walking test (6MWT). The study received ethical approval from the scientific committee of the Aristotle University of Thessaloniki and patients provided written informed consent.

Pulmonary function testing

Pulmonary involvement was assessed by high sensitivity computed tomography. Lung function tests measured FVC, forced expiratory volume in 1 second (FEV 1 ), and diffusing capacity for carbon monoxide (DLCO).

Echocardiography

Complete two-dimensional and Doppler echocardiographic examinations were performed in all patients at the same day of the UA evaluation. We used a commercially available system (Vivid 7, Vingmed, GE, Norway). From the apical 4chamber view, RV inflow view, and parasternal short axis or subcostal view, continuous-wave Doppler echocardiography was used to assess the peak tricuspid regurgitant velocity [18]. Systolic pulmonary artery pressure was estimated by calculating the maximal velocity of the tricuspid regurgitant jet and then, further using the Bernoulli equation, adding to this value an estimated right atrial pressure based on both the size of the inferior vena cava and the change in diameter of this vessel during respiration [START_REF] Nagueh | Relation of mean right atrial pressure to echocardiographic and Doppler parameters of right atrial and right ventricular function[END_REF]. PH was considered present when the systolic pulmonary artery pressure value ≥ 40 mmHg [START_REF] Mukerjee | Echocardiography and pulmonary function as screening tests for pulmonary arterialhypertension in systemic sclerosis[END_REF].

Exercise capacity

The 6MWT was performed using a standardized protocol in accordance to guidelines in the American Thoracic Society statement 2002 [START_REF]ATS statement: guidelines for the six-minute walk test[END_REF]. Patients walked unassisted an enclosed-level corridor; length to first turn around was 30 metres. All patients were told to use their own pace, but to cover as much ground as possible in 6 min.

Blood sampling for UA measurement

For the measurement of the serum UA, creatinine, triglyceride, total cholesterol, ESR and CRP levels, venous blood sample was drawn after an overnight fast, at the study entry. Serum UA levels were determined in close temporal automated chemistry analysers. This test is a colometric assay based on the oxidation of UA by the specific enzyme uricase.

Statistical analysis

Variables are presented as mean ± S.D, except from the duration of the disease, which is expressed as median [25th-75th percentile]. UA and 6MWT values were normally distributed (Kolmogorov-Smirnov and P-P plot). Categorical data are presented as absolute value, and comparisons were tested by Fisher's exact test.

Clinical and biological characteristics between the two groups of SSc patients were compared using the Student's t-test. The relationship between UA and patient's functional state as measured by the 6MWT and other laboratory or functional pulmonary indices was assessed by the Pearson's correlation coefficient. A value of P<0.05 was considered significant. Statistical analysis was performed using SPSS 13 for Windows (SPSS Inc., Chicago, USA).

Results

In 24 (36%) patients the diagnosis of PH was established (22 females, mean age 60 ± 10 years; 19 with limited cutaneous and 6 with diffuse cutaneous SSc). The clinical and biologic characteristics of SSc patients with and without elevated pulmonary pressure are summarized in Table 1. SSc patients with PH had higher UA values (5.1 ± 2.1 mg/dl) compared to patients with normal pulmonary pressure (4.2 ± 0.9 mg/dl, p = 0.04) (Figure 1). Amongst PH patients, UA values were inversely correlated with the SMWT distance (r=-0.51, p=0.01) (Figure 2). No significant correlations were found between UA and serum creatinine, FVC, FEV1 and DLCO values in the same population. Moreover, UA was not related with sex, subtype of the disease, autoantibodies and the clinical features of the disease. In the subgroup of patients without PH, no correlation between UA levels and exercise capacity as assessed by the 6MWT was detected.

Discussion

The mechanisms responsible for the elevated serum UA concentrations in SSc patients with PH may be multifactorial. Increased UA levels have been suggested to reflect raised xanthine oxidase activityan important source of oxygen free radicalsin patients with heart failure [START_REF] Leyva | Uric acid in chronic heart failure: a marker of chronic inflammation[END_REF]. The latter provides a pathophysiological link of UA with a variety of detrimental processes, including increased cytokine production, inflammation and endothelial dysfunction, which are considered to be among the primary events during the progression of the scleroderma induced vasculopathy and fibrosis. Superoxide dismutase represents the major antioxidant defense system against oxygen free radicals. Elevated levels of UA are correlated with decreased superoxide dismutase activity in patients with chronic heart failure [START_REF] Alcaino | Serum uric acid correlates with extracellular superoxide dismutase activity in patients with chronic heart failure[END_REF]. A recent study showed diminished superoxide dismutase activity in SSc patients [START_REF] Frent-Cornateanu | Antioxidant defense capacity in scleroderma patients[END_REF], which also supports the hypothesis that oxidative stress has an important role in the progression of vascular disease, leading to an increase in UA levels in SSc patients with PH.

In the current study we demonstrated inverse correlation between UA levels and functional status in patients with PH. The 6MWT is commonly used for the objective assessment of functional capacity and is used in the evaluation, management, and follow-up of patients with PH [START_REF]ATS statement: guidelines for the six-minute walk test[END_REF]. As SSc patients usually present with cutaneous leg ulcers or pulmonary interstitial involvement that may invalidate the 6MWT, we excluded patients with the aforementioned clinical entities. UA concentration has emerged as a significant predictor of cardiopulmonary exercise parameters and measures of diminished functional capacity in patients with chronic heart failure [START_REF] Leyva | Serum uric acid as an index of impaired oxidative metabolism in chronic heart failure[END_REF]26].

Measurement of serum UA levels is simple, inexpensive and routinely available. Although other biochemical markers including brain natriuretic peptide [27-30], asymmetric dimethylarginine [31] and vascular endothelial growth factor [32] have been proposed as biochemical markers for the detection of PH in SSc, UA may also serve as an additional indicator for the assessment of exercise capacity in this population.

The findings have to be considered against the background of the limits of the study. This is a pilot single centre cross-sectional study and it is limited by the lack of a controlled healthy population. The small number of patients may have contributed to the borderline difference of UA and SMWT values noticed between the two groups of patients (with and without PH). Another important point is the fact that the patients did not undergo right heart catheterization and the lack of correlation between UA levels and hemodynamic indices. Additionally, UA levels are dependent on age and sex, and are affected by the renal impairment and diuretic therapy; hence their interpretation may be difficult in some patients. We took particular care to exclude patients under diuretics, but additional investigations in larger cohorts of patients are warranted to confirm the findings of the present study.

In conclusion, we established that serum UA values increase in proportion to the clinical severity of PH in patients with SSc. Further investigations in prospective studies will unfold in detail the pathophysiological significance of UA in PH in this group, and determine its role as a biochemical marker in the assessment and monitoring of the disease.

  study in chronic obstructive pulmonary disease and obstructive sleep apnea before and after treatment with nasal continuous positive airway pressure. Am Rev Respir Dis 148:173-8 10. Sari I, Akar S, Pakoz B (2009) Hyperuricemia and its related factors in an urban population, Izmir, Turkey. Rheumatol Int. 29:869-74 11. Yang Z, Liang Y, Xi W, et al (2010) Association of serum uric acid with lupus nephritis in systemic lupus erythematosus. Rheumatol Int (In press) 12. De Leonardis F, Govoni M, Colina M, et al (2007) Elderly-onset gout: a review. Rheumatol Int 28:1-6 13. Nagaya N, Uematsu M, Satoh T, et al (1999) Serum uric acid levels correlate with the severity and the mortality of primary pulmonary hypertension. Am J Respir Crit Care Med 160:487-92 14. Bendayan D, Shitrit D, Ygla M, Huerta M, Fink G, Kramer MR (2003) Hyperuricemia as a prognostic factor in pulmonary arterial hypertension. Respir Med 97:130-3 15. Voelkel MA, Wynne KM, Badesch DB, Groves BM, Voelkel NF (2000) Hyperuricemia in severe pulmonary hypertension. Chest 117:19-24 16. Subcommittee for Scleroderma Criteria of the American Rheumatism Association Diagnostic and Therapeutic Criteria Committee (1980) Preliminary criteria ror the classification of the systemic sclerosis (scleroderma). Arthritis Rheum 23:581-90 17. Quadrelli SA, Molinari L, Ciallella LM, et al (2009) Patterns of pulmonary function in smoking and nonsmoking patients with progressive systemic sclerosis. Rheumatol Int 29:995-9 26. Bergamini C, Cicoira M, Rossi A, Vassanelli C (2009) Oxidative stress and hyperuricaemia: pathophysiology, clinical relevance, and therapeutic implications in chronic heart failure. Eur J Heart Fail 11:444-52 27. Dimitroulas T, Giannakoulas G, Karvounis H, Gatzoulis MA, Settas L (2010) Natriuretic peptides in systemic sclerosis-related pulmonary arterial hypertension. Semin Arthritis Rheum 39:278-84 28. Mukerjee D, Yap LB, Holmes AM, et al (2003) Significance of plasma Nterminal pro-brain natriuretic peptide in patients with systemic sclerosisrelated pulmonary arterial hypertension. Respir Med 97:1230-6 29. Allanore Y, Borderie D, Meune C, et al (2003) N-terminal pro-brain natriuretic peptide as a diagnostic marker of early pulmonary artery hypertension in patients with systemic sclerosis and effects of calcium-channel blockers. Arthritis Rheum 48:3503-8 30. Dimitroulas T, Giannakoulas G, Karvounis H, et al (2007) Neurohormonal activation in patients with systemic sclerosis-related pulmonary arterial hypertension. Int J Cardiol 121:135-7 31. Dimitroulas T, Giannakoulas G, Sfetsios T, et al (2008) Asymmetrical dimethylarginine in systemic sclerosis-related pulmonary arterial hypertension. Rheumatology 47:1682-5 32. Papaioannou AI, Zakynthinos E, Kostikas K, et al (2009) et al Serum VEGF levels are related to the presence of pulmonary arterial hypertension in systemic sclerosis. BMC Pulm Med 9:18

Figure 1 :

 1 Figure 1:

Figure 2 :

 2 Figure 2:

Table 1 :

 1 TABLESClinical and biochemical data of SSc patients serum creatinine, ESR: erythrocyte sedimentation rate, CRP: C-reactive protein, ANA: antinuclear antibodies, Scl-70 : anti-topoisomerase I antibodies, ACA: anticentromere antibodies, 6MWT: six-minute walk test

		All SSc patients	SSc patients without	SSc patients with	p
		(n = 66)	PH (n = 42)	PH (n = 24)
	Age (years)	57.7 ± 12.1	55.3 ± 13.1	60 ± 10	NS
	Females/Males	61/3	41/1	22/2	NS
	Duration of the disease (years)	10[0.2-31]	9.7[0.2-31]	11[2-25]	NS
	Limited/Diffuse SSc	47/19	28/14	19/5	NS
	Raynaud's	47 (71.2%)	30 (71.4%)	17 (70%)	NS
	FVC (% predicted)	86.8 ± 22.6	88 ± 25.4	83 18	NS
	FEV 1 (% predicted)	89.4 ± 23.9	89.8 26.8	88 20	NS
	DLCO (% predicted)	65.1 ± 22.1	65 22	62 21	NS
	EF (%)	65.5 ± 4.2	65.5 4.4	65.5 4.8	NS
	Creatinine (mg/dl)	0.73 ± 0.1	0.8 0.1	0.73 0.1	NS
	ESR	32.3 ± 2	30.5 21	34.5 20	NS
	CRP (mg/ml)	0.46 ± 0.1	0.5 0.3	0.4 0.2	NS
	ANA	58 (87.3%)	37 (88.1%)	19 (79.1%)	NS
	SCL-70	21 (31.8%)	13 (30.9%)	9 (37.5%)	NS
	ACA	18 (27.2%)	10 (23.8%)	8 (33.3%)	NS
	Uric Acid (mg/dl)	4.68 ± 1.84	4.27 ± 0.98	5.14 ± 2.18	0.04
	6MWT (metres)	405.4 ±106.4	423 98	383 113	0.04
	Categorical data are presented as absolute value (percentages).

FVC: forced vital capacity, FEV 1 : forced expiratory volume in 1 second, DLCO: diffusing capacity for carbon monoxide, EF: ejection fraction of left ventricular, sPAP: systolic pulmonary artery pressure, Cre:

Acknowledgements

Dr Giannakoulas has received funding from the Hellenic Atherosclerosis Society.

The authors have to declare no conflicts of interest.


MSc Nancy Broos
email: n.broos@lareb.nl

Eugène P Puijenbroek

MD Eugène P Van Puijenbroek

Interaction between topical miconazole and coumarins

Keywords: coumarin, cytochrome P450, interaction, INR, topical miconazole, cutaneous, topical

Interaction between topical miconazole and coumarins

The potential of miconazole to influence coumarin metabolism, which could result in overanticoagulation, has been well described for the oral and vaginal administration of this antifungal agent [1][2][3][4][5][6]. However, topical use of miconazole is thought not to be able to cause clinically relevant drug interactions because of limited systemic absorption [7]. We describe a series of cases where loss of coagulation control has occurred with simultaneous use of topical miconazole and a coumarin. These cases were reported to the Netherlands Pharmacovigilance Centre Lareb between April 2001 and July 2010. Lareb is responsible for the collection and analysis of reports of suspected adverse drug reactions in the Netherlands and covers a population of 16.5 million people [8]. In 2008, the defined daily doses (DDDs) in this population were 64,448,438 for coumarins and 230,070 for topical miconazole [9].

Patient A was an 80-year-old female who was treated with phenprocoumon for atrial fibrillation.

She had a history of congestive heart failure and diabetes mellitus. She was prescribed miconazole 20 mg/g cream for a fungal infection around the anal region. Thirteen days later her INR increased to 10.8. In the previous six months her INR varied from 1.7 to 4.8. The miconazole cream was withdrawn and she was treated with a single dose of 5 mg vitamin K after which her INR dropped back again.

Patient B was a 77-year-old female with congestive heart failure and diabetes mellitus. She was treated with phenprocoumon, losartan, simvastatin, metformin, glimepiride, carvedilol and furosemide. For an intertrigo under and on her breasts, she was prescribed miconazole 20 mg/g cream. Normally her INR ranged between 2.5 and 4. Twelve days after the start of the miconazole cream her INR increased to 20. She recovered after withdrawal of the miconazole and treatment with a total dose of 15 mg vitamin K.

Patient C was a 69-year-old female with a history of dementia, cerebrovascular accident and diabetes mellitus. She was treated with acenocoumarol, metoprolol, bisacodyl suppositories, metformin, glimepiride, fosinopril and chlorothalidone. She was prescribed miconazole 20 mg/g cream for an extensive mycosis. One week later her INR was increased to 12.1.

Patient D was a 59-year-old male with cardiac arrhythmia who was treated with acenocoumarol, sodium valproate, flecainide, sotalol and losartan. He was prescribed miconazole 20 mg/g cream for a fungal infection in his groins and armpits. Two weeks after initiation of miconazole cream his INR increased from his normal values, ranging between 3 and 4, to 5.9. His INR recovered after withdrawal of the miconazole and acenocoumarol dose reduction. In literature, we found one case report describing an interaction between topical miconazole and a coumarin. This report concerns the loss of coagulation control in a patient taking miconazole cream for flexural intertrigo. The patient had been stable for months on warfarin with an INR ranging between 2.2 and 3.1. After two weeks of applying topical miconazole in his right groin his INR increased to 21.4 [10]. Besides, Alexandra et al [11] described six cases of overanticoagulation with coumarin therapy in patients treated with topical econazole or bifonazole. These patients applied the drugs to the vulva under a disposable diaper or to their buttocks, groin or trunk.

The pharmacological mechanism by which miconazole can interfere with the metabolism of coumarins, the inhibition of CYP2C9, is widely known. Systemic absorption of miconazole is necessary to initiate this interaction. However, in clinical trials in healthy people the biological availability of topical miconazole was less than one percent [7]. Higher systemic absorption has only been proven after repeated administration of the drug in children with diaper rash [7,12].

Nevertheless, application of miconazole cream under occlusion, on large surfaces and close to mucous membranes could increase systemic absorption and therefore make patients more vulnerable for a systemic effect of miconazole [11]. As a matter of fact, all cases described in this letter were at risk for increased systemic absorption. Furthermore, in dermatomycosis and intertrigo skin lesions may occur, which may also enhance systemic concentrations. Moreover, other explanations that could explain the loss of anticoagulant control in these patients were absent.

In addition, genetic polymorphisms in the CYP2C9 gene might also increase the susceptibility to this interaction. CYP2C9 poor metabolizers have a lower activity of this enzyme, wherefore it can be expected that inhibition of CYP2C9 could occur at low miconazole plasma levels [13]. The potential of miconazole to influence coumarin metabolism, which could result in overanticoagulation, has been well described for the oral and vaginal administration of this antifungal agent [1][2][3][4][5][6]. However, topical use of miconazole is thought not to be able to cause clinically relevant drug interactions because of limited systemic absorption [7]. We describe a series of cases where loss of coagulation control has occurred with simultaneous use of topical miconazole and a coumarin. These cases were reported to the Netherlands Pharmacovigilance Centre Lareb between April 2001 and July 2010. Lareb is responsible for the collection and analysis of reports of suspected adverse drug reactions in the Netherlands and covers a population of 16.5 million people [8]. In 2008, the defined daily doses (DDDs) in this population were 64,448,438 for coumarins and 230,070 for topical miconazole [9].

Patient A was an 80-year-old female who was treated with phenprocoumon for atrial fibrillation.

She had a history of congestive heart failure and diabetes mellitus. She was prescribed miconazole 20 mg/g cream for a fungal infection around the anal region. Thirteen days later her INR increased to 10.8. In the previous six months her INR varied from 1.7 to 4.8. The miconazole cream was withdrawn and she was treated with a single dose of 5 mg vitamin K after which her INR dropped back again.

Patient B was a 77-year-old female with congestive heart failure and diabetes mellitus. She was treated with phenprocoumon, losartan, simvastatin, metformin, glimepiride, carvedilol and furosemide. For an intertrigo under and on her breasts, she was prescribed miconazole 20 mg/g cream. Normally her INR ranged between 2.5 and 4. Twelve days after the start of the miconazole cream her INR increased to 20. She recovered after withdrawal of the miconazole and treatment with a total dose of 15 mg vitamin K.

Patient C was a 69-year-old female with a history of dementia, cerebrovascular accident and diabetes mellitus. She was treated with acenocoumarol, metoprolol, bisacodyl suppositories, metformin, glimepiride, fosinopril and chlorothalidone. She was prescribed miconazole 20 mg/g cream for an extensive mycosis. One week later her INR was increased to 12.1.

Patient D was a 59-year-old male with cardiac arrhythmia who was treated with acenocoumarol, sodium valproate, flecainide, sotalol and losartan. He was prescribed miconazole 20 mg/g cream for a fungal infection in his groins and armpits. Two weeks after initiation of miconazole cream his INR increased from his normal values, ranging between 3 and 4, to 5.9. His INR recovered after withdrawal of the miconazole and acenocoumarol dose reduction. In literature, we found one case report describing an interaction between topical miconazole and a coumarin. This report concerns the loss of coagulation control in a patient taking miconazole cream for flexural intertrigo. The patient had been stable for months on warfarin with an INR ranging between 2.2 and 3.1. After two weeks of applying topical miconazole in his right groin his INR increased to 21.4 [10]. Besides, Alexandra et al [11] described six cases of overanticoagulation with coumarin therapy in patients treated with topical econazole or bifonazole. These patients applied the drugs to the vulva under a disposable diaper or to their buttocks, groin or trunk.

The pharmacological mechanism by which miconazole can interfere with the metabolism of coumarins, the inhibition of CYP2C9, is widely known. Systemic absorption of miconazole is necessary to initiate this interaction. However, in clinical trials in healthy people the biological availability of topical miconazole was less than one percent [7]. Higher systemic absorption has only been proven after repeated administration of the drug in children with diaper rash [7,12].

Nevertheless, application of miconazole cream under occlusion, on large surfaces and close to mucous membranes could increase systemic absorption and therefore make patients more vulnerable for a systemic effect of miconazole [11]. As a matter of fact, all cases described in this letter were at risk for increased systemic absorption. Furthermore, in dermatomycosis and intertrigo skin lesions may occur, which may also enhance systemic concentrations. Moreover, other explanations that could explain the loss of anticoagulant control in these patients were absent.

In addition, genetic polymorphisms in the CYP2C9 gene might also increase the susceptibility to this interaction. CYP2C9 poor metabolizers have a lower activity of this enzyme, wherefore it can be expected that inhibition of CYP2C9 could occur at low miconazole plasma levels [13].

The described cases indicate that topical miconazole use can cause loss of coagulation control when simultaneously used with coumarins. Clinicians should be aware of this possible interaction to prevent serious consequences in vulnerable patients. When simultaneous use is necessary, close monitoring of the INR is indicated during topical miconazole therapy.

1

 1 The described cases indicate that topical miconazole use can cause loss of coagulation control when simultaneously used with coumarins. Clinicians should be aware of this possible interaction to prevent serious consequences in vulnerable patients. When simultaneous use is necessary, close monitoring of the INR is indicated during topical miconazole therapy.Interaction between topical miconazole and coumarinsNancy Broos MSc 1 , Eugène P. van Puijenbroek MD, PhD1 1. Netherlands Pharmacovigilance Centre Lareb, Goudsbloemvallei 7, 's-Hertogenbosch, the NetherlandsCorresponding authorNancy Broos, telephone: 0031-736469700, fax: 0031-736426136, n.broos@lareb.nlInformative titleInteraction topical miconazole and coumarinsWord count (excluding references)

Formatted

  Where possible miconazole cream should be switched to a noninteracting antifungal agent or,

European Journal of Clinical Pharmacology

Acknowledgements

Nothing to declare. The cases described in this study have been reported to the Netherlands Pharmacovigilance Centre Lareb (patient A: NL-LRB-43273; patient B: NL-LRB-44960; patient C: NL-LRB-59743; patient D: NL-LRB-70136). Acknowledgements Nothing to declare. The cases described in this study have been reported to the Netherlands Pharmacovigilance Centre Lareb (patient A: NL-LRB-43273; patient B: NL-LRB-44960; patient C: NL-LRB-59743; patient D: NL-LRB-70136).

Conflicts of interest

Nothing to declare Keywords coumarin, cytochrome P450, interaction, INR, topical miconazole, cutaneous, topical Deleted: 684

Conflicts of interest

Nothing to declare