算法模型---时间序列模型

文章来源

1、时间序列

时间序列是时间间隔不变的情况下收集的不同时间点数据集合，这些集合被分析用来了解长期发展趋势及为了预测未来。
时间序列与常见的回归问题的不同点在于:
1、时间序列是跟时间有关的;而线性回归模型的假设：观察结果是独立的在这种情况下是不成立的。
2、随着上升或者下降的趋势，更多的时间序列出现季节性趋势的形式；

常用的时间序列模型有AR模型(Autoregressive model:自回归模型)、MA模型(moving average model:滑动平均模型)、ARMA模型(Auto-Regressive and Moving Average Model:自回归滑动平均模型)和ARIMA模型(Autoregressive Integrated Moving Average Model:自回归积分滑动平均模型)等。

2、时间序列的预处理(使数据平稳化和随机化)

拿到一个观察值序列之后，首先要对它的平稳性和纯随机性进行检验，这两个重要的检验称为序列的预处理。根据检验的结果可以将序列分为不同的类型，对不同类型的序列我们会采用不同的分析方法。

2.1 平稳化处理

平稳就是围绕着一个常数上下波动且波动范围有限，即有常数均值和常数方差。如果有明显的趋势或周期性，那它通常不是平稳序列。序列平稳不平稳，一般采用三种方法检验：

2.1.1 时序图检验

看看上面这个图，很明显的增长趋势，不平稳。

2.1.2 利用自相关系数和偏相关系数

自相关系数和偏相关系数的概念可参考《算法模型— 概率论基础—相关系数相关》
下面是不平稳数据的自相关和偏相关系数的一种情形。

左边第一个为自相关图（Autocorrelation），第二个偏相关图(Partial Correlation)。
平稳的序列的自相关图和偏相关图要么拖尾，要么是截尾。截尾就是在某阶之后，系数都为 0 。怎么理解呢，看上面偏相关的图，当阶数为 1 的时候，系数值还是很大， 0.914;二阶长的时候突然就变成了 0.050. 后面的值都很小，认为是趋于 0 ，这种状况就是截尾。什么是拖尾，拖尾就是有一个缓慢衰减的趋势，但是不都为 0 。
自相关图既不是拖尾也不是截尾。以上的图的自相关是一个三角对称的形式，这种趋势是单调趋势的典型图形，说明这个序列不是平稳序列。
平稳序列的自相关系数会快速衰减。

2.1.3 单位根检验

单位根检验是指检验序列中是否存在单位根，如果存在单位根就是非平稳时间序列。
单位根检验：ADF是一种常用的单位根检验方法，他的原假设为序列具有单位根，即非平稳，对于一个平稳的时序数据，就需要在给定的置信水平上显著，拒绝原假设。ADF只是单位根检验的方法之一，如果想采用其他检验方法，可以安装第三方包arch，里面提供了更加全面的单位根检验方法，个人还是比较钟情ADF检验。以下为检验结果，其p值大于0.99，说明并不能拒绝原假设。

2.1.4 利用差分将序列数据转换为平衡序列

差分可以将数据转换为平稳序列。
一阶差分指原序列值相距一期的两个序列值之间的减法运算；k阶差分就是相距k期的两个序列值之间相减。如果一个时间序列经过差分运算后具有平稳性，则该序列为差分平稳序列，可以使用ARIMA模型进行分析。
确定不平稳后，依次进行1阶、2阶、3阶…差分，直到平稳为止。

2.2 随机化处理

对于纯随机序列，又称白噪声序列，序列的各项数值之间没有任何相关关系，序列在进行完全无序的随机波动，可以终止对该序列的分析。白噪声序列是没有信息可提取的平稳序列。对于平稳非白噪声序列，它的均值和方差是常数。通常是建立一个线性模型来拟合该序的发展，借此提取该序列的有用信息。ARMA模型是最常用的平稳序列拟合模型。

3、平稳时间序列建模

