Hello, I was trying to implement a simple local linear trend model, defined as follow:
y_t = x_t + \epsilon_t
\mu_{t+1} = \mu_t + v_t + \xi_t
v_{t+1} = v_t + \zeta_t
I used the same code provided by (RPubs - 状態空間時系列分析入門 第3章 (1/3))
which I copy below:
data {
int<lower=1> n;
vector[n] y;
}
parameters {
vector[n] mu;
vector[n-1] v;
real<lower=0> sigma_level;
real<lower=0> sigma_drift;
real<lower=0> sigma_irreg;
}
transformed parameters {
vector[n] yhat;
for(t in 1:n) {
yhat[t] <- mu[t];
}
}
model {
for(t in 2:n-1)
v[t] ~ normal(v[t-1], sigma_drift);
for(t in 2:n)
mu[t] ~ normal(mu[t-1] + v[t-1], sigma_level);
for(t in 1:n)
y[t] ~ normal(yhat[t], sigma_irreg);
}
I fitted the model on the Nile data with RStan, but I get quite poor fitting:
llt <- stan_model("fig03_01.stan")
fit <- sampling(llt, data=list(y=Nile, n=length(Nile)), chains=3, cores=3)
Warning messages:
1: There were 4 divergent transitions after warmup. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.
2: There were 2000 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
http://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded
3: There were 3 chains where the estimated Bayesian Fraction of Missing Information was low. See
http://mc-stan.org/misc/warnings.html#bfmi-low
4: Examine the pairs() plot to diagnose sampling problems
5: The largest R-hat is 1.83, indicating chains have not mixed.
Running the chains for more iterations may help. See
http://mc-stan.org/misc/warnings.html#r-hat
6: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
Running the chains for more iterations may help. See
http://mc-stan.org/misc/warnings.html#bulk-ess
7: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
Running the chains for more iterations may help. See
http://mc-stan.org/misc/warnings.html#tail-ess
The sigma for the varying slope is especially hard to estimate:
print(fit, pars=grep("sigma.*", names(fit), value = T), digits=5)
Inference for Stan model: fig03_01.
3 chains, each with iter=2000; warmup=1000; thin=1;
post-warmup draws per chain=1000, total post-warmup draws=3000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
sigma[1] 0.43861 0.25419 0.55280 0.03073 0.10549 0.26646 0.51332 2.20960 5 1.53057
sigma[2] 84.86377 0.32654 7.39183 69.03075 80.43926 84.85241 89.78151 98.74073 512 1.01249
sigma[3] 96.38806 0.39164 9.57975 81.67447 89.60247 95.25685 101.21278 118.55008 598 1.00298
This usually replicate on any dataset I fit this model on.
any ideas on how to improve this?
thank you
