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
 A: 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
\begin{equation}
y_{t+1} = (1-\rho)\,\mu + \rho\,y_t + \varepsilon_{t+1} ,
\end{equation}
where $\mu \sim \textsf{N}(m,s^2)$, where $s$ is the standard deviation. I suggest expressing this instead as follows:
\begin{equation}
y_{t+1} = \alpha + \rho\,y_t + \varepsilon_{t+1} ,
\end{equation} 
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:
\begin{equation}
\mu \sim \textsf{N}\left(\frac{m}{1-\rho}, \frac{s^2}{(1-\rho)^2}\right) .
\end{equation}
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$.  
