其他分享
首页 > 其他分享> > 【时间序列】-航空数据预测

【时间序列】-航空数据预测

作者：互联网

学习预测时序数据，如有侵权，请联系删除。

主要参考：

A comprehensive beginner’s guide to create a Time Series Forecast (with Codes in Python)

 

目的：

 

介绍

时间序列（从现在开始称为TS）被认为是数据科学领域中鲜为人知的技能之一。我为自己设定了解决时间序列问题的基本步骤，并在此与大家分享。这些肯定会帮助您在未来的任何项目中获得一个不错的模型！

在阅读本文之前，我强烈建议您阅读A Complete Tutorial on Time Series Modeling in R并参加free Time Series Forecasting course。它侧重于基本概念，我将专注于使用这些概念来解决端到端问题以及Python中的代码。R中有许多时间序列的资源，但是很少有Python的资源，所以我将在本文中使用Python。

 

我们的旅程将经历以下步骤：

  1. 什么使时间序列特别？
  2. 在Pandas中加载和处理时间序列
  3. 如何检查时间序列的平稳性？
  4. 如何制作时间序列固定？
  5. 预测时间序列

1.什么使时间序列特别？

顾名思义，TS是以恒定时间间隔收集的数据点的集合。对这些进行分析以确定长期趋势，以便预测未来或执行其他形式的分析。但是什么使TS不同于说常规回归问题？有两件事：

  1. 这是时间依赖的。因此，在这种情况下，观察是独立的线性回归模型的基本假设不成立。
  2. 随着趋势的增加或减少，大多数TS具有某种形式的季节性趋势，即特定时间范围的特定变化。例如，如果你看到一件羊毛外套随着时间的推移而销售，那么你将在冬季找到更高的销售额。

由于TS的固有属性，分析它涉及各种步骤。这些将在下面详细讨论。让我们从在Python中加载TS对象开始。我们将使用流行的AirPassengers数据集，可在此处

下载。

请注意，本文的目的是让您熟悉TS一般使用的各种技术。这里考虑的例子仅用于说明，我将重点关注广泛的主题，而不是做出非常准确的预测。

 

2.在Pandas中加载和处理时间序列

Pandas有专门的库来处理TS对象，特别是datatime64 [ns]类，它存储时间信息并允许我们快速执行某些操作。让我们从启动所需的库开始：

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()
print '\n Data Types:'
print data.dtypes

该数据包含特定月份和该月旅行的乘客数量。但是这仍然不是作为TS对象读取的，因为数据类型是'object'和'int'。为了将数据作为时间序列读取，我们必须将特殊参数传递给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()

让我们逐一理解这些论点：

现在我们可以看到数据有时间对象作为索引，＃Passengers作为列。我们可以使用以下命令交叉检查索引的数据类型：

data.index

注意dtype ='datetime [ns]'确认它是一个datetime对象。作为个人偏好，我会将列转换为Series对象，以防止每次使用TS时引用列名。请随意使用作为数据帧，这对您更有效。

ts = data[‘#Passengers’] 
ts.head(10)

在进一步讨论之前，我将讨论TS数据的一些索引技术。让我们首先选择Series对象中的特定值。这可以通过以下两种方式完成：

#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)]

两者都将返回值'112'，这也可以从先前的输出中确认。假设我们想要所有数据到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']

两者都会产生以下输出：

这里有两件事需要注意：

考虑另一个需要1949年所有值的实例。这可以通过以下方式完成：

ts['1949']

月份部分被省略了。同样，如果您在特定月份的所有日期，可以省略日期部分。

现在，让我们转到分析TS。

 

3.如何检查时间序列的平稳性？

如果TS的统计特性（例如均值，方差）随时间保持不变，则称TS是静止的。但为什么这很重要？大多数TS模型的工作假设TS是静止的。直觉上，我们可以认为，如果TS随着时间的推移具有特定的行为，那么将来它很可能会遵循相同的行为。此外，与非平稳序列相比，与静止序列相关的理论更成熟，更容易实现。

使用非常严格的标准来定义平稳性。然而，出于实际目的，如果系列随时间具有恒定的统计特性，即我们可以假设该系列是静止的。下列：

我将跳过细节，因为它在本文中已经非常明确地定义了。让我们进入测试平稳性的方法。首要的是简单绘制数据并进行视觉分析。可以使用以下命令绘制数据：

plt.plot(ts)

很明显，数据总体上呈增长趋势，并伴有一些季节性变化。然而，可能并不总是可以做出这样的视觉推断（我们稍后会看到这样的情况）。因此，更正式地说，我们可以使用以下方法检查平稳性：

回到检查平稳性，我们将使用滚动统计图和Dickey-Fuller测试结果，因此我定义了一个函数，它将TS作为输入并为我们生成它们。请注意，我绘制了标准偏差而不是方差，以使单位与平均值相似。

from statsmodels.tsa.stattools import adfuller
def test_stationarity(timeseries):
    
    #Determing rolling statistics
    #rolmean = pd.rolling_mean(timeseries, window=12)
    #rolstd = pd.rolling_std(timeseries, window=12)
    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)

