0、前言

什么是时间序列
时间序列简单的说就是各时间点上形成的数值序列，时间序列分析就是通过观察历史数据预测未来的值。在这里需要强调一点的是，时间序列分析并不是关于时间的回归，它主要是研究自身的变化规律的（这里不考虑含外生变量的时间序列）。
时间序列的分析的步骤是先对数据进行平稳性和非白噪声检验（如不满足需对数据进行平滑或差分等预处理），然后才是模型调参跟预测，因此本文分为2大部分介绍，答题思路见脑图。

一、数据准备&探索

1、平稳性

序列平稳性是进行时间序列分析的前提条件为什么要满足平稳性的要求呢？在大数定理和中心定理中要求样本同分布（这里同分布等价于时间序列中的平稳性），而我们的建模过程中有很多都是建立在大数定理和中心极限定理的前提条件下的，如果它不满足，得到的许多结论都是不可靠的。以虚假回归为例，当响应变量和输入变量都平稳时，我们用t统计量检验标准化系数的显著性。而当响应变量和输入变量不平稳时，其标准化系数不在满足t分布，这时再用t检验来进行显著性分析，导致拒绝原假设的概率增加，即容易犯第一类错误，从而得出错误的结论。
判断一个序列是不是平稳序列有三个评判标准：

1. 均值 ，是与时间t 无关的常数。
2. 方差 ，是与时间t 无关的常数。这个特性叫做方差齐性。
3. 协方差 ，只与时期间隔k有关，与时间t 无关的常数。

如下所示分别为稳定的数据和均值不稳定、方差不稳定、自协方差随时间变化

1.1 平稳性检验

时间序列平稳性检验主要是观察法和单位根检验法，观察法
观察法，通俗的说就是通过观察序列的趋势图与相关图是否随着时间的变化呈现出某种规律。所谓的规律就是时间序列经常提到的周期性因素，现实中遇到得比较多的是线性周期成分，这类周期成分可以采用差分或者移动平均来解决，而对于非线性周期成分的处理相对比较复杂，需要采用某些分解的方法。

一般都是使用单位根检验法 Dickey-Fuller Test 进行判断，即在一定置信水平下，假设时序数据Null hypothesis: 非稳定、序列具有单位根。因此对于一个平稳的时序数据，就需要在给定的置信水平上显著，拒绝原假设。如下通过statsmodels包中的adfuller来进行检测：
数据下载

import os
import sys
import pandas as pd
import pandas_datareader.data as web
import numpy as np
import statsmodels.formula.api as smf
import statsmodels.tsa.api as smt
import statsmodels.api as sm
import scipy.stats as scs
from arch import arch_model
import matplotlib.pyplot as plt
import matplotlib as mpl
from statsmodels.stats.diagnostic import acorr_ljungbox #白噪声检验
#######读取数据等
data_ts = pd.read_csv("test0923.csv")
data_ts.columns = ['week_no','send_date','day_week','send_date_pre','export_type','depart_thr_letter_code','send_weight','real_weight_qty','meterage_weight_qty','send_weight_pre','real_weight_qty_pre','meterage_weight_qty_pre']data_ts.columns
data_ts.head()
data_ts['send_date'] = pd.to_datetime(data_ts['send_date'])
data_ts.dtypes
data_ts.index=data_ts['send_date']
df_data_ts=data_ts[data_ts['export_type']=='HKG']
#ads=df_data_ts['send_weight'][-100:]
ads=df_data_ts['send_weight'][-65:]###########ADF检验adftest = adfuller(ads, autolag='AIC')
print (adftest)
adfuller(np.diff(ads), autolag='AIC')#检查不通过用一阶差分后的数据看看，即np.diff(ads)

结果如下：

如何确定该序列能否平稳呢？主要看：
（1）1%、%5、%10不同程度拒绝原假设的统计值和ADF Test result的比较，ADF Test result同时小于1%、5%、10%即说明非常好地拒绝该假设，本数据中，adf结果为-10.6， 小于三个level的统计值，可以拒绝原假设，即数据是平稳的。
（2）P-value是否非常接近0.本数据中，P-value 为 4.8e-19，接近0。

