有什么好的模型可以做高精度的时间序列预测呢？

近期整理了一下 Facebook 的 Prophet，个人感觉这是一个非常不错的时间序列预测工具。

张戎：Facebook 时间序列预测算法 Prophet 的研究

Facebook 时间序列预测算法 Prophet 的研究

Prophet 简介

Facebook 去年开源了一个时间序列预测的算法，叫做 fbprophet，它的官方网址与基本介绍来自于以下几个网站：

1. Github：https://github.com/facebook/prophet
2. 官方网址：https://facebook.github.io/prophet/
3. 论文名字与网址：Forecasting at scale，https://peerj.com/preprints/3190/

从官网的介绍来看，Facebook 所提供的 prophet 算法不仅可以处理时间序列存在一些异常值的情况，也可以处理部分缺失值的情形，还能够几乎全自动地预测时间序列未来的走势。从论文上的描述来看，这个 prophet 算法是基于时间序列分解和机器学习的拟合来做的，其中在拟合模型的时候使用了 pyStan 这个开源工具，因此能够在较快的时间内得到需要预测的结果。除此之外，为了方便统计学家，机器学习从业者等人群的使用，prophet 同时提供了 R 语言和 Python 语言的接口。从整体的介绍来看，如果是一般的商业分析或者数据分析的需求，都可以尝试使用这个开源算法来预测未来时间序列的走势。

Prophet 的算法原理

Prophet 数据的输入和输出

首先让我们来看一个常见的时间序列场景，黑色表示原始的时间序列离散点，深蓝色的线表示使用时间序列来拟合所得到的取值，而浅蓝色的线表示时间序列的一个置信区间，也就是所谓的合理的上界和下界。prophet 所做的事情就是：

1. 输入已知的时间序列的时间戳和相应的值；
2. 输入需要预测的时间序列的长度；
3. 输出未来的时间序列走势。
4. 输出结果可以提供必要的统计指标，包括拟合曲线，上界和下界等。

就一般情况而言，时间序列的离线存储格式为时间戳和值这种格式，更多的话可以提供时间序列的 ID，标签等内容。因此，离线存储的时间序列通常都是以下的形式。其中 date 指的是具体的时间戳，category 指的是某条特定的时间序列 id，value 指的是在 date 下这个 category 时间序列的取值，label 指的是人工标记的标签（'0' 表示异常，'1‘ 表示正常，'unknown' 表示没有标记或者人工判断不清）。

而 fbprophet 所需要的时间序列也是这种格式的，根据官网的描述，只要用 csv 文件存储两列即可，第一列的名字是 'ds', 第二列的名称是 'y'。第一列表示时间序列的时间戳，第二列表示时间序列的取值。通过 prophet 的计算，可以计算出 yhat，yhat_lower，yhat_upper，分别表示时间序列的预测值，预测值的下界，预测值的上界。两份表格如下面的两幅图表示。

Prophet 的算法实现

在时间序列分析领域，有一种常见的分析方法叫做时间序列的分解（Decomposition of Time Series），它把时间序列 y_{t} 分成几个部分，分别是季节项 S_{t} ，趋势项 T_{t} ，剩余项 R_{t} 。也就是说对所有的 t\geq 0 ，都有

y_{t} = S_{t} + T_{t} + R_{t}.

除了加法的形式，还有乘法的形式，也就是：

y_{t} = S_{t} \times T_{t} \times R_{t}.

以上式子等价于 \ln y_{t} = \ln S_{t} + \ln T_{t} + \ln R_{t} 。所以，有的时候在预测模型的时候，会先取对数，然后再进行时间序列的分解，就能得到乘法的形式。在 fbprophet 算法中，作者们基于这种方法进行了必要的改进和优化。

一般来说，在实际生活和生产环节中，除了季节项，趋势项，剩余项之外，通常还有节假日的效应。所以，在 prophet 算法里面，作者同时考虑了以上四项，也就是：

y(t) = g(t) + s(t) + h(t) + \epsilon_{t}.

其中 g(t) 表示趋势项，它表示时间序列在非周期上面的变化趋势； s(t) 表示周期项，或者称为季节项，一般来说是以周或者年为单位； h(t) 表示节假日项，表示在当天是否存在节假日；\epsilon_{t} 表示误差项或者称为剩余项。Prophet 算法就是通过拟合这几项，然后最后把它们累加起来就得到了时间序列的预测值。

趋势项模型 g(t)

在 Prophet 算法里面，趋势项有两个重要的函数，一个是基于逻辑回归函数（logistic function）的，另一个是基于分段线性函数（piecewise linear function）的。

首先，我们来介绍一下基于逻辑回归的趋势项是怎么做的。

如果回顾逻辑回归函数的话，一般都会想起这样的形式： \sigma(x) = 1/(1+e^{-x}), 它的导数是\sigma'(x) = \sigma(x) \cdot(1-\sigma(x)), 并且 \lim_{x\rightarrow +\infty} \sigma(x) = 1,\lim_{x\rightarrow -\infty} \sigma(x) = 0. 如果增加一些参数的话，那么逻辑回归就可以改写成：

f(x) = C / (1 + e^{-k(x-m)}),

这里的 C 称为曲线的最大渐近值， k 表示曲线的增长率， m 表示曲线的中点。当 C=1, k = 1, m =0 时，恰好就是大家常见的 sigmoid 函数的形式。从 sigmoid 的函数表达式来看，它满足以下的微分方程： y'=y(1-y) 。

那么，如果使用分离变量法来求解微分方程 y'=y(1-y) 就可以得到：

\frac{y'}{y} + \frac{y'}{1-y} = 1 \Rightarrow \ln\frac{y}{1-y} = 1 \Rightarrow y = 1/(1+K e^{-x}).

但是在现实环境中，函数 f(x) = C / (1+e^{-k(x-m)}) 的三个参数 C, k, m 不可能都是常数，而很有可能是随着时间的迁移而变化的，因此，在 Prophet 里面，作者考虑把这三个参数全部换成了随着时间而变化的函数，也就是 C = C(t), k = k(t), m = m(t).

除此之外，在现实的时间序列中，曲线的走势肯定不会一直保持不变，在某些特定的时候或者有着某种潜在的周期曲线会发生变化，这种时候，就有学者会去研究变点检测，也就是所谓 change point detection。例如下面的这幅图的 t_{1}^{*}, t_{2}^{*} 就是时间序列的两个变点。

