1# Copyright 2017 The TensorFlow Authors. All Rights Reserved.
2#
3# Licensed under the Apache License, Version 2.0 (the "License");
4# you may not use this file except in compliance with the License.
5# You may obtain a copy of the License at
6#
7#     http://www.apache.org/licenses/LICENSE-2.0
8#
9# Unless required by applicable law or agreed to in writing, software
10# distributed under the License is distributed on an "AS IS" BASIS,
11# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
12# See the License for the specific language governing permissions and
13# limitations under the License.
14# ==============================================================================
15"""Implements a time series model with seasonality, trends, and transients."""
16
17from __future__ import absolute_import
18from __future__ import division
19from __future__ import print_function
20
21from tensorflow.contrib.timeseries.python.timeseries.state_space_models import level_trend
22from tensorflow.contrib.timeseries.python.timeseries.state_space_models import periodic
23from tensorflow.contrib.timeseries.python.timeseries.state_space_models import state_space_model
24from tensorflow.contrib.timeseries.python.timeseries.state_space_models import varma
25
26from tensorflow.python.ops import variable_scope
27from tensorflow.python.util import nest
28
29
30def _replicate_level_trend_models(multivariate_configuration,
31                                  univariate_configuration):
32  """Helper function to construct a multivariate level/trend component."""
33  with variable_scope.variable_scope("adder"):
34    # Construct a level and trend model for each feature, with correlated
35    # transition noise.
36    adder_features = []
37    for feature in range(multivariate_configuration.num_features):
38      with variable_scope.variable_scope("feature{}".format(feature)):
39        adder_features.append(level_trend.AdderStateSpaceModel(
40            configuration=univariate_configuration))
41    adder_part = state_space_model.StateSpaceCorrelatedFeaturesEnsemble(
42        ensemble_members=adder_features,
43        configuration=multivariate_configuration)
44  return adder_part
45
46
47class StructuralEnsemble(state_space_model.StateSpaceIndependentEnsemble):
48  r"""A structural state space time series model.
49
50  In the spirit of:
51
52  Scott, Steven L., and Hal R. Varian. "Predicting the present with bayesian
53    structural time series." International Journal of Mathematical Modelling and
54    Numerical Optimisation 5.1-2 (2014): 4-23.
55
56  Without the spike-and-slab prior, and with point estimates of parameters
57  instead of sampling.
58
59  The model includes level, trend, seasonality, and a transient moving average.
60
61  An observation at time t is drawn according to:
62    observation_t = level_t + seasonality_t + moving_average_t
63        + observation_noise_t
64    level_t = level_{t-1} + trend_{t-1} + level_noise_t
65    trend_t = trend_{t-1} + trend_noise_t
66    seasonality_t = -\sum_{n=1}^{num_seasons-1} seasonality_{t-n} +
67        seasonality_noise_t
68    moving_average_t = transient_t
69        + \sum_{j=1}^{moving_average_order} ma_coefs_j * transient_{t - j}
70
71  `observation_noise`, `level_noise`, `trend noise`, `seasonality_noise`, and
72  `transient` are (typically scalar) Gaussian random variables whose variance is
73  learned from data, and that variance is not time dependent in this
74  implementation. Level noise is optional due to its similarity with observation
75  noise in some cases. Seasonality is enforced by constraining a full cycle of
76  seasonal variables to have zero expectation, allowing seasonality to adapt
77  over time. The moving average coefficients `ma_coefs` are learned.
78
79  When presented with a multivariate series (more than one "feature", here
80  referring to endogenous features of the series), the model is replicated
81  across these features (one copy per feature of each periodic component, and
82  one level/trend model per feature), and correlations in transition noise are
83  learned between these replicated components (see
84  StateSpaceCorrelatedFeaturesEnsemble). This is in addition to the learned
85  correlations in observation noise between features. While this is often the
86  most expressive thing to do with multiple features, it does mean that the
87  model grows quite quickly, creating and computing with square matrices with
88  each dimension equal to num_features * (sum(periodicities) +
89  moving_average_order + 3), meaning that some operations are approximately
90  cubic in this value.
91  """
92  # TODO(allenl): Implement partial model replication/sharing for multivariate
93  # series (to save time/memory when the series presented can be modeled as a
94  # smaller number of underlying series). Likely just a modification of the
95  # observation model so that each feature of the series is a learned linear
96  # combination of the replicated models.
97
98  def __init__(self,
99               periodicities,
100               moving_average_order,
101               autoregressive_order,
102               use_level_noise=True,
103               configuration=state_space_model.StateSpaceModelConfiguration()):
104    """Initialize the Basic Structural Time Series model.
105
106    Args:
107      periodicities: Number of time steps for cyclic behavior. May be a list, in
108          which case one periodic component is created for each element.
109      moving_average_order: The number of moving average coefficients to use,
110          which also defines the number of steps after which transient
111          deviations revert to the mean defined by periodic and level/trend
112          components.
113      autoregressive_order: The number of steps back for autoregression.
114      use_level_noise: Whether to model the time series as having level
115          noise. See level_noise in the model description above.
116      configuration: A StateSpaceModelConfiguration object.
117    """
118    component_model_configuration = configuration._replace(
119        use_observation_noise=False)
120    univariate_component_model_configuration = (
121        component_model_configuration._replace(
122            num_features=1))
123
124    adder_part = _replicate_level_trend_models(
125        multivariate_configuration=component_model_configuration,
126        univariate_configuration=univariate_component_model_configuration)
127    with variable_scope.variable_scope("varma"):
128      varma_part = varma.VARMA(
129          autoregressive_order=autoregressive_order,
130          moving_average_order=moving_average_order,
131          configuration=component_model_configuration)
132
133    cycle_parts = []
134    periodicity_list = nest.flatten(periodicities)
135    for cycle_number, cycle_periodicity in enumerate(periodicity_list):
136      # For each specified periodicity, construct models for each feature with
137      # correlated noise.
138      with variable_scope.variable_scope("cycle{}".format(cycle_number)):
139        cycle_features = []
140        for feature in range(configuration.num_features):
141          with variable_scope.variable_scope("feature{}".format(feature)):
142            cycle_features.append(periodic.CycleStateSpaceModel(
143                periodicity=cycle_periodicity,
144                configuration=univariate_component_model_configuration))
145        cycle_parts.append(
146            state_space_model.StateSpaceCorrelatedFeaturesEnsemble(
147                ensemble_members=cycle_features,
148                configuration=component_model_configuration))
149
150    super(StructuralEnsemble, self).__init__(
151        ensemble_members=[adder_part, varma_part] + cycle_parts,
152        configuration=configuration)
153
154
155# TODO(allenl): Implement a multi-resolution moving average component to
156# decouple model size from the length of transient deviations.
157class MultiResolutionStructuralEnsemble(
158    state_space_model.StateSpaceIndependentEnsemble):
159  """A structural ensemble modeling arbitrary periods with a fixed model size.
160
161  See periodic.ResolutionCycleModel, which allows a fixed number of latent
162  values to cycle at multiple/variable resolutions, for more details on the
163  difference between MultiResolutionStructuralEnsemble and
164  StructuralEnsemble. With `cycle_num_latent_values` (controlling model size)
165  equal to `periodicities` (controlling the time over which these values
166  complete a full cycle), the models are
167  equivalent. MultiResolutionStructuralEnsemble allows `periodicities` to vary
168  while the model size remains fixed. Note that high `periodicities` without a
169  correspondingly high `cycle_num_latent_values` means that the modeled series
170  must have a relatively smooth periodic component.
171
172  Multiple features are handled the same way as in StructuralEnsemble (one
173  replication per feature, with correlations learned between the replicated
174  models). This strategy produces a very flexible model, but means that series
175  with many features may be slow to train.
176
177  Model size (the state dimension) is:
178    num_features * (sum(cycle_num_latent_values)
179      + max(moving_average_order + 1, autoregressive_order) + 2)
180  """
181
182  def __init__(self,
183               cycle_num_latent_values,
184               moving_average_order,
185               autoregressive_order,
186               periodicities,
187               use_level_noise=True,
188               configuration=state_space_model.StateSpaceModelConfiguration()):
189    """Initialize the multi-resolution structural ensemble.
190
191    Args:
192      cycle_num_latent_values: Controls the model size and the number of latent
193          values cycled between (but not the periods over which they cycle).
194          Reducing this parameter can save significant amounts of memory, but
195          the tradeoff is with resolution: cycling between a smaller number of
196          latent values means that only smoother functions can be modeled. For
197          multivariate series, may either be a scalar integer (in which case it
198          is applied to all periodic components) or a list with length matching
199          `periodicities`.
200      moving_average_order: The number of moving average coefficients to use,
201          which also defines the number of steps after which transient
202          deviations revert to the mean defined by periodic and level/trend
203          components. Adds to model size.
204      autoregressive_order: The number of steps back for
205          autoregression. Learning autoregressive coefficients typically
206          requires more steps and a smaller step size than other components.
207      periodicities: Same meaning as for StructuralEnsemble: number of steps for
208          cyclic behavior. Floating point and Tensor values are supported. May
209          be a list of values, in which case one component is created for each
210          periodicity. If `periodicities` is a list while
211          `cycle_num_latent_values` is a scalar, its value is broadcast to each
212          periodic component. Otherwise they should be lists of the same length,
213          in which case they are paired.
214      use_level_noise: See StructuralEnsemble.
215      configuration: A StateSpaceModelConfiguration object.
216    Raises:
217      ValueError: If `cycle_num_latent_values` is neither a scalar nor agrees in
218          size with `periodicities`.
219    """
220    component_model_configuration = configuration._replace(
221        use_observation_noise=False)
222    univariate_component_model_configuration = (
223        component_model_configuration._replace(
224            num_features=1))
225
226    adder_part = _replicate_level_trend_models(
227        multivariate_configuration=component_model_configuration,
228        univariate_configuration=univariate_component_model_configuration)
229    with variable_scope.variable_scope("varma"):
230      varma_part = varma.VARMA(
231          autoregressive_order=autoregressive_order,
232          moving_average_order=moving_average_order,
233          configuration=component_model_configuration)
234
235    cycle_parts = []
236    if periodicities is None:
237      periodicities = []
238    periodicity_list = nest.flatten(periodicities)
239    latent_values_list = nest.flatten(cycle_num_latent_values)
240    if len(periodicity_list) != len(latent_values_list):
241      if len(latent_values_list) != 1:
242        raise ValueError(
243            ("`cycle_num_latent_values` must either be a list with the same "
244             "size as `periodicity` or a scalar. Received length {} "
245             "`cycle_num_latent_values`, while `periodicities` has length {}.")
246            .format(len(latent_values_list), len(periodicity_list)))
247      latent_values_list *= len(periodicity_list)
248    for cycle_number, (cycle_periodicity, num_latent_values) in enumerate(
249        zip(periodicity_list, latent_values_list)):
250      with variable_scope.variable_scope("cycle{}".format(cycle_number)):
251        cycle_features = []
252        for feature in range(configuration.num_features):
253          with variable_scope.variable_scope("feature{}".format(feature)):
254            cycle_features.append(
255                periodic.ResolutionCycleModel(
256                    num_latent_values=num_latent_values,
257                    periodicity=cycle_periodicity,
258                    configuration=univariate_component_model_configuration))
259        cycle_parts.append(
260            state_space_model.StateSpaceCorrelatedFeaturesEnsemble(
261                ensemble_members=cycle_features,
262                configuration=component_model_configuration))
263
264    super(MultiResolutionStructuralEnsemble, self).__init__(
265        ensemble_members=[adder_part, varma_part] + cycle_parts,
266        configuration=configuration)
267