ARIMA时间序列模型建立与求解——加利福尼亚州每天的女性出生人数预测

ARIMA时间序列模型

ARIMA（AutoRegressive Integrated Moving Average，自回归差分移动平均）模型是一种广泛用于时间序列分析的统计模型。它将自回归 (AR)、差分 (I) 和移动平均 (MA) 三个部分结合起来，适用于平稳和非平稳时间序列的数据建模与预测。

模型构成

ARIMA 模型由三个主要成分组成：

• 自回归 (AR): 通过时间序列自身的历史值预测当前值。
• 差分 (I): 用于消除时间序列中的趋势或季节性，使其成为平稳序列。
• 移动平均 (MA): 通过过去的预测误差来调整当前值的预测。

ARIMA 模型通常表示为 ARIMA(p, d, q)，其中：

• p: 自回归项的阶数，即使用多少个过去的值来预测当前值。
• d: 差分阶数，即需要几次差分操作才能使时间序列平稳。
• q: 移动平均项的阶数，即使用多少个过去的误差项来预测当前值。

模型建立

差分过程

对于非平稳的时间序列 $y_t$，通过对序列进行差分操作来消除趋势，使其成为平稳序列。差分操作公式为：
$$y_t^{\prime }=y_t-y_{t-1}$$
如果序列在第一次差分后仍不平稳，则可以进行二次差分：
$$y_t^{\prime \prime }=y_t^{\prime }-y_{t-1}^{\prime }=(y_t-y_{t-1})-(y_{t-1}-y_{t-2})=y_t-2y_{t-1}+y_{t-2}$$
一般情况下，经过 $d$次差分后，序列可以转化为平稳序列,

自回归模型 AR($p$)

自回归模型表示当前值与前 $p$个时刻值之间的线性关系，其数学表达式为：
$$y_t={\varphi }_1{y}_{t-1}+{\varphi }_2{y}_{t-2}+\cdots +{\varphi }_p{y}_{t-p}+{ϵ}_t$$
其中：

• ${\varphi }_1,{\varphi }_2,\dots ,{\varphi }_p$为自回归系数。
• ${ϵ}_t$为白噪声。

滑动平均模型MA($q$)

滑动平均模型表示当前值与前 $q$个时刻的误差项之间的关系，其数学表达式为：
$$y_t={ϵ}_t+{\theta }_1{ϵ}_{t-1}+{\theta }_2{ϵ}_{t-2}+\cdots +{\theta }_q{ϵ}_{t-q}$$
其中：

• ${\theta }_1,{\theta }_2,\dots ,{\theta }_q$ 为滑动平均系数。
• ${ϵ}_t$为白噪声。

综合模型 ARIMA($p,d,q$)

将自回归部分和滑动平均部分结合，并考虑差分操作后，ARIMA模型的表达式为：
$${\Delta }^d{y}_t={\varphi }_1{\Delta }^d{y}_{t-1}+{\varphi }_2{\Delta }^d{y}_{t-2}+\cdots +{\varphi }_p{\Delta }^d{y}_{t-p}+{ϵ}_t+{\theta }_1{ϵ}_{t-1}+{\theta }_2{ϵ}_{t-2}+\cdots +{\theta }_q{ϵ}_{t-q}$$
其中：

• ${\Delta }^d{y}_t$表示对序列$y_t$进行$d$次差分后的结果。

模型求解

最大似然估计：
$$\text{Log-Likelihood Function}\ell (\theta )=-\frac{n}{2}\log (2\pi {\sigma }^2)-\frac{1}{2{\sigma }^2}\sum _{t=1}^n{ϵ}_t^2$$
其中，${ϵ}_t=X_t-{\hat{X}}_t$ 是预测误差，${\sigma }^2$ 是误差的方差，$\theta$ 表示模型参数的集合，包括 ${\varphi }_i$ 和 ${\theta }_i$。