在 Prophet 里面，是需要设置变点的位置的，而每一段的趋势和走势也是会根据变点的情况而改变的。在程序里面有两种方法，一种是通过人工指定的方式指定变点的位置；另外一种是通过算法来自动选择。在默认的函数里面，Prophet 会选择 n_changepoints = 25 个变点，然后设置变点的范围是前 80%（changepoint_range），也就是在时间序列的前 80% 的区间内会设置变点。通过 forecaster.py 里面的 set_changepoints 函数可以知道，首先要看一些边界条件是否合理，例如时间序列的点数是否少于 n_changepoints 等内容；其次如果边界条件符合，那变点的位置就是均匀分布的，这一点可以通过 np.linspace 这个函数看出来。

下面假设已经放置了 S 个变点了，并且变点的位置是在时间戳 s_{j}, 1\leq j\leq S 上，那么在这些时间戳上，我们就需要给出增长率的变化，也就是在时间戳 s_{j} 上发生的 change in rate。可以假设有这样一个向量： \boldsymbol{\delta}\in\mathbb{R}^{S}, 其中 \delta_{j} 表示在时间戳 s_{j} 上的增长率的变化量。如果一开始的增长率我们使用 k 来代替的话，那么在时间戳 t 上的增长率就是 k + \sum_{j:t>s_{j}} \delta_{j} ，通过一个指示函数 \mathbf{a}(t)\in \{0,1\}^{S} 就是

a_{j}(t) = \begin{cases} 1, \text{ if } t\geq s_{j},\\ 0, \text{ otherwise.} \end{cases}

那么在时间戳 t 上面的增长率就是 k + \mathbf{a}^{T}\boldsymbol{\delta}. 一旦变化量 k 确定了，另外一个参数 m 也要随之确定。在这里需要把线段的边界处理好，因此通过数学计算可以得到：

\gamma_{j} = \bigg(s_{j} - m - \sum_{\ell <j} \gamma_{\ell} \bigg) \cdot \bigg( 1- \frac{k + \sum_{\ell < j} \delta_{\ell}}{k + \sum_{\ell\leq j}\delta_{\ell}} \bigg).

所以，分段的逻辑回归增长模型就是：

