Python | 单变量时间序列分析与预测

时间序列数据是机器学习中最具挑战性的任务之一，也是与数据相关的现实问题之一，因为数据实体不仅取决于物理因素，而且主要取决于它们发生的时间顺序。我们可以基于一个单变量特征和两个双变量或多变量特征来预测时间序列中的目标值，但在本文中，我们将学习如何对来自Kaggle的Rainfall数据集执行单变量预测。

什么是单变量预测？

单变量预测通常用于您想要预测单个变量的值时，特别是当该变量存在可用的历史数据点时。它是经济、金融、天气预报和供应链管理中的需求预测等领域的基础和广泛应用的技术。

对于多个变量或外部因素可能产生影响的更复杂的预测任务，使用多变量预测技术。这些模型考虑到多个变量及其相互作用进行预测。

单变量预测的重要概念

• 趋势：时间序列的长期运动或方向由其趋势表示。它显示数据中的基本模式，例如值随时间的推移而上升或下降。确定和模拟趋势对于理解变量的总体轨迹和做出精确的预测至关重要。
• 季节性：在时间序列中，季节性是用于描述定期出现的周期性波动或模式的术语。例如，由于天气相关因素或假期，零售销售中经常出现季节性模式。考虑季节性对于识别重复发生的趋势和适当修改预测至关重要。
• 平稳性：当一个时间序列的统计特征，如其均值和方差，不随时间变化时，它被认为是平稳的。由于非平稳数据可能产生不准确的预测，平稳性是许多预测模型的重要前提。
• 时间序列数据：时间序列数据，或一系列定期进行的观察或测量，是单变量预测的主题。销售数字、温度读数、GDP增长率和股票价格就是几个例子。

单变量预测技术

在单变量时间序列分析中，有几种方法用于理解、建模和预测单个变量随时间的行为。在单变量时间序列分析中，经常使用以下方法：

• 自回归（AR）：它利用观测值与预定数量的滞后观测值（较早的时间间隔）之间的相关性。
• 移动平均（MA）：它使用应用于滞后观测值的移动平均模型来模拟观测值和残差之间的关系。
• 自回归移动平均（ARIMA）：它使时间序列平稳（即，该数据具有一致的平均值和方差随时间的推移），通过结合AR和MA方法，并考虑到原始观测值的差异。
• 季节性自回归移动平均（SARIMA）：它扩展了ARIMA，以考虑时间序列数据的季节性成分。
• 指数平滑（ETS）：它使用历史观测值的加权平均值来预测下一个时间点，为最近的观测值提供更多权重。
• 长短期记忆（LSTM）：一种RNN，专门用于识别长时间内时间序列数据中的模式。

单变量预测的实现

1.导入所需库

import pandas as pd
import numpy as np
import seaborn as sb
import matplotlib.pyplot as plt
import statsmodels.api as sm

import warnings
warnings.filterwarnings('ignore')

Python库使我们能够非常容易地处理数据，并通过一行代码执行典型和复杂的任务。

• Pandas -此库有助于以2D数组格式加载数据框，并具有多个功能，可以一次性执行分析任务。
• Numpy - Numpy数组非常快，可以在很短的时间内执行大型计算。
• Matplotlib/Seaborn -此库用于绘制可视化。
• Sklearn -该模块包含多个库，这些库具有预实现的功能，可以执行从数据预处理到模型开发和评估的任务。

2.加载数据集

df = pd.read_csv('Rainfall_data.csv')
df.head()

输出

   Year  Month  Day  Specific Humidity  Relative Humidity  Temperature  \
0  2000      1    1               8.06              48.25        23.93   
1  2000      2    1               8.73              50.81        25.83   
2  2000      3    1               8.48              42.88        26.68   
3  2000      4    1              13.79              55.69        22.49   
4  2000      5    1              17.40              70.88        19.07   
   Precipitation  
0           0.00  
1           0.11  
2           0.01  
3           0.02  
4         271.14  

