文献阅读+ARIMA模型学习

摘要

本周阅读的文献中提出了一种基于趋势分量和随机增量特征的ARIMA模型与K-means聚类模型相结合的水质预测方法，将ARIMA模型与聚类模型相结合，可以弥补单一ARIMA模型不能考虑具有随机变化特征数据的不足。通过聚类分析发现可能存在的影响因素，提高对于具有随机特征的数据预测精度。文献中采用自AR模型、MA模型和差分法结合的ARIMA模型能够对数据本身的历史信息来预测未来信息和k-means聚类算法、能够把相似的数据样本聚到一组。同时本文利用阳朔站点检测的水流样本数据使用该方法进行模拟预测。

Abstract

The literature read this week proposes a water quality prediction method that combines an ARIMA model based on trend components and random increment characteristics with a K-means clustering model. By integrating the ARIMA model with the clustering model, it can compensate for the insufficiency of a single ARIMA model in considering data with trend characteristics. Through cluster analysis, potential influencing factors can be identified, which improves the prediction accuracy for data with random characteristics. The literature employs an ARIMA model that combines autoregressive (AR) models, moving average (MA) models, and differencing, which can predict future information based on the historical information of the data itself. The K-means clustering algorithm is capable of grouping similar data samples together. Finally, the paper uses the water flow sample data detected at the Yangshuo station to simulate predictions using this method.

1 文献阅读

文献题目：Application of Time Serial Model in Water Quality Predicting

该文献提出了一种结合ARIMA和聚类模型的水质预测方法。以某流域水质监测数据为样本，选取水质TP指数作为预测对象。首先对样本数据进行清洗、平稳性分析、白噪声分析。其次，根据BIC原理选取合适的参数，利用ARIMA方法获得趋势分量特征，进行水质预测；再采用k -means聚类方法分析监测水域降水与TP指数的关系，计算降水对水质变化的随机增量特征。最后，结合趋势分量特征和随机增量特征，计算水质预测结果。

现有问题： 传统的时间序列ARIMA模型对具有趋势特征的数据预测效果较好，但对具有随机特征的数据预测效果较差。
提出办法： 针对水质数据变化的趋势性和随机性，将ARIMA
模型与聚类模型相结合，提出了一种基于趋势分量和随机增量特征的ARIMA模型与聚类模型相结合的水质预测方法。

1.1 相关工作

• 时间序列
  用ARIMA模型拟合时间序列数据有几个基本步骤，主要包括构造数据的时间序列图、执行适当的数据转换、模型排序、参数估计、模型诊断和模型选择

• k-means 聚类
  k -means算法是一种简单的迭代方法，用于将给定的数据集划分为用户指定数量的k个聚类；每个聚类中的数据特征在某种意义上是相似的。
  k -means将整个样本集分成k个组，同一组中样本之间的距离最小；常用的计算距离的方法是欧拉距离：
  假设$x=(x_1,x_2,.....x_n)和y=(y_1,y_2,...y_n)$为n维向量。
  $$d=\sqrt{\sum _{i=1}^n|x_i-y_i{|}^2}$$

1. 首先需要选择k个质心，这些质心作为每组的中心点。
2. 将样本集中的每个特征值分类到离该特征值距离最小的质心所在的组中。计算并分类完所有特征值后，完成k个聚类分组。然后再计算每一组的质心，也就是计算该组中所有特征值数据的平均值。
3. 当每一组的质心发生变化时，重新计算质心，再对特征值进行计算和分类
4. 重复第2步第3步，直至每个组的质心都不发生改变

• 评价指标
  使用MAE、MSE和MAPE检验模型的预测精度;其中$z$为原始序列值，$z_j$为预测值，n为样本数

1. MAE是绝对误差的平均值，能较好地反映预测误差值的实际情况。
   $$MAE=\frac{1}{n}\sum _{j=1}^n|z-z_j|$$

2. MSE是均方误差，是指参数预测值与实际值之差的平方的期望值。MSE评估数据的波动。MSE越小，模型的预测精度越高。

$$MSE=\frac{1}{n}\sum _{j=1}^n(z-z_j{)}^2$$

3. MAPE是平均绝对百分比误差。MAPE值越小，模型的预测精度越高。
   $$MAPE=\frac{100\%}{n}\sum _{j=1}^n|\frac{z-z_j}{z}|$$

1.2 数据集