g(t) = \frac{C(t)}{1+\exp(-(k+\boldsymbol{a}(t)^{t}\boldsymbol{\delta}) \cdot (t - (m+\boldsymbol{a}(t)^{T}\boldsymbol{\gamma})},

其中，

\boldsymbol{a}(t) = (a_{1}(t),\cdots,a_{S}(t))^{T}, \boldsymbol{\delta} = (\delta_{1},\cdots,\delta_{S})^{T}, \boldsymbol{\gamma} = (\gamma_{1},\cdots,\gamma_{S})^{T}.

在逻辑回归函数里面，有一个参数是需要提前设置的，那就是 Capacity，也就是所谓的 C(t) ，在使用 Prophet 的 growth = 'logistic' 的时候，需要提前设置好 C(t) 的取值才行。

再次，我们来介绍一下基于分段线性函数的趋势项是怎么做的。众所周知，线性函数指的是 y=kx+b, 而分段线性函数指的是在每一个子区间上，函数都是线性函数，但是在整段区间上，函数并不完全是线性的。正如下图所示，分段线性函数就是一个折线的形状。

因此，基于分段线性函数的模型形如：

g(t)=(k+\boldsymbol{a}(t)\boldsymbol{\delta})\cdot t+(m+\boldsymbol{a}(t)^{T}\boldsymbol{\gamma}),

其中 k 表示增长率（growth rate）， \boldsymbol{\delta} 表示增长率的变化量， m 表示 offset parameter。而这两种方法（分段线性函数与逻辑回归函数）最大的区别就是 \boldsymbol{\gamma} 的设置不一样，在分段线性函数中， \boldsymbol{\gamma}=(\gamma_{1},\cdots,\gamma_{S})^{T},\gamma_{j}=-s_{j}\delta_{j}. 注意：这与之前逻辑回归函数中的设置是不一样的。

在 Prophet 的源代码中，forecast.py 这个函数里面包含了最关键的步骤，其中 piecewise_logistic 函数表示了前面所说的基于逻辑回归的增长函数，它的输入包含了 cap 这个指标，因此需要用户事先指定 capacity。而在 piecewise_linear 这个函数中，是不需要 capacity 这个指标的，因此 m = Prophet() 这个函数默认的使用 growth = 'linear' 这个增长函数，也可以写作 m = Prophet(growth = 'linear')；如果想用 growth = 'logistic'，就要这样写：

m = Prophet(growth='logistic') 
df['cap'] = 6 
m.fit(df) 
future = m.make_future_dataframe(periods=prediction_length, freq='min') 
future['cap'] = 6

变点的选择（Changepoint Selection）

在介绍变点之前，先要介绍一下 Laplace 分布，它的概率密度函数为：

f(x|\mu, b) = \exp\bigg(-|x-\mu|/b\bigg)/2b,

其中 \mu 表示位置参数， b>0 表示尺度参数。Laplace 分布与正态分布有一定的差异。

在 Prophet 算法中，是需要给出变点的位置，个数，以及增长的变化率的。因此，有三个比较重要的指标，那就是

1. changepoint_range，
2. n_changepoint，
3. changepoint_prior_scale。

changepoint_range 指的是百分比，需要在前 changepoint_range 那么长的时间序列中设置变点，在默认的函数中是 changepoint_range = 0.8。n_changepoint 表示变点的个数，在默认的函数中是 n_changepoint = 25。changepoint_prior_scale 表示变点增长率的分布情况，在论文中， \delta_{j} \sim Laplace(0,\tau) ，这里的 \tau 就是 change_point_scale。

在整个开源框架里面，在默认的场景下，变点的选择是基于时间序列的前 80% 的历史数据，然后通过等分的方法找到 25 个变点（change points），而变点的增长率是满足 Laplace 分布\delta_{j} \sim Laplace (0,0.05) 的。因此，当 \tau 趋近于零的时候， \delta_{j} 也是趋向于零的，此时的增长函数将变成全段的逻辑回归函数或者线性函数。这一点从 g(t) 的定义可以轻易地看出。

对未来的预估（Trend Forecast Uncertainty）

从历史上长度为 T 的数据中，我们可以选择出 S 个变点，它们所对应的增长率的变化量是 \delta_{j} \sim Laplace(0,\tau) 。此时我们需要预测未来，因此也需要设置相应的变点的位置，从代码中看，在 forecaster.py 的 sample_predictive_trend 函数中，通过 Poisson 分布等概率分布方法找到新增的 changepoint_ts_new 的位置，然后与 changepoint_t 拼接在一起就得到了整段序列的 changepoint_ts。

changepoint_ts_new = 1 + np.random.rand(n_changes) * (T - 1)
changepoint_ts = np.concatenate((self.changepoints_t, changepoint_ts_new))

第一行代码的 1 保证了 changepoint_ts_new 里面的元素都大于 change_ts 里面的元素。除了变点的位置之外，也需要考虑 \delta 的情况。这里令 \lambda = \sum_{j=1}^{S}|\delta_{j}|/S ，于是新的增长率的变化量就是按照下面的规则来选择的：当 j>T 时，

\delta_{j}=\begin{cases} 0 \text{, with probability } (T-S)/T \\ \sim Laplace(0,\lambda) \text{, with probability } S/T \end{cases}.

季节性趋势

几乎所有的时间序列预测模型都会考虑这个因素，因为时间序列通常会随着天，周，月，年等季节性的变化而呈现季节性的变化，也称为周期性的变化。对于周期函数而言，大家能够马上联想到的就是正弦余弦函数。而在数学分析中，区间内的周期性函数是可以通过正弦和余弦的函数来表示的：假设 f(x) 是以 2\pi 为周期的函数，那么它的傅立叶级数就是a_{0} + \sum_{n=1}^{\infty}(a_{n}\cos(nx) + b_{n}\sin(nx)) 。

在论文中，作者使用傅立叶级数来模拟时间序列的周期性。假设 P 表示时间序列的周期， P = 365.25 表示以年为周期， P = 7 表示以周为周期。它的傅立叶级数的形式都是：

s(t) = \sum_{n=1}^{N}\bigg( a_{n}\cos\bigg(\frac{2\pi n t}{P}\bigg) + b_{n}\sin\bigg(\frac{2\pi n t}{P}\bigg)\bigg).

就作者的经验而言，对于以年为周期的序列（ P = 365.25 ）而言， N = 10 ；对于以周为周期的序列（P = 7 ）而言， N = 3 。这里的参数可以形成列向量：

\boldsymbol{\beta} = (a_{1},b_{1},\cdots,a_{N},b_{N})^{T}.

当 N = 10 时，

X(t) = \bigg[\cos\bigg(\frac{2\pi(1)t}{365.25}\bigg),\cdots,\sin\bigg(\frac{2\pi(10)t}{365.25}\bigg)\bigg]

当 N = 3 时，

X(t) = \bigg[\cos\bigg(\frac{2\pi(1)t}{7}\bigg),\cdots,\sin\bigg(\frac{2\pi(3)t}{7}\bigg)\bigg].

因此，时间序列的季节项就是： s(t) = X(t) \boldsymbol{\beta}, 而 \boldsymbol{\beta} 的初始化是 \boldsymbol{\beta} \sim Normal(0,\sigma^{2})。这里的 \sigma 是通过 seasonality_prior_scale 来控制的，也就是说 \sigma= seasonality_prior_scale。这个值越大，表示季节的效应越明显；这个值越小，表示季节的效应越不明显。同时，在代码里面，seasonality_mode 也对应着两种模式，分别是加法和乘法，默认是加法的形式。在开源代码中， X(t) 函数是通过 fourier_series 来构建的。

节假日效应（holidays and events）

在现实环境中，除了周末，同样有很多节假日，而且不同的国家有着不同的假期。在 Prophet 里面，通过维基百科里面对各个国家的节假日的描述，hdays.py 收集了各个国家的特殊节假日。除了节假日之外，用户还可以根据自身的情况来设置必要的假期，例如 The Super Bowl，双十一等。

由于每个节假日对时间序列的影响程度不一样，例如春节，国庆节则是七天的假期，对于劳动节等假期来说则假日较短。因此，不同的节假日可以看成相互独立的模型，并且可以为不同的节假日设置不同的前后窗口值，表示该节假日会影响前后一段时间的时间序列。用数学语言来说，对与第 i个节假日来说， D_{i} 表示该节假日的前后一段时间。为了表示节假日效应，我们需要一个相应的指示函数（indicator function），同时需要一个参数 \kappa_{i} 来表示节假日的影响范围。假设我们有 L 个节假日，那么

h(t)=Z(t) \boldsymbol{\kappa}=\sum_{i=1}^{L} \kappa_{i}\cdot 1_{\{t\in D_{i}\}},

其中 Z(t)=(1_{\{t\in D_{1}\}},\cdots,1_{\{t\in D_{L}\}}), \text{ } \boldsymbol{\kappa}=(\kappa_{1},\cdots,\kappa_{L})^{T}.

其中 \boldsymbol{\kappa}\sim Normal(0,v^{2}) 并且该正态分布是受到 v = holidays_prior_scale 这个指标影响的。默认值是 10，当值越大时，表示节假日对模型的影响越大；当值越小时，表示节假日对模型的效果越小。用户可以根据自己的情况自行调整。

模型拟合（Model Fitting）

按照以上的解释，我们的时间序列已经可以通过增长项，季节项，节假日项来构建了，i.e.

y(t)=g(t)+s(t)+h(t)+\epsilon

下一步我们只需要拟合函数就可以了，在 Prophet 里面，作者使用了 pyStan 这个开源工具中的 L-BFGS 方法来进行函数的拟合。具体可以参考 forecast.py 里面的 stan_init 函数。

Prophet 中可以设置的参数

在 Prophet 中，用户一般可以设置以下四种参数：

1. Capacity：在增量函数是逻辑回归函数的时候，需要设置的容量值。
2. Change Points：可以通过 n_changepoints 和 changepoint_range 来进行等距的变点设置，也可以通过人工设置的方式来指定时间序列的变点。
3. 季节性和节假日：可以根据实际的业务需求来指定相应的节假日。
4. 光滑参数： \tau= changepoint_prior_scale 可以用来控制趋势的灵活度， \sigma=seasonality_prior_scale 用来控制季节项的灵活度， v= holidays prior scale 用来控制节假日的灵活度。

如果不想设置的话，使用 Prophet 默认的参数即可。

Prophet 的实际使用

Prophet 的简单使用

因为 Prophet 所需要的两列名称是 'ds' 和 'y'，其中，'ds' 表示时间戳，'y' 表示时间序列的值，因此通常来说都需要修改 pd.dataframe 的列名字。如果原来的两列名字是 'timestamp' 和 'value' 的话，只需要这样写：

df = df.rename(columns={'timestamp':'ds', 'value':'y'})

如果 'timestamp' 是使用 unixtime 来记录的，需要修改成 YYYY-MM-DD hh:mm:ss 的形式：

df['ds'] = pd.to_datetime(df['ds'],unit='s')

在一般情况下，时间序列需要进行归一化的操作，而 pd.dataframe 的归一化操作也十分简单：

df['y'] = (df['y'] - df['y'].mean()) / (df['y'].std())

然后就可以初始化模型，然后拟合模型，并且进行时间序列的预测了。

初始化模型：m = Prophet()
拟合模型：m.fit(df)
计算预测值：periods 表示需要预测的点数，freq 表示时间序列的频率。
future = m.make_future_dataframe(periods=30, freq='min')
future.tail()
forecast = m.predict(future)

而 freq 指的是 pd.dataframe 里面的一个指标，'min' 表示按分钟来收集的时间序列。具体参见文档：http://pandas.pydata.org/pandas-docs/stable/timeseries.html#offset-aliases

在进行了预测操作之后，通常都希望把时间序列的预测趋势画出来：

画出预测图：
m.plot(forecast)
画出时间序列的分量：
m.plot_components(forecast)

如果要画出更详细的指标，例如中间线，上下界，那么可以这样写：

x1 = forecast['ds']
y1 = forecast['yhat']
y2 = forecast['yhat_lower']
y3 = forecast['yhat_upper']
plt.plot(x1,y1)
plt.plot(x1,y2)
plt.plot(x1,y3)
plt.show()

其实 Prophet 预测的结果都放在了变量 forecast 里面，打印结果的话可以这样写：第一行是打印所有时间戳的预测结果，第二行是打印最后五个时间戳的预测结果。

print(forecast[['ds', 'yhat', 'yhat_lower', 'yhat_upper']])
print(forecast[['ds', 'yhat', 'yhat_lower', 'yhat_upper']].tail())

Prophet 的参数设置

Prophet 的默认参数可以在 forecaster.py 中看到：

def __init__(
    self,
    growth='linear',
    changepoints=None,
    n_changepoints=25, 
    changepoint_range=0.8,
    yearly_seasonality='auto',
    weekly_seasonality='auto',
    daily_seasonality='auto',
    holidays=None,
    seasonality_mode='additive',
    seasonality_prior_scale=10.0,
    holidays_prior_scale=10.0,
    changepoint_prior_scale=0.05,
    mcmc_samples=0,
    interval_width=0.80,
    uncertainty_samples=1000,
):

增长函数的设置

在 Prophet 里面，有两个增长函数，分别是分段线性函数（linear）和逻辑回归函数（logistic）。而 m = Prophet() 默认使用的是分段线性函数（linear），并且如果要是用逻辑回归函数的时候，需要设置 capacity 的值，i.e. df['cap'] = 100，否则会出错。

m = Prophet()
m = Prophet(growth='linear')
m = Prophet(growth='logistic')

变点的设置

在 Prophet 里面，变点默认的选择方法是前 80% 的点中等距选择 25 个点作为变点，也可以通过以下方法来自行设置变点，甚至可以人为设置某些点。

m = Prophet(n_changepoints=25)
m = Prophet(changepoint_range=0.8)
m = Prophet(changepoint_prior_scale=0.05)
m = Prophet(changepoints=['2014-01-01'])

而变点的作图可以使用：

from fbprophet.plot import add_changepoints_to_plot
fig = m.plot(forecast)
a = add_changepoints_to_plot(fig.gca(), m, forecast)

周期性的设置

通常来说，可以在 Prophet 里面设置周期性，无论是按月还是周其实都是可以设置的，例如：

m = Prophet(weekly_seasonality=False)
m.add_seasonality(name='monthly', period=30.5, fourier_order=5)
m = Prophet(weekly_seasonality=True)
m.add_seasonality(name='weekly', period=7, fourier_order=3, prior_scale=0.1)

节假日的设置

有的时候，由于双十一或者一些特殊节假日，我们可以设置某些天数是节假日，并且设置它的前后影响范围，也就是 lower_window 和 upper_window。

playoffs = pd.DataFrame({
  'holiday': 'playoff',
  'ds': pd.to_datetime(['2008-01-13', '2009-01-03', '2010-01-16',
                        '2010-01-24', '2010-02-07', '2011-01-08',
                        '2013-01-12', '2014-01-12', '2014-01-19',
                        '2014-02-02', '2015-01-11', '2016-01-17',
                        '2016-01-24', '2016-02-07']),
  'lower_window': 0,
  'upper_window': 1,
})
superbowls = pd.DataFrame({
  'holiday': 'superbowl',
  'ds': pd.to_datetime(['2010-02-07', '2014-02-02', '2016-02-07']),
  'lower_window': 0,
  'upper_window': 1,
})
holidays = pd.concat((playoffs, superbowls))

m = Prophet(holidays=holidays, holidays_prior_scale=10.0)

结束语

对于商业分析等领域的时间序列，Prophet 可以进行很好的拟合和预测，但是对于一些周期性或者趋势性不是很强的时间序列，用 Prophet 可能就不合适了。但是，Prophet 提供了一种时序预测的方法，在用户不是很懂时间序列的前提下都可以使用这个工具得到一个能接受的结果。具体是否用 Prophet 则需要根据具体的时间序列来确定。

参考文献：

1. https://otexts.org/fpp2/components.html
2. https://en.wikipedia.org/wiki/Decomposition_of_time_series
3. A review of change point detection methods, CTruong, L. Oudre, N.Vayatis
4. https://github.com/facebook/prophet
5. https://facebook.github.io/prophet/

时间序列预测是我们实际项目场景中经常碰到的一类主题。在这篇文章里简单介绍一下我们观远在时序问题上的一些探索和心得体会。

时序问题的定义和分类

顾名思义，时间序列指的是按照时间顺序先后发生的数据序列，在此基础上，我们会对这个序列做各种任务，如分类，聚类，异常检测，预测等。本文主要的关注点会放在时间序列的预测类任务上。

时序预测的用途非常的广泛，如气象降雨，交通流量，金融，商业公司的销售经营，医学上的药物反应，各类系统的运作负荷，等等方面都可以见到相关的应用场景。麦肯锡的一个关于 AI 场景价值的研究中也把时序类问题的价值排在了各类数据形式中的第二位：

下图是 Google 内部的各种时序预测的场景，看起来是不是应用范围一点也不亚于搜广推呢 :)

在很多研究中，会把时间序列进一步区分成单变量时间序列和多变量时间序列。例如单支股票的价格走势是一个单变量时序问题，但某支股票的价格波动除了自身因素，也会受到其它股票价格波动的影响，甚至是股票市场外的一些其它信息影响，所以我们需要通盘考虑来做预测，就形成了所谓的多变量时间序列问题。

简单来设想，我们只要把输入和输出的维度从一维变到高维就能处理多变量的时序问题了。但实际问题中，情况往往更为复杂。例如股票市场中可能会有新股上市，也可能会有退市，这个多维本身可能也会有数量的波动。如果考虑所有曾经出现过的股票信息，那么维度可能就会很高，但序列本身长度却相对受到局限，不像 NLP 场景中有海量的序列数据可以做各种数据增强，自监督学习，迁移学习的玩法，一般很难达到比较好的预测效果。

因此，业界在实际做时序问题时，通常采用的手段还是对每一个序列来做训练学习（但不一定是单独模型），把相关的动态，静态信息通过特征的方式输入到模型中。比如股票代码就不再是维度上的区别，而是做成一个类别变量，输入到模型中进行训练。这也是目前较为主流的问题形式和对应的实践方式，后面我们会再展开详述。

时序预测的效果检验

在介绍具体预测方法之前，我们先来看下如何检验时序预测的效果。

Metric

在指标方面，作为一个回归问题，我们可以使用 MAE，MSE 等方式来计算。但这类 metric 受到具体预测数值区间范围不同，展现出来的具体误差值区间也会波动很大。比如预测销量可能是几万到百万，而预测车流量可能是几十到几百的范围，那么这两者预测问题的 MAE 可能就差距很大，我们很难做多个任务间的横向比较。

所以实际问题中，我们经常会使用对数值量纲不敏感的一些 metric，尤其是 SMAPE 和 WMAPE 这两种：

这类误差计算方法在各类不同的问题上都会落在 0~1 的区间范围内，方便我们来进行跨序列的横向比较，十分方便。

在实际项目中我们还会经常发现，很多真实世界的时序预测目标，如销量，客流等，都会形成一个类似 tweedie 或 poisson 分布的情况。如果我们用 WMAPE 作为指标，模型优化目标基本可以等价为 MAE（优化目标为中位数），则整体的预测就会比平均值小（偏保守）。在很多业务问题中，预测偏少跟预测偏多造成的影响是不同的，所以实际在做优化时，我们可能还会考察整体的预测偏差（总量偏大或偏小），进而使用一些非对称 loss 来进行具体的优化。

交叉验证

我们在一开始接触机器学习的交叉验证概念时，教科书上一般都是对数据做随机的 split。不过在时序问题上，需要特别注意不能做随机 split，而需要在时间维度上做前后的 split，以保证与实际预测应用时的情况一致。比如用 1-6 月的数据做训练，在 7 月数据上做验证，用 2-7 月的数据做训练，在 8 月数据上做验证，以此类推。

后面在介绍机器学习方法时，我们也会再次提到类似的思想，即使用一段时间窗口的历史数据作为模型输入，对未来的一个时间窗口做预测输出和效果验证。这是时序问题中极为重要的一点。

稳定性

时序问题整体来说是个难度很大的问题，纵观 Kaggle 相关比赛，能够稳定在时序问题上拿金的大佬几乎没有，但其它像 CV，分类等问题上排名靠前的 GM，熟悉的名字比例明显就会高很多。这其中很大一部分原因来自于真实世界中的时序问题基本都是不符合 i.i.d 假设的，世界变幻莫测，如果我们能准确预知未来，那还做什么算法工程师 :)

