# 1) class option
model = LinearTrend() + FourierSeasonality()
model = HierarchicalLinearTrend(group='subject') + FourierSeasonality()
model = GroupedLinearTrend(group='subject') + FourierSeasonality()

# 2) metaclass option
model = LinearTrend() + FourierSeasonality()
model = Hierarchical(LinearTrend(), on='subject') + FourierSeasonality()
model = Distinct(LinearTrend(), on='subject') + FourierSeasonality()

# 3) seperate hierarchy object 1
partial_subject = Share(type='partially_shared', on='subject')
distinct_subject = Share(type='distinct', on='subject')

model = LinearTrend() + FourierSeasonality()
model = LinearTrend().share(partial_subject) + FourierSeasonality()
model = LinearTrend().share(distinct_subject) + FourierSeasonality()

# 4) seperate hierarchy object 2
shared_subject = Share(type='shared')  # default
partial_subject = Share(type='partially_shared', on='subject')
distinct_subject = Share(type='distinct', on='subject')

model = LinearTrend(share=shared_subject) + FourierSeasonality(share=shared_subject)
model = LinearTrend(share=partial_subject) + FourierSeasonality()
model = LinearTrend(share=distinct_subject) + FourierSeasonality()

# 5) keyword option
model = LinearTrend(shared='full') + FourierSeasonality()  #  shared='full' is default
model = LinearTrend(group='subject', shared='partial') + FourierSeasonality()
model = LinearTrend(group='subject', shared='none') + FourierSeasonality()

# Cannibalism in this option
model = Constant(range=(0, None), shared='none') + LinearTrend(shared='full') + ((FourierSeasonality(perios=365) + FourierSeasonality(period=7)) * Constant(range=[-1, 1], group='subject', shared='none'))
model = LinearTrend() + (FourierSeasonality(perios=365) * Constant(range=[-1, 1], group='subject', shared='none')) + FourierSeasonality(period=7)

# 6) keyword option
model = LinearTrend() + FourierSeasonality()
model = LinearTrend(partially_shared_on='subject') + FourierSeasonality()
model = LinearTrend(distinct_on='subject') + FourierSeasonality()

