从0开始机器学习--4.4机器学习模型-回归（线性回归模型原理、假设、衍生，ARIMA（SARIMA），灰色预测模型，train_test_split 函数，含代码和场景实战）

从误差最小化角度：线性回归使用最小二乘法（见下）来最小化预测值与真实值之间的平方误差。均方误差是回归任务中常用的性能度量（具体可见本专栏《7.模型评价方式》）
从概率最大化角度：假设误差服从正态分布，线性回归通过最大化似然函数（见下）来选择最优参数，即最大化观测数据在模型下的出现概率。（更多用于分类任务，见本专栏《4.4分类》）

误差

误差是独立并且具有相同分布，并且服从均值为0方差为 $\theta^{2}$ 的高斯（正态）分布。

均方误差（MSE）： $E=\frac{1}{n}\sum_{i=1}^n(y_i-y_i^*)^2$
平均绝对误差（MAE）： $E=\frac{1}{n}\sum_{i=1}^n|y_i-y_i^*|$
均方根误差（RMSE）： $E=\sqrt{\frac{1}{n}\sum_{i=1}^n(y_i-y_i^*)^2}$

似然函数

似然函数是基于概率论的概念。在统计学中，似然函数衡量某一组参数下数据出现的可能性。在线性回归中，如果假设误差项（残差）服从正态分布，则似然函数表示给定模型参数下数据的观测概率。

对所有观测值 {y1,y2,…,yn}，联合似然函数是所有观测数据的条件概率（边缘概率）（概率论与数理统计基础知识）的乘积：

似然函数： $L(\theta)=\prod_{i=1}^n\frac{1}{\sqrt{2\pi\sigma^2}}\exp\left(-\frac{(y_i-\theta_0-\theta_1x_i)^2}{2\sigma^2}\right)$

计算累乘的联合概率时，由于大量的乘法计算，随着数据量增加，乘积很快会变得非常复杂，并且数值上可能出现溢出问题。为了简化计算，我们通常对联合概率取对数，将乘法转化为加法：

对数似然函数： $l(\theta)=\sum_{i=1}^n\left[-\frac{1}{2}\log(2\pi\sigma^2)-\frac{(y_i-\theta_0-\theta_1x_i)^2}{2\sigma^2}\right]$

最小二乘法（Ordinary Least Squares, OLS）

是线性回归中最常用的一种参数估计方法，其目标是通过找到一组回归系数，使得模型预测值与实际观察值之间的误差平方和最小化。

$OLS=min\sum_{i=1}^n(y_i-y_i^*)^2$

通过最小化这个目标函数，最小二乘法可以找到最佳的回归系数（权重和截距），使得模型的预测值与实际值尽可能接近。即：试图找到一条直线，使得样本到直线上的欧氏距离之和最小。

似然函数与最小二乘法的关系：

线性回归模型的最大似然估计与最小二乘法的目标一致，即最大化似然函数相当于最小化平方误差。对数似然函数的最大化与最小化误差平方和（即最小二乘法）的优化方向是一致的。
当假设误差服从正态分布时，最大化似然函数等价于最小化最小二乘误差。换句话说，最小二乘法是最大似然估计在正态分布假设下的特例。

最小二乘法依赖于几个经典的基本假设:

1. 线性关系假设

回归模型假设因变量𝑌与自变量𝑋之间存在线性关系。

2. 误差项期望为0假设

3. 误差项同方差性假设

误差项𝜖的方差应当是恒定的，不随自变量𝑋的变化而变化。

异方差性

误差项的方差随𝑋的变化而变化，则OLS估计量仍然无偏（见下），但不再是最有效的。（异方差性对估计结果的偏差没有直接影响，但会影响标准误差的计算，从而影响假设检验的可靠性）

异方差性检测

残差图：绘制残差与拟合值的散点图，如果散点图呈现漏斗形状或其他系统性的模式，可能存在异方差性。
Breusch-Pagan检验：通过检验残差的方差是否与自变量相关来检测异方差。

异方差性解决

加权最小二乘法（Weighted Least Squares, WLS）

加权最小二乘法是对普通最小二乘法的一种扩展，用于处理异方差问题，即数据点的误差不具有相同的方差。它通过给每个数据点分配不同的权重，使方差较大的点权重较小，方差较小的点权重大。核心思想是通过加权(w)来调整模型的拟合，使得拟合结果对低方差的数据点更敏感。权重通常与误差的方差成反比：

