0
votes

I am using JAGS to run Bayesian analyses for ARMA models. My data is simulated data, so I know the results.

So far, if I estimate (for example) a stationary AR(1) process, I get good results for the autoregressive parameter.

Now, I have a time series which has a unit root after half of the observations. So 1:500 stationary AR(1) (with autoregressive parameter 0.3), 500:1000 unit root. My goal is to get a density which has mass on both the value for the stationary autoregressive parameter (which would be for example 0.3) and mass on the unit root value (so around 1). I want to show that a unit root is in a part of the time series...

My expectation was that if I use a noninformative uniform prior for the autoregressive parameter like rho1~dunif(-10,10), I should get mass on both values. What really happens is that I just get mass on a value in between (at around 0.6).

  • Should I use a different prior for the autoregressive term? What other non-centered possibilities do I have?
  • How is it possible that the GIBBS sampler runs through the stationary and non-stationary part but a histogram plots no observations around 0.3 (stationary ar-parameter) and around the unit root?

*edit:

It's a bit difficult to post the code, because it is both in R and JAGS. The following is the JAGS model. I use this JAGS model to estimate the following time series with 1000 observations: 1:500 AR(1) process: y= alpha + rho1*y[i-1] with rho1=0.2 and alpha=0. For 501:1000 the time series has a unit root (random walk).

    model  {
    for (i in 2:N)  
{   
y[i]~dnorm(f[i],tau)
f[i] <- alpha + rho1*y[i-1]
    }
rho1~dunif(-10,10)
tau~dgamma(0.001,0.001)
alpha~dnorm(0,0.001)
}
1
It would be helpful if you posted your model code, or better yet, code and toy data to reproduce your problem (i.e. a minimal working example).Jacob Socolar

1 Answers

1
votes

You defined one parameter, but you want it to take either of two separate values, depending on the situation. This isn't how parameters work. The Gibbs sampler doesn't "run through" the stationary and non-stationary parts separately. It looks at the whole model, and realizes that a value of 0.3 would be exceedingly unlikely given the nonstationary part, while a value of 1 would be similarly unlikely given the stationary part. So it samples from the posterior distribution given all the data, and the best it can do is something in between 0.3 and 1.

You need two separate parameters, one to use in the stationary part and one in the nonstationary part. Since you know a priori which parts are which, this should be trivial to code:

model  {
for (i in 2:500)  
{   
y[i]~dnorm(f[i],tau)
f[i] <- alpha + rho1*y[i-1]
}
for (i in 501:N)  
{   
y[i]~dnorm(f[i],tau)
f[i] <- alpha + rho2*y[i-1]
}
rho1~dunif(-10,10)
rho2~dunif(-10,10)
tau~dgamma(0.001,0.001)
alpha~dnorm(0,0.001)
}