PYMC3 Bayesian Prediction Cones

Question

I'm still learning PYMC3, but I cannot find anything on the following problem in the docs. Consider the Bayesian Structure Time Series (BSTS) model from this question with no seasonality. This can be modeled in PYMC3 as follows:

import pymc3, numpy, matplotlib.pyplot
# generate some test data
t = numpy.linspace(0,2*numpy.pi,100)
y_full = numpy.cos(5*t)
y_train = y_full[:90]
y_test = y_full[90:]

# specify the model
with pymc3.Model() as model:
  grw = pymc3.GaussianRandomWalk('grw',mu=0,sd=1,shape=y_train.size)
  y = pymc3.Normal('y',mu=grw,sd=1,observed=y_train)
  trace = pymc3.sample(1000)
  y_mean_pred = pymc3.sample_ppc(trace,samples=1000,model=model)['y'].mean(axis=0)

  fig = matplotlib.pyplot.figure(dpi=100)
  ax = fig.add_subplot(111)
  ax.plot(t,y_full,c='b')
  ax.plot(t[:90],y_mean_pred,c='r')
  matplotlib.pyplot.show()

Now I would like to predict the behavior for the next 10 time steps, i.e., y_test. I would also like to include credible regions over this area produce a Bayesian cone, e.g., see here. Unfortunately the mechanism for producing the cones in the aforementioned link is a little vague. In a more conventional AR model one could learn the mean regression coefficients and manually extend the mean curve. However, in this BSTS model there is no obvious way to do this. Alternatively, if there were regressors, then I could use a theano.shared and update it with a finer/extended grid to impute and extrapolate with sample_ppc, but thats not really an option in this setting. Perhaps sample_ppc is a red herring here, but its unclear how else to proceed. Any help would be welcome.

Answer 1

I think the following work. However, its super clunky and requires that I know how far in advance I want to predict before I train (in particular it percludes streaming usage or simple EDA). I suspect there is a better way and I would much rather accept a better solution by someone with more Pymc3 experience

import numpy, pymc3, matplotlib.pyplot, seaborn

# generate some data
t = numpy.linspace(0,2*numpy.pi,100)
y_full = numpy.cos(5*t)
# mask the data that I want to predict (requires knowledge 
#   that one might not always have at training time).
cutoff_idx = 80
y_obs = numpy.ma.MaskedArray(y_full,numpy.arange(t.size)>cutoff_idx)

# specify and train the model, used the masked array to supply only 
#   the observed data
with pymc3.Model() as model:
  grw = pymc3.GaussianRandomWalk('grw',mu=0,sd=1,shape=y_obs.size)
  y = pymc3.Normal('y',mu=grw,sd=1,observed=y_obs)
  trace = pymc3.sample(5000)
  y_pred = pymc3.sample_ppc(trace,samples=20000,model=model)['y']
  y_pred_mean = y_pred.mean(axis=0)

  # compute percentiles
  dfp = numpy.percentile(y_pred,[2.5,25,50,70,97.5],axis=0)

  # plot actual data and summary posterior information
  pal = seaborn.color_palette('Purples')
  fig = matplotlib.pyplot.figure(dpi=100)
  ax = fig.add_subplot(111)
  ax.plot(t,y_full,c='g',label='true value',alpha=0.5)
  ax.plot(t,y_pred_mean,c=pal[5],label='posterior mean',alpha=0.5)
  ax.plot(t,dfp[2,:],alpha=0.75,color=pal[3],label='posterior median')
  ax.fill_between(t,dfp[0,:],dfp[4,:],alpha=0.5,color=pal[1],label='CR 95%')
  ax.fill_between(t,dfp[1,:],dfp[3,:],alpha=0.4,color=pal[2],label='CR 50%')
  ax.axvline(x=t[cutoff_idx],linestyle='--',color='r',alpha=0.25)
  ax.legend()
  matplotlib.pyplot.show()

This outputs the following which seems like a really bad prediction, but at least the code is supplying out of sample values.