1.2 数据处理（平滑、变换、差分、分解）

如果，序列不平稳，或者数据预测不是很好，就要考虑进行数据处理。处理方法主要有以下四种：

1.2.1 对数变换

对数变换主要是为了减小数据的振动幅度，使其线性规律更加明显。对数变换相当于增加了一个惩罚机制，数据越大其惩罚越大，数据越小惩罚越小。这里强调一下，变换的序列需要满足大于0，小于0的数据不存在对数变换。

ads_log = np.log(ads)
adfuller(ads_log, autolag='AIC')

1.2.2 平滑法（移动平均&指数平均）

移动平均即利用一定时间间隔内的平均值作为某一期的估计值，也有升级版的加权平均。参考代码，如下：

############ 移动平均图----------
def draw_trend(timeSeries, size):f = plt.figure(facecolor='white')# 对size个数据进行移动平均,min_periods=1 则前面size数不为空rol_mean = timeSeries.rolling(window=size,min_periods=1).mean()# 对size个数据进行加权移动平均# 加权移动平均是越近是越远的2倍#rol_weighted_mean = pd.ewma(timeSeries, span=size)rol_weighted_mean=timeSeries.ewm(span=size,min_periods=1).mean()timeSeries.plot(color='blue', label='Original')rol_mean.plot(color='red', label='Rolling Mean')rol_weighted_mean.plot(color='black', label='Weighted Rolling Mean')plt.legend(loc='best')plt.title('Rolling Mean')plt.show()return rol_meanads_rol_mean=draw_trend(ads,3)    

用平均法的数据带入模型后预测结果需要还原

######移动平均还原
ads_rec=ads_rol_mean*3 - ads.rolling(window=2,min_periods=1).sum().shift(1)
##前size个数值需专门计算 ##or 另一种还原方法（直接计算肯定更不容易出错）
ads_rec=ads_rol_mean
for i in range(len(ads_rol_mean)) :print(i)if i<=1 :ads_rec[0]=ads_rol_mean[0]ads_rec[1]=ads_rol_mean[1]*2-ads_rol_mean[0]else :ads_rec[i]=ads_rol_mean[i]*3-ads_rec[i-1]-ads_rec[i-2]

指数平均则是用变权的方法来计算均值。具体怎么做的呢？一个简单的数学式为：
$$y_t^s=\mu \ast y_t+(1-\mu )\ast y_{t-1}^s,$$
式子中的 $\mu$表示平滑因子，它定义我们“遗忘”当前真实观测值的速度有多快。 $\mu$越小，表示当前真实观测值的影响力越小，而前一个模型预测值的影响力越大，最终得到的时间序列将会越平滑。
那么指数体现在哪呢？指数就隐藏在递归函数之中，我们上面的函数，每次都要用( 1 − $\mu$)乘以模型的上一个预测值。
代码可参考如下（自己项目未用的到，代码来自这篇博客）：

def exponential_smoothing(series, alpha):"""series - dataset with timestampsalpha - float [0.0, 1.0], smoothing parameter"""result = [series[0]] # first value is same as seriesfor n in range(1, len(series)):result.append(alpha * series[n] + (1 - alpha) * result[n-1])return resultdef plotExponentialSmoothing(series, alphas):with plt.style.context('seaborn-white'):    plt.figure(figsize=(15, 7))for alpha in alphas:plt.plot(exponential_smoothing(series, alpha), label="Alpha {}".format(alpha))plt.plot(series.values, "c", label = "Actual")plt.legend(loc="best")plt.axis('tight')plt.title("Exponential Smoothing")plt.grid(True);plotExponentialSmoothing(ads.Ads, [0.3, 0.05])
plotExponentialSmoothing(currency.GEMS_GEMS_SPENT, [0.3, 0.05])

上述是但指数平滑法，还有双指数平滑法等，可参考上面的博客。

1.2.3 差分法

时间序列最常用来剔除周期性因素的方法当属差分了，它主要是对等周期间隔的数据进行线性求减。
参考代码如下：