某个时间序列经过预处理，被判定为平稳非白噪声序列，就可以进行时间序列建模。
建模步骤：
（1）计算出该序列的自相关系数（ACF）和偏相关系数（PACF）；
（2）模型识别，也称模型定阶。根据系数情况从AR(p)模型、MA(q)模型、ARMA(p，q)模型、ARIMA（p，d，q）模型中选择合适模型，其中p为自回归项，d为差分阶数，q为移动平均项数。
若平稳序列的偏相关函数是截尾的，而自相关函数是拖尾的，可断定序列适合AR模型；若平稳序列的偏相关函数是拖尾的，而自相关函数是截尾的，则可断定序列适合MA模型；若平稳序列的偏相关函数和自相关函数均是拖尾的，则序列适合ARMA模型。（截尾是指时间序列的自相关函数（ACF）或偏自相关函数（PACF）在某阶后均为0的性质（比如AR的PACF）；拖尾是ACF或PACF并不在某阶后均为0的性质（比如AR的ACF）。）
（3）估计模型中的未知参数的值并对参数进行检验；
（4）模型检验；
（5）模型优化；
（6）模型应用：进行短期预测。

4、混合自回归移动平均过程（公式表达）

差分方程
所谓差分方程即将变量 yt$y_t$与它的滞后期联系起来的表达式。
研究变量在第t期的值记为yt$y_t$。假定给出的动态方程将变量y第t期的值与另外的变量wt$w_t$以及y的前一期联系起来:

yt=ϕyt−1+wt$$y_t=\varphi y_{t-1}+w_t$$

上述称为一阶差分方程是因为仅仅只有变量的一阶滞后(yt−1)$(y_t-1)$出现在方程中。
移动平均
“移动平均”的含义源于Yt$Y_t$是最近两期的ϵ$ϵ$的加权平均。
令ϵt${ϵ}_t$是一个白噪声序列。

Yt=μ+ϵt+θϵt−1$$Y_t=\mu +{ϵ}_t+\theta {ϵ}_{t-1}$$

其中μ和θ可以是任意的常数。这个时间序列称为一阶移动平均过程，记为MA(1)$MA(1)$。
自回归过程
一阶自回归过程
一个一阶自回归，记作AR(1)$AR(1)$满足下面的差分方程

Yt=c+ϕYt−1+ϵt$$Y_t=c+\varphi Y_{t-1}+{ϵ}_t$$

ϵt${ϵ}_t$是一个白噪声序列。
p阶自回归过程
一个p阶自回归，记作AR(p)$AR(p)$满足下式

Yt=c+ϕ1Yt−1+ϕ2Yt−2+...+ϕpYt−p+ϵt$$Y_t=c+{\varphi }_1{Y}_{t-1}+{\varphi }_2{Y}_{t-2}+...+{\varphi }_{p}Yt-p+{ϵ}_t$$

混合自回归移动平均过程
一个ARMA(p,q)$ARMA(p,q)$过程包括自回归和移动平均项：

Yt=c+ϕ1Yt−1+ϕ2Yt−2+...+ϕpYt−p+ϵt+θ1ϵt−1+θ2ϵt−2+...+θqϵt−q$$Y_t=c+{\varphi }_1{Y}_{t-1}+{\varphi }_2{Y}_{t-2}+...+{\varphi }_{p}Yt-p+{ϵ}_t+{\theta }_1{ϵ}_{t-1}+{\theta }_2ϵt-2+...+{\theta }_q{ϵ}_{t-q}$$

5、利用statsmodels库实现时间序列的分析处理

Time Series analysis tsa
statsmodels.tsa.arima_model.ARIMA
AR(I)MA时间序列建模过程——步骤和python代码
Seasonal ARIMA with Python
Python 3中使用ARIMA进行时间序列预测的指南
Python_Statsmodels包_时间序列分析_ARIMA模型
时间序列实战（一）
How to Make Predictions for Time Series Forecasting with Python
How to Use and Remove Trend Information from Time Series Data in Python

python中有专门的库statsmodels可以用来作时间序列分析，但spark中则没有现成的。

#-*- coding: utf-8 -*-
import matplotlib.pyplot as plt
import pandas as pd
data=pd.read_csv('./HIS_MONTH.csv')
b=data['fData']
plt.plot(data['fData'])
plt.show()