在这里，在这段代码中，我们将数据集加载到pandas数据框架中，以便我们可以探索数据集的不同方面。

df.shape

输出

(252, 7)

df.info()

输出

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 252 entries, 0 to 251
Data columns (total 7 columns):
 #   Column             Non-Null Count  Dtype  
---  ------             --------------  -----  
 0   Year               252 non-null    int64  
 1   Month              252 non-null    int64  
 2   Day                252 non-null    int64  
 3   Specific Humidity  252 non-null    float64
 4   Relative Humidity  252 non-null    float64
 5   Temperature        252 non-null    float64
 6   Precipitation      252 non-null    float64
dtypes: float64(4), int64(3)
memory usage: 13.9 KB

查看数据统计量

print(df.describe().T)

输出

                   count         mean         std      min        25%  \
Year               252.0  2010.000000    6.067351  2000.00  2005.0000   
Month              252.0     6.500000    3.458922     1.00     3.7500   
Day                252.0     1.000000    0.000000     1.00     1.0000   
Specific Humidity  252.0    14.416746    4.382599     5.74    10.0100   
Relative Humidity  252.0    67.259524   17.307101    34.69    51.8450   
Temperature        252.0    16.317262    6.584842     4.73    10.8650   
Precipitation      252.0   206.798929  318.093091     0.00     0.4025   
                        50%       75%      max  
Year               2010.000  2015.000  2020.00  
Month                 6.500     9.250    12.00  
Day                   1.000     1.000     1.00  
Specific Humidity    15.200    18.875    20.57  
Relative Humidity    66.655    84.610    92.31  
Temperature          16.915    22.115    29.34  
Precipitation        11.495   353.200  1307.43  

3.探索性数据分析

EDA是一种使用可视化技术分析数据的方法。它用于发现趋势和模式，或在统计摘要和图形表示的帮助下检查假设。在执行该数据集的EDA时，我们将尝试查看独立特征之间的关系，即一个特征如何影响另一个特征。

# Feature engineering
df['Date'] = pd.to_datetime(df['Year'].astype(str) + '-' + df['Month'].astype(str) + '-' + df['Day'].astype(str),
							format='%Y-%m-%d')
df['DayOfWeek'] = df['Date'].dt.dayofweek
df['Month'] = df['Date'].dt.month
df['Quarter'] = df['Date'].dt.quarter
df['Year'] = df['Date'].dt.year

# EDA
cat_cols = ['DayOfWeek', 'Month', 'Quarter', 'Year']

for col in cat_cols:
	df[['Precipitation', col]].groupby(col).mean().plot.bar()
	plt.title(f'Mean Precipitation by {col}')
	plt.show()

输出

释义：

• x轴表示一周中的每一天（0到6，其中0是星期一，6是星期日）。
• y轴代表一周中每一天的平均降水量。
• 根据条形图的高度，分析一周中哪一天的平均降水量会更高或更低。因此，在一周的第5天，降水率高于其他日子。

释义：

• x轴代表一年中的每个月（1到12）。
• y轴代表每个月的平均降水量。
• 确定不同月份的降水模式。某些月份可能会经历比其他更多的降水。7月降水量较多。

释义：

• x轴代表一年中的每个季度（1到4）。
• y轴代表每个季度的平均降水量。
• 通过检查各季度的平均降水量如何不同来探索季节变化。这说明第三季度的降水量比其他三个季度（第一、第二和第四）多。

释义：

• x轴代表每年。
• y轴代表各年份的平均降水量。
• 评估多年来平均降水量的任何趋势或变化。此图显示2019年的降水量较其他年份多。

此代码使用带有降雨量数据的时间序列数据集进行特征工程。它从“Date”列中获取数据，并添加额外的时间元素，如“DayOfWeek”，“Month”，“Quarter”和“Year”。然后进行探索性数据分析（EDA），重点是新设计的时间特征的每个独特值的平均降水量。为了直观地评估一周、一个月、一个季度和一年中的平均降水模式，代码将遍历指定的类别列并创建条形图。通过协助发现降雨数据中可能的时间趋势和模式，这些可视化使人们能够进一步了解数据集的时间特征。

