Does the ftp.sts module support ARMA modeling?

[boxu08]

tfp.sts.Autoregressive provides AR modeling. Does tfp.sts have a class or function that provides MA or ARMA modeling?

[brianwa84]

AIUI ARMA would be equivalent to tfp.sts.Sum([tfp.sts.Autogregressive(order=1), tfp.sts.LocalLevel()])

[brianwa84]

@davmre fyi

[davmre]

I believe the model Brian posted is not exactly equivalent to ARMA, although they're quite similar. In ARMA, the autoregressive term is in terms of the past k observed values (which, in particular, include the contribution of the MA component), vs a Sum([Autoregressive, LocalLevel]) model where the autoregression is an independent latent process.

There are state-space representations of ARMA processes, e.g., as discussed in these lecture notes http://www-stat.wharton.upenn.edu/~stine/stat910/lectures/14_state_space.pdf and it should be straightforward to implement one of them as an STS component (I'd guess a minor change to the plain AR model). We'd love to see that as a contribution!

Dave

On Tue, Dec 10, 2019 at 10:53 AM Brian Patton notifications@github.com wrote:

Assigned #677 https://github.com/tensorflow/probability/issues/677 to @davmre https://github.com/davmre.

— You are receiving this because you were assigned. Reply to this email directly, view it on GitHub https://github.com/tensorflow/probability/issues/677?email_source=notifications&email_token=AAHSFCQGRC257IA7GX36KALQX7QSTA5CNFSM4JYSWYS2YY3PNVWWK3TUL52HS4DFWZEXG43VMVCXMZLOORHG65DJMZUWGYLUNFXW5KTDN5WW2ZLOORPWSZGOVMVFJDI#event-2871678093, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAHSFCRTZAQESA75U5TRUALQX7QSTANCNFSM4JYSWYSQ .

[jakee417]

@davmre @jburnim Following these notes on Hamilton ARMA SSMs (which it seems autoregressive.py was already following), as @davmre suggests, it is easy to modify the H' matrix to arrive at a ARMA(p, p-1) process. I would be happy to contribute this, if you think this is of value?

[davmre]

Thanks for the pointer, Jake! I'd known that a construction existed but hadn't realized how simple the change actually is. :-)

We've also been looking at some applications recently that touch on ARMA/ARIMA models, so it seems like the stars are aligning to add a proper STS component. I've put together an internal prototype that implements ARMA with the representation from the Hamilton chapter (as in your code), plus an additional wrapper to support integrated processes (the 'I' in ARIMA), so it seems quite feasible to get something done soon.

If you're interested in creating a PR and adding a lightweight test, I think we'd be able to check your code in essentially as-is, so you'd get the credit for contributing to TFP (and then I'd follow on in a separate change to add ARIMA and expose a StructuralTimeSeries component). Is that appealing? If you'd prefer not to do that work, we'll in any case likely add a similar object as part of a larger ARIMA change, so it's up to you whether you'd like to go through the contribution process. :-)

[jakee417]

@davmre, okay, I'd love to contribute (seeing that we use the library quite extensively and really love it!). Here is a pull request.. I also added a quick check against statsmodel's version, not sure if that is the kind of thing that can also be added to the unit tests. Still working on the explicit log_prob calculation, but can update later this week! Also, looking forward to the integrated differences wrapper, we had wondered about this on our team as well (seemed like it is quite natural in the Harvey formulation mentioned in chapter 3.4 of this)

[jakee417]

By the way, now that we have a super cool ARMA(p, p - 1) sts model, some users (like me) might be interested in prediction intervals from the marginal distribution of y (the observations). I believe these can be evaluated analytically (without sampling) since this is a Multivariate Normal Distribution. Any talk about adding something like this R function to tfd.MultivariateNormalTriL which might be a nice way to handle this? Or maybe you know of other hacks that I could get pretty error bars on the ARMA distribution :)?

[davmre]

Thanks for the PR! I just left a few comments, but overall it looks great. I'll look forward to the updated test if/when you get a chance. (I think we could include the statsmodels test, if you like, though it's also okay as-is so I wouldn't spend too much time on it).

A couple of thoughts re prediction intervals:

1. It looks like a prereq for qvmnorm is pvmnorm, i.e., CDFs for multivariate normal distributions. There's been some previous discussion here: https://github.com/tensorflow/probability/issues/422 but no current TFP implementation as far as I know. We'd certainly be open to a contribution.

2. On the other hand, since STS works with univariate time series I don't think multivariate cdfs are strictly necessary for what you're asking---at each time step the marginal predictive distribution is a mixture of univariate Gaussians (or even just a univariate Gaussian if you don't care about integrating over parameters), so you just need the CDF of a MixtureSameFamily of Normal distributions, which is implemented.

3. We do have univariate root search methods that can invert a CDF into a quantile function. For example, here's code using root search to get intervals from a tfp.sts.one_step_predictive distribution: https://github.com/tensorflow/probability/blob/main/tensorflow_probability/python/sts/anomaly_detection/anomaly_detection_lib.py#L304 . This should also work for STS forecast distributions, but you'd need to convert them to a similar mixture-of-univariate-normals form, something like

forecast_dist = tfp.sts.forecast(...)  # type: MixtureSameFamily of LinearGaussianStateSpaceModels
forecast_ssms = forecast_dist.components_distribution  # Extract the batch of state space models.
forecast_marginals = tfd.MixtureSameFamily(
  mixture_distribution=forecast_dist.mixture_distribution,
  components_distribution=tfd.Normal(loc=forecast_ssms.mean(), scale=forecast_ssms.stddev()))

so that forecast_marginals.log_cdf(y) returns a batch of CDFs (one for each time step) that can be plugged into a root search. An STS utility that does something like this would also be a great potential contribution. :-)

One other note if you're trying this is that both root search algorithms and state space models are dramatically faster inside of a tf.function context (especially with jit_compile=True)---the eager overhead from rerunning Python code at every loop iteration is pretty immense, so it's worth taking advantage of tf.function tracing if you're not already.

[jakee417]

Okay, yes this is super helpful! We are more in the consideration phase of importing the sts module outright (we have heavy use of just the vanilla distributions module), but if we go that route 2) and then 3) certainly seem like ideal options.

Regardless, I think we will use the ARMA model independently of sts so I would be happy to make more contributions since we try to avoid sampling at scale in production.

Also - thanks for the note about tf.function, I will keep that in mind as we convert our development code into production :)