本篇文章的数据来源是网络上的航班日期和对应的乘客人数的数据集AirPassengers.csv。
初步的探索性分析
我们先看看数据集的大致情况。
数据集的大致情况
可以看得出来，数据集只有两列，第一列是日期，后面要转化成pandas的TimeStamp格式，第二列是对应的乘客数。
接下来直接读取数据。
import pandas as pd
#定义一个返回日期的函数
dates=lambda date: pd.datetime.strptime(date, '%m/%d/%Y')
data=pd.read_csv('AirPassengers.csv', parse_dates=['TravelDate'], date_parser=dates, index_col='TravelDate', dayfirst=False)
#index_col可以指定哪一列为index
输出结果
输出的结果如上。
接下来，我们用plot来先绘制一下整体的趋势。
import matplotlib.pyplot as plt
plt.plot_date(x=data.index, y='Passengers', data=data, xdate=True)
plt.xlabel('TravelDate')
plt.ylabel('Passengers')
plt.show()
日期与人数的关系
以天计的折线图
从图上看出来，人数和年数似乎成正相关关系。
接下来我们单独分析一下1949年的情况。
plt.plot_date(x=data['1949-01-01':'1950-01-01'].index, y='Passengers', data=data['1949-01-01':'1950-01-01'], xdate=True)
plt.xlabel('TravelDate')
plt.ylabel('Passengers')
plt.show()
又或者我们做另外的处理，直接先导入数据
data=pd.read_csv('AirPassengers.csv')
data.index=pd.to_datetime(data['TravelDate'], format='%m/%d/%Y')
然后把多出的一列舍弃
data=data.drop(['TravelDate'], axis=1)
data.head()
输出结果
plt.plot_date(x=data['1949'].index, y='Passengers', data=data['1949'], xdate=True, fmt='bo-')
plt.xlabel('1949')
plt.ylabel('Passengers')
plt.show()
1949年的变化
重新采样，看看年的情况
data_year=data.resample('1Y').sum()
data_year=data.resample('12M').sum()
plt.plot_date(x=data_year.index, y='Passengers', data=data_year, xdate=True, fmt='bo-')
plt.show()
全部年份汇总以后的变化
上面这个图是把一年的总人数合计之后得出的图，对比一下把全部点绘制出来的图，可以看得到把差异给抹平了，更加凸显出一个上升的趋势。
如何检查一个时序分析的不变性
时序分析依赖于两者：一个是时间另外一个是季节，就好像貂皮大衣在冬天销量要比夏天要高很多一样。
时序分析目标是要在变动性比较大的数据当中（比如以天计的数据变动性就很大，参考第二张图）看有没有一个比较稳定的规律，我们就可以用这个规律来预测将来的趋势。在经过把天转换成年的这段时间间隔后，我们就可以发现一个光滑得多的一个明显趋势，然而这样趋势只能帮助我们分析以年为单位的规律，实际用处并不是很大。
实际上，我们如果要看一组数据随着时间变化是否具有某种不变的统计学特性，要看三个统计学变量：
- 固定的平均数
- 固定的方差（接下来为了作图的单位统一，用标准差代表）
- 固定的自协方差
关于时序分析的进一步内容，可以参考：时序分析的完整介绍
以天计的折线图
回过头来看这张图，确实有部分的随着季节变换的规律，但我们不能只是做个图，通过肉眼来观察就得出结论。我们必须要有更加确切的量化指标，这里有两个指标来评估统计学意义的不变性：
-
计算滚动统计量。
滚动统计量的定义可以参考wikipedia，简单说就是有一个总体，然后我固定样本大小是n，第一个平均值就是前面所有的平均值，接下来就是把样本里面的第一个值除掉，从样本的队尾进入总体的下一个值，又可以计算出这个样本的平均值作为第二个均值，就想像成排队就行了。这样子我们可以计算出随着时间的变动，滚动均值rolling mean和滚动方差rolling variance，然后可以求出的变化是否随着时间的变动而改变。
-
进行Dickey-Fuller Test
类似假设检验，零假设就是我们认为这个时间序列TS并不是有统计学上的不变性。输入一个时间序列TS，就可以得出对应的t value，然后把t value与对应的critcal value相比，看是否落在置信区间内，如果低于了critical value，也就是置信水平，那么我们就有理由拒绝零假设，接受备择假设， 也就是认为这个TS确实是存在着统计学上的不变性，这个不变性可以理解为这个Series具有的规律，可以帮助我们预测，这一点规律是不变的。
关于实现Dickey-Fuller Test的模块官方文档请查看这里。
接下来我们就定义一个函数能够绘制出第一个变量的图以及返回第二个统计学参数的值的函数。
我们先看一下adfuller这个接口传入的内容和返回的参数。
adfuller接口
传入的是单个的一维series，regression表示拟合的方程里面具有的形式（有没有常数项，二次项等等），其他的有待进一步学习。
返回
返回的是一个元组，分别对应这些信息。
我先执行一个简单的adfuller test来看看。
adf=adfuller(data['Passengers'], autolag='AIC')
print adf
#导入进行假设检验的模块，这个方法会返回对应的p-value值还有置信水平对应的值。
from statsmodels.tsa.stattools import adfuller
def test_stationary(timeseries):
rolling_mean=pd.rolling_mean(timeseries, window=12)
rolling_std=pd.rolling_std(timeseries, window=12)
#plot the plot
plt.plot(timeseries, 'b-', label='Origin')
plt.plot(rolling_mean, 'r-',label='Mean')
plt.plot(rolling_std, 'k-', label='standard deviation')
plt.legend(loc='best')
plt.title('rolling mean and rolling standard deviation')
plt.xlabel('Date')
plt.ylabel('Passengers')
plt.show()
#conduct the adfuller test and return the value
adf=adfuller(x=timeseries, regression='ctt', autolag='AIC')
#we packed the result into a pandas series
adf_test_output=pd.Series(adf[:4], index=['Test statistic', 'p-value', 'number of lags used', 'number of observation used'])
for key, value in adf[4].iteritems():
adf_test_output['%s'%key]=value
print (adf_test_output)
接下来我们执行结果看看。
test_stationary(data['Passengers'])
滚动统计学变量
打印结果
从初步结果看来，p-value是比较大的，有19.59%，初步估计没有统计学平稳性。但是从图上看来，rolling mean确实是逐步增大。关于p-value和alpha水平的关系可以参考这篇文章。
平稳性处理
因为图线不平稳，所以要进行平稳性处理。 平稳性处理有多种方法。其中一种是转变。我们试一下用对数转变。
import numpy as np
log_data=np.log(data['Passengers'])
plt.plot(log_data)
plt.title('log-transformed variation')
plt.xlabel('Date')
plt.ylabel('Log Passengers')
plt.show()
经过对数转换的数据
我们跟上面那张图对比，确实可以感觉到上下波动的幅度是减小了，也就是稳定性变好了。但是这个稳定性还是远远不达标。
我们还可以通过绘制一条更加光滑的曲线来标示出这个明显逐步上升的趋势先看看。有三种常用方法：
1.聚合的方法。也就是可以把时间单位从每日转变到每个月，每一年取平均或者是汇总都可以，这一步在以上做过了。
2.取滚动平均值。也就是用滑动窗不断取样本的平均值，参考滚动统计量曲线变光滑。
3.对曲线进行多次方的方程拟合。
我们试一下按照年份进行汇总之后用多次方的方程进行拟合，seaborn有对曲线进行拟合的模块regplot。
import seaborn as sns
#对数据进行重采样，然后统计12个月的总和
data_year=data.resample('12M').sum()
data_year['id']=np.arange(len(data_year))
sns.regplot('id', 'Passengers', data=data_year, x_ci='ci', marker='o', scatter_kws={'s':5}, truncate=True, order=2)
一次方程的情况
可以看得出来，一次方拟合的情况很不错，周围透明部分是置信区间。
二次方程的情况
二次方方程拟合情况也不错，但是在首尾可以看得出来偏离比较大。
滚动的平均值
在这里我们试一下取一下滚动的log对数转变后的平均值。
test_stationary(log_data)
经过对数转换的滚动统计学变量的图
adfuller的输出结果
我们尝试一下求把随着时间增长的值与滚动平均值的差值，然后作图看看。
diff_log_data=log_data-pd.rolling_mean(log_data, window=12)
print (diff_log_data)
输出结果
从输出结果来看，由于是以12为一个采样集，前面11个不够12个，所以没有rolling_mean，所以返回的数值是NaN。先把前面11个数值去除掉。
diff_log_data.dropna(axis=0, inplace=True, how='any')
test_stationary(diff_log_data)
去除平均值影响后
输出结果
消除掉平均值随着时间的增长的影响之后，可以看到上升的趋势被消除掉，只剩下随着季节变化的波动，输出结果的p-value也变小了，说明稳定性的趋势增强了。
接下来对数据进行分解，用seasonal模块分解后返回的是季节波动，总体趋势以及残差。
from statsmodels.tsa.seasonal import seasonal_decompose
season_model=seasonal_decompose(log_data, model='additive')
seasonal=season_model.seasonal
trend=season_model.trend
resid=season_model.resid
#Next we will plot the figure
plt.subplot(411)
plt.plot(log_data, label='log_data')
plt.legend(loc='best')
plt.subplot(412)
plt.plot(seasonal, label='seasonal')
plt.legend(loc='best')
plt.subplot(413)
plt.plot(trend, label='trend')
plt.legend(loc='best')
plt.subplot(414)
plt.plot(resid, label='resid')
plt.legend(loc='best')
plt.show()
分解
可以分别绘制出趋势，季节性变化和残差。在这里就把趋势和季节性波动从数据当中分离开了。
接下来我们看一下残差的变化。
resid.dropna(inplace=True, axis=0)
test_stationary(resid)
残差的滚动统计学数值随着时间的变化
结论
从上面的结果可以看出，用转换和合并的方法可以减小时序数据的变动性，呈现出比较稳定的趋势；可以用分解的办法来将趋势和季节性波动从数据当中分开。
后期优化可以考虑用ARIMA模型来预测整个模型的趋势以及预测值，效果应该会很不错。参考这篇文章。