$WLS=min\sum_{i=1}^nw_i(y_i-y_i^*)^2$

广义最小二乘法（Generalized Least Squares, GLS）

广义最小二乘法用于处理自相关或异方差的情况。在普通最小二乘法中，假设误差项是独立同分布的（即具有相同的方差且互相不相关），而 GLS 放宽了这一假设，允许误差项之间存在相关性或不同方差。

无偏

估计量的期望值应该与真实值相等，在大量重复抽样的情况下，估计量不会有系统性偏差。

内生变量问题（Endogeneity Problem）

指的是在回归分析中，模型中的某个解释变量（自变量）与误差项之间存在相关性。换句话说，内生变量是指那些在模型中被解释但同时也受到其他因素影响的变量，这种影响会导致估计结果偏误（vs 异方差性-见上）。可能源于遗漏变量、测量误差、同时性。

内生变量问题检测

1. 统计检验：

Durbin-Wu-Hausman检验

Durbin-Wu-Hausman（DWH）检验是一种用于检验模型中内生性问题的方法，特别是用于判断模型中某个解释变量是否与误差项相关。如果解释变量与误差项相关，使用普通最小二乘法（OLS）的估计将会是有偏的，DWH 检验能够帮助判断是否需要使用工具变量法（IV）进行修正。如果检验结果表明解释变量是内生的，则应当使用工具变量法（IV）（见下）来修正偏差。如果检验结果表明解释变量是外生的，则普通最小二乘法（OLS）的估计是有效的，可以继续使用 OLS。

Wu-Hausman 检验

Wu-Hausman 检验是 DWH 检验的一种变体，但其核心思想一致。

2. 可视化分析：

通过绘制散点图观察自变量与因变量之间的关系。如果出现非随机分布，可能表明存在内生性问题。

内生变量问题解决

工具变量法（Instrumental Variable, IV）：通过引入一个或多个与内生变量相关但不与误差项相关的外生变量（仪器变量）来替代内生变量。使用两阶段最小二乘法（2SLS）进行估计。

第一阶段回归：用工具变量对有问题的解释变量进行回归，得到解释变量的估计值。
第二阶段回归：用第一阶段估计的解释变量去进行回归，得到回归模型的系数估计值。

内生变量问题关注的是解释变量（自变量）与误差项之间的相关性，而异方差性则涉及误差项本身的变异性。

4. 误差项的独立性假设

5. 误差项的正态性假设

6. 无多重共线性假设

多重共线性

指的是回归分析中自变量之间（vs 内生变量问题--见上）存在高度相关性，导致它们彼此之间存在显著的线性依赖性。简单来说，当一个自变量可以通过其他自变量的线性组合近似表示时，就存在多重共线性。

多重共线性问题

回归系数估计不稳定
模型解释困难
拟合效果看似良好

多重共线性检测

相关系数矩阵（Correlation Matrix）
方差膨胀因子（Variance Inflation Factor, VIF）
条件数（Condition Number）

多重共线性解决

删除高度相关的自变量/合并自变量
主成分分析（PCA）（在上一篇《4.2降维》中有过具体涉及）
岭回归（Ridge Regression）：正则化（会在《5.正则化技术》中具体涉及）
偏最小二乘回归（Partial Least Squares, PLS）（降维，类似PCA）
增加样本量

工作原理（步骤）

建立模型：假设线性模型的形式。
确定损失函数：通常使用最小二乘法。
求解参数：通过正规方程或梯度下降法（运筹学基础概念）估计模型参数。
评估模型：使用 R 方、MSE 等指标评估模型拟合效果。（会在《7.模型评估》中具体涉及）
预测：利用训练好的模型进行新数据的预测。
诊断模型：检查假设条件和模型问题。

代码

线性回归模型本身没有传统意义上的超参数。

以上都是线性回归背后的原理，看起来很复杂，但是实际运用起来是相对来说很轻松的，sklearn库帮我们封装好了一切背后的操作，我们只要调用它就可以了，以下是一个简单的例子：

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib import rcParams
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error