正因如此，除了验证准确率的 metric，我们还需要考察模型的稳定性。例如在不同的序列之间的精度波动，不同预测 step 之间的精度波动，以及不同时间点的精度波动等。综合表现稳定的模型往往才是业务上的更优选择。

传统时序方法

从这一节开始，我们介绍具体的预测方法。

时间序列问题有着悠久的研究历史，如果去看这方面的相关资料，会找到很多经典的时间序列检测，预处理，和预测方法。例如很多时序模型都要求序列本身是平稳的，如果没有达到检测要求的话就需要通过差分操作来转换为平稳序列。这方面的资料比较多，理论也略显复杂，我们在这里不做太多展开，只列举几种比较常用的方法。

移动平均

在业界做过几个时序预测场景的同学，经常会发出这样的感慨：“MA 好强啊”。移动平均虽然是一个很简单的方法，但往往在很多实际问题中有着非常良好的表现，是一个不容易打败的 baseline。

最简单的移动平均，就是使用过去 n 个时间点的观测值的平均，作为下一个点的预测。

可以看到移动平均给出的预测非常的“靠谱”，如果把 MA 的参考范围缩小，比如只考虑前一个时间点，我们很容易就能得出一个看起来非常准确，但只是“滞后”了一点的预测效果。但是这个滞后并不是那么好修复的，这也是很多初学者经常有的一个疑问。