4.季节性分解

在时间序列分析中使用的一种统计技术，用于分离数据集的组成部分，称为季节分解。确定了时间序列的三个基本组成部分：趋势，季节性和残差。长期运动或方向由趋势表示，定期重复的模式由季节性捕获，随机波动由残差捕获。通过将季节性的影响与更广泛的趋势和异常分离开来，分解时间序列有助于理解各个组成部分的具体贡献，从而实现更准确的分析和预测。

# Seasonal decomposition
ts = df.set_index('Date')['Precipitation'] + 0.01 # Add a small constant
result = seasonal_decompose(ts, model='multiplicative', period=12)
result.plot()
plt.show()

对于季节性组成：

• 图的上半部分表示季节性分量。
• x轴对应于时间，通常以月为单位，给定指定的周期=12。
• y轴代表季节性变化的幅度。

对于趋势组成：

• 图的中间部分表示趋势分量。
• x轴对应于时间，反映了整个时间序列的总体趋势。
• y轴表示趋势的大小。

对于残差组成：

• 图的底部表示残差分量（也称为余数）。
• x轴对应于时间。
• y轴表示观测值与季节和趋势分量之和之间的差值。

该函数使用表示降水数据的时间序列（“ts”）执行季节分解。添加一个小常数以解决零值或负值可能出现的问题。该代码由三个部分组成：趋势，季节性和残差，它使用了一个12个月的季节性乘法模型。由此产生的图形通过提供这些元素的视觉表示，有助于确定降水数据中的长期趋势和经常性模式。

5.自相关和偏自相关图

自相关：时间序列与其滞后值的关联通过自相关来衡量。每个滞后都是相关的，并且自相关图中的峰值在特定延迟处显示出高相关性。通过揭示时间序列数据中的重复模式或季节性，这有助于理解其时间结构，并支持为时间序列分析选择合适的模型参数。

偏自相关：当测量变量与其滞后的直接相关性时，偏自相关消除了中间延迟的影响。部分自相关函数（PACF）图中的显著峰值表示特定滞后对当前观测值有直接影响。它有助于捕捉每个滞后的不同贡献，通过协助在时间序列建模中的自回归分量的适当排序。

# Autocorrelation and Partial Autocorrelation
plot_acf(ts, lags=30)
plot_pacf(ts, lags=30)
plt.show()

释义：
ACF测量时间序列与其在不同时间间隔的滞后值之间的相关性。
在ACF图中：

• x轴表示滞后或时间间隔的数量。

• y轴表示相关系数。

• 蓝色阴影区域上方的点被认为具有统计学显著性。

• 正滞后表示当前观测值与该滞后时的过去观测值之间的正相关性。

• 负滞后表示负相关。

释义：
PACF测量时间序列与其滞后值之间的相关性，控制其他滞后的影响。
在PACF图中：

• x轴表示滞后或时间间隔的数量。

• y轴表示偏相关系数。

• 蓝色阴影区域上方的点被认为具有统计学显著性。

• 特定滞后的偏自相关表示当前观测值与该滞后的过去观测值之间的相关性，不包括中间滞后的影响。

• 它有助于识别当前观测值与特定滞后之间的直接关系。

对于表示降水数据的时间序列“ts”，提供的代码创建自相关函数（ACF）和部分自相关函数（PACF）的图。时间序列与其滞后值之间的相关性显示在这些图表中，这些图表限制为30个滞后。ACF图的峰值表明在特定的滞后处可能是季节性的。通过PACF图可以更容易地为时间序列建模找到适当的自回归项，该图说明了每个滞后对当前观测值的直接影响。这些图作为选择模型参数和识别时间模式的一般参考。

6.滞后图

在时间序列分析中，滞后图是一种图形工具，用于显示变量与其滞后值之间的关系。它通过将每个数据点与之前的观察结果进行比较，帮助发现数据中的模式、趋势或不可预测性。如果点与其滞后之间存在实质相关性，则该图可显示自相关的存在。这有助于理解时间依赖性，并指导选择最佳模型进行时间序列分析。