# 设置字体为 SimHei 以显示中文
rcParams['font.sans-serif'] = ['SimHei']
rcParams['axes.unicode_minus'] = False

# 假设有一个房屋价格数据集
# data = pd.read_csv('house_prices.csv')
# 创建一个简单的示例数据（无实际意义）
data = pd.DataFrame({'size': [800, 900, 1000, 1100, 1200], 'price': [200000, 230000, 260000, 300000, 320000]})

# 分割数据为训练集和测试集
house_X = data[['size']]
house_y = data['price']
X_train, X_test, y_train, y_test = train_test_split(house_X, house_y, test_size=0.2, random_state=42)

# 创建线性回归模型并训练
model = LinearRegression()
model.fit(X_train, y_train)
# 预测测试集
y_pred = model.predict(X_test)
# 计算均方误差
mse = mean_squared_error(y_test, y_pred)
print("线性回归均方误差:", mse)  # (predictd-real)^2
# 打印模型参数
print("线性回归斜率:", model.coef_)
print("线性回归截距:", model.intercept_)

# 对新的数据点进行预测
new_X = np.array([[1500]])
predicted_y = model.predict(new_X)
print("线性回归预测值:", predicted_y)

# 绘图部分
plt.figure(figsize=(8, 6))
# 绘制散点图 - 房屋大小与价格
plt.scatter(house_X, house_y, color='blue', label='数据点')
# 绘制拟合的回归线
plt.plot(house_X, model.predict(house_X), color='red', linewidth=2, label='回归线')
# 为新数据点预测值绘制点
plt.scatter(new_X, predicted_y, color='green', marker='x', s=100, label=f'预测值: {predicted_y[0]:.2f}')
# 添加图形标题与标签
plt.title('房屋大小与价格的线性回归模型')
plt.xlabel('房屋大小 (平方英尺)')
plt.ylabel('价格 (美元)')
plt.legend()
# 显示图形
plt.show()

输出：

线性回归均方误差: 2040816.3265305886
线性回归斜率: [308.57142857]
线性回归截距: -46285.71428571432
线性回归预测值: [416571.42857143]

train_test_split 函数

这里我们还要介绍和记住一个函数，这个函数在以后所有的代码中几乎都会出现，及其有用。

train_test_split 是 scikit-learn 库中的一个函数，用于将数据集按比例随机分成训练集和测试集。它的作用是为了在模型训练时，可以使用一部分数据来训练模型，而使用另一部分数据来评估模型的性能，避免过拟合或欠拟合的情况。

参数

下面详细介绍train_test_split(house_X, house_y, test_size=0.2, random_state=42) 中每个参数的作用和意义：

1. house_X 和 house_y

house_X：特征数据（输入变量），是自变量（独立变量），即我们用来进行预测的数据集。在房屋价格预测的示例中，house_X 是房屋的大小。
house_y：目标数据（输出变量），是因变量（依赖变量），即我们想要预测的值。在示例中，house_y 是房屋的价格。

2. test_size=0.2

test_size：这个参数控制测试集的比例。0.2 表示将 20% 的数据作为测试集，80% 的数据作为训练集。
- 训练集用于训练模型，帮助模型学习输入特征和输出之间的关系。
- 测试集用于评估模型的泛化能力，检查模型在未见过的数据上表现如何，帮助我们了解模型的实际预测性能。
比如数据集有 100 条数据，test_size=0.2 意味着 80 条数据会进入训练集，20 条数据会进入测试集。

3. random_state=42

random_state：这个参数控制数据拆分的随机性，使得每次执行时都能得到相同的分割结果。random_state=42 是一种常用的约定值，任何整数都可以，42 只是一个有名的“魔术数字”。
- 如果不设置 random_state，每次运行时分割的数据集都会不同。
- 设置 random_state=42 后，数据集划分的结果是可复现的，这对调试和比较模型效果时很有帮助。

4. X_train, X_test, y_train, y_test

X_train：训练集中的特征数据，用于模型的训练。
X_test：测试集中的特征数据，用于模型的评估。
y_train：训练集中的目标数据，与 X_train 对应，用于模型训练时的真实输出。
y_test：测试集中的目标数据，与 X_test 对应，用于模型评估时的真实输出。

函数的整体意义

