ARIMA进行时间序列预测-python实现

用ARIMA进行时间序列预测

本文翻译于Kaggle，原文链接时间序列预测教程。中文论坛很少有对整个过程进行描述，所以想转载过来供大家一起学习。数据在原文也有，我也放了云盘天气数据。
英文水平和学术水平都比较低，所以翻译问题和理论问题在所难免，如果不能理解，请查看原文。

• 我们将使用最常见的方法ARIMA
• ARIMA：差分整合移动平均自回归模型。我将在下一部分详细解释。
• 接下来我们会看到：
  • 什么是时间序列（Time Series）？
  • 时间序列的平稳性（Stationarity )
  • 使一个时间序列平稳？
  • 预测一个时间序列

什么是时间序列？

• 时间序列是以国定的时间间隔收集的数据点的集合。
• 它依赖于时间。
• 大多数的时间序列都有某种形式的季节性趋势。例如，如果我们销售冰淇淋，很可能在夏季会有更高的销售额。因此，该事件序列具有季节性趋势。
• 再举一个例子，假如我们在一年内每天投掷一次骰子。正如你想的那样，不会出现数字6主要出现在夏季或者数字5出现在1月份这样的情景。因此，该时间序列不具有季节性趋势。

时间序列的平稳性

• 一个时间序列是平稳序列还是非平稳序列，有三个基本的判断标准。
  • 时间序列的统计特征比如均值、方差应随时间保持恒定，那么我们说这个时间序列是平稳的。
    • 常数平均值（constant mean）
    • 常数方差（constant variance）
    • 不依赖于时间的自协方差。自协方差是指时间序列与滞后时间序列之间的协方差。

# Mean temperature of Bindikuri area
plt.figure(figsize=(22,10))
plt.plot(weather_bin.Date,weather_bin.MeanTemp)
plt.title("Mean Temperature of Bindukuri Area")
plt.xlabel("Date")
plt.ylabel("Mean Temperature")
plt.show()

# lets create time series from weather 
timeSeries = weather_bin.loc[:, ["Date","MeanTemp"]]
timeSeries.index = timeSeries.Date
ts = timeSeries.drop("Date",axis=1)

• 从上面的图中可以看出，我们的时间序列具有季节性变化。每年夏季平均气温较高，冬季平均气温较低。
• 现在我们来检验一下时间序列的平稳性。我们可以用以下方法检验平稳性:
  • 绘制滚动数据：我们有一个窗口假设窗口大小为6然后通过滚动均值和方差来检验是否平稳。
  • 迪基-福勒检验（Dickey-Fuller Test）：测试结果包括一个测试统计量和一些差异置信水平的临界值。如果检验统计量小于临界值，我们可以说时间序列是平稳的。

# adfuller library 
from statsmodels.tsa.stattools import adfuller
# check_adfuller
def check_adfuller(ts):
    # Dickey-Fuller test
    result = adfuller(ts, autolag='AIC')
    print('Test statistic: ' , result[0])
    print('p-value: '  ,result[1])
    print('Critical Values:' ,result[4])
# check_mean_std
def check_mean_std(ts):
    #Rolling statistics
    rolmean = pd.rolling_mean(ts, window=6)
    rolstd = pd.rolling_std(ts, window=6)
    plt.figure(figsize=(22,10))   
    orig = plt.plot(ts, color='red',label='Original')
    mean = plt.plot(rolmean, color='black', label='Rolling Mean')
    std = plt.plot(rolstd, color='green', label = 'Rolling Std')
    plt.xlabel("Date")
    plt.ylabel("Mean Temperature")
    plt.title('Rolling Mean & Standard Deviation')
    plt.legend()
    plt.show()
    
# check stationary: mean, variance(std)and adfuller test
check_mean_std(ts)
check_adfuller(ts.MeanTemp)

Test statistic:  -1.409596674588769
p-value:  0.5776668028526388
Critical Values: {'1%': -3.439229783394421, '5%': -2.86545894814762, '10%': -2.5688568756191392}

• 平稳的第一个标准是常数均值。因为平均值不是恒定的所以它不是平稳序列。你可以从上面的图（黑线）中看到。
• 第二个是常数方差。它看起来像是一个常数。
• 第三，如果检验统计量小于临界值，我们可以说时间序列是平稳的。我们来看一看：
  • test statistic = -1.4 and critical values = {‘1%’: -3.439229783394421, ‘5%’: -2.86545894814762, ‘10%’: -2.5688568756191392}.检验统计值大于临界值。（no stationary)
• 因此，我们确信我们的时间序列不是平稳的。
• 接下来让我们把时间序列变得平稳。

使一个时间序列平稳？

• 如前文所述，时间序列的非平稳性有两个原因：
  • 趋势（Trend）:随时间变化的平均值。对于平稳的时间序列，我们需要常数均值。
  • 季节性（Seasonality）:在特定时间发生变化。我们需要常数变化来表示时间序列的平稳。
• 首先解决趋势（常数均值）的问题
  • 最常用的方法是移动平均线。
    • 移动平均（Moving average）:我们有一个窗口，取过去n个样本的平均值。'n’是窗口大小。