# Lag Plots
for i in range(1, 3):
	lag_plot(ts, lag=i)
	plt.title(f'Lag Plot with lag = {i}')
	plt.show()

释义：

• x轴表示时间序列在时间t的值。

• y轴表示时间序列在时间t+1（滞后1）的值。

• 如果图中的点遵循明确定义的模式或趋势，则表示滞后1处的自相关。

• 如果这些点是随机分散的，则表明在滞后1处缺乏自相关。

释义：

• 检查时间t的值与时间t+2的值之间的关系。
• 帮助识别滞后2处的自相关。

将滞后设置为1和2时，代码将为时间序列“ts”创建滞后图。在每次迭代中绘制每个数据点与其指定延迟处的先行观测之间的关联。通过揭示自相关模式，这些表示可以更容易地发现数据中可能的趋势和时间关系。为了使每个滞后图更容易理解和评估自相关的强度，标题提供了滞后值。

7.平稳性检查

为了确保数据集的统计特征不随时间变化，平稳性检查是时间序列分析中的一个重要步骤。当时间序列是平稳的时，模型预测更容易，因为它的均值、方差和自相关保持不变。视觉检查，如滚动统计图，和正式的统计测试，如增强的Dickey-Fuller测试，是常见的方法。通过减少非常数模式的影响，平稳性确保了可靠的建模，并有助于更精确的预测和时间序列数据的趋势分析。

# Stationarity check
adf_result = adfuller(ts)
print('ADF Statistic:', adf_result[0])
print('p-value:', adf_result[1])

输出

ADF Statistic: -2.46632490177327
p-value: 0.12388427626757847

代码使用增强的Dickey-Fuller检验来验证时间序列“ts”的平稳性。将显示p值和ADF统计量。p值表示ADF统计量估计时间序列是否存在单位根的重要性。对于更平稳的时间序列，较低的p值和更负的ADF统计量指向更多反对平稳性的证据，这有助于确定是否需要差分。

8.滚动和汇总

滚动：滚动是一种用于时间序列分析的统计方法，它计算数据集连续子集的汇总统计量，例如移动平均值。固定大小的窗口在滚动统计上下文中遍历数据，并且在每个步骤中基于该窗口内的观测值计算新值。这有助于通过减少波动性和突出时间序列中的趋势或模式来可视化和理解基本动态。

汇总：汇总是时间序列分析中确定大趋势的常用技术。它们需要将多个数据点合并和汇总为一个值。在这个意义上，计算特定时间间隔内的统计措施，如平均数或总数，可以被认为是汇总数据。汇总通过结合相关信息简化了解释复杂时间序列数据的过程。这使得分析师更容易发现趋势，模式或异常，并促进更有效的预测和决策。

# Rolling and Aggregations
rolling_mean = ts.rolling(window=12).mean()
rolling_std = ts.rolling(window=12).std()

plt.plot(ts, label='Actual Data')
plt.plot(rolling_mean, label='Rolling Mean')
plt.plot(rolling_std, label='Rolling Std')
plt.legend()
plt.show()

释义：
蓝线表示原始时间序列数据。

滚动平均值：

• 橙线表示时间序列的滚动平均值。

• 在12个数据点的窗口（在这种情况下为月）内计算滚动平均值。

• 如果滚动平均数平滑了实际数据中的波动，它有助于确定趋势。

• 上升或下降趋势可以通过比较滚动平均数和实际数据来观察。

滚动标准差：

• 绿色线表示时间序列的滚动标准差。

• 在12个数据点的窗口内计算滚动标准差。

• 表示时间序列的波动性或可变性。

• 滚动标准差中的峰值可以表示可变性增加的时期。

