# Hierarchical time series / grouped time series coherence

Hello! Thank you for maintaining Pyro.

I am working on a hierarchical time series problem where sub-time series aggregate into larger ones, as described in this book.

Thanks to the nice Pyro tutorials I have already fit a simple multivariate time series model, but I do not see how to enforce coherence. I do not understand how to incorporate the S matrix (using the linked book’s notation) into the model, since all the time series (including the ‘summed up’ one) are treated the same in this formulation. I would appreciate any advice on the matter!

Minimal example:

``````import torch
import pyro
from pyro import distributions as dist
from pyro import poutine
from pyro.contrib.forecast import ForecastingModel, Forecaster, eval_crps
from pyro.infer.reparam import LocScaleReparam, SymmetricStableReparam
from pyro.ops.tensor_utils import periodic_repeat
from pyro.ops.stats import quantile
torch.random.manual_seed(123)

t = torch.arange(50)
w1 = torch.tensor([1, 2, 3, 5, 2, 1, 1])
w2 = torch.tensor([2, 1, 2, 3, 2, 2, 2])

y2 = 25 + 0.3 * t + 0.5*torch.randn(len(t)) + w2[t % 7]
y1 = 10 + 0.15 * t + 0.5*torch.randn(len(t)) + w1[t % 7]
y = y1 + y2
Y = torch.stack((y, y1, y2), 1)

class Model(ForecastingModel):
def model(self, zero_data, covariates):
duration, data_dim = zero_data.shape

drift_stability = pyro.sample("drift_stability", dist.Uniform(1, 2))
drift_scale = pyro.sample("drift_scale", dist.HalfNormal(10))
with pyro.plate("origin", data_dim, dim=-2):
with self.time_plate:
with poutine.reparam(config={"drift": LocScaleReparam()}):
with poutine.reparam(config={"drift": SymmetricStableReparam()}):
drift = pyro.sample("drift", dist.Stable(drift_stability, 0, drift_scale))

with pyro.plate("day_of_week", 7, dim=-1):
seasonal = pyro.sample("seasonal", dist.Normal(0, 10))

seasonal = periodic_repeat(seasonal, duration, dim=-1)
motion = drift.cumsum(dim=-1)
prediction = motion + seasonal

assert prediction.shape[-2:] == (data_dim, duration)
prediction = prediction.unsqueeze(-1).transpose(-1, -3)
assert prediction.shape[-3:] == (1, duration, data_dim), prediction.shape

obs_scale = pyro.sample("obs_scale", dist.HalfNormal(1))
noise_dist = dist.Normal(0, obs_scale.unsqueeze(-1))
self.predict(noise_dist, prediction)

pyro.set_rng_seed(1)
pyro.clear_param_store()
covariates = torch.zeros(len(t), 0, dtype=torch.float)
forecaster = Forecaster(Model(), Y, covariates, learning_rate=0.1, num_steps=1000)
for name, value in forecaster.guide.median().items():
if value.numel() == 1:
print("{} = {:0.4g}".format(name, value.item()))
``````

So I thought a bit more about this and I think it’s actually very simple if we can subset on the time series that are “unaggregated”. Then we run the model as before on the subset of series, and apply the S matrix to sum them together. Perhaps this should be done outside of the model though, since the summation is completely deterministic. I have not thought enough about the pros and cons yet but wanted to share this for posterity.

``````
class CoherentModel(ForecastingModel):
def model(self, zero_data, covariates):
duration, full_data_dim = zero_data.shape
data_dim = full_data_dim - 1 # remove the aggregated dimension

drift_stability = pyro.sample("drift_stability", dist.Uniform(1, 2))
drift_scale = pyro.sample("drift_scale", dist.HalfNormal(10))
with pyro.plate("origin", data_dim, dim=-2):
with self.time_plate:
with poutine.reparam(config={"drift": LocScaleReparam()}):
with poutine.reparam(config={"drift": SymmetricStableReparam()}):
drift = pyro.sample("drift", dist.Stable(drift_stability, 0, drift_scale))

with pyro.plate("day_of_week", 7, dim=-1):
seasonal = pyro.sample("seasonal", dist.Normal(0, 10))

seasonal = periodic_repeat(seasonal, duration, dim=-1)
motion = drift.cumsum(dim=-1)
part_prediction = motion + seasonal

# construct and apply S matrix
S = torch.tensor([
[1, 1],
[1, 0],
[0, 1]
], dtype=torch.float)
prediction = S@part_prediction

assert prediction.shape[-2:] == (full_data_dim, duration)
prediction = prediction.unsqueeze(-1).transpose(-1, -3)
assert prediction.shape[-3:] == (1, duration, full_data_dim), prediction.shape

obs_scale = pyro.sample("obs_scale", dist.HalfNormal(1))
noise_dist = dist.Normal(0, obs_scale.unsqueeze(-1))
self.predict(noise_dist, prediction)

``````