train_test_split 的主要作用是将数据集分为两部分，分别用于模型的训练和测试。

训练集（Training Set）：用于拟合模型，帮助模型学习特征与目标之间的关系。我们希望模型在训练集上表现良好，但我们并不关心模型在训练集上的准确度是否极高，因为训练集准确度高可能意味着模型过拟合。
测试集（Testing Set）：用于评估模型的泛化能力，模拟模型在实际生产环境中遇到的全新数据的表现。如果模型在测试集上的表现良好，说明模型的泛化能力强，能较好地处理新的、未见过的数据。此处其实test集合为“验证集”-帮助模型训练得到反馈预测误差的集合。

为什么要分离数据集？

防止过拟合：如果我们用同一个数据集来训练和评估模型，模型可能会过拟合，也就是它会学习到数据中的噪音和细节，而非真正的规律。这样会导致模型在新数据上表现不佳。分离数据集能确保模型在没有见过的数据上也能表现好。
模型评估：分离出测试集能帮助我们独立评估模型的泛化能力。通过比较模型在训练集和测试集上的表现，我们可以判断模型是否存在过拟合或欠拟合问题。

因此，train_test_split 函数为我们提供了一个简单、有效的方法，将数据集分为训练集和测试集，确保我们可以公正地评估模型的实际性能。

对数线性回归

线性模型虽简单，却有丰富的变化。可否逼近y的衍生物？
例如：输出标记在指数尺度上变化，那就可将输出标记的对数作为线性模型逼近的目标，即:

lny= $W^{T}x+b$

在形式上仍然是线性回归，但实质上是在求取输入空间到输出空间的非线性函数映射。

更一般的，可以考虑单调可微函数g(.)--广义线性模型。

多元非线性回归（多项式）

原理

多项式回归是线性回归的一种扩展，它可以对非线性关系建模。通过将特征转换为多项式形式，它允许模型在输入数据中捕捉更复杂的关系。

目标函数

多项式回归的目标是最小化预测值与真实值之间的误差，常用的误差衡量方式是最小二乘法。

假设

数据之间的关系可以通过多项式函数表示。
数据点是独立的。

公式

假设我们有特征 x，其多项式回归的形式为： $y = \beta_0 + \beta_1 x + \beta_2 x^2 + \beta_3 x^3 + \cdots + \beta_n x^n$

其中， $\beta_0, \beta_1, \ldots, \beta_n$ 为多项式系数。

常用场景

适用于数据表现出非线性关系的情况，如经济数据预测、环境数据建模等。

代码

超参数

多项式的阶数。这里要注意：

阶数太高：模型容易过拟合（会在下一篇《4.4分类》中具体涉及），复杂度增加。
阶数太低：模型可能欠拟合，无法捕捉数据中的非线性关系。
该超参数通常通过交叉验证（见《2.ai模型概念+基础》）或经验判断来确定。

import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import PolynomialFeatures
# 准备数据
poly_X = np.array([24,25,26,27,28,29,30])
y = np.array([25,26,20,20,19,23,17])
# 转换特征矩阵为多项式特征
poly = PolynomialFeatures(degree=8) #degree=8时，在构建多项式特征时会生成包含输入特征的所有组合，包括一次、二次和三次方特征，以及它们的交叉项。
#degree参数是可以根据需求进行自行修改的。通过调整degree的值，你可以控制生成多项式特征的阶数，从而改变拟合曲线的灵活性和复杂性。
X_poly = poly.fit_transform(poly_X[:, np.newaxis])
# 创建线性回归模型，并拟合多项式特征
model = LinearRegression()
model.fit(X_poly, y)
y_pred = model.predict(X_poly)
# 计算准确率
# 计算决定系数
r2 = model.score(X_poly, y)
print("多项式回归degree=8决定系数(score):", r2)
# 绘制原始数据与拟合曲线
plt.scatter(poly_X, y, color='blue', label='Data')
x_range = np.linspace(poly_X.min(), poly_X.max(), 100)
x_poly_range = poly.transform(x_range[:, np.newaxis])
y_pred = model.predict(x_poly_range)
plt.plot(x_range, y_pred, color='red', label='Polynomial Regression')
plt.xlabel('X:date')
plt.ylabel('y:highest_temperture')
plt.legend()
plt.show()
print( )