数据通过水质监测平台获取某流域2019年1月1日至2020年12月31日的水质数据。共有16876行数据，主要有水温、pH、溶解氧、电导率、浊度、高锰酸盐指数、氨氮、TP和总氮9项指标。其中，高锰酸盐指标、氨氮、总磷、总氮每4 h采集一次，其他5项指标每小时采集一次。

以某流域水质监测数据为样本，选取水质总磷(TP)指数作为预测对象。

1.3 实验步骤

1. 首先对数据进行清洗，平稳性检验、白噪声检验分析
2. 根据贝叶斯信息准则（BIC）热图确定合适的参数p,q值；利用ARIMA方法获得趋势分量特征，进行水质预测；
3. 采用k均值聚类方法分析监测水域降水与TP指数的关系，计算降水对水质变化的随机增量特征。
4. 结合趋势分量特征和随机增量特征，计算水质预测结果，比较ARIMA模型与ARIMA结合聚类模型两种方法的预测结果。

1.4 实验结果

对TP数据进行聚类分析如图所示，整体数据范围在0.02-0.035之间，有少量数据分布在这个范围之外，有少量数据分布在这个范围之外。故分析监测水域附近降雨对TP指数的影响。通过相关分析可以确定降雨量增量与TP增量之间存在很强的相关性。

只使用ARIMA模型进行预测结果
ARIMA模型结合聚类模型的预测效果

从评价指标也可以看出，利用ARIMA模型预测TP指数，其预测精度不高。这是因为ARIMA模型只提取了水质变化的趋势分量特征，而不能提取随机变化特征。与ARIMA水质预测方法相比，ARIMA结合聚类的方法具有更高的精度，其平均绝对误差(MAE)、均方误差(MSE)和平均绝对百分比误差(MAPE)分别降低了44.6%、56.8%和45.8%

1.5 ARIMA模型python实现

数据集：阳朔站点水流数据2020-01-01 00_2020-12-31。

平稳性检验

from statsmodels.tsa.stattools import adfuller


def test_stationarity(timeseries):
    # Determing rolling statistics
    rolmean = timeseries.rolling(12).mean()  #计算时间序列的12个月滚动平均值
    rolstd = timeseries.rolling(12).std()  #计算时间序列的12个月滚动方差
    # 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)
    # 保存图像到文件
    plt.savefig('E:\用户文件\Desktop\water_predict\my_pic\stationarity_test.png', bbox_inches='tight')
    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_3mon)

绘制ACF图和PACF图

import matplotlib.pyplot as plt
from statsmodels.graphics.tsaplots import plot_acf

# 假设 ts 是您的时间序列数据
# ts = ...
ts_3mon=ts.iloc[:763]
print(ts_3mon)
# 绘制自相关图
plot_acf(ts_3mon, lags=40)  # lags 参数表示您希望计算的滞后项数
plt.xlabel('delay order')  # 设置x轴标签
plt.ylabel('autocorrelation coeffieient')  # 设置y轴标签
plt.title('Autocorrelation Function')
plt.savefig('E:\用户文件\Desktop\water_predict\my_pic\ACF.png', bbox_inches='tight')
plt.show(block=False)


# 绘制偏自相关图
plot_pacf(ts_3mon, lags=40)  # lags 参数表示您希望计算的滞后项数
plt.xlabel('delay order')  # 设置x轴标签
plt.ylabel('autocorrelation coeffieient')  # 设置y轴标签
plt.title('Partial Autocorrelation Function')
plt.savefig('E:\用户文件\Desktop\water_predict\my_pic\PACF.png', bbox_inches='tight')
plt.show(block=False)

ACF图

PACF图

根据AIC热图确定p,q值为p=6,q=0;

拟合效果图

#拟合模型
import pandas as pd
from statsmodels.tsa.arima.model import ARIMA


# 假设 ts 是您的时间序列数据
obj=data['总磷mg/L']
print(obj)
ts = pd.Series(obj, index=data.index)
ts = pd.to_numeric(ts, errors='coerce')# 尝试将 timeseries 转换为数值类型，无法转换的值将被设置为 NaN
ts = ts.dropna()# 删除含有 NaN 的行
print('以下是ts:')
print(ts)
ts_3mon=ts.iloc[:763]

ts1 = pd.read_csv('E:\用户文件\Desktop\water_predict\my_new_data_dropped.csv', index_col='时间', parse_dates=True)
ts1=ts1['总磷mg/L']
print('以下是ts1:')
print(ts1)   #ts1不能拟合，因为数据类型不对
# 选择模型参数
p, d, q = 6, 1, 0  # 这些值应该根据您的数据进行选择