使用估计出的模型参数进行预测：
$${\hat{X}}_{t+1}={\varphi }_1{X}_t+{\varphi }_2{X}_{t-1}+\cdots +{\varphi }_p{X}_{t+1-p}+{ϵ}_{t+1}+{\theta }_1{ϵ}_t+\cdots +{\theta }_q{ϵ}_{t+1-q}$$
对未来的时间点进行预测。

使用ARIMA预测1959年加利福尼亚州每天的女性出生人数。（Python语言实现）

数据集 链接：https://pan.baidu.com/s/1Njfx5YlGj07rKcqfGZGFbw 提取码：1111

问题描述

我们手头有一个数据集，记录了1959年加利福尼亚州每天女性出生人数。这类数据是时间序列数据，因为它包含了按时间顺序排列的观测值。在这个项目中，我们的目标是根据已有的历史数据预测未来某个时间点的女性出生人数。

数据准备与初步分析

导入所需的库

# 导入库
import pandas as pd 
import matplotlib.pyplot as plt 
from sklearn.metrics import mean_squared_error 
import warnings 

from statsmodels.tsa.stattools import adfuller 
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf 
from statsmodels.tsa.arima.model import ARIMA 

# 关闭警告
warnings.filterwarnings( 'ignore' )

数据预处理

数据集的加载

# 加载数据集
data = pd.read_csv( "daily-total-female-births-CA.csv" ) 
data.head()

检查数据集的基本信息

data.info()

统计缺失值的数目

data.isnull().sum()

时间格式转换：将日期列转换为日期时间格式，以便能够按时间顺序对数据进行处理和分析。

data["Date"] = pd.to_datetime(data["Date"])
data.set_index("Date", inplace=True)

• 第一行数据框 data 中的“Date”列转换为 datetime 对象。pd.to_datetime() 函数可以将字符串格式的日期转换为 pandas 的 datetime 对象，从而便于之后的时间序列分析操作。

• 第二行将“Date”列设置为数据框 data 的索引（index）。使用 set_index() 方法可以将指定的列作为数据框的索引，使得每一行的数据可以根据日期进行访问和操作。

时间序列可视化：通过绘制出生人数随时间变化的图形，观察数据的总体趋势、季节性、周期性等特征。

# plot the time series
plt.plot(data["Births"])
plt.title('Daily Total Female Births in California')
plt.xlabel('Date')
plt.ylabel('Number of births')
plt.show()

模型构成

ARIMA模型是一种用于时间序列预测的常用模型，它由以下三个部分组成：

• 自回归部分 (AR)：表示当前值与之前若干时刻的值之间的关系。模型中的阶数由参数$p$表示。
• 差分部分 (I)：用于处理时间序列的非平稳性，通过计算时间序列的差分来使其变为平稳序列。差分的阶数由参数$d$表示。
• 移动平均部分 (MA)：表示当前值与之前若干时刻的误差之间的关系。模型中的阶数由参数$q$表示。

因此，ARIMA模型可以表示为$ARIMA(p,d,q)$，其形式为：
$$y_t=c+{\varphi }_1{y}_{t-1}+{\varphi }_2{y}_{t-2}+\cdots +{\varphi }_p{y}_{t-p}+{ϵ}_t+{\theta }_1{ϵ}_{t-1}+{\theta }_2{ϵ}_{t-2}+\cdots +{\theta }_q{ϵ}_{t-q}$$
其中：

• $y_t$为当前时刻的观测值。
• ${\varphi }_1,{\varphi }_2,\dots ,{\varphi }_p$为AR部分的系数。
• ${\theta }_1,{\theta }_2,\dots ,{\theta }_q$为MA部分的系数。
• $c$为常数项（模型截距）。
• ${ϵ}_t$为白噪声误差项。

模型建立

确定差分阶数$d$：通过单位根检验（如ADF检验）判断时间序列是否平稳，若非平稳则进行差分处理，直至序列平稳为止。差分次数即为参数$d$。

adf_test = adfuller(data["Births"])
print('ADF Statistic: %f' % adf_test[0])
print('p-value: %f' % adf_test[1])