在最基础的移动平均基础上，我们还可以有加权移动平均，指数移动平均等方法。当年 M4 比赛的冠军方法是由 Uber 提出的 ES-RNN，其中 ES 就是指数平滑的缩写 :)

移动平均在具体实现时也比较方便，例如在 pandas 里就有 rolling, ewm 等方法，可以直接进行计算。在 sql 中也可以用 window function 的方法来快速实现。相对其它方法来说，移动平均的计算速度很快，在海量序列场景下也能够适用。

不过如果要做多步预测，移动平均方法就会显得有些困难了。

ARIMA

时序预测领域最知名的算法，应该没有之一。其中 AR 部分可以类比为高级的加权移动平均，而 MA 虽然是移动平均的缩写，但其实是对 AR 部分的残差的移动平均。

ARIMA 相关的理论基础非常多，这里就略过了（我也不懂）。实际在使用时还需要做平稳性检验，确定 p, d, q 参数的值，使用起来有点麻烦。好在我们有 Auto ARIMA（原版为 R 中的 forecast 包，Python 中也有类似的库如 pmdarima），可以通过 AutoML 的手段来自动搜寻最佳的参数组合。

与 ARIMA 类似的还有很多，例如改进版的 SARIMA，ARIMAX，还有虽然听过，但从来没用过的 ARCH，GARCH 模型等……这类模型相比简单的移动平均，拟合能力明显会更强一些，但缺点是运行时间也明显变长了。通常来说，这类传统模型我们都需要对每一条时间序列都单独拟合和预测。如果我们要对淘宝所有的商品做销量预测，可以预见到序列的数量会非常之大，这类方法在执行时就需要花费很长的时间，而且需要用户自己来开发相应的并发执行机制。