# 进行差分
ads=pd.Series(ads) ###也可不用pd.Series，用np.diff(ads)
ads.head()
ads_diff=ads.diff(1)
ads_diff.head()
ads_diff=ads_diff.dropna()############## 差分还原
diff_shift=ads.shift(1)
ads_recov=ads_diff.add(diff_shift)
ads_recov[0]=ads[0]
# or第二种方法
ads_recov=ads_diff  ###看0 is nan?
for i in range(len(ads_diff)) :print(i)if i<=0 :ads_recov[0]=ads[0]else :ads_recov[i]=ads_recov[i]+ads_recov[i-1]
############

1.2.4 分解

分解就是将时序数据分离成不同的成分。statsmodels使用的X-11分解过程，它主要将时序数据分离成长期趋势、季节趋势和随机成分。与其它统计软件一样，statsmodels也支持两类分解模型，加法模型和乘法模型，这里我只实现加法，乘法只需将model的参数设置为"multiplicative"即可。本人未用该方法，代码参考于这篇文章

from statsmodels.tsa.seasonal import seasonal_decompose
decomposition = seasonal_decompose(ads, model="additive")trend = decomposition.trend
seasonal = decomposition.seasonal
residual = decomposition.resid

2、非白噪声检验

非白噪声检验也可以用初步观察的方法，可用如下代码观察残差图等：

def tsplot(y, lags=None, figsize=(10, 8), style='bmh'):if not isinstance(y, pd.Series):y = pd.Series(y)with plt.style.context(style):    fig = plt.figure(figsize=figsize)#mpl.rcParams['font.family'] = 'Ubuntu Mono'layout = (3, 2)ts_ax = plt.subplot2grid(layout, (0, 0), colspan=2)acf_ax = plt.subplot2grid(layout, (1, 0))pacf_ax = plt.subplot2grid(layout, (1, 1))qq_ax = plt.subplot2grid(layout, (2, 0))pp_ax = plt.subplot2grid(layout, (2, 1))y.plot(ax=ts_ax)ts_ax.set_title('Time Series Analysis Plots')smt.graphics.plot_acf(y, lags=lags, ax=acf_ax, alpha=0.5)smt.graphics.plot_pacf(y, lags=lags, ax=pacf_ax, alpha=0.5)sm.qqplot(y, line='s', ax=qq_ax)qq_ax.set_title('QQ Plot')        scs.probplot(y, sparams=(y.mean(), y.std()), plot=pp_ax)plt.tight_layout()return 

非白噪声检验代码如下：

###########纯随机性检验（白噪声检验）
p_value = acorr_ljungbox(ads, lags=1) 
print (p_value)
acorr_ljungbox(np.diff(ads), lags=1)

结果如图：

统计量的P值小于显著性水平0.05，则可以以95%的置信水平拒绝原假设，认为序列为非白噪声序列（否则，接受原假设，认为序列为纯随机序列。）

由于P值为3.9e-5远大于0.05所以拒绝原假设，认为时间序列是非白噪声的，即非是随机产生的序列，具有时间上的相关性。
从严格意义上说，只有通过平稳性检验和非白噪声检验才能用时间序列模型。不过很多博客都没有非白噪声检验这一步，可能这跟数据分析挖掘最终只看效果有关吧。

二、模型（ARIMA）

1、 自回归移动平均模型（ARMA）

ARIMA比ARMA模型多了I这个差分，其他一样。自回归移动平均模型(ARMA(p，q))是时间序列中最为重要的模型之一，它主要由两部分组成： AR代表p阶自回归过程，MA代表q阶移动平均过程，其公式特征如下：

1. AR§：autoregression model(自回归模型)，即时间序列对自身的回归。基本假设是当前的序列值取决于它之前的值，并且存在一定的滞后。 模型最大的滞后值称为 p ，要确定初始的 p 值，我们需要在 PACF 图中找到最大的滞后点。
2. MA(q)：moving average model(移动平均值模型)，这里不讨论它的细节，它基于当前的误差依赖于先前的误差，有一定的滞后性（滞后值记为q） 这一假设，对时间序列进行建模。和上面一样，可以在 ACF 图找到滞后值 q qq 的初始值；

详细例子和说明可参考该文章