输出：

多项式回归degree=8决定系数(score): 0.9999999566380657

决定系数等其他评价指标是什么会在第七篇《评估》中具体涉及。

ARIMA

ARIMA（AutoRegressive Integrated Moving Average，自回归积分滑动平均模型）模型是一种广泛使用的时间序列分析方法，它可以用于对未来的数据进行预测。特别适用于那些具有一定趋势或季节性但又相对平稳的数据。

组成

ARIMA模型由自回归模型（AR模型）、差分整合模型（I模型）和移动平均模型（MA模型）组成，因此也被称为ARIMA(p,d,q)模型。

p: 自回归（AR）部分的阶数，表示模型中滞后观测值的数量。
d: 差分阶数，表示数据需要几次差分操作才能达到平稳。
q: 移动平均（MA）部分的阶数，表示模型中滞后预测误差项的数量。

SARIMA模型的参数比ARIMA还多了seasonal_order 参数，通常写作 (P, D, Q, s)：

P：季节性自回归（SAR）阶数，用于捕捉季节性周期中自回归的部分。
D：季节性差分次数，为了去除季节性趋势而进行的差分次数。
Q：季节性移动平均（SMA）阶数，用于捕捉季节性周期中移动平均的部分。
s：季节性周期的长度。例如，月度数据可能有12个月的季节周期，s 常取12；如果是季度数据且一年4个季度，s 常取4。

参数选择

差分操作和自相关函数 (ACF)、偏自相关函数 (PACF) 是非常重要的工具，帮助我们识别数据的趋势性、平稳性以及选择模型的合适阶数。

一阶差分 (First-order Difference)

一阶差分是将当前时间点的值与前一个时间点的值做差，用来消除数据中的线性趋势，使时间序列数据变得更加平稳。

作用：一阶差分通常用于去除时间序列中的趋势成分（即长期的线性变化），帮助数据变得平稳。特别是在时间序列分析或建模（如PCA或ARIMA）中，如果数据存在上升或下降的趋势，应用一阶差分可以消除这种趋势，确保数据更适合平稳模型。

二阶差分 (Second-order Difference)

二阶差分是对一阶差分后的数据再次做差，主要用于处理更加复杂的趋势，比如非线性或二次趋势。

作用：二阶差分通常用于消除数据中的非线性趋势或加速变化的趋势。如果一阶差分不足以使数据变平稳，应用二阶差分可以进一步消除更复杂的趋势。

自相关函数 (ACF, AutoCorrelation Function)

自相关函数 (ACF) 是用来衡量一个时间序列与其滞后值之间的相关性。它展示了当前值与其之前多个时间步的值的线性关系。

作用：ACF图可以帮助识别时间序列中存在的依赖性或周期性。如果 ACF 图在特定的滞后处显著非零（即远离零线），说明该滞后处存在相关性。ACF 还可以用来判断是否需要做差分处理以及估计 ARIMA 模型中的 MA（滑动平均）部分的阶数。

偏自相关函数 (PACF, Partial AutoCorrelation Function)

偏自相关函数 (PACF) 是在考虑中间滞后值的影响后，衡量当前值与滞后 kkk 个时间步值之间的相关性。它排除了较小滞后的影响，直接反映当前值与滞后值之间的关系。

作用：PACF图用于识别AR部分的阶数。在AR模型中，PACF会在特定的滞后阶数后急剧减小，表明该滞后阶数之后的自相关性不再显著。PACF特别适用于选择ARIMA模型中的AR（自回归）部分的阶数。

ARIMA模型的应用

ARIMA 模型在金融市场、经济分析、销售预测等领域应用广泛，适合处理趋势性、周期性或季节性较弱的时间序列。
例如，ARIMA可以用于预测股票价格、宏观经济指标（如GDP增长率）以及商品销售数据等。

扩展模型

SARIMA (Seasonal ARIMA)：在标准ARIMA的基础上，增加季节性成分，用于处理季节性时间序列。

from statsmodels.tsa.statespace.sarimax import SARIMAX
import statsmodels.api as sm
from statsmodels.tsa.stattools import adfuller

# 执行ADF检验
adf_test = adfuller(train_data)
# 打印检验结果
print(f'ADF Statistic: {adf_test[0]}')
print(f'p-value: {adf_test[1]}')