# Moving average method
window_size = 6
moving_avg = pd.rolling_mean(ts,window_size)
plt.figure(figsize=(22,10))
plt.plot(ts, color = "red",label = "Original")
plt.plot(moving_avg, color='black', label = "moving_avg_mean")
plt.title("Mean Temperature of Bindukuri Area")
plt.xlabel("Date")
plt.ylabel("Mean Temperature")
plt.legend()
plt.show()

ts_moving_avg_diff = ts - moving_avg
ts_moving_avg_diff.dropna(inplace=True) # first 6 is nan value due to window size

# check stationary: mean, variance(std)and adfuller test
check_mean_std(ts_moving_avg_diff)
check_adfuller(ts_moving_avg_diff.MeanTemp)

Test statistic:  -11.138514335138474
p-value:  3.150868563164652e-20
Critical Values: {'1%': -3.4392539652094154, '5%': -2.86546960465041, '10%': -2.5688625527782327}

• 常数平均标准:平均值看起来像常数，你可以从上面的图(黑线)中看到。
• 第二个是常数方差。它看起来像个常数。
• 检验统计量小于1%的临界值所以我们可以说有99%的信心说这是一个平稳序列。
• 这样我们可以得到一个平稳时间序列。然后让我们一下另一种避免趋势和季节性的方法。
  • 差分法（Diffrencing method）:是最常用的方法之一。方法是在时间序列和移位（shifted）的时间序列之间取差。

# differencing method
ts_diff = ts - ts.shift()
plt.figure(figsize=(22,10))
plt.plot(ts_diff)
plt.title("Differencing method") 
plt.xlabel("Date")
plt.ylabel("Differencing Mean Temperature")
plt.show()

ts_diff.dropna(inplace=True) # due to shifting there is nan values
# check stationary: mean, variance(std)and adfuller test
check_mean_std(ts_diff)
check_adfuller(ts_diff.MeanTemp)

Test statistic:  -11.678955575105382
p-value:  1.7602075693558453e-21
Critical Values: {'1%': -3.439229783394421, '5%': -2.86545894814762, '10%': -2.5688568756191392}

• 常数平均标准:平均值看起来像常数，你可以从上面的图(黑线)中看到。
• 第二个是常数方差。它看起来像个常数。
• 检验统计量小于1%的临界值所以我们可以说有99%的信心说这是一个平稳序列。

预测一个时间序列

• 我们学习了两种不同的方法，即移动平均和差分法来避免趋势和季节性问题。
• 对于预测(prediction、forecasting），我们将使用ts_diff时间序列，它是差分法的结果。
• 预测方法为ARIMA。
  • AR：auto-Regressive（p）：AR项是因变量的滞后。举个例子p = 3，我们将用x(t-1),x(t-2),x(t-3)来预测x(t)
  • I:Intergraed(d): 这些事非季节性差异的数目。举个例子，在这里我们取一阶差分。所以我们传递这个变量并让d = 0
  • MA：Moving Averages(q)：MA项是预测方程中的滞后预测误差。
• (p,d,q)是ARIMA模型的三个参数。
• 为了确定p,d,q三个参数的值，我们将使用两个不痛的图。
  • 自相关函数 (ACF):时间序列与滞后时间序列之间相关性的度量。
  • 部分自相关函数(PACF):

# ACF and PACF 
from statsmodels.tsa.stattools import acf, pacf
lag_acf = acf(ts_diff, nlags=20)
lag_pacf = pacf(ts_diff, nlags=20, method='ols')
# ACF
plt.figure(figsize=(22,10))

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

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

• 两条虚线是置信区间。我们使用这些线来确定p和q的值
  • p: PACF图第一次跨越上置信区间时的滞后值。p = 1。
  • q: ACF图第一次跨越上置信区间时的滞后值。q = 1。
• 现在使用（1，0，1）作为ARIMA模型的三个参数并进行预测
  • ARIMA：来自statsmodels库
  • datetime: 我们将使用它作为预测方法的起始和结束指标

# ARIMA LİBRARY
from statsmodels.tsa.arima_model import ARIMA
from pandas import datetime

# fit model
model = ARIMA(ts, order=(1,0,1)) # (ARMA) = (1,0,1)
model_fit = model.fit(disp=0)

# predict
start_index = datetime(1944, 6, 25)
end_index = datetime(1945, 5, 31)
forecast = model_fit.predict(start=start_index, end=end_index)

# visualization
plt.figure(figsize=(22,10))
plt.plot(weather_bin.Date,weather_bin.MeanTemp,label = "original")
plt.plot(forecast,label = "predicted")
plt.title("Time Series Forecast")
plt.xlabel("Date")
plt.ylabel("Mean Temperature")
plt.legend()
plt.show()

# predict all path
from sklearn.metrics import mean_squared_error
# fit model
model2 = ARIMA(ts, order=(1,0,1)) # (ARMA) = (1,0,1)
model_fit2 = model2.fit(disp=0)
forecast2 = model_fit2.predict()
error = mean_squared_error(ts, forecast2)
print("error: " ,error)
# visualization
plt.figure(figsize=(22,10))
plt.plot(weather_bin.Date,weather_bin.MeanTemp,label = "original")
plt.plot(forecast2,label = "predicted")
plt.title("Time Series Forecast")
plt.xlabel("Date")
plt.ylabel("Mean Temperature")
plt.legend()
plt.savefig('graph.png')

plt.show()

error:  1.8625819952314109

结论

• 我们学习如何用pyplot绘制地图。
• 我们学习如何进行时间序列预测。