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

时间：2024-02-23 13:08:01

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

主要参考：

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

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

  • parse_dates：指定包含日期时间信息的列。如上所述，列名称为“Month”。
  • index_col：将Pandas用于TS数据背后的一个关键思想是索引必须是描述日期时间信息的变量。所以这个论点告诉pandas使用\'Month\'列作为索引。
  • date_parser：这指定一个将输入字符串转换为datetime变量的函数。默认情况下，Pandas以“YYYY-MM-DD HH：MM：SS”格式读取数据。如果数据不是这种格式，则必须手动定义格式。类似于此处定义的数据分析函数的东西可用于此目的。

现在我们可以看到数据有时间对象作为索引，＃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\']

两者都会产生以下输出：

这里有两件事需要注意：

  • 与数字索引不同，此处包含结束索引。例如，如果我们将列表索引为[：5]，那么它将返回索引处的值 - [0,1,2,3,4]。但这里的指数\'1949-05-01\'已包含在输出中。
  • 索引必须按范围排序才能工作。如果你随机乱改索引，这将无法正常工作。

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

ts[\'1949\']

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

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

 

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

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

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

  • 恒定的均值
  • 恒定的方差
  • 一种不依赖于时间的自协方差

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

plt.plot(ts)

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

  • 绘制滚动统计：我们可以绘制移动平均值或移动方差，并查看它是否随时间变化。通过移动平均值/方差，我的意思是在任何时刻“t”，我们将采用去年的平均值/方差，即过去12个月。但这又是一种视觉技术。
  • Dickey-Fuller测试：这是检查平稳性的统计测试之一。这里的零假设是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对静止时间序列的预测只不过是一个线性（如线性回归）方程。预测变量取决于ARIMA模型的参数（p，d，q）：

  • AR（自回归）项的数量（p）：AR项只是因变量的滞后。例如，如果p是5，则x（t）的预测值将是x（t-1）... .x（t-5）。
  • MA（移动平均线）项的数量（q）：MA项是预测方程中的滞后预测误差。例如，如果q为5，则x（t）的预测值将为e（t-1）... .e（t-5），其中e（i）是第i个瞬时的移动平均值与实际值之间的差值。
  • 差异数量（d）：这些是非季节性差异的数量，即在这种情况下我们采用了第一阶差异。因此，我们可以传递该变量并将d = 0或传递原始变量并放入d = 1。两者都会产生相同的结果。

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

  • 自相关函数（ACF）：它是TS与其自身滞后版本之间相关性的度量。例如，在滞后5处，ACF将在时刻\'t1\'...\'t2\'的系列与在时刻\'t1-5\'...\'t2-5\'（t1-5和t2是终点）的系列进行比较。
  • 部分自相关函数（PACF）：它测量TS与其自身的滞后版本之间的相关性，但是在消除了已经通过干预比较解释的变化之后。例如，在滞后5处，它将检查相关性，但是除去已经由滞后1到4解释的效果。

差分后的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。

现在，让我们考虑个人和组合效果制作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)))

 

相关文章