How to represent this structural time series model efficiently?


I’m new to PyMC3 and trying to write a fairly simple time series model without having to resort to loops. I’m trying to implement the structural time series model from “Machine Learning: A Probabilistic Perspective” by Murphy:


where y_t is the observed value, a_t is “local level”, and b_t is “local linear trend.”

To my understanding, loops are inefficient in PyMC3, so i’m trying to avoid them. I know I can write b_t as a GaussianRandomWalk, but then I’m at a loss for how to proceed from there. I can write

a_t - a_{t-1} ~ N(b_{t-1}, Q)

so that I can take a cumulative sum and obtain the time series, but the problem is that it doesn’t look like PyMC3 supports the cumulative sum operation according to its math docs. Is there a way to write this efficiently? Thanks!


You can model this using GaussianRandomWalk.

e.g. this very partial sketch:

b = GaussianRandomWalk(‘b’, sd=epsilon_b)
a = GaussianRandomWalk(‘a’, mu=b, sd=epislon_a)
y = pm.Normal(‘y’, mu=a, sd=epsilon_y)