# 拟合ARIMA模型
model = ARIMA(ts_3mon, order=(p, d, q))
results = model.fit()

#绘制原始数据和拟合数据
import matplotlib.pyplot as plt

# 绘制原始时间序列
plt.figure(figsize=(12, 6))
plt.plot(ts_3mon, label='TP')

# 绘制拟合数据
fitted_values = results.fittedvalues
plt.plot(fitted_values, label='forecast', color='red')
plt.title('Fitting Effect of ARIMA Model')
plt.xlabel('Time')
plt.ylabel('TP')
plt.legend()
plt.savefig('E:\用户文件\Desktop\water_predict\my_pic\/fitting_effect.png', bbox_inches='tight')
plt.show(block=False)

2 ARIMA 模型

自回归差分移动平均模型（Autoregressive Integrated Moving Average Model，ARIMA）
ARIMA模型的基本构成：

1. 自回归（AR）部分
2. 差分（I）部分
3. 移动平均（MA）部分

1. 自回归模型(AR)
   描述当前值和历史值之间的关系
   用变量自身的历史数据对自身进行预测
   其必须要满足平稳性要求。*p为自回归项

$$Y_t=c+{\phi }_1{Y}_{t-1}+{\phi }_2{Y}_{t-2}+......+{\phi }_p{Y}_{t-p}+{ϵ}_t$$
或

$$y_t=\mu +\sum _{i=1}^P{\gamma }_i{y}_{t-1}+e_t$$

2. 移动平均模型(MA)
   描述自回归模型中的误差项累加。
   有效地消除预测中的随机波动。*q为移动平均项数
   $$Y_t=c+{ϵ}_t+{\theta }_1{ϵ}_{t-1}+{\theta }_2{ϵ}_{t-2}+......+{\theta }_p{ϵ}_{t-p}$$
   或
   $$y_t=\mu +e_t+\sum _{i=1}^P{\theta }_i{e}_{t-1}$$

3. 自回归移动平均模型(ARMA)
   自回归和移动平均模型的结合。
   $$Y_t=c+{\phi }_1{Y}_{t-1}+{\phi }_2{Y}_{t-2}+......+{\phi }_p{Y}_{t-p}+{\theta }_1{ϵ}_{t-1}$$
   $$+{\theta }_2{ϵ}_{t-2}+......+{\theta }_p{ϵ}_{t-p}+{ϵ}_t$$

或

$$y_t=\mu +\sum _{i=1}^P{\gamma }_i{y}_{t-1}+e_t+\sum _{i=1}^P{\theta }_i{e}_{t-1}$$
4. 自回归差分移动平均模型(ARIMA)
将自回归模型（AR）、移动平均模型（MA）和差分法结合
对非平稳时间序列进行d阶差分运算得到平稳序列，进而，对其通过ARMA(p，q)模型预测。因此我们要确定(p，q，d)
*d为时间序列成为平稳时所做的差分次数

$y_i$为当期值，$\mu$为常数项，P为阶数，${\gamma }_i$为自相关系数，$e_t$为误差项

差分法实现
在时间序列分析中，差分法用于将非平稳时间序列转换为平稳时间序列。差分操作可以帮助消除数据中的趋势和季节性成分，使其更适合进行统计分析和预测

1. 一阶差分
   一阶差分是指将时间序列中的每一项与其前一项相减，得到的新序列称为一阶差分序列。
   若$Y_i$是原始时间序列，那么一阶差分$\Delta Y_t=Y_t-Y_{t-1}$

2. 二阶差分
   二阶差分是在一阶差分的基础上再次进行差分操作。如果$\Delta Y_t$是一阶差分序列，那么二阶差分${\Delta }^2{Y}_t=\Delta Y_t-\Delta Y_{t-1}=(Y_t-Y_{t-1})-(Y_{t-1}-Y_{t-2})=Y_t-2Y_{t-1}+Y_{t-2}$

示例：

假设有一个简单的时间序列数据：
$Y_t=[10,12,13,15,18,20]$
一阶差分$\Delta Y_t=[NaN,12-10,13-12,15-13,18-15,20-18]=[NaN,2,1,2,3,2]$
二阶差分${\Delta }^2{Y}_t=[NaN,1-2,2-1,3-2,2-3]=[NaN,-1,1,1,-1]$

