有两种有趣的时间序列预测方法称为BATS和TBATS ，它们能够对具有多个季节性的时间序列进行建模。

名称是模型关键特征的首字母缩略词：Trigonometric seasonality, Box-Cox transformation, ARMA errors, Trend and Seasonal components。

TBATS模型以指数平滑方法为基础，可以通过以下公式描述：

每个季节性都是基于傅里叶级数的三角表示。这种方法的一个主要优点是无论周期长短，它只需要2个seed状态。另一个优点是能够模拟非整数长度的季节性影响。例如，给定日常观察，可以模拟闰年，长度为365.25。

BATS与TBATS的的区别仅仅在于它模拟季节效应的方式。在BATS中，我们有一种更传统的方法，每种季节性都通过以下方式建模：

这意味着BATS只能模拟整数周期长度。在BATS中采取的方法需要第i季的m_i seed状态，如果这个季节很长，模型可能会变得难以处理。

TBATS如何选择最终模型

TBATS将考虑各种替代方案，它会考虑模型：

• 有没有Box-Cox转换
• 有没有趋势
• 有和没有趋势阻尼
• 使用和不使用ARMA（p，q）过程来模拟残差
• 非季节性模型
• 用于模拟季节性影响的各种谐波

最终的模型将使用Akaike information criterion (AIC)进行选择。

特别是auto ARMA用于决定残差是否需要建模，以及哪些p和q值是合适的。

有关详细信息，我们邀请您阅读原始论文https://www.tandfonline.com/doi/abs/10.1198/jasa.2011.tm09771。

现有实现

到目前为止，唯一的实现已经在 forecast包中以R语言提供。使用Python的数据科学家要么不得不放弃测试这些模型，或者被迫使用Rpy2等R wrapper 来运行它们。

新的实现

Python中的TBATS实现，可以在GitHub上找到（https://github.com/intive-DataScience/tbats）。在本文的其余部分，我们将提供示例用法，并将此实现的性能与其他方法进行比较。

时间序列数据集

为了测试预测方法，我们需要一些时间序列数据。让我们使用Kaggle Store Item Demand Forecasting Challenge时间序列（https://www.kaggle.com/c/demand-forecasting-kernels-only）。这里的数据包含5年期间10个商店中每天50个商品的销售额（总共500个不同的时间序列）。出于我们的目的，我们只需要一个时间序列，因此我将在商店1任意取第1项的销售额。

import pandas as pd
df = pd.read_csv('kaggle_sales.csv')
df = df[(df['store'] == 1) & (df['item'] == 1)] # item 1 in store 1
df = df.set_index('date')
y = df['sales']
y_to_train = y.iloc[:(len(y)-365)]
y_to_test = y.iloc[(len(y)-365):] # last year for testing

图1：商店1的第1项的每日销售额

销售数据包含每日观察。它展示了每周和每年的季节性模式。这意味着我们正在处理包含多个季节性影响的时间序列。其中一个季节性很长，包含365个（闰年366个）观测值。这是TBATS的设计目标。

TBATS模型

为了开始预测，我们需要安装tbats包并拟合模型。我们必须手动提供的模型是季节长度，Python代码如下：

from tbats import TBATS, BATS
# Fit the model
estimator = TBATS(seasonal_periods=(7, 365.25))
model = estimator.fit(y_to_train)
# Forecast 365 days ahead
y_forecast = model.forecast(steps=365)

您可能已经注意到，每年的季节长度不是整数。它等于365.25以说明闰年，这是TBATS能够处理的特征。

TBATS似乎在模拟两种季节性效应方面做得非常出色：

图2：3年的销售数据。TBATS正在模拟年度季节性影响

图3：12周的数据。TBATS也在模拟每周季节性影响

如果我们仔细研究并查看模型参数，我们将发现3个季节性谐波用于模拟每周模式，11个谐波用于模拟年度模式。TBATS选择使用了λ为0.234955的box - cox变换。没有对趋势进行建模，也没有使用ARMA对残差进行建模，因为p, q为0。

Use Box-Cox: True
Use trend: False
Use damped trend: False
Seasonal periods: [ 7. 365.25]
Seasonal harmonics [ 3 11]
ARMA errors (p, q): (0, 0)
Box-Cox Lambda 0.234955
Smoothing (Alpha): 0.015789

具有每周季节性的SARIMA模型

让我们将TBATS与广泛使用和广为人知的另一种方法进行比较：SARIMA。事实证明，SARIMA为时间序列预测提供了最先进的解决方案。不幸的是，它有两个主要缺点：（1）只能模拟一个季节性效果，（2）季节长度不应太长。

让我们使用pmdarima包中的auto_arima构建SARIMA模型。我们将忽略每年的季节性，并专注于模拟每周季节性模式，Python代码如下：

from pmdarima import auto_arima
arima_model = auto_arima(y_to_train, seasonal=True, m=7)
y_arima_forecast = arima_model.predict(n_periods=365)

Auto arima选择了SARIMA（0,1,1）x（1,0,1,7）模型。正如预期的那样，年度模式没有建模（见图4）。

图4：SARIMA模型仅为每周模式。与图2相比

SARIMAX与傅里叶项

我们可以应用技巧 [https://content.pivotal.io/blog/forecasting-time-series-data-with-multiple-seasonal-periods]来利用SARIMAX中的外生变量来模拟傅立叶项的其他季节性。

我们将继续使用SARIMA的季节性部分对每周模式进行建模。对于每年的季节模式，我们将使用上述技巧。我比较了傅里叶级数的多项选择，其中2项的预测最准确。因此我们将使用2个傅里叶项作为外生变量。带有傅里叶项的SARIMAX模拟了每周和每年的模式的Python代码如下：

# prepare Fourier terms
exog = pd.DataFrame({'date': y.index})
exog = exog.set_index(pd.PeriodIndex(exog['date'], freq='D'))
exog['sin365'] = np.sin(2 * np.pi * exog.index.dayofyear / 365.25)
exog['cos365'] = np.cos(2 * np.pi * exog.index.dayofyear / 365.25)
exog['sin365_2'] = np.sin(4 * np.pi * exog.index.dayofyear / 365.25)
exog['cos365_2'] = np.cos(4 * np.pi * exog.index.dayofyear / 365.25)
exog = exog.drop(columns=['date'])
exog_to_train = exog.iloc[:(len(y)-365)]
exog_to_test = exog.iloc[(len(y)-365):]
# Fit model
arima_exog_model = auto_arima(y=y_to_train, exogenous=exog_to_train, seasonal=True, m=7)
# Forecast
y_arima_exog_forecast = arima_exog_model.predict(n_periods=365, exogenous=exog_to_test)

模型比较

让我们使用365天的预测来比较模型的性能。我们将使用平均绝对误差作为我们的指标：

• TBATS：3.8577
• SARIMA：7.2249
• SARIMAX有2个傅立叶项：3.9045

正如预期的那样，SARIMA提供了一个糟糕的模型，因为它无法模拟年度季节性。具有傅立叶项的TBATS和SARIMAX提供了更好的模型。

缺点

不幸的是，BATS和TBATS并非免费提供。在底层，它构建并评估了许多候选模型。这导致计算的缓慢。当需要训练大量并行时间序列的模型时，这可能是至关重要的。

与SARIMAX不同，BATS和TBATS不允许将外生变量添加到模型中以改进预测。