用pandas处理时序数据
导入时序数据，设置时间为索引
# 用pandas导入和处理时序数据
import numpy as np
import pandas as pd
from datetime import datetime
import matplotlib.pyplot as plt
from statsmodels.tsa.stattools import adfuller
# 读取数据，pd.read_csv默认生成DataFrame对象，需将其转换成Series对象
data = pd.read_csv('./arima_data_pems04.csv', encoding='utf-8', index_col='Date')
df = data[0:288]
df.index = pd.to_datetime(df.index) # 将字符串索引转换成时间索引
ts = df['flow'] # 生成pd.Series对象
ts.head()
检验序数据的稳定性
平稳性
平稳性就是要求经由样本时间序列所得到的拟合曲线在未来一段时间内仍能顺着现有的形态惯性地延续下去。平稳性要求序列的均值和方差不发生明显变化。
检验方法
平稳性检验一般采用观察法和单位根检验法。
观察法
需计算每个时间段内的平均的数据均值和标准差。
查看原始数据的均值和方差：
# 移动平均图
def draw_trend(timeseries, size):
f = plt.figure(facecolor='white',figsize=(8,4))
# 对size个数据进行移动平均
rol_mean = timeseries.rolling(window=size).mean()
# 对size个数据移动平均的方差
rol_std = timeseries.rolling(window=size).std()
timeseries.plot(color='blue', label='Original')
rol_mean.plot(color='red', label='Rolling Mean')
rol_std.plot(color='black', label='Rolling standard deviation')
plt.legend(loc='best')
plt.title('Rolling Mean & Standard Deviation')
plt.show()
def draw_ts(timeseries):
f = plt.figure(facecolor='white')
timeseries.plot(color='blue')
plt.show()
#Dickey-Fuller test:
def teststationarity(ts):
dftest = adfuller(ts)
# 对上述函数求得的值进行语义描述
dfoutput = pd.Series(dftest[0:4], index=['Test Statistic','p-value','#Lags Used','Number of Observations Used'])
for key,value in dftest[4].items():
dfoutput['Critical Value (%s)'%key] = value
return dfoutput
# 查看原始数据的均值和方差
draw_trend(ts,12)
通过上图，我们可以发现数据的移动平均值/标准差有越来越大的趋势，是不稳定的。接下来我们再看Dickey-Fuller的结果
单位根检验法
通过Dickey-Fuller Test 进行判断，大致意思就是在一定置信水平下，对于时序数据假设 Null hypothesis: 非稳定。这是一种常用的单位根检验方法，它的原假设为序列具有单位根，即非平稳，对于一个平稳的时序数据，就需要在给定的置信水平上显著，拒绝原假设。
# 接下来我们再看Dickey-Fuller的结果
teststationarity(ts)
# 此时p值为0.794490，说明并不能拒绝原假设。
此时p值为0.794490，说明并不能拒绝原假设。通过DF的数据可以明确的看出，在任何置信度下，数据都不是稳定的。
处理时序数据变成稳定数据
数据不稳定的原因主要有以下两点：
趋势（trend）：数据随着时间变化。比如说升高或者降低。
季节性（seasonality）：数据在特定的时间段内变动。比如说节假日，或者活动导致数据的异常。
由前面的分析可知，该序列是不平稳的，然而平稳性是时间序列分析的前提条件，故我们需要对不平稳的序列进行处理将其转换成平稳的序列。
对数变换
对数变换主要是为了减小数据的振动幅度，使其线性规律更加明显，同时保留其他信息。这里强调一下，变换的序列需要满足大于0，小于0的数据不存在对数变换。
ts_log = np.log(ts)
绘制ts_log图像可以看出经过对数变换后，数据值域范围缩小了，振幅也没那么大了。
平滑法
根据平滑技术的不同，平滑法具体分为移动平均法和指数平均法。
移动平均法
指数平均法
移动平均即利用一定时间间隔内的平均值作为某一期的估计值，而指数平均则是用变权的方法来计算均值。
def draw_moving(timeSeries, size):
f = plt.figure(facecolor='white')
# 对size个数据进行移动平均
rol_mean = timeSeries.rolling(window=size).mean()
# 对size个数据进行加权移动平均
rol_weighted_mean = pd.DataFrame.ewm(timeSeries, span=size).mean()
#rol_weighted_mean=timeSeries.ewm(halflife=size,min_periods=0,adjust=True,ignore_na=False).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()
draw_moving(ts_log,12)
从上图可以发现窗口为12的移动平均能较好的剔除周期性因素，而指数平均法是对周期内的数据进行了加权，能在一定程度上减小周期因素，但并不能完全剔除，如要完全剔除可以进一步进行差分操作。
差分
时间序列最常用来剔除周期性因素的方法当属差分了，它主要是对等周期间隔的数据进行线性求减。前面我们说过，ARIMA模型相对ARMA模型，仅多了差分操作，ARIMA模型几乎是所有时间序列软件都支持的，差分的实现与还原都非常方便。
diff_12 = ts_log.diff(12)
diff_12.dropna(inplace=True)
diff_12_1 = diff_12.diff(1)
diff_12_1.dropna(inplace=True)
teststationarity(diff_12_1)
结果如下：
从上面的统计检验结果可以看出，经过12阶差分和1阶差分后，该序列满足平稳性的要求了。
分解
所谓分解就是将时序数据分离成不同的成分。statsmodels使用的X-11分解过程，它主要将时序数据分离成长期趋势、季节趋势和随机成分。与其它统计软件一样，statsmodels也支持两类分解模型，加法模型和乘法模型，这里只实现加法，乘法只需将model的参数设置为"multiplicative"即可。
from statsmodels.tsa.seasonal import seasonal_decompose
def decompose(timeseries):
# 返回包含三个部分 trend（趋势部分） ， seasonal（季节性部分） 和residual (残留部分)
decomposition = seasonal_decompose(timeseries,model='addictive',period=12)
trend = decomposition.trend
seasonal = decomposition.seasonal
residual = decomposition.resid
plt.subplot(411)
plt.plot(ts_log, label='Original')
plt.legend(loc='best')
plt.subplot(412)
plt.plot(trend, label='Trend')
plt.legend(loc='best')
plt.subplot(413)
plt.plot(seasonal,label='Seasonality')
plt.legend(loc='best')
plt.subplot(414)
plt.plot(residual, label='Residuals')
plt.legend(loc='best')
plt.tight_layout()
plt.show()
return trend , seasonal, residual
trend , seasonal, residual = decompose(ts_log)
residual.dropna(inplace=True)
draw_trend(residual,12)
teststationarity(residual)
如图所示，数据的均值和方差趋于常数，几乎无波动(看上去比之前的陡峭，但是要注意他的值域只有[-0.05,0.05]之间)，所以直观上可以认为是稳定的数据。另外DFtest的结果显示，Statistic值原小于1%时的Critical value，所以在99%的置信度下，数据是稳定的。
时序数据的预测
在前面的分析可知，该序列具有明显的周期与长期成分。对于周期成分我们使用窗口为12的移动平均进行处理，对于长期趋势成分我们采用1阶差分来进行处理。
rol_mean = ts_log.rolling(window=12).mean()
rol_mean.dropna(inplace=True)
ts_diff_1 = rol_mean.diff(1)
ts_diff_1.dropna(inplace=True)
teststationarity(ts_diff_1)
观察其统计量发现该序列在置信水平为95%的区间下并不显著，我们对其进行再次一阶差分。再次差分后的序列其自相关具有快速衰减的特点，t统计量在99%的置信水平下是显著的，这里不再做详细说明。
ts_diff_2 = ts_diff_1.diff(1)
ts_diff_2.dropna(inplace=True)
teststationarity(ts_diff_2)
数据平稳后，需要对模型定阶，即确定p、q的阶数。先画出ACF,PACF的图像,代码如下：
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
def draw_acf_pacf(ts,lags):
f = plt.figure(facecolor='white')
ax1 = f.add_subplot(211)
plot_acf(ts,ax=ax1,lags=lags)
ax2 = f.add_subplot(212)
plot_pacf(ts,ax=ax2,lags=lags)
plt.subplots_adjust(hspace=0.5)
plt.show()
draw_acf_pacf(ts_diff_2,30)
观察上图，发现自相关和偏相系数都存在拖尾的特点，并且他们都具有明显的一阶相关性，所以我们设定p=1, q=1。下面就可以使用ARMA模型进行数据拟合了。（Ps.PACF是判定AR模型阶数的，也就是p。ACF是判断MA阶数的，也就是q）
from statsmodels.tsa.arima.model import ARIMA
model = ARIMA(ts_diff_1, order=(1,1,1))
result_arima = model.fit()
模型拟合完后，我们就可以对其进行预测了。由于ARMA拟合的是经过相关预处理后的数据，故其预测值需要通过相关逆变换进行还原。
predict_ts = result_arima.predict()
# 一阶差分还原
diff_shift_ts = ts_diff_1.shift(1)
diff_recover_1 = predict_ts.add(diff_shift_ts)
# 再次一阶差分还原
rol_shift_ts = rol_mean.shift(1)
diff_recover = diff_recover_1.add(rol_shift_ts)
# 移动平均还原
rol_sum = ts_log.rolling(window=11).sum()
rol_recover = diff_recover*12 - rol_sum.shift(1)
# 对数还原
log_recover = np.exp(rol_recover)
log_recover.dropna(inplace=True)
评估模型样本内拟合的好坏：
ts = ts[log_recover.index] # 过滤没有预测的记录plt.figure(facecolor='white')
import math
from sklearn.metrics import mean_squared_error as mse, mean_absolute_error as mae
def mape(y_true, y_pred):
y_true, y_pred = np.array(y_true), np.array(y_pred)
non_zero_index = (y_true > 0)
y_true = y_true[non_zero_index]
y_pred = y_pred[non_zero_index]
mape = np.abs((y_true - y_pred) / y_true)
mape[np.isinf(mape)] = 0
return np.mean(mape) * 100
rmse_score = math.sqrt(mse(log_recover,ts))
mae_score = mae(log_recover,ts)
mape_score = mape(log_recover,ts)
#计算指定置信水平下的预测准确率
#level为小数
def comput_acc(real,predict,level):
num_error=0
for i in range(len(real)):
if abs(real[i]-predict[i])/real[i]>level:
num_error+=1
return 1-num_error/len(real)
acc = comput_acc(ts,log_recover,0.2)
log_recover.plot(color='blue', label='Predict')
ts.plot(color='red', label='Original')
plt.legend(loc='best')
plt.title(f'rmse:{rmse_score:.4f} mae:{mae_score:.4f} mape:{mape_score:.4f} acc:{acc:.4f}')
plt.show()
