I’ve been playing around with a Stan version of the semilocal linear trend model, as found in the R package bsts. The model both automatically scales and matches bsts’s fit fairly well. Still, I’m hoping to get some feedback about:

- starting points
`mu[1]`

and`nu[1]`

, - choice of priors (I offer my thoughts in line), and
- possible ways to eliminate the few and sporadic divergences, as the simple solution of increase adapt_delta is already used.

```
data {
int <lower=1> T;
vector[T] y;
}
transformed data {
real sd_y = sd(y);
}
parameters {
real<lower=0> sigma_y;
vector[T] gamma;
real<lower=0> sigma_gamma;
vector[T] zeta;
real<lower=0> sigma_zeta;
real eta;
real<lower=-1, upper=1> phi;
}
transformed parameters {
vector[T] mu;
vector[T] nu;
mu[1] = y[1] + sigma_gamma * gamma[1];
nu[1] = zeta[1];
for (t in 2:T) {
mu[t] = mu[t-1] + nu[t-1] + sigma_gamma * gamma[t];
nu[t] = eta + phi * (nu[t-1] - eta) + sigma_zeta * zeta[t];
}
}
model {
// priors
//sigma_y ~ exponential(1 / sd_y);
sigma_y ~ normal(0, 2.5 * sd_y);
gamma ~ normal(0, 1);
// sigma_gamma ~ gamma(2, 1 / sd_y);
sigma_gamma ~ normal(0, 2.5 * sd_y);
zeta ~ normal(0, 1);
// sigma_zeta ~ gamma(2, 1 / sd_y);
sigma_zeta ~ normal(0, 2.5 * sd_y);
eta ~ normal(0, sd_y);
phi ~ normal(0, 0.5);
// likelihood
y ~ normal(mu, sigma_y);
}
generated quantities {
vector[T] y_pred;
for (t in 1:T)
y_pred[t] = normal_rng(mu[t], sigma_y);
}
```

If you’re interested, here’s some R code to get you started. The code below compares some features from Stan and bsts fits across two data sets, one from the `datasets`

package and one from `bsts`

.

```
library(rstan)
options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)
library(bsts)
sm <- stan_model('semilocal_linear_trend.stan')
## Nile
y <- as.vector(Nile)
T <- length(y)
x <- 1:T
## rsxfs
data(rsxfs)
y <- as.vector(rsxfs)
T <- length(y)
x <- 1:T
### Stan
sfit <- sampling(sm, data=list(y = y, T = T),
control=list(adapt_delta=0.99,
max_treedepth=15))
### bsts
m <- list()
m <- AddSemilocalLinearTrend(m, y=y)
bfit <- bsts(y, m, niter=2000)
### plot
plot(x, y, type='l')
lines(x, colMeans(bfit$state.contributions), col='red')
lines(x, colMeans(extract(sfit, "y_pred")$y_pred), col='blue')
plot(density(bfit$trend.slope.mean), col='red')
lines(density(extract(sfit, "eta")$eta), col='blue')
summary(bfit$trend.slope.mean)
summary(extract(sfit, "eta")$eta)
```

Thanks in advance for any and all comments.