Prophet

由 Facebook 开源的 Prophet 是另一个非常知名的时序预测模型。因为 API 设计比较友好，还附带一系列可视化和模型解释，在广大人民群众之中迅速的流行了起来。

Prophet 背后使用的是加性模型模式，将时间序列分解为趋势，季节，节假日等外部变量这三类模型之和，且利用了概率建模方式，在预测时可以输出预测值的概率分布情况。具体可以参考我们组鼎彦同学的这篇优秀的 介绍 Prophet 原理的文章。

Prophet 相比原版 ARIMA，在非线性趋势，季节性，外部变量方面都具有优势，做多步预测也会更加自然一些。但同样，Prophet 的训练预测也需要在每一条序列维度来进行，大规模序列的性能会是一个挑战。

最近 Uber 也推出了一个有些类似的时序预测库 Orbit，据称效果比 Prophet 更好。另外我们还尝试过基于 PyTorch 的 NeuralProphet，API 跟 Prophet 非常接近，不过实测下来预测的稳定性没有 Prophet 好，可能神经网络比较容易跑飞……

问题

这里我们来总结一下传统时序预测方法的一些问题：

1. 对于时序本身有一些性质上的要求，需要结合预处理来做拟合，不是端到端的优化；
2. 需要对每条序列做拟合预测，性能开销大，数据利用率和泛化能力堪忧，无法做模型复用；
3. 较难引入外部变量，例如影响销量的除了历史销量，还可能有价格，促销，业绩目标，天气等等；
4. 通常来说多步预测能力比较差。

正因为这些问题，实际项目中一般只会用传统方法来做一些 baseline，主流的应用还是属于下面要介绍的机器学习方法。

机器学习方法

如果我们去翻阅一下 Kaggle 或其它数据科学竞赛平台上的相关时序预测比赛，会发现绝大多数的获胜方案使用的是传统机器学习的方式，更具体地来说，一般就是 xgboost 和 lightgbm 这类梯度提升树模型。其中有个有意思的例外是当年的 Web Traffic Forecasting，我当时看了这个比赛也很激动，尝试了 N 多深度学习的方法来做时序问题，可惜大都没有很好的结果。砍手豪大佬的这篇文章 也对相关原因做了一些分析。下面我们对这类方法做个简单介绍。

建模方式

机器学习方法处理时序问题的基本思路跟前面提到的时序验证划分一致，就是把时序切分成一段历史训练窗口和未来的预测窗口，对于预测窗口中的每一条样本，基于训练窗口的信息来构建特征，转化为一个表格类预测问题来求解。

如上图中，浅蓝色的部分即为我们构建特征的窗口，我们利用这部分的信息输入构建特征后，再去预测深蓝色预测窗口中的值，计算误差，再不断迭代改进。这个窗口可以不断往前滑动，就形成了多个预测窗口的样本，一定程度上可以提高我们的数据利用率。

实际场景中，一般我们需要确定几个参数：

1. 历史窗口的大小，即我们预测未来时，要参考过去多少时间的信息作为输入。太少可能信息量不充分，太多则会引入早期不相关的信息（比如疫情前的信息可能目前就不太适用了）。
2. 预测点 gap 的大小，即预测未来时，我们是从 T+1 开始预测，还是 T+2，T+3？这与现实的业务场景有关，例如像补货场景，预测 T+1 的销量，可能已经来不及下单补货了，所以我们需要扩大这个提前量，做 T+3 甚至更多提前时间的预测。
3. 预测窗口的大小，即我们需要连续预测多长的未来值。比如从 T+1 开始一直到 T+14 都需要预测输出。这一点也跟实际的业务应用场景有关。

另外值得一提的是，上图中画的是一条时间序列，实际上如果我们有成百上千个序列，是可以把这些数据放在一起做训练的。这也是机器学习方法对于传统时序方法的一个较为明显的优势。

在看一些文章的时候，我们也会看到一些额外加入时序预处理步骤的方法，比如先做 STL 分解再做建模预测。我们尝试下来这类方法总体来说效果并不明显，但对于整个 pipeline 的复杂度有较大的增加，对于 AutoML，模型解释等工作都造成了一定的困扰，所以实际项目中应用的也比较少。

特征工程

这张图更明确的指出了我们构建特征和建模的方式。为了便于理解，我们可以假设预测的 horizon 长度仅为 1 天，而历史的特征 window 长度为 7 天，那么我们可以构建的最基础的特征即为过去 7 天的每天的历史值，来预测第 8 天的值。这个历史 7 天的值，跟之前提到的移动平均，AR（自回归）模型里所使用的值是一样的，在机器学习类方法中，一般被称为 lag 特征。

对于时间本身，我们也可以做各类日期衍生特征，例如我们以天为粒度做预测，我们可以添加这天是星期几，是一个月的第几天，是哪个月份，是否是工作日等等特征输入。

另外一类最常见的基础特征，就是区分不同序列的类别特征，例如不同的门店，商品，或者不同的股票代码等。通过加入这个类别特征，我们就可以把不同的时间序列数据放在一张大表中统一训练了。模型理论上来说可以自动学习到这些类别之间的相似性，提升泛化能力。

类别属性实际上可以归类为静态特征，即随着时间的变化，不会发生变化的信息。除了最细粒度的唯一键，还可以加入其它形式的静态特征。例如商品属于的大类，中类，小类，门店的地理位置特性，股票所属的行业等等。除了类别型，静态特征也可能是数值型，例如商品的重量，规格，一般是保持不变的。

