Hierarchical time series / grouped time series coherence

Misc.
1

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()))
2

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)