# 7) bitwise or Operator option
model = LinearTrend()                                         # single estimation for all subjects
model = (LinearTrend() | 'subject') + FourierSeasonality()    # hierarchical estimation
model = (LinearTrend() || ('subject') + FourierSeasonality()  # independent estimation

#return day well spent

I have been working on this feature, and had chance to play with it on some highly seasonal series. If you think this would be a nice feature, I can wrap it up, open a PR and hopefully share some plots.

Sometimes there is an unknown special effect that results in very high values in the time series. (let's say a hurricane). I might be okay not getting accurate predictions for those unknown special days. But the neighboring points might also be affected due to the Fourier seasonality. I think it could be nice to provide a model that is robust to outliers. One approach could be adding an option to use pm.StudentT for the observations.

PyMC3 docs has a nice example that includes other (more advanced) methods as well. GLM: Robust Regression example

Here is a motivating example:

from timeseers.utils import seasonal_data
from timeseers import FourierSeasonality, LinearTrend
import numpy as np
import pandas as pd

np.random.seed(42)

df, _ = seasonal_data(5)

df.loc[630, "value"] = 100.
df.loc[631, "value"] = 100.
df.loc[632, "value"] = 100.
df.loc[633, "value"] = 100.

model = FourierSeasonality(n=5, period=pd.Timedelta(days=365)) + LinearTrend(n_changepoints=1)
model.fit(df[["t"]], df["value"], tune=2000, cores=1, chains=1)
model.plot_components(X_true=df, y_true=df['value'])

Using normal distribution I get the following predictions.

Changing to StudentT (two line diff), I get:

Do you think this would be a useful addition to the library?

Afterthought: There is a similar discussion on fbprophet too. See: #1500.

It might be nice if the users could define the likelihood in TimeSeriesModel as they please.

Sorry for writing this question as an issue. One of the applications I deal with is forecasting short-mid life sku (2-10months). With prophet, you cannot model the trend for the product lifecycle with introduction-growth-maturity-decline phases. However, we can model trend with some kind of scikit/random forest type model based on other SKUs that have completed their lifecycle.

Is there any way to insert an external model for the trend component, for an application such as above?

Closes #39

Add RBFSeasonality class and its tests. One might directly specify the peaks by passing peaks argument. Otherwise default peaks are obtained by evenly placing 20 peaks in a single period. I have also added rbf_seasonal_data utility for testing purposes.

I have tried running it with seasonal_data and rbf_seasonal_data. Both seems fine. I have also tested it some real sales data (both with partial pooling and complete pooling). It seems to perform well especially with highly seasonal series (such as fruit sales). I hope I can share such plots in the future.

Here is a sample of peaks (features):

Code to reproduce:

import pandas as pd
import numpy as np
from matplotlib import pyplot as plt
from timeseers.utils import MinMaxScaler


def _X_t(t, periods, p, sigma, year):
    year = year / p
    mod = (t % year)[:, None]
    left_difference = np.sqrt( (mod - periods[None, :]) **2 )
    right_difference = np.abs(year - left_difference)
    return  np.exp(- ((np.minimum(left_difference, right_difference)) ** 2) / (2 * sigma**2))

t = pd.Series(pd.date_range("2010-01-01", "2014-01-01"))
scaler = MinMaxScaler()
scaled_t = scaler.fit_transform(t)
periods = [pd.Timedelta(days=i * 18.2625) for i in range(20)]
periods_ = [p / scaler.scale_factor_ for p in periods]
periods_ = np.array(periods_)
scale_factor = t.max() - t.min()
res = _X_t(scaled_t, periods_, scaler.scale_factor_, sigma=0.01, year=pd.Timedelta(days=365.25))
plt.figure(figsize=(24,10))
plt.plot(res);

Here is a fit when the data is obtained with seasonal_data utility:

Code to reproduce:

from timeseers import RBFSeasonality
from timeseers.utils import seasonal_data, MinMaxScaler
import pandas as pd
import numpy as np

df, _ = seasonal_data(5)
periods = [pd.Timedelta(days=i * 18.2625) for i in range(20)]

pers = np.array(periods)
model = RBFSeasonality(period=pd.Timedelta(days=365.25), sigma=0.01)
model.fit(df[['t']], df['value'], tune=2000, y_scaler=MinMaxScaler)
model.plot_components(X_true=df, y_true=df['value']);

Here is a sample fit with rbf_seasonal_data:

Code to reproduce:

import pandas as pd

from timeseers.utils import MinMaxScaler, get_periodic_peaks, rbf_seasonal_data
from timeseers import RBFSeasonality

df, _ = rbf_seasonal_data(5)
ps = get_periodic_peaks(10)
period = pd.Timedelta(days=365.25)
model = RBFSeasonality(peaks=ps, period=period, sigma=0.015)
model.fit(df[['t']], df['value'], tune=2000, y_scaler=MinMaxScaler)
model.plot_components(X_true=df, y_true=df['value']);

Not sure about the API, maybe you can point me towards some improvements.

Here's my take on this.

• Add base Likehood class
• Allow users to pass arbitrary likelihood terms as instances of Likelihood
• Define Gaussian likelihood and use it as the default
• Add StudentT likelihood. We might remove this one and expect the users to extend as they please. Or we can create a pool of likelihoods (i.e. negative binomial).

I'm not sure about the variable naming.

Closes #30

Closes #2 Added piecewise logistic growth and tested with some data. I have tested with the example data at prophet docs

• I had to add a MaxScaler to get the model working properly. I am not sure if it's okay. Fbprophet does something like this:

        if self.growth == 'logistic' and 'floor' in df:
            self.logistic_floor = True
            floor = df['floor']
        else:
            floor = 0.
        self.y_scale = (df['y'] - floor).abs().max()

        if 'y' in df:
            df['y_scaled'] = (df['y'] - df['floor']) / self.y_scale

• I have added a test case following other examples.
• Currently we have a constant capacity and no floor.

Hey,

I was browsing through the codebase and found it very interesting, considering playing a bit more with it.

I saw then, that in some modules, there seems to be an unnecessary dependency on the theano.

https://github.com/MBrouns/timeseers/blob/50275a4f50eae9b7af0083b8db91d57b63684692/src/timeseers/constant.py#L6

No big issue, but more accessible if this would be removed while being quite a quick fix. What do you think?

I went ahead and created a PR. Feel free to reject if you have other plans here.

• y_scale_factor changed to std_
• parameter names are updated
• I forgot to remove cores=1 and chains=1 from logistic growth test. That might be slowing down.

Closes #35

The weighted regression is implemented via pm.Potential and the logp method for each pymc3 distribution. To implement, added define_priors abstract method to generalize the previous abstract observed method.

The sample_weight argument is added to the fit method of TimeSeriesModel in the same format as sklearn arguments.

Currently working through tests.

Hi, Matthijs!

You’ve done a good job with your timeseers package! Thanks! I’m a data scientist in a big retail company and we are now working on a new product in which we want to forecast sales. It seems like timeseers could be a great decision because we have lots of SKUs and it is really difficult to builld a model for each of them individually. In our model we would like to use some exogenous regressors. It seems like timeseers already has a “Regressor” module. However, it is not very clear if it works and how it works. Could you please give some comments about it or some toy example how to use it? Does it work only with categorical features or continuous features can also be used? What “ON” argument is used for?

Thanks in advance for your reply!

Kate

When using the LogisticGrowth model the _y_scale attribute of the base class is not found. This seem to be because the method that actually sets this attribute is never called.

For example, line 38 - 44 in linear_trend with n_groups=1 makes delta a matrix of 1 by n_changepoints. It might be faster if this is just a vector in that case?

We noticed that if we run the a model with a LinearTrend without any arguments it errors with an unclear error message: " unsupported operand type(s) for +: 'NoneType' and 'int' "