# State space models with brms

**URL:** <https://discourse.mc-stan.org/t/state-space-models-with-brms/4087>\
**Category:** brms\
**Created:** [May 3, 2018, 9:15pm UTC](https://discourse.mc-stan.org/t/state-space-models-with-brms/4087 "2018-05-03T21:15:32Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![jebyrnes](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/jebyrnes/32/2301_2.png) [@jebyrnes](https://discourse.mc-stan.org/u/jebyrnes)\
**Post date:** [May 3, 2018, 9:15pm UTC](https://discourse.mc-stan.org/t/state-space-models-with-brms/4087/1 "2018-05-03T21:15:32Z")

</div>

Hello, all! I’ve been futzing around with state space models recently and have been trying to figure out their proper implementation.

So, consider, for example, the Nile dataset.

```
library(fpp2)
library(forecast)

#from base
fitNile <- StructTS(Nile, "level")

cat("Transitional variance:", fitNile$coef["level"],
    "\n", "Observational variance:", fitNile$coef["epsilon"],
    "\n", "Initial level:", fitNile$model0$a, "\n")

autoplot(Nile) +
  geom_line(y=fitted(fitNile), col="red", lty=2) +
  ggtitle("Fit model from State Space Model")

```

Works great! Now, I’m trying to figure out the proper implementation in brms. Using

```

library(brms)
library(tidybayes)
#n.b. tidybayes is amazing - from
#https://github.com/mjskay/tidybayes

Nile_df <- as.data.frame(Nile) %>%
  mutate(t = 1871:1970) %>%
  mutate(x_cent = x-mean(x)) #centering helps

#not sure if this works - trying
#to find a better example....
nile_brms <- brm(x_cent ~ 1,
                 autocor = cor_bsts(~ t),
                 data = Nile_df)

```

seems to produce OK parameter estimates (although still not quite at good convergence). But, I only get a fitted mean

```

#plot
nile_results <- Nile_df %>%
  add_fitted_samples(nile_brms) 

#show the fit and CIs
ggplot(nile_results, aes(x=t)) +
  stat_lineribbon(aes(y = estimate+mean(x)), .prob = c(.99, .95, .8, .5)) +
  geom_line(aes(y = x)) +
  scale_fill_brewer()

```

shows this quite well.

Thoughts - have any of you done something similar\>

---

<div class="post-metadata">

**Author:** ![paul.buerkner](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/paul.buerkner/32/3303_2.png) [@paul.buerkner](https://discourse.mc-stan.org/u/paul.buerkner)\
**Post date:** [May 4, 2018, 6:16am UTC](https://discourse.mc-stan.org/t/state-space-models-with-brms/4087/2 "2018-05-04T06:16:27Z")

</div>

The problem is that the `cor_bsts()` is ignored when predicting with newdata currently (as indicated by the warning you see). Because of this and considerable convergence problems, I recently deprecated this correlation structure. See the end of [https://github.com/paul-buerkner/brms/issues/148](https://github.com/paul-buerkner/brms/issues/148) for more details.

To get predictions including the BSTS parameters you have to run them on the existing data via `fitted(nile_brms)`. Not sure if tidybayes has an option for this.

---

<div class="post-metadata">

**Author:** ![paul.buerkner](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/paul.buerkner/32/3303_2.png) [@paul.buerkner](https://discourse.mc-stan.org/u/paul.buerkner)\
**Post date:** [May 4, 2018, 6:22am UTC](https://discourse.mc-stan.org/t/state-space-models-with-brms/4087/3 "2018-05-04T06:22:07Z")

</div>

Here is how you can get an appropriate plot

```
library(tidyverse)
pred <- cbind(Nile_df, fitted(nile_brms))
names(pred) <- make.names(names(pred))
pred %>%
  ggplot(aes(t, Estimate)) +
  geom_smooth(
    aes(ymin = X2.5.ile, ymax = X97.5.ile),
    stat = "identity"
  ) + 
  geom_point(aes(t, x_cent), inherit.aes = FALSE)

```

---

<div class="post-metadata">

**Author:** ![paul.buerkner](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/paul.buerkner/32/3303_2.png) [@paul.buerkner](https://discourse.mc-stan.org/u/paul.buerkner)\
**Post date:** [May 4, 2018, 6:37am UTC](https://discourse.mc-stan.org/t/state-space-models-with-brms/4087/4 "2018-05-04T06:37:22Z")

</div>

You could alternatively use splines or Gaussian processes:

```
nile_brms2 <- brm(x_cent ~ s(t), data = Nile_df)
nile_brms3 <- brm(x_cent ~ gp(t), data = Nile_df, chains = 1) # takes very long

```

They should work with tidybayes out of the box.

---

<div class="post-metadata">

**Author:** ![jebyrnes](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/jebyrnes/32/2301_2.png) [@jebyrnes](https://discourse.mc-stan.org/u/jebyrnes)\
**Post date:** [May 9, 2018, 1:18pm UTC](https://discourse.mc-stan.org/t/state-space-models-with-brms/4087/5 "2018-05-09T13:18:49Z")

</div>

Thanks - those work, but produce fairly different answers as compared to likelihood fits of state space models - and look at the data conceptually in a different way - than state space models. I’ll try bsts - [http://www.unofficialgoogledatascience.com/2017/07/fitting-bayesian-structural-time-series.html](http://www.unofficialgoogledatascience.com/2017/07/fitting-bayesian-structural-time-series.html) - or other solution. Hrm. I feel like there should be a way to incorporate the kalman filter somehow (e.g. [https://jrnold.github.io/ssmodels-in-stan/](https://jrnold.github.io/ssmodels-in-stan/) ) but I am not a STAN wizard. It’s why I use brms! Thanks, though!

---

<div class="post-metadata">

**Author:** ![roualdes](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/roualdes/32/1594_2.png) [@roualdes](https://discourse.mc-stan.org/u/roualdes)\
**Post date:** [May 9, 2018, 3:38pm UTC](https://discourse.mc-stan.org/t/state-space-models-with-brms/4087/6 "2018-05-09T15:38:43Z")

</div>

Here’s the in-progress Stan code I’ve been playing with to fit Bayesian structural time series, at least the local level model:

```
data {                                                                                                
  int<lower=0> T;                                                                                     
  vector[T] y;                                                                                        
}                                                                                                     
transformed data {                                                                                    
  real sd_y;                                                                                          
}                                                                                                     
parameters {                                                                                          
  vector[T] u_err;                                                                                    
  real<lower=0> u_tau;                                                                                
                                                                                                      
  real<lower=0> y_std;                                                                                
  real<lower=0> y_tau;                                                                                
}                                                                                                     
transformed parameters {                                                                              
  vector[T] u;                                                                                        
  real y_err;                                                                                         
                                                                                                      
  u[1] = y[1] + u_tau * u_err[1];                                                                     
  for (t in 2:T) {                                                                                    
    u[t] = u[t-1] + u_tau * u_err[t-1];                                                               
  }                                                                                                   
  y_err = y_tau * y_std;                                                                              
}                                                                                                     
model {                                                                                               
  // priors                                                                                           
  u_err ~ normal(0, 1);                                                                               
  u_tau ~ exponential(1); // divide by sd_?                                                                  
  y_std ~ normal(0, 1);                                                                               
  y_tau ~ exponential(1 / sd_y);                                                                      
                                                                                                      
  // likelihood                                                                                       
  y ~ normal(u, y_err);                                                                               
}                                                                                                     
generated quantities {                                                                                
  vector[T] y_pred;                                                                                   
  for (t in 1:T)                                                                                      
    y_pred[t] = normal_rng(u[t], y_err);                                                              
}  

```

In my opinion, the above model is still missing a reasonable default scaling for the level, `u_tau`. I’d appreciate some thoughts here, if anybody has them. Thanks.