Scaling GaussianRandomWalk to multiple timeseries

Questions
1

I have a collection of timeseries that I am trying to model using GaussianRandomWalk. I had success modeling a single timeseries, as follows:

with pm.Model() as model:
    
    alpha = .1
    sigma = pm.HalfCauchy('sigma', 10)
    mu = pm.GaussianRandomWalk('mu', 
                              sigma=sigma * (1. - alpha), 
                              shape=len(y)
                              )
    likelihood = pm.Normal('timeseries', 
                      mu=mu, 
                      sigma=sigma * alpha, 
                      observed=y
                          )

However, when I try to scale it to the entire dataset by building a hierarchical model, the sampling fails:

y = data # (n series x D time periods)
s = (y.shape[0], 1)

with pm.Model() as model:
    
    # hyperparameter
    alpha = .1
    
    # sigma hyperprior
    sigma_sd = pm.HalfCauchy('sigma_sd', beta=10, shape=s)
    
    # sigma 
    sigma = pm.HalfNormal('sigma', sd=sigma_sd, shape=s)
    
    mu = pm.GaussianRandomWalk('mu', 
                              sigma=sigma * (1. - alpha), 
                              shape=y.shape
                              )

    likelihood = pm.Normal('timeseries', 
                      mu=mu, 
                      sigma=sigma * alpha, 
                      observed=y
                          )

Is the specification of the hierarchical model wrong? Should I try something else here?

A lot of my timeseries are correlated, so I am inclined to think that I can fit a single multivariate model, but I am not sure how. Any advice would be highly appreciated. Thank you.

2

There is a bug right now in pm.GaussianRandomWalk, take a look at the discussion here might help you: Multiple independent random walks with differing innovations

1 Like
3

Thank you @junpenglao. The MvGRW code you posted in that thread works like a charm!

1 Like