# 差分并确定参数d
d = 0
while adfuller(train_data)[1] > 0.05:
    train_data = train_data[1:] - train_data[:-1]  # 一阶差分
    d += 1
print(f"进行第 {d} 次差分，数据已经平稳")
    
adf_test = adfuller(train_data)
print(f'ADF Statistic: {adf_test[0]}')
print(f'p-value: {adf_test[1]}')

# 使用 AIC 和 BIC 选择最佳参数 (p, d, q)
train_result = sm.tsa.arma_order_select_ic(monthly_reports, ic=['aic', 'bic'], trend='c', max_ar=4, max_ma=4)
print('AIC 最优参数:', train_result.aic_min_order)
print('BIC 最优参数:', train_result.bic_min_order)

# 使用SARIMA模型预测未来的传播数量
seasonal_order = (1, 1, 0, 12)  # 设置季节性部分的参数
model = SARIMAX(monthly_reports, order=(6,1,0), seasonal_order=seasonal_order)
model_fit = model.fit()

ARIMAX：在ARIMA模型中加入外生变量，用于考虑外部影响因素。

优缺点

优点：对平稳性较好的时间序列，ARIMA可以捕捉时间序列的自相关结构，预测效果较好。
缺点：对季节性强的数据，普通ARIMA模型可能不适用，需扩展为SARIMA；另外，ARIMA模型假设时间序列是线性的，因此对非线性关系的捕捉较弱。

模型建立过程

平稳性检测：首先对时间序列进行平稳性分析，常用方法如ADF（Augmented Dickey-Fuller）检验。如果数据不平稳，需要通过差分（d）处理使其平稳。
参数选择：确定 ARIMA(p, d, q) 中的 p 和 q 的值，通常使用自相关函数（ACF）和偏自相关函数（PACF）图来帮助选择合适的阶数。

使用 AIC（Akaike 信息准则）或 BIC（贝叶斯信息准则）等指标来选择最优的 (p,d,q) 及 (P,D,Q,s)参数组合。
模型估计：基于已确定的参数，使用最大似然估计或其他方法对模型进行拟合。
模型诊断：检查残差是否为白噪声，是否满足独立性假设。如果残差仍然存在结构性信息，则需要调整模型。
预测：使用拟合好的模型进行未来值的预测。

代码

这是一段预测未来股票价格的ARIMA代码：（源代码参考自Python3实现基于ARIMA模型来预测茅台股票价格趋势_股票预测python-CSDN博客）

import yfinance as yf #Python和相关的金融数据API
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
import statsmodels.api as sm
import warnings
warnings.filterwarnings('ignore')
sns.set_theme(font='SimHei')
plt.rcParams['font.sans-serif'] = ['SimHei'] #解决中文显示
plt.rcParams['axes.unicode_minus'] = False   #解决符号无法显示

# 设置股票代码和时间范围
ticker = "600519.SS"
start_date = "1990-01-01"
end_date = "2024-05-08"

# # 下载股票数据
# stock_data = yf.download(ticker, start=start_date, end=end_date)#雅虎财经下载股票的历史数据。
# # 将数据保存到指定路径的 CSV 文件中
# save_path = r"文件路径.csv"
# stock_data.to_csv(save_path)
# print(stock_data.head())

stock_file=r"文件路径.csv"

# 数据预处理
    # 重点分析收盘价并预测，对原始数据进行重采样，以周且指定周一为单位求平均值
    # 导入数据集并将其转换为时间序列
df = pd.read_csv(stock_file, index_col='Date', parse_dates=True)
stock_week = df['Close'].resample('W-MON').mean()
# 取出2001-1月到2024-2月的数据作为训练数据
stock_train = stock_week['2001-9':'2024-2']
# 做出折线图
stock_train.plot(figsize=(12,7))
plt.legend()
plt.title('Stock Close')
sns.despine() #去除图表边框，更美观
plt.show()

#差分并确定参数d
    # 对数据进行拆分的目的是保证数据的平稳性，
    # 因为通过上图我们发现原始数据波动的幅度很大，需要进行拆分操作。
    #1. 将时间序列进行差分，直到其成为平稳序列
