How to structure hierarchical linear regression with multiple time series

[GallowayJ]

I ended up using the model below to solve my problems. Main things I did were:

• Restructure my dataset so that each time series in its own column
• Create a shared variable from my dataframe using pm.Data (follow this example in the docs)
• When wanting to multiple a vector by another vector with different first dimensions I have used np.asarray(1, my_vector) Based on this talk by Meenal Jhajharia, specifically the section regarding broadcasting

The final code and model graph:

coords = {
    'calendar_features': np.arange(seasonal_indicator_df.shape[1]),
    'time': new_data.index.to_numpy(),
    'time_series': (new_data.columns)
}

with pm.Model(coords=coords) as hier_linear_w_seasonality:
    data = pm.Data('time series', new_data, dims=('time', 'time_series'))
    mu_a = pm.Normal('mu_a', mu=.30, sigma=0.2)
    sigma_a = pm.HalfNormal('sigma_a', 0.2)
    mu_b = pm.Normal('mu_b', mu=.15, sigma=0.2)
    sigma_b = pm.HalfNormal('sigma_b', 0.2)
   
    a = pm.Normal('a', mu=mu_a, sigma=sigma_a, dims='time_series')
    b = pm.Normal('b', mu=mu_b, sigma=sigma_b, dims='time_series')
    
    b_calendar = pm.Normal('b_cal', mu=0, sigma=0.3, dims=('calendar_features', 'time_series'))
    seasonality = pm.Deterministic('seasonality', pm.math.dot(seasonal_indicator_df.to_numpy().T, b_calendar))
    
    trend = pm.Deterministic('trend', a + (b * np.asarray(1, new_data.index)))
    y_estimate = trend + seasonality
    error = pm.HalfNormal('error', 0.15)
    
    data_like = pm.Normal('y_expected', mu=y_estimate, sigma=error, observed=data)