Lag 特征，日期特征这类，则属于动态特征，随着时间变化会发生改变。这其中又可以分成两类，一类是在预测时无法提前获取到的信息，例如预测值本身，跟预测值相关的不可知信息，如未来的客流量，点击量等。对于这类信息，我们只能严格在历史窗口范围内做各种特征构建的处理，一般以 lag 为主。另一类则是可以提前获取到的信息，例如我们有明确的定价计划，可以预知在 T+1 时计划售卖的商品价格是多少。对于这类特征，我们则可以直接像静态特征那样直接加入对应时间点的信息进去。

以上提到的基本属于直接输入的信息，基于这些信息，我们还可以进一步做各种复杂的衍生特征。例如在 lag 的基础上，我们可以做各种窗口内的统计特征，比如过去 n 个时间点的平均值，最大值，最小值，标准差等。进一步，我们还可以跟之前的各种维度信息结合起来来计算，比如某类商品的历史均值，某类门店的历史均值等。也可以根据自己的理解，做更复杂计算的衍生，例如过去 7 天中，销量连续上涨的天数，过去 7 天中最大销量与最低销量之差等等。很多数据科学比赛的获胜方案中都会有大量篇幅来讲解这方面的衍生特征如何来构建。

最后值得一提的是还有很多将各类特征工程手段自动化的工具，在时间序列领域最有名的莫过于 tsfresh 了。除了前面提到的一些基础操作，tsfresh 还能够支持 wavelet 等高深操作，但缺点就是运行时间相对有点长，且需要结合特征选择来达到更好的效果。

模型选择

模型这块，基本上没有什么花样，大家的主流选择基本都是 GBDT 和 NN。个人最常使用的选择是 LightGBM 和 fastai，然后选择好时序验证方式，做自动参数优化就可以了（比如使用 Optuna 或 FLAML）。Lgb 的训练速度快，而且在某些业务特征比较重要的情况下，往往能达到比神经网络更好更稳定的效果。而 NN 的主要优势在类别变量的表达学习上，理论上可以达到更好的 embedding 表示。此外 NN 的 loss 设计上也会比较灵活，相对来说 lgb 的 loss 或者多目标学习限制条件就比较多了。更多的讨论也可以参考我的这篇 表格数据模型对比 的文章。总体来说，目前最常见的选择仍然是树模型一族。

有一个值得注意的考量点在于 local 模型与 global 模型的取舍。前面提到的经典时序方法中都属于 local 模型，即每一个序列都要构建一个单独的模型来训练预测；而我们提到的把所有数据都放在一起训练则是 global 模型的方式。实际场景中，可能需要预测的时序天然就会有很不一样的规律表现，比如科技类股票，跟石油能源类股票的走势，波动都非常不一样，直接放在一起训练反而可能导致整体效果下降。所以很多时候我们要综合权衡这两种方式，在适当的层级做模型的拆分训练。深度学习领域有一些工作如 DeepFactor 和 FastPoint 也在自动适配方面做了些尝试。

深度学习方法

前面有提到过在 Kaggle 2018 年的 Web Traffic Forecasting 比赛中，冠军选手采用了深度学习的方案，当时年幼的我看到这个分享大受震撼，感觉深度学习统治时间序列领域的时代就要到来了！后面也花了不少时间调研和尝试各类针对时序问题设计的 NN 模型（有不少是从 NLP 领域借鉴过来的）。不过几年尝试下来，发现绝大多数论文中实验数字看起来很漂亮的模型，在真实世界场景中应用的效果都不太理想，包括后来的很多比赛也仍然是树模型占据获胜方案的主流地位。这里有一个原因可能跟前面介绍传统时序方法中的问题类似，很多学术研究的数据集（参考 papers with code）都是比较单一的时间序列（比如气象信息，医学记录），没有包含什么其它输入信息和业务目标。而现实应用中的时序场景很多时候都是海量序列，包含了很多层级维度，促销，气候，外部事件等异常丰富的业务输入信息，其预测场景也更加丰富多样。

总体来说，深度学习的思路是尽量只使用原始的序列和其它相关输入信息，基本不做特征工程，希望通过各类模型结构自动学习到时序的隐含表达，进而做端到端的预测输出。所以我把特征工程 + NN 的方案归类到了上面机器学习方法中。这一节里简要介绍一下这几年我们做过的深度学习相关尝试，可以供同学们参考。

RNN 系列

直接借鉴 NLP 领域中经典的 RNN, GRU, LSTM 模型来做时间序列的预测应该是最直观的一种思路了。使用这个方法，甚至可以直接做任意步的预测窗口输出。但实际场景中，一般时序问题的输入输出窗口大小带有比较重要的业务含义，也需要针对性进行训练，评估，优化，所以往往不会直接使用这类原始 RNN 的方式来做训练预测。

Seq2Seq

这就是前面提到的 Web Traffic Forecasting 比赛冠军方案中主要采用的模型。基本是借鉴了 NLP 里的经典架构，使用 RNN/GRU/LSTM 作为基本单元，encoder 中做训练窗口中的信息提取，然后在 decoder 中做预测 horizon 的多步输出。作者在方案里还尝试了在 decoder 时同时引入注意力机制，但发现效果并不稳定，最后直接改成了 lag 特征来捕捉固定周期的历史信息。

在训练预测方面，作者也花了不少功夫，例如使用 SMAC3 进行自动调参，使用 COCOB 作为优化器，通过一系列 SGD averaging，多模型，多 checkpoint 输出的均值来提升模型的稳定性等，具体可以参考作者的这篇 总结文档。

我们当时也把这套架构迁移到了很多我们内部的项目中，但整体用下来发现调参的计算开销要比跑树模型大的多得多，训练稳定性却远不如树模型，很难调整到一个稳定预测输出的状态。再加上整体的误差分析和模型解释也比较难做，所以后来也并没有推广使用。砍手豪大佬之后也分析过这次比赛之所以是神经网络模型获胜，跟使用的 SMAPE 指标也有很大关系。

WaveNet

这也是当年非常火的一个模型，主要是 RNN 系列模型不好并行，所以突然发现谷歌提出的这个空洞因果卷积感觉很高级，性能理论上也比 RNN 之类的好很多，它的结构大致长这样：

除了使用一维 CNN 来做序列预测外，WaveNet 里还加入了 residual connection 和 skip connection，以及一系列复杂的“门机制”：

不过我们实际使用下来，感觉 CNN 整体对于序列问题的预测效果还是不如 RNN 系列。事后来看可能跟缺乏位置编码这类信息有关。