数据结构如下

#自相关图
from statsmodels.graphics.tsaplots import plot_acf
plot_acf(data['fData']).show()

如下所示

从图上可以自相关系数基本是拖尾的，稍做调整就可以使用。

#偏自相关图
from statsmodels.graphics.tsaplots import plot_pacf
plot_pacf(data['fData']).show()

如下图所示

偏相关系数也基本符合拖尾。

#平稳性检测
from statsmodels.tsa.stattools import adfuller as ADF
print('原始序列的ADF检验结果为：')
print(ADF(data['fData']))
#返回值依次为adf、pvalue、usedlag、nobs、critical values、icbest、regresults、resstore

如下所示

原始序列的ADF检验结果为：
(-6.561077625309946, 8.378411469638636e-09, 0L, 47L, {'5%': -2.925338105429433, '1%': -3.5778480370438146, '10%': -2.6007735310095064}, 273.5266547319719)

p值远小于0.05，拒绝原假设（原假设认为存在单位根），所以该序列为平衡时间序列。

#白噪声检验
from statsmodels.stats.diagnostic import acorr_ljungbox
print('差分序列的白噪声检验结果为：')
print(acorr_ljungbox(data['fValueData'], lags=1))

结果如下

差分序列的白噪声检验结果为：
(array([0.05198092]), array([0.81965149]))

P=0.81965149，统计量的P值大于显著性水平0.05，则接受原假设（原假设认为为白噪声序列）。

下面对原序列做一阶差分运算

#差分后的结果
D_data = data['fValueData'].diff().dropna()
D_data.columns = ['fValueData']
#时序图
D_data.plot()

结果如下：

从图上看，一阶差分后，感觉数据分布更对称为随机了。

#自相关图
plot_acf(D_data).show()
plt.show()

结果如下

从图上看，一阶差分对数据有改善，振荡衰减更明显

#偏自相关图
plot_pacf(D_data).show()

结果如下

从图上看，偏相关系数虽然没有自相关系数好，但基本能满足平稳序列的要求。

#平稳性检测
print('差分序列的ADF检验结果为：')
print( ADF(D_data))

结果如下

差分序列的ADF检验结果为：
(-4.595500765524432, 0.0001316587309452837, 10L, 36L, {'5%': -2.9459512825788754, '1%': -3.626651907578875, '10%': -2.6116707716049383}, 270.35268975914374)

p值远小于0.05，拒绝原假设（原假设认为存在单位根），所以该序列为平衡时间序列。

#白噪声检验
from statsmodels.stats.diagnostic import acorr_ljungbox
#返回统计量和p值
print('差分序列的白噪声检验结果为：')
print(acorr_ljungbox(D_data, lags=1))

结果如下：

差分序列的白噪声检验结果为：
(array([12.35619393]), array([0.00043953]))

P=0.00043953，统计量的P值小于显著性水平0.05，则拒绝原假设（原假设认为为白噪声序列），所以一阶差分后该序列为非白噪声序列，有进一步分析的价值。

6、Holt-Winters模型原理及python实现

本部分来源，也可以参考该文

6.1、原理分析

公式前面也写过，这里再表述一下，有时间再来整理
假设随机时间序列为：
xt(t=1,2,…,N,N=n)$x_t(t=1,2,\dots ,N,N=n)$
时间序列数据一般有以下几种特点：1.趋势(Trend) 2. 季节性(Seasonality)。
趋势描述的是时间序列的整体走势，比如总体上升或者总体下降。下图所示的时间序列是总体上升的：

季节性描述的是数据的周期性波动，比如以年或者周为周期，如下图：

三次指数平滑算法可以对同时含有趋势和季节性的时间序列进行预测，该算法是基于一次指数平滑和二次指数平滑算法的。

6.1.1、 移动平均(The simple moving average (MA))

直观上,最简单的平滑时间序列的方法是实现一个无权重的移动平均,目前已知的方法是用窗口函数,平滑统计量 St$S_t$就是最近k个观察值的均值。公式如下:

