面向初学者的指南：创建时间序列预测 (使用Python）

import pandas as pd
import numpy as np
import matplotlib.pylab as plt
%matplotlib inline
from matplotlib.pylab import rcParams
rcParams['figure.figsize'] = 15, 6

data = pd.read_csv('AirPassengers.csv')
print(data.head(n = 10))
# 默认 n为5)
print('\n Data Types:')
print(data.dtypes)

#数据类型是“object”和“int”，所以它仍然不是作为TS对象读取的。
#为了将数据读取为时间序列，我们必须将特殊参数传递给read_csv命令:
dateparse = lambda dates: pd.datetime.strptime(dates, '%Y-%m')
data = pd.read_csv('AirPassengers.csv', parse_dates=['Month'], index_col='Month',date_parser=dateparse)
print(data.head())

#parse_dates:指定包含日期-时间信息的列。正如我们上面所说，列名是’ Month '。
#index_col:将panda用于TS数据背后的一个关键思想是，索引必须是描述日期-时间信息的变量。所以这个参数告诉熊猫使用“Month”列作为索引。
#date_parser:指定一个函数，该函数将输入字符串转换为datetime变量。默认情况下，熊猫读取的数据格式为“YYYY-MM-DD HH:MM:SS”。如果数据不是这种格式，则必须手动定义格式。类似于这里定义的dataparse函数可以用于此目的。

#可以看到数据的索引是time对象，列是#passenger。我们可以用以下命令交叉检查索引的数据类型
data.index

#dtype= 'datetime[ns] ’ 确认它是一个datetime对象。作为个人偏好，我将把列转换成一个Series对象，以防止每次使用TS时引用列名。
ts = data['#Passengers'] 
ts.head(10)

#1. Specific the index as a string constant:
ts['1949-01-01']

#2. Import the datetime library and use 'datetime' function:
from datetime import datetime
ts[datetime(1949,1,1)]

#假设我们想要到1949年5月为止的所有数据。这可以通过两种方式实现:
#1. Specify the entire range:
ts['1949-01-01':'1949-05-01']

#2. Use ':' if one of the indices is at ends:
ts[:'1949-05-01']

#与数字索引不同，这里 包含了结束索引。例如，如果我们将列表索引为[:5]，那么它将返回索引处的值-[0,1,2,3,4]。但是这里的输出中包含索引 ‘1949-05-01’。
#索引必须按照范围排序。如果随机打乱索引，这将不起作用

#1949年所有值的例子。可以这样做:

ts['1949']

#简单地绘制数据图并进行可视化分析。数据可以使用以下命令绘制:
plt.plot(ts)

#回到检查平稳性，我们将经常使用滚动统计图和Dickey-Fuller测试结果，所以我定义了一个函数，它接受TS作为输入并为我们生成它们。请注意，我画的是标准差而不是方差来保持单位接近均值。

from statsmodels.tsa.stattools import adfuller
def test_stationarity(timeseries):
    
    #Determing rolling statistics
    rolmean = timeseries.rolling(12).mean()
    rolstd = timeseries.rolling(12).std()

    #Plot rolling statistics:
    orig = plt.plot(timeseries, color='blue',label='Original')
    mean = plt.plot(rolmean, color='red', label='Rolling Mean')
    std = plt.plot(rolstd, color='black', label = 'Rolling Std')
    plt.legend(loc='best')
    plt.title('Rolling Mean & Standard Deviation')
    plt.show(block=False)
    
    #Perform Dickey-Fuller test:
    print('Results of Dickey-Fuller Test:')
    dftest = adfuller(timeseries, autolag='AIC')
    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
    print(dfoutput)
test_stationarity(ts)

#估计和消除趋势

#减少趋势的首要技巧之一是转变。
#例如，在这种情况下，我们可以清楚地看到，有一个显著的积极趋势。
#所以我们可以用变换来惩罚更大的值而不是更小的值。可以取对数，平方根，立方根，等等。为了简单起见，我们在这里做一个对数变换:
ts_log = np.log(ts)
plt.plot(ts_log)

#移动平均数

#在这种方法中，我们根据时间序列的频率取“k”连续值的平均值。这里我们可以取过去一年的平均值，也就是最近12个值。panda定义了用于确定滚动统计数据的特定函数。

moving_avg = ts_log.rolling(12).mean()
plt.plot(ts_log)
plt.plot(moving_avg, color='red')

#红线表示滚动平均值。从原级数中减去这个。注意，由于我们取最后12个值的平均值，所以前11个值没有定义滚动平均值。这可以观察到:

ts_log_moving_avg_diff = ts_log - moving_avg
ts_log_moving_avg_diff.head(12)

#注意前11位是Nan。让我们删除这些NaN值，并检查这些图以测试平稳性

ts_log_moving_avg_diff.dropna(inplace=True)
test_stationarity(ts_log_moving_avg_diff)

#一种流行的方法是 指数加权移动平均法，它将权重分配给所有具有衰减因子的先前值。在这里找到详细信息。这可以在panda中实现为:
expwighted_avg = pd.DataFrame.ewm(ts_log, halflife=12).mean()
plt.plot(ts_log)
plt.plot(expwighted_avg, color='red')
#注意，这里使用参数“半衰期”来定义指数衰减的数量。这只是一个假设，主要取决于业务领域。其他参数，如跨度和质心也可以用来定义衰减，这是讨论在上面共享的链接。
    

#现在，让我们把这个从级数中移除，并检查平稳性:

ts_log_ewma_diff = ts_log - expwighted_avg
test_stationarity(ts_log_ewma_diff)

#这个TS在平均值和标准差上的变化更小。同时，测试统计量 小于1%的临界值，优于前一种情况。注意，在这种情况下，不会缺少值，因为从一开始所有的值都是给定权重的。所以即使没有之前的值，它也能工作。

#这似乎大大降低了趋势。让我们用我们的图来验证:

ts_log_diff.dropna(inplace=True)
test_stationarity(ts_log_diff)
#我们可以看到，平均值和std随时间的变化有很小的变化。此外，Dickey-Fuller检验统计量小于10%的临界值，因此TS是稳定的，有90%的置信度。我们也可以采取二阶或三阶的差异，这可能在某些应用中得到更好的结果。我把它们留给你去尝试。

#分解：
#在这种方法中，趋势和季节性分别建模，并返回系列的其余部分。我将跳过统计数据，来看一下结果:

from statsmodels.tsa.seasonal import seasonal_decompose
decomposition = seasonal_decompose(ts_log)

#趋势
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()

#这里我们可以看到趋势，季节性从数据中分离出来我们可以对残差进行建模。检验残差的平稳性:

ts_log_decompose = residual
ts_log_decompose.dropna(inplace=True)
test_stationarity(ts_log_decompose)
#Dickey-Fuller检验统计量显著低于1%的临界值(p值<0.05)。所以这个TS非常接近于静止。您还可以尝试高级分解技术，这些技术可以生成更好的结果。另外，您应该注意，在这种情况下，将残差转换为未来数据的原始值并不十分直观。

#经差分后的TS的ACF和PACF图可以绘制为:

#ACF and PACF plots:
from statsmodels.tsa.stattools import acf, pacf

lag_acf = acf(ts_log_diff, nlags=20)
lag_pacf = pacf(ts_log_diff, nlags=20, method='ols')

#Plot ACF: 
plt.subplot(121) 
plt.plot(lag_acf)
plt.axhline(y=0,linestyle='--',color='gray')
plt.axhline(y=-1.96/np.sqrt(len(ts_log_diff)),linestyle='--',color='gray')
plt.axhline(y=1.96/np.sqrt(len(ts_log_diff)),linestyle='--',color='gray')
plt.title('Autocorrelation Function')

#Plot PACF:
plt.subplot(122)
plt.plot(lag_pacf)
plt.axhline(y=0,linestyle='--',color='gray')
plt.axhline(y=-1.96/np.sqrt(len(ts_log_diff)),linestyle='--',color='gray')
plt.axhline(y=1.96/np.sqrt(len(ts_log_diff)),linestyle='--',color='gray')
plt.title('Partial Autocorrelation Function')
plt.tight_layout()

#在这个图中，0两边的虚线是置信区间。这些可以用来确定“p”和“q”的值如下:
#p - PACF图第一次越过上置信区间时的滞后值。如果你仔细观察，这里p=2。
#q - ACF图第一次越过上置信区间时的滞后值。如果你仔细观察，这里q=2。

from statsmodels.tsa.arima_model import ARIMA

#可以使用ARIMA的order参数指定p、d、q值，该参数接受一个元组(p、d、q)。让我们对这3种情况建模:

#AR MODEL

model = ARIMA(ts_log, order=(2, 1, 0))  
results_AR = model.fit(disp=-1)  
plt.plot(ts_log_diff)
plt.plot(results_AR.fittedvalues, color='red')
plt.title('RSS: %.4f'% sum((results_AR.fittedvalues-ts_log_diff)**2))

#MA Model

model = ARIMA(ts_log, order=(0, 1, 2))  
results_MA = model.fit(disp=-1)  
plt.plot(ts_log_diff)
plt.plot(results_MA.fittedvalues, color='red')
plt.title('RSS: %.4f'% sum((results_MA.fittedvalues-ts_log_diff)**2))

#Combined Model

model = ARIMA(ts_log, order=(2, 1, 2))  
results_ARIMA = model.fit(disp=-1)  
plt.plot(ts_log_diff)
plt.plot(results_ARIMA.fittedvalues, color='red')
plt.title('RSS: %.4f'% sum((results_ARIMA.fittedvalues-ts_log_diff)**2))

#回到原来的比例

#由于组合模型给出了最佳结果，因此我们可以将其缩放回原始值并查看其在那里的表现。第一步是将预测结果存储为单独的系列并观察它。

predictions_ARIMA_diff = pd.Series(results_ARIMA.fittedvalues, copy=True)
print(predictions_ARIMA_diff.head())

predictions_ARIMA_diff_cumsum = predictions_ARIMA_diff.cumsum()
print(predictions_ARIMA_diff_cumsum.head())

#您可以使用以前的输出快速进行一些回溯计算，以检查这些计算是否正确。接下来我们要把它们加到底数上。为此，让我们创建一个以所有值为基数的系列，并将差异添加到其中。可以这样做:

predictions_ARIMA_log = pd.Series(ts_log.ix[0], index=ts_log.index)
predictions_ARIMA_log = predictions_ARIMA_log.add(predictions_ARIMA_diff_cumsum,fill_value=0)
predictions_ARIMA_log.head()

#这里第一个元素是基数本身，并且从那里累加的值。最后一步是取指数并与原始系列进行比较。

predictions_ARIMA = np.exp(predictions_ARIMA_log)
plt.plot(ts)
plt.plot(predictions_ARIMA)
plt.title('RMSE: %.4f'% np.sqrt(sum((predictions_ARIMA-ts)**2)/len(ts)))