ts = df['Close']
d = 0
while not sm.tsa.stattools.adfuller(ts)[1] < 0.05:
    ts = ts.diff().dropna() #dropna()：去除差分后产生的 NaN 值。
    d += 1
print('参数d为：',d)

    #2. 将时间序列进行差分并确定参数d
    # 一阶差分
stock_diff_1 = stock_train.diff()
stock_diff_1.dropna(inplace=True)
    # 二阶差分
stock_diff_2 = stock_diff_1.diff()
stock_diff_2.dropna(inplace=True)
plt.figure(figsize=(12,6))
plt.subplot(2,1,1) #分别创建两行子图。
plt.plot(stock_diff_1)
plt.title('一阶差分')
plt.subplot(2,1,2) #分别创建两行子图。
plt.plot(stock_diff_2)
plt.title('二阶差分')
plt.show()

#确定参数q和p
    #1. 做出ACF、PACF图确定
    # 做出ACF（自相关函数）图确定参数q
    #ACF图显示了数据点与其滞后值之间的相关性。可以判断数据的滞后效应。
sm.graphics.tsa.plot_acf(stock_diff_1)
plt.title('ACF')
plt.show()
    # 做出PACF（偏自相关函数）图并确定参数p
    #显示了在考虑其他滞后值影响的情况下数据点与其特定滞后值之间的相关性。对确定自回归AR阶数很有用。
sm.graphics.tsa.plot_pacf(stock_diff_1)
plt.title('PACF')
plt.show()
    #2. 使用BIC和AIC的值确定q和p
    # AIC（Akaike 信息准则）: 衡量模型质量的准则，考虑了模型的复杂度。较小的 AIC 值表示模型更好。
    # BIC（贝叶斯信息准则）: 类似于 AIC，但对模型复杂度的惩罚更严格。较小的 BIC 值表示模型更好。
train_result = sm.tsa.arma_order_select_ic(stock_train,ic=['aic','bic'],trend='c',max_ar=4,max_ma=4)
print('AIC',train_result.aic_min_order)
print('BIC',train_result.bic_min_order)
#若BIC和AIC的值不一样，两个结果都试看哪个参数组合训练的模型效果最好。这里取p=q=1的参数结果。

#训练模型并预测
# 拟合ARIMA模型
model = sm.tsa.ARIMA(stock_train, order=(4, 1, 4),freq='W-MON')
    #order（p,d,q）,将前面确定数值填进去即可，freq是为了和前面重采样保持一致。
result = model.fit()
# 使用该模型进行预测
forecast = result.predict(start='2022-01-03', end='2025-02-01')
#预测的时候需要填写起始时间和终止时间，起始时间必须在训练数据中出现

#预测的结果和真实值可视化
plt.figure(figsize=(12,6))
plt.xticks(rotation=45)
plt.plot(forecast,label='预测值')
plt.plot(stock_train,label='真实值')
plt.legend()
plt.show()

输出：

参数d为： 1
AIC (3, 4)
BIC (3, 3)

ACF 和 PACF 图的解读

ACF 图：通常用于确定 MA 部分的阶数。ACF 图的缓慢衰减或拖尾现象表明数据有 MA 结构，显著的滞后值可能暗示需要使用特定阶数的 MA 模型。
PACF 图：用于确定 AR 部分的阶数。如果 PACF 在某个滞后阶数之后突然接近零，这通常意味着可以用该滞后阶数作为 AR 模型的阶数。

在运行结果的最后一幅图中我们可以看到，预测值的最后是一条水平的直线，在这里是正确的，就是模型最终拟合认为未来会不变。下面介绍一种另外的方法以打破这种明显的错误结果：

滚动预测（Rolling Forecasting）

是一种时间序列分析方法，旨在通过不断更新模型的输入数据来进行短期的递归预测。每次预测完一个时间点后，模型会将最新的预测结果作为下一次预测的输入，以便进行后续时间点的预测。这种方法可以提高预测的准确性，尤其适用于非平稳时间序列或长期预测。

滚动预测的步骤

模型初始化：基于历史数据训练模型，通常是用一段较早的时间序列数据作为训练集。
单步预测：使用训练好的模型预测未来某个时间点的值（例如，下一个时间点）。
更新模型：将新预测的结果加入原始数据中，重新拟合模型，以反映最新的信息。
重复预测：递归地对下一个时间点进行预测，直到达到预定的预测时间段。