st=1k∑n=0k−1xt−n=xt+xt−1+…+xt−k+1k=st−1+xt−xt−kk$$s_t=\frac{1}{k}\sum _{n=0}^{k-1}{x}_{t-n}=\frac{x_t+x_{t-1}+\dots +x_{t-k+1}}{k}=s_{t-1}+\frac{x_t-x_{t-k}}{k}$$

这样的方法存在明显的缺陷,当k比较小时,预测的数据平滑效果不明显,而且突出反映了数据最近的变化;当k较大时,虽然有较好的平滑效果,但是预测的数据存在延迟。而且最少需要k个值(窗口有限)。

6.1.2、 加权移动平均

一种稍微复杂的方法是先选择一组权重因子来计算加权移动平均

ω1,ω2,…,ωk,满足∑n=1kωn=1$${\omega }_1,{\omega }_2,\dots ,{\omega }_k,满足\sum _{n=1}^k{\omega }_n=1$$

然后用这些权重来计算这些平滑统计量:

st=∑n=1kωnxt−n+1=ω1xt+ω2xt−1+…+ωkxt−k+1$$s_t=\sum _{n=1}^k{\omega }_n{x}_{t-n+1}={\omega }_1{x}_t+{\omega }_2{x}_{t-1}+\dots +{\omega }_k{x}_{t-k+1}$$

在实践中,通常在选择权重因子时,赋予时间序列中的最新数据更大的权重,并减少对旧数据的权重。这个方法也需要最少k个值,并且计算复杂。

6.1.3、 简单指数平滑法

幸运地是有一种方法可以避免上述问题,它叫做指数平滑法。最简单的指数平滑法如下:

si=αxi+(1−α)si−1$$s_i=\alpha x_i+(1-\alpha )s_{i-1}$$

其中α$\alpha$是平滑因子,0<α<1$0<\alpha <1$。si$s_i$是当前时刻的平滑值。换句话说,平滑统计值Si$S_i$是当前统计值Xt$X_t$与上一时间平滑值Si−1$S_{i-1}$加权平均。这个简单指数平滑是很容易被应用的,因为只要有两个观察值就能计算了。这里α的选取,我们可以采用最小二乘来决定α$\alpha$(最小化(st−xt)2$(s_t-x_t{)}^2$)。
一次指数平滑算法进行预测的公式为：

xi+h=si$$x_{i+h}=s_i$$

其中i为当前最后的一个数据记录的坐标，亦即预测的时间序列为一条直线，不能反映时间序列的趋势和季节性。

为什么被称为“指数”平滑法
从它的递推公式就能发现:

简单指数平滑法适用于没有总体趋势的时间序列。如果用来处理有总体趋势的序列,平滑值将往往滞后于原始数据,除非α的值接近1,但这样一来就会造成不够平滑。

6.1.4、 二次指数平滑（HoltWinters-无季节趋势）

为了解决上述问题,于是引出了二次指数平滑,能够保留总体趋势信息。因为将指数平滑应用了两次,所以被称为二次指数平滑。与简单指数平滑相比,二次指数平滑加入了时间趋势统计量ti$t_i$,公式如下:

si=αxi+(1−α)(si−1+ti−1)ti=β(si−si−1)+(1−β)ti−1$$\begin{gathered} s_i=\alpha x_i+(1-\alpha )(s_{i-1}+t_{i-1}) \\ {t}_i=\beta (s_i-s_{i-1})+(1-\beta )t_{i-1} \end{gathered}$$

si$s_i$添加了与一次指数平滑相比，添加了一个新的趋势项ti$t_i$；而ti$t_i$的表达式与一次指数平滑形式上是一样的，只不过用将xi$x_i$换成(si−si−1)$(s_i-s_{i-1})$,趋势可以想像成增量，就是在此刻与上一刻的差值，趋势当然应该用平滑后的数据的差值；同时将一次平滑中的si−1$s_{i-1}$换成ti−1$t_{i-1}$
二次指数平滑的预测公式为

xi+h=si+h ti$$x_{i+h}=s_i+h{t}_i$$

二次指数平滑的预测结果是一条斜的直线

6.1.5、 三次指数平滑

