我这里补充一篇之前我做的时间序列部分的工作记录，留个纪念怕忘。
这里用的数据说我们公司自己的数据，分析过程是随着统计学的书来的，当时对时间序列的理解不是很充分，如果有错误，欢迎指出。
1,时间序列分析的流程
首先，在数据分析方面数据预处理是极其的重要，甚至比真正的数据分析还重要，因为数据预处理关乎着你的数据分析过程能否有正确的结果。数据预处理包括缺失数据不全，异常值处理等。想得到正确的模型，就需要对数据中缺失的数据进行补全，对异常增大或者异常减小的数据进行处理，常见的有移动平均。
其次，在拿到正确且完整的数据之后，我们要做的就是先画时序图，看数据的分布，有明显上升或者下降趋势或者有明显周期性的数据是非平稳的。
第三，当时序图不能显然的得出结论的时候就需要通过自相关图和偏自相关图进行判断。如果数据的统计学性质比较好，通过自相关和偏自相关图不但能判断平稳还能直接得出模型的阶数也就是ARMA(P,Q)中的PQ值，进而确定了模型。
| 自相关 | 偏自相关 | 模型 |
|---|---|---|
| 拖尾 | p阶截尾 | AR(P) |
| q阶截尾 | 拖尾 | MA(Q) |
| 拖尾 | 拖尾 | ARMA(P,Q) |
第四，当然通过图形判断还是比较主观的，还有一个从数值上判断平稳性的方法就是单位根检验（即ADF检验），这里给出结论，ADF检验的p_value值大于0.05就证明有单位根，不平稳。反之，当p_value值小于0.05就证明无单位根，平稳。这个p_value值的意思就是我们有多大的把握认为他有单位根，0.05就是统计学上的95%置信区间。落在置信区间之内就认为接受，否则就是拒绝。
第五，如果我们的数据是平稳序列，我们还要进行白噪声检验，即Q统计量检验或者LB统计量检验。同上，p_value小于0.05我们认为不是白噪声，p_value大于0.05我们认识是白噪声。如果证明是白噪声就停止分析了，因为白噪声是随机的不可分析的。
第六，如果我们的序列是非平稳序列，那就不用进行白噪声检验，因为肯定不会白噪声的。非平稳序列还要考虑是什么样的非平稳。带有明显上升下降趋势的，可以通过低阶差分（1,2,3阶差分），以消除趋势达到平稳。带有周期的数据可以通过以周期为步长的差分消除周期性，达到平稳。所有的非平稳序列都要处理成平稳序列才能进行建模。上面属于书中的非平稳序列的随机分析，趋势和周期只有一种。
第七，非平稳序列的确定性分析。X-11季节模型。可以将数据分为趋势部分，周期性部分（就是季节部分），还有残差部分。我们在处理的过程中。只需要将趋势部分做模型，然后将周期部分和残差部分相叠加就可以了。周期部分是规律的不变的，不用做处理。残差部分应该是零均值白噪声。
第八，当模型确定好之后，我们要对模型的优劣进行选择，因为有可能不止一种模型。这里有一个判断准则，AIC准则，BIC准则（也就是SBC准则），判断每个自变量的重要性。简化模型。
这里对残差还有一些处理，根据残差的一些性质也可以拟合残差的若干种模型。以达到效果的精准。
2,自相关图判断平稳性
from pandas.plotting import autocorrelation_plot
autocorrelation_plot(timeseries)
计算时间序列的自相关图
因为我不会设置自相关的延迟阶数，但是至少可以看到5000以后基本为0.
所以我决定把数据截取前100个看一下自相关图
上图中实线是95%置信区间，虚线是99%置信区间（我看了源码），证明自相关系数基本为0，无相关性？？？？
我不甘心有画了800期延迟的自相关图，又发现有了超过95%置信区间的数据了
顿悟了。重新画了全部数据的自相关图，并且放大：
很多网上的资料都没说清楚，只是直接的给出了结论。其实上面的源码写的很清楚。放大看就可以看出来，图上自带的95%置信区间和99%置信区间了。
从上图可以看出来，一直到两万个点附近，自相关才逐渐收敛到95%置信区间之内。其实书上有这么一个结论，只要延迟的阶数足够大，自相关系数总会收敛到0的。
两万个点就是接近一个月的数据。还可以看出来前5000个点（7天）自相关系数显示相关性最强。而且自相关系数有周期性，可以看出以一天为一个周期（720个点）。
所以此时要做的就是通过差分，消除周期性。
3,单位根（ADF）检验平稳性
检验是否平稳还可以通过单位根检验。
首先先看数据的分布图。
看样子是在某一数字周围震荡，基本平稳，先用单位根检验一下平稳性。
from statsmodels.tsa.stattools import adfuller
dftest = adfuller(timeseries, autolag='AIC')
# 这里打印的还是adf检验的p_value，很小（小于0.05），证明平稳
# TestStatistic
print(dftest[0])
# p_value
print(dftest[1])
# 使用的延迟阶数
print(dftest[2])
print(dftest[3])
print(dftest[4])
print(dftest[5])
-10.3959989166 #adf统计量的值
1.96457899653e-18 #p值
46 #延迟阶数
24799
{'1%': -3.4306137193843012, '5%': -2.8616565559416833, '10%': -2.566832039327017}
214670.53643 #99%，95%，90%置信区间的临界值（adf统计量）对照
[ 9.10487879e-273] #LB统计量的p值
[ 1245.21359197] #LB统计量的值
这是打印的结果，根据p值可以看出，p值比0.05（95%置信区间）小得多，因此应该不存在单位根，序列应该是平稳的才对。
可是对照看一下adf统计量的值，计算的adf统计量的值比三种置信区间得到的临界统计量的值都小得多，因此，应该是不平稳的才对。
源码里有一段note是这么说的，如果p值接近于0.05，那么应该依据统计量和临界值判断。但是这个p值明显比0.05小的多。网上的数据都是实验数据，统计学性质表现很好。而且网上都是根据p值判断的。参照之前的自相关图，我虽然知道了，我的数据不可能是平稳的，（现实生活中数据几乎不可能是平稳的）但是我也不知道这个p值和adf统计量的值怎么看，而且adf统计量的值网上都是正值，书上有负值，但是也没说用统计量怎么判断，书上也是p值判断的。我也不懂这里应该是用绝对值判断还是数值判断。（留个疑问在这里，以后搞懂了再来补充）。
搞懂了，参见11，ADF检验补充。
为了搞清楚，我决定再用tps数据做一下试验，计算一下adf检验的值。
这是tps数据的分布图，显然是不平稳的。
-0.66255326589
0.856212167775
18
581
{'1%': -3.4416553818946145, '5%': -2.8665274458710064, '10%': -2.5694261699959413}
4709.61042801
[ 2.68092399e-100] #LB统计量的p值
[ 451.97508501] #LB统计量的值
这里的p值就是网上的很显然的那种数字，显然比0.05大，就是说不平稳。但是这里的adf统计量的值-0.66就比后面给的对照的临界值的-3.44，-2.866，-2.56数值上大，但是绝对值上小。也就是说各方各面都跟baseline数据相反，我也不懂了。
4,白噪声判断是否是纯随机序列（LB）
如果是平稳序列那么就要判断是否是白噪声序列，要进行纯随机性检验（LB统计量）。
from statsmodels.stats.diagnostic import acorr_ljungbox
p_value = acorr_ljungbox(timeseries, lags=1)[1]
这里打印的是LB统计量的p值，可以看出，p值都很小（小于0.05），根据书上的理论，毫无疑问不是白噪声，可以分析。说毫无疑问是因为，（可以跟进去看一下源码）这里的参数很简单，返回值也很简单。
参数，时间序列，延迟阶数。
返回值，LB统计量值，LB统计量的p值，Q统计量值，Q统计量的p值。所以基本不会出错。
5,确定ARMA的阶数
import statsmodels.tsa.stattools as st
from statsmodels.tsa.arima_model import ARMA
# timeseries是时间序列,max_ar、max_ma是p、q值的最大备选值。
# 通过ar、ma的最大迭代次数反复迭代
# 在'aic', 'bic', 'hqic'在这三个准则之下
# 返回值是许多的ARMA模型
order = st.arma_order_select_ic(timeseries, max_ar=5, max_ma=5, ic=['aic', 'bic', 'hqic'])
# 这句话的意思是，返回以BIC准则(BIC最小，BIC就是SBC)确定的阶数,是一个tuple类型
order = order.bic_min_order
# 把时间序列和确定的阶数赋值给ARMA，得到一个模型
model = ARMA(timeseries, order)
print(order)
# 用模型来拟合，不过这个fit为啥没有输入数据
result_arma = model.fit(disp=-1, method='css')
print(result_arma)
# 预测，这个预测也没有输入数据，困惑
predict_ts = result_arma.predict()
print(predict_ts)
为什么我之前说，我的数据应该不平稳，首先是因为现实数据不太可能平稳。其次就是，上面的语句可以根据BIC准则，反复迭代，自动找到最优的ARMA的pq值，但是我在运行的时候，总是会提示“最大似然函数无法收敛”。当我手动给出一组低阶参数时，拟合效果极差。尝试了几组，基本都是下面这种效果，和之前用机器学习算法差不多，无法拟合突变，总是想着回归出一个趋势。所以，现在要研究非平稳序列消除周期性。
下图是使用机器学习算法中效果相对最好的，随机森林回归。因为数据变化很大，不平滑，大家都知道数据发生突变对回归算法来说是非常不利的。
6,以周期为步长，做一阶差分
网上关于周期数据该如何处理的东西很少。所以要自己研究了。
下图是做了一阶查分之后的图。（并没有以步长为周期，就是普通的一阶差分）
仔细看还是有很多点在95%置信区间之外。
还有，在5000以内，可以看到有一个点非常的突出，根据数据知识判断这个点应该是以每天为周期的第720个点左右。于是放大坐标看。
结论得到验证。因此一阶差分后的数据还是有很多点在95%置信区间之外，因此，还是不可用。
现在进行以周期为步长的一阶差分
可以看出，在720个点左右，自相关仍旧很大，还是有很多点都在95%置信区间的边界之外。仍旧不能认为是平稳的。由于偏自相关的计算时间比较缓慢，因此，只截取前1000个点计算一下偏自相关。可以看出，偏自相关的性质能稍好一点，但是在720个点位置，还是出现了超出置信界的现象。
下图是进行了两次周期差分的自相关图：
可以看出来，经过两次差分效果依然非常差劲，依然非常的不平稳。
7,缺失数据的处理
虽然我没有进行专业的数据补全的处理，但是仍旧在这里绕了一大圈。
首先一开始我觉得检查数据太麻烦了，因为数据每天720个点特别多，如果缺失不多，也不影响预测结果。一直是这么想的，结果，今天我走投无路，实在没办法了，就想着还是老老实实先补数据吧。我把数据以一天为单位检查了一下。结果发现前7天的数据根本要不得，直接就缺失太多了，每天都是400-500个点缺失太多了肯定是没办法补的，就算是勉强补全了也会和真实值差距太远影响结果。
不过好在后面的数据还是比较整齐的。后面将近30天的数据，基本都是720个点，只有两天有缺失，一个缺6个，一个缺1个，这种量级是可以补全的。因为比较少，我手动补全了，就是用前一天和后一天同一时刻的数据平均了一下，很简单。现在再用正确的数据用机器学习算法预测一下。留最后一天数据做test。前面数据做train。
现在用新的数据重新最分析。
自相关图：没啥变化，仍旧是周期性，非常的不平稳
感觉一切结果跟之前都差不多。单位根检验依旧有疑惑。白噪声检验仍旧证明不是白噪声。
进行一阶差分（非周期）
自相关图显示，仍旧有很多点在95%置信区间之外，不平稳。
因为偏自相关的计算速度比较慢，所以，考虑到包含了数据性质，选择了前800个点，偏自相关图显示出明显的截尾性质。
下图是周期差分的自相关图
周期差分的偏自相关
这是周期差分的单位根检验和白噪声检验，到底平稳不？
从自相关图上看应该是不平稳的。
8,季节分解
虽然把数据补全了，但是效果不如料想的那么好，仍旧困难重重，不知道如何处理。显然我的数据是周期性数据，处理周期数据一定要先去除周期性。去除周期性的方法，首先是以周期为步长的差分，我做了效果不好。然后现在考虑的就是对数据进行季节分解。
from pandas import read_csv
import numpy as np
from test_stationarity import *
from statsmodels.tsa.seasonal import seasonal_decompose
# 读训练数据
ts = read_csv('E:\\workfile\\data\\baseline\\time\\train.csv')
def qushifenjie():
# 将数据分解成趋势、季节、残差
# statsmodels使用的是X-11分解，因此seasonal_decompose也有两种分解方法
# 加法additive，乘法multiplicative
decomposition = seasonal_decompose(np.array(ts), model="additive", freq=720, two_sided=False)
trend = decomposition.trend
seasonal = decomposition.seasonal
residual = decomposition.resid
# 画图
decomposition.plot()
plt.show()
# 打印保存趋势序列
# print(trend)
# np.savetxt('E:\\workfile\\data\\baseline\\time\\trend', trend, fmt='%s')
图片是通过这些语句得到的。
第一部分是时序图，数据本身的分布图。
第二部分是趋势部分，trend。
第三部分是周期部分，seasonal。
第四部分是残差部分，residual。
因为我不是很会使用这些API，因此决定。将trend，seasonal，residual单独保存起来，单独分析，要用的时候在读文件。我虽然知道麻烦，但是目前只好这样了。
分析趋势的统计属性。
这是趋势的自相关图，显示出拖尾的性质，逐渐减少，收敛到两倍标准差范围之内，效果明显比之前好多了。长舒一口气。现在看偏自相关。
因为偏自相关图计算的比较慢。不失一般性的画了前1440个点，也就是说两天的数据。
明显呈现出截尾的性质。从这两个图上看出来统计学的性质就非常的明显了。自相关拖尾，偏自相关截尾，明显的AR(P)模型
给一个放大图，可以看出来，偏自相关2阶截尾。因此，这是典型的AR(2)模型。
单位根检验和白噪声检验的结果也非常的明确，显然不平稳，显然不是白噪声。
按道理来说，通过上述语句提取趋势，季节，残差之后。趋势应该是平稳非白噪声序列。季节是周期序列。残差是零均值白噪声序列。网上的例子也是这么做的，直接就给趋势部分的数据求ARMA模型了。我也这么做了，代码提示报错，错误是，你的数据不是平稳序列，请差分后再来求解。我也没办法了，只好差分。但是这里有个问题，提取了趋势之后，趋势和残差部分数据，由于计算的原因，已经缺失了前一个周期的数据了，即720个点。如果再进行周期差分，差分一次少一个周期的点。这样对数据的完整性来说是非常不利的。
走到这里又进了死胡同，由于趋势部分也是不平稳的，完全没办法进行，于是尝试对趋势数据进行一阶差分（普通的一阶差分），将差分后的数据做自相关图，发现，自相关图在720个点附近还是超过了两倍的标准差范围。。。
这是两次差分的结果：还是在720个点附近超过了2倍标准差。
三次差分的结果：
下图是残差的检验结果。平稳的非白噪声序列。
我不甘心，对trend数据又进行了一次周期为步长的差分，显示自相关效果差不多
这个偏自相关是啥情况？？？完全没有了。
如果数据随着时间的推移变化幅度越来越大，那就要用乘法，如果变化幅度没有发生变化就要用加法，这里的数据变化幅度并没有发生变化，理论上应该是加法。
我还是不甘心，季节分解的方法有四种，这个语句支持两种即加法和乘法。
又把季节分解的加法，变成了乘法，没有任何作用。下面是分解的图，一模一样的。是的，我还是不甘心，又把trend文件保存了一份，打开看了一下，是的，一模一样的。加法分解和乘法分解生成的trend文件是一模一样的。
-------------------------中间省略了我用随机森林的部分----------------------
11,ADF检验补充
之前的检验结果显示序列是平稳的，但是序列实际上是带有周期性的，我终于搞懂问题在哪了。
跟进去看一下，有一个regression的参数，在王燕的书里，6.3.2ADF检验部分，提到了，ADF检验可以用于三种类型的单位根检验（代码里分了四种）
regression : {'c','ct','ctt','nc'}
Constant and trend order to include in regression
* 'c' : constant only (default)
* 'ct' : constant and trend
* 'ctt' : constant, and linear and quadratic trend
* 'nc' : no constant, no trend
第一种是默认的，仅有常数均值。
第二种是有常数均值和趋势。
第三种是有常数均值，有线性和二次趋势。
第四种是无常数均值，无趋势的。
这里我的数据是第四种，而我不知道这个参数是干什么的，一直用的是默认的。选第四种以后。这是检验结果：
可见p_value值比0.05大，这时候的τ统计量值在5%-10%之间。不平稳。
得出结果使用的延迟阶数为44。
AIC准则的值为20115。
12,季节分解后的残差分析
昨天看到书上有一个例题是这么做的，季节差分后的结果并不是对趋势进行分析，而是对残差进行分析，于是我今天也尝试这么做。
这是残差的自相关图，仔细看，有很多数据是在2倍标准差之外的。
这是自相关图的放大图。可以清楚的看到，有很多数据是在2倍标准差之外的，而且绵延期数很长，显而易见是不平稳的。
这是延迟800期的偏自相关。
这是偏自相关的局部放大图。非常的不平稳。