• adfuller 函数是 statsmodels 库中的一个函数，用于执行 Augmented Dickey-Fuller (ADF) 检验。该函数的输入是一个时间序列数据，这里是 data["Births"]，即每天的女性出生人数数据。
• adf_test[0]表示ADF统计量。ADF统计量是用于检验时间序列是否具有单位根的核心统计量。数值越小，表明时间序列越有可能是平稳的。
• adf_test[1]表示p值。p值用于判断检验结果的显著性。通常，如果p值小于设定的显著性水平（如0.05），则拒绝原假设，认为时间序列是平稳的；如果p值大于显著性水平，则无法拒绝原假设，时间序列可能具有单位根，即可能是非平稳的。

输出的结果如图所示：

• 这个值是Augmented Dickey-Fuller (ADF) 检验的统计量。ADF统计量衡量的是时间序列的单位根假设的程度。数值越小（越负），通常意味着时间序列越可能是平稳的。-4.808291 是一个相对较小的（负）值，表明时间序列数据很可能是平稳的。
• 在这个结果中，p值为0.000052，非常小，远远小于常用的显著性水平（如0.05、0.01等）。这表明有足够的证据拒绝原假设，意味着时间序列数据很可能是平稳的
  

所以数据符合标准，不需要进行差分。

PACF和ACF图：确定ARIMA模型的阶数时，可以使用偏自相关函数（PACF）和自相关函数（ACF）图来帮助选择参数$p$和$q$。PACF图有助于识别AR项的阶数，而ACF图有助于识别MA项的阶数。

plot_acf(data['Births'], lags=40)
plot_pacf(data['Births'], lags=40)
plt.show()

• ACF 图展示了自相关系数与滞后期（lag）之间的关系。通常，ACF 图会随着滞后期的增加而迅速下降。

• 显著自相关: 如果在某些滞后期的自相关系数超出了置信区间（通常是虚线），这些滞后期被认为是显著的。

• 截尾: 如果自相关系数在某个滞后期后迅速降到零（即接近于零并在置信区间内），则MA模型的阶数通常是这个滞后期的位置。(如下图所示)

本题中的ACF图（上图）就存在截尾现象，在滞后1处有显著的自相关，然后迅速下降到零，说明可以尝试MA(1)模型。

• PACF 图展示了在去除之前所有滞后期的影响后，当前滞后期的偏自相关系数与滞后期（lag）之间的关系。
• 显著偏自相关: 如果在某些滞后期的偏自相关系数超出了置信区间（通常是虚线），这些滞后期被认为是显著的。
• 截尾: 如果偏自相关系数在某个滞后期后迅速降到零（即接近于零并在置信区间内），则AR模型的阶数通常是这个滞后期的位置。

本题中的PACF图（上图）存在截尾现象，在滞后1处有显著的偏自相关，然后迅速下降到零，说明可以尝试AR(1)模型。

综上所述， ARIMA（1, 0, 1）模型可能是一个很好的起点。

model = ARIMA(data['Births'], order=(1, 0, 1))
model_fit = model.fit()

• ARIMA(data['Births'], order=(1, 0, 1)) 创建一个ARIMA模型的实例。1 是自回归部分（AR）的阶数，表示模型中包括一个自回归项。0 是差分部分（I）的阶数，表示数据没有进行差分处理，即数据是平稳的（或者已经在之前的步骤中处理过差分）。1 是移动平均部分（MA）的阶数，表示模型中包括一个移动平均项。
• model.fit() 方法用于拟合模型。它通过最大似然估计来估计模型参数，使得模型的预测最符合给定的数据。

训练和预测

forecast = model_fit.get_forecast(steps=30)

• model_fit 是之前拟合好的ARIMA模型对象，包含了模型的参数和拟合信息。
• get_forecast(steps=30) 方法用于生成未来30天的预测值。steps 参数指定了预测的时间步数，即预测的未来时间段长度。
• forecast 是一个 ForecastResults 对象，包含了未来30天的预测值及相关统计信息。

模型评估

数据分割

# 将数据集拆分为训练集和测试集
train_size = int(len(data) * 0.8)  # 计算训练集的大小，数据总长度的80%
train, test = data[0:train_size], data[train_size:len(data)]  # 切分数据集为训练集和测试集