2、 模型参数确定

根据ARIMA模型理论很容易知道我们需要确定 p(自回归滞后)、d（差分）、q（误差滞后）三个参数。主要有如下三种方法

2.1 PACF和ACF 图 找到

p值可以通过 PACF 图找到，q值可以通过 ACF 图找到。d可以根据差分后的平稳性检验得到。参考代码如下：

def tsplot(y, lags=None, figsize=(10, 8), style='bmh'):if not isinstance(y, pd.Series):y = pd.Series(y)with plt.style.context(style):    fig = plt.figure(figsize=figsize)#mpl.rcParams['font.family'] = 'Ubuntu Mono'layout = (3, 2)ts_ax = plt.subplot2grid(layout, (0, 0), colspan=2)acf_ax = plt.subplot2grid(layout, (1, 0))pacf_ax = plt.subplot2grid(layout, (1, 1))qq_ax = plt.subplot2grid(layout, (2, 0))pp_ax = plt.subplot2grid(layout, (2, 1))y.plot(ax=ts_ax)ts_ax.set_title('Time Series Analysis Plots')smt.graphics.plot_acf(y, lags=lags, ax=acf_ax, alpha=0.5)smt.graphics.plot_pacf(y, lags=lags, ax=pacf_ax, alpha=0.5)sm.qqplot(y, line='s', ax=qq_ax)qq_ax.set_title('QQ Plot')        scs.probplot(y, sparams=(y.mean(), y.std()), plot=pp_ax)plt.tight_layout()return tsplot(df_data_ts['send_weight'], lags=30)

2.2 信息准则确定

就是遍历所有参数值，看信息准则最小的（AIC/BIC/HQIC）,一般选取AIC即可。
信息准则的区别可参考文章

best_aic = np.inf 
best_order = None
best_mdl = None
#一般阶数不超过整体数据的十分之一，因此这里采用穷举法。分别从0~10取p、q。
pq_rng = range(10) # [0,1,2,3,4]
d_rng = range(3) # [0,1] ARMA
for i in pq_rng:for d in d_rng:for j in pq_rng:try:tmp_mdl = smt.ARIMA(ads, order=(i,d,j)).fit(method='mle', trend='nc')tmp_aic = tmp_mdl.aicif tmp_aic < best_aic:best_aic = tmp_aicbest_order = (i, d, j)best_mdl = tmp_mdlexcept: continueprint('aic: {:6.5f} | order: {}'.format(best_aic, best_order))
_ = tsplot(best_mdl.resid, lags=8)
print(best_mdl.summary())## 德宾-沃森（Durbin-Watson）检验
print(sm.stats.durbin_watson(best_mdl.resid.values))

【注意】 代码中专门加了德宾-沃森（Durbin-Watson）检验。对于一阶差分，当D-W检验值接近于2时，不存在自相关性，说明模型较好。

2.3 热力图确定

其实热力图定阶的方式和信息准则定阶的方式类似，只是用热力图的方式呈现了。

#设置遍历循环的初始条件，以热力图的形式展示，跟AIC定阶作用一样p_min = 0q_min = 0p_max = 5q_max = 5d_min = 0d_max = 5# 创建Dataframe,以BIC准则results_aic = pd.DataFrame(index=['AR{}'.format(i) \for i in range(p_min,p_max+1)],\columns=['MA{}'.format(i) for i in range(q_min,q_max+1)])# itertools.product 返回p,q中的元素的笛卡尔积的元组for p,d,q in itertools.product(range(p_min,p_max+1),\range(d_min,d_max+1),range(q_min,q_max+1)):if p==0 and q==0:results_aic.loc['AR{}'.format(p), 'MA{}'.format(q)] = np.nancontinuetry:model = sm.tsa.ARIMA(ads, order=(p, d, q))results = model.fit()#返回不同pq下的model的BIC值results_aic.loc['AR{}'.format(p), 'MA{}'.format(q)] = results.aicexcept:continueresults_aic = results_aic[results_aic.columns].astype(float)#print(results_bic)fig, ax = plt.subplots(figsize=(10, 8))ax = sns.heatmap(results_aic,#mask=results_aic.isnull(),ax=ax,annot=True, #将数字显示在热力图上fmt='.2f',)ax.set_title('AIC')plt.show()

