# Dynamax for Multiple Time Series problems

**URL:** <https://community.intuitivebayes.com/t/dynamax-for-multiple-time-series-problems/1056>\
**Category:** Dynamax\
**Created:** [April 21, 2023, 1:42am UTC](https://community.intuitivebayes.com/t/dynamax-for-multiple-time-series-problems/1056 "2023-04-21T01:42:00Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![jroberayalas](https://avatars.discourse-cdn.com/v4/letter/j/7feea3/32.png) [@jroberayalas](https://community.intuitivebayes.com/u/jroberayalas)\
**Post date:** [April 21, 2023, 1:42am UTC](https://community.intuitivebayes.com/t/dynamax-for-multiple-time-series-problems/1056/1 "2023-04-21T01:42:00Z")

</div>

Hey folks,

Just wondering if any of you has worked with multiple time series problems using the Dynamax library. I’m thinking on a problem where there are 100 time series with some trend and seasonality components (weekly, yearly, holidays) and you want to model them using a SSM while accounting for any correlation between them. The closest example I’ve seen is in the Time Series Forecasting [notebook](https://num.pyro.ai/en/stable/tutorials/time_series_forecasting.html) in the documentation of the numpyro library, but this is for a univariate time series problem. Any ideas / resources / examples are more than welcome!

---

<div class="post-metadata">

**Author:** ![RavinKumar](https://yyz1.discourse-cdn.com/flex007/user_avatar/community.intuitivebayes.com/ravinkumar/32/13_2.png) [@RavinKumar](https://community.intuitivebayes.com/u/RavinKumar)\
**Post date:** [April 21, 2023, 4:40am UTC](https://community.intuitivebayes.com/t/dynamax-for-multiple-time-series-problems/1056/2 "2023-04-21T04:40:11Z")

</div>

I havent used dynamax specifically but here’s an example with Tensorflow. If all the time series are independent it’d be quite nice to use a library like Dynamax as you could parallelize quite easily

[https://bayesiancomputationbook.com/markdown/chp\_06.html#state-space-models](https://bayesiancomputationbook.com/markdown/chp_06.html#state-space-models)

I’m curious, what are you working on?

---

<div class="post-metadata">

**Author:** ![jroberayalas](https://avatars.discourse-cdn.com/v4/letter/j/7feea3/32.png) [@jroberayalas](https://community.intuitivebayes.com/u/jroberayalas)\
**Post date:** [April 22, 2023, 12:07am UTC](https://community.intuitivebayes.com/t/dynamax-for-multiple-time-series-problems/1056/3 "2023-04-22T00:07:16Z")

</div>

Hey @RavinKumar! I’m working on something similar to the [M5 Forecasting](https://www.kaggle.com/c/m5-forecasting-accuracy) competition where I have hierarchical sales data at the item-store-region-date level. There are thousands of items, hundreds of stores and tens of regions. Forecasting is the main task, but equally important is understanding the correlations between time series, which led me to choose a Bayesian framework.

For this problem, I’ve experimented quite a bit with PyMC after reading @AlexAndorra’s (amazing) [post](https://alexandorra.github.io/pollsposition_blog/popularity/macron/hidden%20markov%20models/polls/2021/05/16/hmm-popularity.html) on estimating latent presidential popularity across time with a Markov chain. I started by aggregating the data to the region level. However, the moment I start going to deeper hierarchical levels, PyMC does not scale well, and most of the time, the model does not sample at all. I checked Sandra’s Scalable Bayesian Modelling [repo](https://github.com/symeneses/SBM) where I tried to leveraged the power of GPUs, but I kept running into Out of Memory errors. I also created a [post](https://discourse.pymc.io/t/pymc-model-on-large-datasets/11891) on the PyMC forum, trying to learn more from the community when it comes to scaling PyMC models to large datasets, but it seems this is an area of active development.

Overall, I don’t have many features (global level, item effect, region effect, weekly seasonality, yearly seasonality, special days effect). After attending the Dynamax book club, I started thinking that maybe I could reframe the problem using SSM. I know I could use a deep learning model, but the interpretable and adaptive nature of SSMs makes them more interesting. My first attempt came from the Time Series Forecasting [notebook](https://num.pyro.ai/en/stable/tutorials/time_series_forecasting.html) I mentioned before, where they implement an SSM model called [SGT](https://cran.r-project.org/web/packages/Rlgt/vignettes/GT_models.html) (Seasonal Global Trend). I managed to extend it to handle multiple correlated time series (check [here](https://forum.pyro.ai/t/expanding-sgt-model-to-multiple-time-series/3925/7)) but this model can only handle one type of seasonality at a time and I wasn’t able to break this constraint. Then I found Pyro’s [tutorials](http://pyro.ai/examples/forecasting_ii.html) on forecasting with SSMs. I’m currently checking them, and I’m surprised to see they can handle multiple time series that are modest in size, with the model running in a matter of minutes. They are using stochastic variational inference (SVI) and it seems this is related to what @ckrapu mentioned in this [post](https://discourse.pymc.io/t/what-are-the-differences-between-nuts-and-advi/4890/2?u=jroberayalas) when he had to train a big model.

> [@ckrapu](#):
>
> used ADVI + GPU to train deep convolutional autoencoders with 10 million+ parameters and Bayesian regularization using the minibatched implementation of ADVI within PyMC3.

So going back to the book club and Dynamax, I’m trying to learn how to leverage this library as well as JAX to tackle this interesting problem using SSMs. Still, there are lots of things to learn and try. I would be happy to hear your thoughts.

---

<div class="post-metadata">

**Author:** ![theorashid](https://yyz1.discourse-cdn.com/flex007/user_avatar/community.intuitivebayes.com/theorashid/32/447_2.png) [@theorashid](https://community.intuitivebayes.com/u/theorashid)\
**Post date:** [April 22, 2023, 7:17pm UTC](https://community.intuitivebayes.com/t/dynamax-for-multiple-time-series-problems/1056/4 "2023-04-22T19:17:15Z")

</div>

Hey @jroberayalas, I’ve been thinking about using dynamax for M5 data too. I think this is better for [sts-jax](https://github.com/probml/sts-jax) rather than pure dynamax. I was planning on contributing this as an example, but I’d be really keen to collaborate if you’re interested. Happy to set up a quick call.

Following the notation of the README, we have

y\_t = H\_t z\_t + u\_t + \epsilon\_t, \qquad \epsilon\_t \sim \mathcal{N}(0, \sigma^2\_t)

z\_{t+1} = F\_t z\_t + R\_t \eta\_t, \qquad \eta\_t \sim \mathcal{N}(0, Q\_t)

We can make some simplifying assumptions:

- \sigma\_t = \sigma\_h static measurement noise
- Q\_t = \sigma\_q^2. This is now a _vector_ of disturbance, the same length as the number of states \text{length}(z\_t) = N\_z, so each state has independent disturbances (no covariance in the AR effects). And as there are the same number of states as disturbances, \eta\_t is the same length as z\_t so R\_t is identity.
- F\_t z\_t= f \cdot z\_t, where f is the diagonal of some matrix F (no cross-correlations between states). This should save us having to invert F in Kalman filter, because inverting a diagonal is easy. We can either set all elements of f to 1 (perfect correlation between timesteps) or estimate the correlation of each effect.
- u\_t = x\_t^T \beta: regression component from external inputs.

Let’s say, following the M5 context, we have N\_s stores and N\_t timesteps. Let’s say each store has it’s own AR effect, and each store has a seasonality effect. Using seasonal dummies, we can have weekly seasonality. So have N\_z = N\_s (1 + 6), because we only need 6 dummies for a seasonality of 7.

Putting this together, and using bold for vectors and matrices because I prefer that, we have

\mathbf{y\_t} = \mathbf{H\_t} \mathbf{z\_t} + \mathbf{u\_t} + \mathbf{\epsilon\_t}, \qquad \mathbf{\epsilon\_t} \sim \mathcal{N}(0, \sigma\_h^2)

\mathbf{z\_{t+1}} = \mathbf{f} \cdot \mathbf{z\_t} + \mathbf{\eta\_t}, \qquad \mathbf{\eta\_t} \sim \mathcal{N}(0, \mathbf{\sigma\_q^2})

\mathbf{u\_t} = \mathbf{X\_t} \mathbf{\beta}

And the form of z\_t is something like

\begin{bmatrix} \text{AR effect, store 1} \\ \text{seasonal effect 1, store 1} \\ ... \\ \text{seasonal effect 6, store 1} \\ ... \\ \text{AR effect, store } N\_s \\ \text{seasonal effect ,1 store } N\_s \\ ... \\ \text{seasonal effect 6, store } N\_s \end{bmatrix}

Then \mathbf{H\_t}, which has shape N\_s \times N\_z, sums the states to the observation dimension \text{length}(y\_t) = N\_s, is something like (for N\_s = 2, and the day of the week is such that the third seasonal component is selected)

\begin{bmatrix} 1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \end{bmatrix}

The 1st column means we select store 1’s AR component. Col 4 is the same as we are selecting the 3rd seasonal effect for store 1. Col 8 means we select store 2’s AR component. Col 11 means we select the 3rd seasonal effect for store 2.

At the moment, this is all independent except for the common \sigma\_h. But we can extend this to share information between stores in a few ways:

- Common AR and seasonal effects between stores
- hierarchical priors on \mathbf{\sigma\_q}, \beta
- Allowing covariance between disturbances (matrix Q), or cross-correlations between AR steps (matrix \textbf{F})

Although I don’t know how easy it is to control priors in dynamax as it isn’t a PPL. Does anyone have experience with this?

Let me know what you all think, but I’d been keen to contribute this.