三次指数平滑将时间序列的季节性这一特征也考虑进去了。
季节性被定义为时间序列数据的趋势,它表现出每一个周期重复自身的行为,就像任何周期函数一样。“季节”这个词用来表示行为每隔时间段L就开始自我重复。在自然界中有不同类型的季节性“累加性”(additive)和“累乘性“(multiplicative),就像加法和乘法是数学的基本运算。
如果每个12月都比每个11月多卖出1000套公寓,我们就说这样的季节趋势是“累加性”的。可以用绝对增⻓来表示。如果我们在夏季比冬季多卖出10%的公寓,那么季节趋势在自然中是“累乘性”的。
累乘性公式如下:
下面的式子中k$k$表示季节频率，这个还需要查资料进行理解

si=αxipi−k+(1−α)(si−1+ti−1)ti=β(si−si−1)+(1−β)ti−1pi=γxisi+(1−γ)pt−k$$\begin{gathered} s_i=\alpha \frac{x_i}{p_{i-k}}+(1-\alpha )(s_{i-1}+t_{i-1}) \\ {t}_i=\beta (s_i-s_{i-1})+(1-\beta )t_{i-1} \\ {p}_i=\gamma \frac{x_i}{s_i}+(1-\gamma )p_{t-k} \end{gathered}$$

上式中k为周期， 累乘三次指数平滑的预测公式为：

xi+h=(si+h ti)pi−k+(h mod k)$$x_{i+h}=(s_i+h{t}_i)p_i-k+(hmodk)$$

或者

xi+h=(si+h ti)pi−k+1+(h−1) mod k$$x_{i+h}=(s_i+h{t}_i)p_{i-k+1+(h-1)modk}$$

两个公式还需要确认下
累乘性公式初始值的计算

sk=1k(s1+s2+⋯+sk)tk=1k[xk+1−x1k+xk+2−x2k+⋯+xk+k−xkk]p1=x1sk,p2=x2sk,⋯,pk=xksk$$\begin{gathered} s_k=\frac{1}{k}(s_1+s_2+\cdots +s_k) \\ {t}_k=\frac{1}{k}\left[\frac{x_{k+1}-x_1}{k}+\frac{x_{k+2}-x_2}{k}+\cdots +\frac{x_{k+k}-x_k}{k}\right] \\ {p}_1=\frac{x_1}{s_k},p_2=\frac{x_2}{s_k},\cdots ,p_k=\frac{x_k}{s_k} \end{gathered}$$

累加性公式如下:

si=α(xi−pi−k)+(1−α)(si−1+ti−1)ti=β(si−si−1)+(1−β)ti−1pi=γ(xi−si)+(1−γ)pi−kxi+h=si+hbi+pi−k+h mod k$$\begin{gathered} s_i=\alpha (x_i-p_{i-k})+(1-\alpha )(s_{i-1}+t_{i-1}) \\ {t}_i=\beta (s_i-s_{i-1})+(1-\beta )t_{i-1} \\ {p}_i=\gamma (x_i-s_i)+(1-\gamma )p_{i-k} \\ {x}_{i+h}=s_i+h{b}_i+p_{i-k+hmodk} \end{gathered}$$

或者

xi+h=si+hbi+pi−k+1+(h−1) mod k$$x_{i+h}=s_i+h{b}_i+p_{i-k+1+(h-1)modk}$$

两个公式还需要确认下
累加性公式初始值的计算

sk=1k(s1+s2+⋯+sk)tk=1k[xk+1−x1k+xk+2−x2k+⋯+xk+k−xkk]p1=x1−sk,p2=x2−sk,⋯,pk=xk−sk$$\begin{gathered} s_k=\frac{1}{k}(s_1+s_2+\cdots +s_k) \\ {t}_k=\frac{1}{k}\left[\frac{x_{k+1}-x_1}{k}+\frac{x_{k+2}-x_2}{k}+\cdots +\frac{x_{k+k}-x_k}{k}\right] \\ {p}_1=x_1-s_k,p_2=x_2-s_k,\cdots ,p_k=x_k-s_k \end{gathered}$$