结果如下：

热力图黑色的位置最好，一般情况下是越小越好。

3、模型预测

模型预测，有2个函数predict()、forecast()，2个函数的各自参数不做介绍，自己查看帮助文档。区别除了模型参与以外，有以下2个重要区别：

1. predict()可以预测样本内和样本外，forecast只能预测样本外，而且第一个数必须是样本最后一个。predict()函数进行预测，其可以传入start、end参数代表预测数据的开始和结束坐标，坐标值不仅可以是int、str，还可以是 datetime时间类型，如果start、end介于原有数据区域内，即为对原数据的预测拟合，如果end超过了原有数据长度，即代表对未来数据进行预测。forecast()函数接收一个steps参数，代表对未来多少个数据进行预测，也可以传入datetime时间，代表从样本数据结束一直预测到该时间。
2. 如果模型参数d=0不需要差分，则2个函数预测结果一样。若d>0需要差分时,则predict()返回的是差分值，而forecast()则返回的是原始值，从这一点上讲建议使用forecast()。

参考代码分别如下：

best_aic = np.inf 
best_order = None
best_mdl = Nonepq_rng = range(10) # [0,1,2,3,4]
d_rng = range(2) # [0,1] ARMA
for i in pq_rng:for d in d_rng:for j in pq_rng:try:tmp_mdl = smt.ARIMA(ads, order=(i,d,j)).fit(method='mle', trend='nc')tmp_aic = tmp_mdl.aicif tmp_aic < best_aic:best_aic = tmp_aicbest_order = (i, d, j)best_mdl = tmp_mdlexcept: continueprint('aic: {:6.5f} | order: {}'.format(best_aic, best_order))
_ = tsplot(best_mdl.resid, lags=8)
print(best_mdl.summary())
## 德宾-沃森（Durbin-Watson）检验print(sm.stats.durbin_watson(best_mdl.resid.values))###最好用forecast，不用predict，predict在需要查分时数据结果需要还原,2种方法结果不同
n_steps = 14
f, err95, ci95 = best_mdl.forecast(steps=n_steps) # 95% CI
_, err99, ci99 = best_mdl.forecast(steps=n_steps, alpha=0.01) # 99% CI
##预测的第一个数是样本最后的日期
idx = pd.date_range(ads.index[-1], periods=n_steps, freq='D')
fc_95 = pd.DataFrame(np.column_stack([f, ci95]), index=idx, columns=['forecast', 'lower_ci_95', 'upper_ci_95'])
fc_99 = pd.DataFrame(np.column_stack([ci99]), index=idx, columns=['lower_ci_99', 'upper_ci_99'])
fc_all = fc_95.combine_first(fc_99)
fc_all.head()############predict
pred = best_mdl.predict(1,66)#不加开始和结束数据则为原始时间索引
########差分还原------
pred= best_mdl.predict()
pred= pred.add(ads.shift(1))
#这种方法简单但不能还原第一个数

3、 其他模型

还有增加季节因素的SARIMA模型，以及ARCH、GARCH模型，可参考文章。

参考文献

1、http://www.blackarbs.com/blog/time-series-analysis-in-python-linear-models-to-garch/11/1/2016
2、https://www.cnblogs.com/foley/p/5582358.html
3、https://blog.csdn.net/theVicTory/article/details/104941483
4、https://blog.csdn.net/Earl211/article/details/50957029
5、https://blog.csdn.net/jh1137921986/article/details/90257764
6、https://blog.csdn.net/foneone/article/details/90141213?utm_medium=distribute.pc_aggpage_search_result.none-task-blog-2_allfirst_rank_v2~rank_v25-2-90141213.nonecase&utm_term=arima%E4%BA%8C%E9%98%B6%E5%B7%AE%E5%88%86%E7%9A%84%E9%A2%84%E6%B5%8B%E5%80%BC%E8%BF%98%E5%8E%9F#%EF%BC%883%EF%BC%89ADF%E6%A3%80%E9%AA%8C