顺带一提 WaveNet 这类 CNN 结构也可以用在 Seq2Seq 框架中的 encoder 部分。

LSTNet

当年应该是在 Papers With Code 上看到 LSTNET 占据了几个 benchmark 的榜首位置，也去简单尝试了一下，模型结构也是愈发花里胡哨：

不过效果嘛，还是不理想，不如特征工程加 fastai 效果来的好。

DeepAR

亚马逊提出的一个网络架构，也是基于 Seq2Seq，不过 DeepAR 的输出是一个概率分布，这也是它与传统 RNN 系列最大的一个区别。

虽然来头很大，但尝试下来仍然是很难稳定收敛，多次训练的精度波动也很大，最终效果也无法与 GBDT 匹敌。不过在尝试 DeepAR 的过程中发现亚马逊开源的一个挺不错的时序预测库 gluon-ts，里面包含了非常多的深度学习时序模型实现，也很方便自己实现相关模型，大大加速了我们的实验尝试，非常值得推荐！

概率分布输出本身是个挺有用的特性，例如用于下游的 Service Level 的满足率计算等。理论上我们也可以用 quantile regression 方法，或者 ngboost ， LightGBMLSS 等库在 GBDT 模型上实现概率预测输出。

N-Beats

这也是一个来头很大的模型，出自 Element AI，Bengio 是其中的 Co-Founder。第一次见到它是来自 M5 比赛亚军的分享，不过他也只是在 top-level 的预测中使用了一下 N-Beats 模型。

从介绍来看，N-Beats 专注于做单变量的时序预测，且可以具有一定的 seasonality，trend 的可解释性，跟 Prophet 很相似。从论文的实验来看，作者使用了非常重的 ensemble，每个序列搞了 180 个模型的 bagging。这感觉有点过于“杀鸡用牛刀”了……我们实测下来也没有取得很好的效果，而且看起来还不好加额外的特征变量，使用场景很受限。

TFT

终于来到了一个尝试下来表现可以与树模型匹敌的深度学习模型了！这就是 Google AI 提出的 Temporal Fusion Transformers。也是本文提到的第一个带 transformer 结构的模型 :)

个人感觉 TFT 的设计里最有意思的是对于特征变量选择网络的考虑，从实现效果上来说跟树模型做特征选择很相似。顺带一提我们在表格类问题上测试下来表现比较好的模型，例如 TabNet，NODE 等，也都有这种模拟决策树行为的设计。具体实验下来在一些场景中 TFT 甚至可以超越特征+GBDT 的建模方案，非常惊人！不过训练计算开销仍然是比较大，这点还是跟树模型有差距。

TFT 还有一点比较有意思的是对于特征输入的设计挺系统化，分成了静态类别/连续变量，动态已知类别/连续变量，和动态未知类别/连续变量。以我们使用的 pytorch-forecasting 库为例，其 dataset 接口大致长这样：

training = TimeSeriesDataSet(
    data[lambda x: x.date <= training_cutoff],
    time_idx= ...,  # column name of time of observation
    target= ...,  # column name of target to predict
    group_ids=[ ... ],  # column name(s) for timeseries IDs
    max_encoder_length=max_encoder_length,  # how much history to use
    max_prediction_length=max_prediction_length,  # how far to predict into future
    # covariates static for a timeseries ID
    static_categoricals=[ ... ],
    static_reals=[ ... ],
    # covariates known and unknown in the future to inform prediction
    time_varying_known_categoricals=[ ... ],
    time_varying_known_reals=[ ... ],
    time_varying_unknown_categoricals=[ ... ],
    time_varying_unknown_reals=[ ... ],
)

这种归类方式非常具有通用性，值得推广。

深度学习总结

总体看来，目前（22 年初）能够大规模推广应用的深度时序模型感觉还基本没有（最新的 Informer, Autoformer 还没尝试）。前阵子有一篇很有意思的论文也讨论了这点，标题就叫 Do We Really Need Deep Learning Models for Time Series Forecasting?，从结果来看基本还是 GBDT 完胜。

这块也有很多综述文章可以参考，例如：

• Learnings from Kaggle’s Forecasting Competitions
• Neural forecasting: Introduction and literature overview

除了模型预测本身，深度学习中的各种强项，例如预训练和迁移学习，表达学习，生成式模型等方面，目前也还很难应用于时序领域。未来想要有所突破，如何能更好的做 数据增强，做 时序问题的预训练，也是很值得深入研究的方向。

最后在整体架构方面，跟机器学习模型 pipeline 来对比看，前者一般模型部分处理会相对简单，但涉及到的预处理，特征工程及后处理方面会比较复杂；而深度学习手段正好相反，在整体 pipeline 层面会相对简单，更提倡 end-to-end 的训练，但模型部分则相对复杂和难以优化。

AutoML

最后我们来看下时序领域的 AutoML，近几年 Github 上也出现了一些针对时序问题的 AutoML 库，例如：

• Auto_TS
• AutoTS

不过总体来说他们面向的还是模型的自动选择和调优，针对时序问题做了一些特定的配置项，如时间粒度，前面有提到过的预测的长度，自动的 validation 等。但从前文的对比来看，目前效果比较好的主流方法还是特征工程+GBDT 模型，尤其是特征工程这块的自动化显得尤为关键，而目前 tsfresh 也并没有在衍生特征的自动化上做多少工作。个人感觉在前面 TFT 的结构化时序数据集接口设计基础上，设计相应的自动化特征工程与模型优化，会是一个能够达到比较好效果的路径。

Update: 22年5月竟然真的出现了一篇文章用了我前面的设想，直接基于 TFT 的架构来做时序预测神经网络架构的搜索，看起来效果还挺好。具体参考：Efficient Automated Deep Learning for Time Series Forecasting，代码实现放在了 Auto-Pytorch 中。

此外，前文中比较少有提到时序数据的各种检测判断，预处理等环节。在之前的 AutoML 比赛中，就有 life-long learning 的设定，即模型的学习预测会随着时间的推移不断有新数据的输入，这也与真实项目的情况非常符合。因此完善的 AutoML 方案中，也需要包含例如 prior shift, covariate shift, concept drift 等方面的检测与处理，以适应复杂的真实预测场景。

可以预见未来会有更多的面向时序问题的 AutoML 框架和产品出现，不断降低使用门槛，扩展相关的应用场景。对这些方向感兴趣的同学也可以多跟我们讨论交流，一起打造行业领先的 AI 产品。