滚动预测的优点

适应性强：随着新数据的加入，模型会不断更新，从而捕捉时间序列中最新的趋势和波动。
提高长期预测精度：递归预测虽然每次只预测一步，但可以通过滚动过程实现多步预测，相较于直接长时间段的预测更为准确。
灵活性：可以根据时间序列的特性，动态调整模型的参数或结构。

滚动预测的实例

假设我们用 ARIMA 模型对股票价格进行预测，通过滚动预测，我们每次用模型预测下周的收盘价，并将预测结果作为下一周预测的基础。如此反复操作，逐步预测多个时间点的未来值。

滚动预测的适用场景

金融时间序列预测：如股票价格、汇率等，市场波动性强且具有不可预见性。
销售量预测：用来预测产品销售量或需求，尤其在长期趋势变化较快的行业。
经济指标预测：如通货膨胀率、GDP 增长等，需根据最新数据进行修正。

滚动预测的注意事项

模型误差累积：递归预测存在误差传播问题，前期预测误差可能会影响后续预测的准确性。因此，对于长时间的预测需要特别注意误差的控制。
模型更新频率：在进行滚动预测时，模型并不需要每次都重新训练，可以通过对部分数据进行调整来提高效率。

灰色预测模型 GM(1,1)

灰色预测模型（Grey Prediction Model），特别是灰色预测模型 GM(1,1)，是一种在小样本或不确定性较高的情况下进行时间序列预测的数学模型。它是灰色系统理论的一部分，主要用于处理具有部分已知信息和部分未知信息的系统。通常可以用于填补缺失值。

"1,1" 表示该模型是一阶单变量模型（即处理单个变量的一阶微分方程）。该模型适用于具有一定规律性的时间序列数据，尤其在样本量较小的情况下仍能有效进行预测。

优势

适用于小样本、短时间序列数据。
不需要数据的严格分布假设。被广泛应用于各个领域，特别是在数据稀缺或不完全的情况下。
计算简单，容易实现。

局限

对数据的光滑性要求较高，噪声大时预测效果较差。
仅适用于具有一定单调性或近似指数规律的时间序列。

步骤

原始数据的累加生成序列：生成一个新的序列，使其波动变小、平稳。
建立微分方程：根据累加序列建立一阶微分方程，描述序列的发展规律。
参数估计：利用最小二乘法估计模型参数。
预测：根据拟合的模型方程进行未来数据点的预测。
还原：对累加序列进行还原，得到原始序列的预测值。

例

使用灰色预测模型 GM(1,1) 对给定的时间序列进行未来数据点的预测。它预测的是该时间序列未来的值，在本例中，预测的是 x0 数组的延续趋势。

代码

import numpy as np
import matplotlib.pyplot as plt

plt.rcParams['font.sans-serif'] = ['SimHei'] #解决中文显示
plt.rcParams['axes.unicode_minus'] = False   #解决符号无法显示

# 原始数据序列
x0 = np.array([322, 299, 311, 328, 357, 392])  # 示例数据

# 1. 累加生成序列
x1 = np.cumsum(x0)

# 2. 构造数据矩阵B和Y
B = np.array([-0.5 * (x1[i] + x1[i+1]) for i in range(len(x0)-1)])
B = np.vstack((B, np.ones(len(B)))).T
Y = x0[1:]

# 3. 计算参数a和b
ab = np.linalg.inv(B.T @ B) @ B.T @ Y
a, b = ab

# 4. 预测模型x^(1)
def gm1_model(a, b, x0, k):
    return (x0[0] - b/a) * np.exp(-a * k) + b/a

# 5. 还原x^(0)预测值
pred_x1 = gm1_model(a, b, x0, np.arange(1, len(x0) + 1))
pred_x0 = np.diff(np.hstack(([x0[0]], pred_x1)))

# 6. 可视化结果
plt.plot(np.arange(1, len(x0) + 1), x0, 'o-', label='真实值')
plt.plot(np.arange(1, len(pred_x0) + 1), pred_x0, 's-', label='预测值')
plt.xlabel('时间')
plt.ylabel('数值')
plt.title('GM(1,1)预测模型')
plt.legend()
plt.show()