Bayesian autoregressive model with second peak at 1 in posterior distirbution of AR parameter

I am trying to run a Bayesian hierarchical AR1 model for a set of fairly short time series. In some of the series I get a second peak around 1 in the posterior distribution of the AR1 parameter.

Here is an example of a non-hierarchical AR1 model on a single series from my data which replicates the behavior

x <- c(47219, 48918, 53911, 56678, 51087, 59203, 78238, 67749, 68439, 63355, 61659, 67342,
65799, 63253, 48673, 66175)
d <- list(N = length(x),
Y = log(x))

cat("
model {
mu ~ dnorm(0, 0.01);
tau.pro ~ dgamma(0.001,0.001);
sd.pro <- 1/sqrt(tau.pro);
phi ~ dnorm(0, 1);

predY[1] <- Y[1];
for(i in 2:N) {
predY[i] <- mu + phi * (Y[i-1] -mu);
Y[i] ~ dnorm(predY[i], tau.pro);
}
}", file = "m0.txt")

library(dclone)
mcmc0 <- jags.fit(data = d, model = "m0.txt", params = c("mu", "phi", "tau.pro"),
n.iter = 1e04, n.chain = 1)
library(ggmcmc)
theta <- ggs(mcmc0)
ggs_density(D = theta)


with traceplots...

ggs_density(D = theta)


Is there anyway to tweak the model to get rid of this behavior? I have tried truncating the prior of mu, but this an impractical approach for the hierarchical model where there are many series with a wide range of means

• Given your parameterization, mu should be the mean of the data if the absolute value of phi is less than one. The sample mean of your data is more than 60,000, but your "estimate" of mu is near 10, so that doesn't make any sense. I believe your prior for mu is centered at zero with a standard deviation of 10. That puts the sample mean 6,000 standard deviations away from the "location" of your prior. If I'm correct, then you have a very strong prior and it's at odds with the data. But even if you fix that, I'm guessing you will still have problems. I'll try to get back to this later. – mef Feb 22 '18 at 15:38
• @ref I log the data before it goes to the model, so the mean looks about right – gjabel Feb 22 '18 at 20:45
• Oops! I see that now. – mef Feb 22 '18 at 23:09

The peak can be eliminated by using a different prior for $\mu$. The simplest way to implement the new prior is to change the parameterization. Currently, you have $$y_{t+1} = (1-\rho)\,\mu + \rho\,y_t + \varepsilon_{t+1} ,$$ where $\mu \sim \textsf{N}(m,s^2)$, where $s$ is the standard deviation. I suggest expressing this instead as follows: $$y_{t+1} = \alpha + \rho\,y_t + \varepsilon_{t+1} ,$$ where $\alpha \sim \textsf{N}(m,s^2)$. This can be estimated using your software with a minimal change in your code. For example, you could set $m = 0$ and $s = 10$.
Given the prior for $\alpha$, the implicit prior for $\mu$ can be computed by a change of variables: $$\mu \sim \textsf{N}\left(\frac{m}{1-\rho}, \frac{s^2}{(1-\rho)^2}\right) .$$ Clearly this prior depends on $\rho$. Trying to use this prior directly with the original parameterization may cause problems due to numerical instability. That's why it's a good idea to switch parameterizations.
By the way, the "meaning" of $\mu$ changes when $|\rho|\ge 1$. For example, when $\rho = 1$, there is no mean and $\mu$ instead represents a time trend. It's not clear to me that you really want or need to allow for $|\rho|\ge 1$. I often impose the restriction $\rho \in [0,1)$. However, this can be considered a separate issue from that of the prior for $\mu$.