Forecasting with Structural Time Series

Start with the linked issue comment and the tfp.sts.forecast call, then compare how build_model_forecast supplies design_matrix values alongside parameter_samples. Done means the forecasting workflow for future regressors is clearly established and the question about rebuilding the model is resolved.

Referencing: https://github.com/tensorflow/probability/issues/343#issuecomment-480952321

Trying to clarify my understanding of best practices for forecasting with STS. So I'm good up until:

q_samples_demand_ = sess.run({k: q.sample(50)
                                for k, q in variational_posteriors.items()})

From this I can get an estimate of the inferred parameters:

Inferred parameters:
observation_noise_scale: 5.289973735809326 +- 0.02318570390343666
hour_of_day_effect/_drift_scale: 5.65079402923584 +- 0.018391167744994164
day_of_week_effect/_drift_scale: 2.2549798488616943 +- 0.02803918346762657
month_of_year_effect/_drift_scale: 0.3897028863430023 +- 0.037868596613407135
x_effect/_weights: [ 5.2341375  -2.3535817   0.51498556] +- [0.17526183 0.05769397 0.78600216]
autoregressive/_coefficients: [0.88036] +- [0.0036582]
autoregressive/_level_scale: 9.465166091918945 +- 0.018752284348011017

Then, before actually building the sts.forecast() object, I should create a new model object with the future regressor values in the design_matrix:

def build_model_forecast(observed_time_series):
    
    hour_of_day_effect = sts.Seasonal(
        num_seasons = 24,
        num_steps_per_season = 1,
        observed_time_series = observed_time_series,
        name = 'hour_of_day_effect')
    
    day_of_week_effect = sts.Seasonal(
        num_seasons = 7,
        num_steps_per_season = 24,
        observed_time_series = observed_time_series,
        name = 'day_of_week_effect')
    
    month_of_year_effect = tfp.sts.Seasonal(
        num_seasons=12,
        num_steps_per_season= 30 * 24,
        drift_scale_prior=tfd.LogNormal(loc=-1., scale=0.1),
        initial_effect_prior=tfd.Normal(loc=0., scale=5.),
        name='month_of_year_effect'
    )
    
    x_effect = sts.LinearRegression(
        design_matrix = tf.stack([load, wind, outage], axis = -1),
        weights_prior = tfd.Normal(loc=0., scale=1.),
        name='x_effect')
    
    autoregressive = sts.Autoregressive(
        order=1,
        observed_time_series=observed_time_series,
        name='autoregressive')

    model = sts.Sum([hour_of_day_effect,
                   day_of_week_effect,
                   month_of_year_effect,
                   x_effect,
                   autoregressive],
                   observed_time_series=observed_time_series)
    
    return model

Then at forecast time I would pass the original model parameter samples, along with the model containing future regressor values:

demand_forecast_dist = tfp.sts.forecast(
    model=build_model_forecast(),
    observed_time_series=target_training_data,
    parameter_samples=q_samples_demand_,
    num_steps_forecast=num_forecast_steps)

Does that seem about right? It feels a bit awkward, but other than appending the original model design matrix directly, I haven't found any alternatives.