import pandas as pd

# 创建一个时间序列
data = [10, 12, 13, 15, 18, 20]
ts = pd.Series(data)

# 计算一阶差分
first_diff = ts.diff()

# 计算二阶差分
second_diff = first_diff.diff()

print("原始时间序列:")
print(ts)
print("\n一阶差分:")
print(first_diff)
print("\n二阶差分:")
print(second_diff)

DataFrame.diff(periods=1, axis=0)

• periods: 整数，默认为1，表示计算相隔多少行或列的差值。
• axis: 整数（0或1），默认为0，表示沿哪个轴进行差分运算。axis=0表示按行计算差值（即每一列的上下元素差分），axis=1表示按列计算差值（即每一行的左右元素差分）

2.1 ARIMA(p,d,q)参数

1. p和q是什么

在 ARIMA(p, d, q) 模型中：

• p 代表 “自回归部分 (Autoregressive)”： 这部分描述了模型中使用的观测值的滞后值（即前面 p 个期的值）。自回归模型的出发点是认为观测值是它前面的 p 个值的线性组合。具体的数学形式如下：
  

• q 代表 “移动平均部分 (Moving Average)”： 这部分描述了模型中使用的错误项的滞后值（即前面 q 个期的值）。移动平均模型是将当前值和过去的白噪声之间建立关系。具体的数学形式如下：

2. d是什么

ARIMA模型中的d代表的就是差分阶数。我们在ARIMA模型中设定d等于一个特定的数值时，我们实际上是在告诉模型，我们应用了多少次滞后运算（也就是进行了多少次差分）来使数据变得平稳。在进行ARIMA模型拟合前，我们需要先通过画图或者ADF检验等方式，确定最小的d使得数据平稳。在确定了d之后，我们就可以将d阶差分后的序列代入模型进行拟合。

2.2 如何确定p,q值

1. 观测ACF图和PACF图确定p,q值

• 自相关系数（ACF）
  有序随机变量与其自身相比较。对于时间序列$x_t$，$x_t$与$x_{t-k}$的相关系数为$x_t$间隔k的自相关系数。

• 偏相关系数（PACF）
  PACF:剔除了其他变量，算出$x_{t-k}$对$x_t$的影响程度假设k=3，那么我们描述的是$y_t$和$y_{t-3}$之间的相关性
  对于ACF，这个相关性还受到$y_{t-1}$和$y_{t-2}$的影响。PACF剔除了这个影响，只考虑$y_{t-3}$对$y_t$的影响

用ACF图确认q值，用PACF图确认p值；
其中：对于ACF图先判断截尾or拖尾，再确定q值，q为最后一个超出二倍标准差的阶数；
对于PACF图先判断截尾or拖尾，再确定p值，p为最后一个超出二倍标准差的阶数

在ACF图和PACF图中，截尾和拖尾

1. 截尾（Cut-off）：如果ACF图或PACF图在某个滞后之后的所有系数都显著为0，那么我们说该图在该滞后处截尾。这意味着时间序列在该滞后之后没有显著的自相关性或偏自相关性。
2. 拖尾（Tail-off）：如果ACF图或PACF图的系数在某个滞后之后逐渐减小并趋向于0，那么我们说该图在该滞后之后拖尾。这意味着时间序列在该滞后之后的自相关性或偏自相关性逐渐减弱。

但是通过此方法求p,q的值，有一定的主观性，所以通过AIC或BIC来确定最优的p和q。

2. 通过AIC或BIC来选择最优的p和q。
赤池信息准则（Akaike Information Criterion，AIC） 和 贝叶斯信息准则（Bayesian Information Criterion，BIC），两个准则都是用于模型选择的，它们考虑了模型的复杂度和拟合优度。AIC和BIC越小，表示模型越好。

• 赤池信息准则（Akaike Information Criterion，AIC）
  公式：$AIC=2k-2ln(L)$
  k 是模型中估计参数的数量，L 是模型拟合的最大对数似然。

• 贝叶斯信息准则（Bayesian Information Criterion，BIC）
  公式：$BIC=ln(n)\ast k-2ln(L)$
  n 是观察的数据数量，k 是模型中估计参数的数量，L 是模型拟合的最大对数似然。

总结

传统的时间序列ARIMA模型对具有趋势特征的数据预测效果较好，但对具有随机特征的数据预测效果较差。，将ARIMA模型与聚类模型相结合，其预测精度明显高于单一ARIMA模型。