虽然标准偏差的变化很小，但平均值随着时间的推移而明显增加，这不是一个固定的系列。此外，测试统计数据远远超过临界值。请注意，应比较有符号值，而不是绝对值。

 

4.如何制作时间序列固定？

让我们了解什么使TS非静止。 TS的非固定性背后有两个主要原因： 

基本原则是模拟或估计系列中的趋势和季节性，并从系列中删除那些以获得固定系列。然后可以在这个系列中实施统计预测技术。最后一步是通过应用趋势和季节性约束将预测值转换为原始比例。

估计和消除趋势

减少趋势的第一个技巧之一可能是转型。例如，在这种情况下，我们可以清楚地看到存在显着的积极趋势。因此，我们可以应用比较小值更多地惩罚更高值的转换。这些可以采用日志，平方根，立方根等。为简单起见，请在此处进行日志转换：

ts_log = np.log(ts)
plt.plot(ts_log)

在这种简单的情况下，很容易看到数据的前向趋势。但它在噪音的存在下不是很直观。因此，我们可以使用一些技术来估计或模拟这种趋势，然后将其从系列中删除。可以有很多方法，最常用的一些方法是：

平滑是指采用滚动估计，即考虑过去的几个实例。可以有各种方式，但我将在这里讨论其中的两个。

移动平均线

在这种方法中，我们根据时间序列的频率取“k”个连续值的平均值。在这里我们可以取过去1年的平均值，即最后12个值。 Pandas具有为确定滚动统计数据而定义的特定功能。 

# moving_avg = pd.rolling_mean(ts_log,12)
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)

这看起来像一个更好的系列。滚动值似乎略有不同，但没有具体趋势。此外，测试统计量小于5％的临界值，因此我们可以95％的置信度说这是一个固定的序列。

然而，这种特定方法的缺点是必须严格定义时间段。在这种情况下，我们可以采取年平均值，但在复杂的情况下，如预测股票价格，很难得出一个数字。因此，我们采用“加权移动平均线”，其中更近期的值被赋予更高的权重。可以有许多分配权重的技术。一种流行的是指数加权移动平均值，其中权重被分配给具有衰减因子的所有先前值。在此查找详情。这可以在Pandas中实现：

#expwighted_avg = pd.ewma(ts_log, halflife=12)
expwighted_avg = ts_log.ewm(halflife=12, min_periods=0,adjust=True,ignore_na=False).mean()
plt.plot(ts_log)
plt.plot(expwighted_avg, color='red')

注意，这里参数'halflife'用于定义指数衰减量。

ts_log_ewma_diff = ts_log - expwighted_avg
test_stationarity(ts_log_ewma_diff)

该TS的平均值和标准偏差的幅度变化甚至更小。此外，检验统计量小于1％临界值，这比前一种情况更好。

 

消除趋势和季节性 

之前讨论的简单趋势减少技术在所有情况下都不起作用，尤其是具有高季节性的技术。让我们讨论两种消除趋势和季节性的方法：

差分

在这种技术中，我们将特定时刻的观察值与前一时刻的观察值进行区分。这主要用于改善平稳性。一阶差分可以在Pandas中完成：

ts_log_diff = ts_log - ts_log.shift()
plt.plot(ts_log_diff)

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

我们可以看到均值和标准变量随时间变化很小。此外，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％临界值。所以这个TS非常接近静止。您也可以尝试高级分解技术，这可以产生更好的结果。此外，您应该注意，在这种情况下，将残差转换为未来数据的原始值并不是非常直观。

 

5.预测时间序列

我们看到了不同的技术，所有这些技术都使得TS固定不动。让我们在差分后在TS上制作模型，因为它是一种非常流行的技术。此外，在这种情况下，相对容易将噪声和季节性添加回预测残差。执行趋势和季节性估计技术后，可能存在两种情况：

让我简单介绍一下ARIMA。我不会详细介绍技术细节，但如果您希望更有效地应用它们，您应该详细了解这些概念。ARIMA代表自动回归综合移动平均线。 ARIMA对静止时间序列的预测只不过是一个线性（如线性回归）方程。预测变量取决于ARIMA模型的参数（p，d，q）：

这里一个重要的问题是如何确定'p'和'q'的值。我们使用两个图来确定这些数字。让我们先讨论一下。

差分后的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'值为：

现在，让我们考虑个人和组合效果制作3种不同的ARIMA模型。我还将为每个打印RSS。请注意，这里的RSS是残差值，而不是实际序列。

ARModel 

from statsmodels.tsa.arima_model import ARIMA
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))

MRModel

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))

在这里我们可以看到AR和MA模型具有几乎相同的RSS，但组合明显更好。现在，我们留下最后一步，即将这些值恢复到原始比例。

 

把它恢复到原始规模

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

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

请注意，这些是从'1949-02-01'而不是第一个月开始的。为什么？这是因为我们将滞后加1，并且第一个元素在它减去之前没有任何东西。将差分转换为对数比例的方法是将这些差异连续添加到基数.一种简单的方法是首先确定索引处的累积和，然后将其添加到基数。累计金额可以是：

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)))

 

标签：plt,log,航空,ARIMA,ts,TS,序列,diff,预测
来源： https://www.cnblogs.com/zheng1076/p/10728818.html