该代码计算滚动平均值和滚动标准差，窗口大小为12，对时间序列“ts”执行滚动统计。滚动平均值、滚动标准差和实际数据叠加在生成的图形上，以帮助可视化模式和变异性。这种方法有助于降噪，增强底层模式的可见性，并提供对数据时间属性的洞察。为便于理解，图例将滚动平均值、滚动标准差和原始数据分开。

9.构建模型

我们将通过使用日期列作为预测的特征来训练单变量预测的SARIMA模型。但首先，我们必须在数据集中创建一个日期列，也在pd中。DateTime格式，因此，我们将使用pandas框架中提供的pd.to_datetime函数。

# Combine 'Year' and 'Month' columns to create a 'Date' column
df['Date'] = pd.to_datetime(df['Year'].astype(str) +\
							'-' + df['Month'].astype(str) +\
							'-' + df['Day'].astype(str),
							format='%Y-%m-%d')
df.head()

输出

   Year  Month  Day  Specific Humidity  Relative Humidity  Temperature  \
0  2000      1    1               8.06              48.25        23.93   
1  2000      2    1               8.73              50.81        25.83   
2  2000      3    1               8.48              42.88        26.68   
3  2000      4    1              13.79              55.69        22.49   
4  2000      5    1              17.40              70.88        19.07   
   Precipitation       Date  
0           0.00 2000-01-01  
1           0.11 2000-02-01  
2           0.01 2000-03-01  
3           0.02 2000-04-01  
4         271.14 2000-05-01  

现在让我们将索引设置为日期列，目标列是降水列，让我们将其与完整的数据集分离。

ts = df.set_index('Date')['Precipitation']
ts

输出

Date
2000-01-01       0.00
2000-02-01       0.11
2000-03-01       0.01
2000-04-01       0.02
2000-05-01     271.14
               ...   
2020-08-01    1203.09
2020-09-01     361.30
2020-10-01     180.18
2020-11-01       0.49
2020-12-01      12.23
Name: Precipitation, Length: 252, dtype: float64

10.训练SARIMA模型

# Fit a SARIMA model
# We need to specify the order (p, d, q) and seasonal order (P, D, Q, S)
# Example order and seasonal order values
p, d, q = 1, 1, 1

# For monthly data with a yearly seasonality
P, D, Q, S = 1, 1, 1, 12

model = sm.tsa.SARIMAX(ts, order=(p, d, q), seasonal_order=(P, D, Q, S))
results = model.fit()

由于模型已经过训练，我们可以使用这个模型来预测下一年的降雨量，并将其与原始数据一起绘制沿着，以了解预测是否遵循以前的趋势。

该代码将时间序列ts拟合到SARIMA模型。模型的阶数由阶数=（p，d，q）参数定义，其中p表示自回归项的数量，d表示差分度，q表示移动平均项的数量。模型的季节分量的阶数由seasonal_order=（P，D，Q，S）参数指定，其中P表示季节自回归项的数量，D表示季节差异的程度，Q表示季节移动平均项的数量，S表示季节性周期的持续时间。

SARIMA模型对象由代码行model = sm.tsa.SARIMAX（ts，order=（p，d，q），seasonal_order=（P，D，Q，S））创建。需要建模的时间序列是ts参数。模型的阶数和季节阶数分别由order=（p，d，q）和seasonal_order=（P，D，Q，S）参数指定。

11.预测

# Make predictions
forecast = results.get_forecast(steps=12)  # Forecast the next 12 periods
forecast_mean = forecast.predicted_mean
 
# Plot the actual data and the forecast
plt.figure(figsize=(12, 6))
plt.plot(ts, label='Actual Data')
plt.plot(forecast_mean, label='SARIMA Forecast')
plt.legend()
plt.show()

此代码绘制实际数据和预测，并使用拟合的SARIMA模型进行预测。最初，使用get_forecast（）函数创建一个包含即将到来的12个时段的预期值的预测对象。为了获得预测的平均值，然后提取预测对象的predicted_mean属性。接下来，制作一个图，并使用plot（）函数绘制预测和实际数据。现在，每行都有一个标签，并显示一个图例。