• train_size：计算训练集的大小，通常是数据总长度的80%。
• train：训练集数据，从数据开始到 train_size 的索引。
• test：测试集数据，从 train_size 的索引到数据结束。

拟合ARIMA模型

# 在训练集上拟合ARIMA模型
model_train = ARIMA(train['Births'], order=(1, 0, 1))  # 定义ARIMA模型，参数为(p=1, d=0, q=1)
model_train_fit = model_train.fit()  # 拟合ARIMA模型

• ARIMA(train['Births'], order=(1, 0, 1))：定义一个ARIMA模型，其中 (p, d, q) 的参数为 (1, 0, 1)。这里的参数通常表示
• model_train_fit：拟合后的ARIMA模型对象。

对测试集进行预测

# 在测试集上进行预测
test_forecast = model_train_fit.get_forecast(steps=len(test))  # 预测测试集长度的数据
test_forecast_series = pd.Series(test_forecast.predicted_mean, index=test.index)  # 将预测结果转换为Series，并设置与测试集相同的索引

• model_train_fit.get_forecast(steps=len(test))：预测未来 len(test) 个时间步的数据，即测试集的长度。
• test_forecast_series：将预测的均值转换为一个新的Series对象，并设置与测试集相同的索引。

计算误差指标

RMSE 是均方误差（MSE, Mean Squared Error）的平方根。均方误差是预测值与实际值之间差异的平方的平均值。具体计算公式为：
$$\mathrm{MSE}=\frac{1}{n}\sum _{i=1}^n(y_i-{\hat{y}}_i{)}^2$$

$$RMSE=\sqrt{MSE}$$

其中：

• $y_i$是实际值。
• ${\hat{y}}_i$是预测值。
• $n$是样本数量。

# 计算均方误差（MSE）和均方根误差（RMSE）
mse = mean_squared_error(test['Births'], test_forecast_series)  # 计算均方误差
rmse = mse**0.5  # 计算均方根误差

• mean_squared_error(test['Births'], test_forecast_series)：计算均方误差（MSE），衡量预测值与实际值之间的平均平方差异。

• rmse：均方根误差（RMSE），是MSE的平方根，更容易解释为实际单位的误差。

均方根误差（RMSE）： 6.970853444715778

表明预测效果良好。

绘制结果图

# 绘制训练数据、实际数据和预测数据的比较图
plt.figure(figsize=(14,7))  # 设置图形的大小
plt.plot(train['Births'], label='训练数据')  # 绘制训练数据
plt.plot(test['Births'], label='实际数据', color='orange')  # 绘制测试集中的实际数据
plt.plot(test_forecast_series, label='预测数据', color='green')  # 绘制预测的数据
plt.fill_between(test.index, 
                 test_forecast.conf_int().iloc[:, 0], 
                 test_forecast.conf_int().iloc[:, 1], 
                 color='k', alpha=.15)  # 绘制预测的置信区间，表示预测的不确定性
plt.title('ARIMA模型评估')  # 设置图标题
plt.xlabel('日期')  # 设置x轴标签
plt.ylabel('出生人数')  # 设置y轴标签
plt.legend()  # 显示图例
plt.show()  # 显示图形

# 输出均方根误差（RMSE）
print('均方根误差（RMSE）：', rmse)  # 打印RMSE值，以评估预测的准确性

• plt.plot(train['Births'], label='Training Data')：绘制训练数据。
• plt.plot(test['Births'], label='Actual Data', color='orange')：绘制测试集中的实际数据。
• plt.plot(test_forecast_series, label='Forecasted Data', color='green')：绘制预测的数据。
• plt.fill_between(test.index, test_forecast.conf_int().iloc[:, 0], test_forecast.conf_int().iloc[:, 1], color='k', alpha=.15)：绘制预测的置信区间，以表示预测的不确定性。

ARIMA 模型是时间序列预测的强大工具。通过了解时间序列数据中的潜在模式并使用统计测试和图表来确定模型参数，我们可以构建一个提供准确可靠预测的模型。