其中 α$\alpha$是数据平滑因子, 0<α<1$0<\alpha <1$;β$\beta$是趋势平滑因子,0<β<1$0<\beta <1$; γ$\gamma$是季节改变平滑因子0<γ<1$0<\gamma <1$。
α，β，γ$\alpha ，\beta ，\gamma$的值都位于[0,1]之间，可以多试验几次以达到最佳效果。

s,t,p$s,t,p$初始值的选取对于算法整体的影响不是特别大，通常的取值为s0=x0,t0=x1−x0$s_0=x_0,t_0=x_1-x_0$,累加时p=0$p=0$,累乘时p=1$p=1$.
对三次指数平滑法而言,我们必须初始化一个完整的“季节”pi$p_i$的值,不过我们可以简单地设置为全1(针对累乘式)或全0(针对累加式)。只有当序列的⻓度较短时,我们才需要慎重考虑初始值的选取。
我们这里讲的Holt-Winters模型就是三次指数平滑法。哇,终于切入正题了。
所有的指数平滑法都要更新上一时间步⻓的计算结果,并使用当前时间步⻓的数据中包含的新信息。它们通过“混合”新信息和旧信息来实现,而相关的新旧信息的权重由一个可调整的拌和参数来控制。各种方法的不同之处在于它们跟踪的量的个数和对应的拌和参数的个数。三次指数平滑法,功能最强大,既能体现趋势性又能体现季节性,所以三次指数平滑法的参数最多,有三个。
下图为使用累加三次指数平滑进行预测的效果：其中红色为源时间序列，蓝色为预测的时间序列，α，β，γ$\alpha ，\beta ，\gamma$的取值为0.45，0.2，0.95：

下图为累乘三次指数平滑进行预测的效果，α，β，γ$\alpha ，\beta ，\gamma$的取值为0.4，0.05，0.9：

可以看到三次指数平滑算法可以很好的保存时间序列数据的趋势和季节性信息，在International Airline Passengers数据集上累乘平滑指数算法的效果更好。

6.1.6、 python代码实现

我们知道HoltWinters模型有三个可调参数,我们的目的就是训练出有效的α,β, γ
。我们有两种方法,一种就是自己取值来试试,一种就是采用数值优化的思想,比如
前面我们提到的最小二乘来最小化误差来求参数(注意不一定能全局收敛!这个问题
实在是让人头痛。。。)我们就采用最小二乘法(L-BFGS)。

RMSE的实现

7、基于SVM作短期时间序列的预测

传统的做法是提取1、2、3、4、5、7、9、13个单位时间的数据作为特征进行预测；
举个例子进行分析，比如每天都有口香糖的销量，那么如何通过几周的数据预测明天的数据，
就可以选择前1、2、3、4、5、7、14天的数据作为特征，从而预测明天的数据，
通过构建特征，再选择核函数进行预测，其中调参的参数尽量要进行最优化，
参考方法：如果选择RBF核函数，那么其中就会有三个参数，固定两个，然后不停的优化另外一个，直到得到最优解。

具体应用的例子：

（1）SVM预测风场：http://wenku.baidu.com/link?url=SCCIJJe8tXLbTjLMZ81x5Qy6elsceAKIOwtkZ0QxfSCQQ4KaWKwo8Biepjs3Ss2LJ2ewhisNR0ixrDY4kV1Rd7BcqWRenuTaG85K80E-30y
（2）SVM预测股票指数：基于SVM修正的模糊时间序列模型在沪指预测中的应用
（3）SVM预测时间序列其他方面：http://www.docin.com/p-233353900.html

8、LSTM模型分析及对时序数据预测的具体实现（python实现）

来源：2017年09月30日 10:28:08

9、时间序列的并行实现

spark里面的库是没有时间序列算法的，但是国外有人已经写好了相应的算法。其github网址是:https://github.com/sryza/spark-timeseries
sryza/spark-timeseries
Spark-TimeSeries使用方法
A New Library for Analyzing Time-Series Data with Apache Spark
【Spark Summit East 2017】使用Spark进行时间序列分析

其他资料

时间序列挖掘-预测算法-三次指数平滑法(Holt-Winters)
Holt-Winters原理和初始值的确定