期待今晚的莎莎和大头😎！
🏆🏆🏆
愿巴黎登顶💜

正文开始
主要是跟着实战：时间序列模型(五)：时间序列案例_实现销售额预测

💜本人是我们组的编程手，本博客只为记录自己学习的过程，大家可以参考学习！
💜我用的数据中下面的截图日期一列应加入/符号分开年月日（2022/03/15），不然出错。

一.前期准备工作

• 为了更好的理解ARIMA时间序列预测模型的相关定义和建模过程，可以先移步以下四篇文章：
  时间序列模型(一)：时序的相关概念
  时间序列模型(二)：AR模型
  时间序列模型(三)：MA模型
  时间序列模型(四)：ARIMA模型
• 获取文件地址：
  我一般保存在桌面，右击-复制文件地址即可，（或选中使用快捷键，Ctrl+Shift+C）
  如何把Excel的数据导入python?

二.数据预处理

2.1去除重复项

函数pandas.DataFrame.drop_duplicates(subset=None, keep=‘first’, inplace=False, ignore_index= False)主要用来去除重复项，返回DataFrame类型的数据。

subset：默认为None
去除重复项时要考虑的标签，当subset=None时所有标签都相同才认为是重复项

keep： {‘first’, ‘last’, False}，默认为‘first’。
keep=‘first’时保留该标签第一次出现时的样本，之后重复出现的全部丢弃。
keep=‘last’表示保留该标签最后一次出现的样本。
keep=False时该标签重复的样本全部丢弃。

inplace:bool，默认为False
inplace=False时返回去除重复项后的DataFrame,原来的DataFrame不改变。
inplace=True时返回空值，原来DataFrame被改变

ignore_index:bool，默认为False
ignore_index=False时，丢弃重复值之后的DataFrame的index不改变
ignore_index=True时，丢弃重复值之后的DataFrame的index重新变为0, 1, …, n-1

这是我的数据

#删除重复项
import pandas as pd
import xlsxwriter
import xlrd
io = r'D:\Users\25705\Desktop\乌方每半月兵力损失.xlsx'
df = pd.read_excel(io,sheet_name=0)
df.head()
df_deduplicated = df.drop_duplicates(subset=None,inplace=False,ignore_index=True)  # 去除重复项  subset='日期'
output_file_path = 'D:\\Users\\25705\\Desktop\\去重后数据.xlsx'  # 替换为你想要保存的文件路径
df_deduplicated.to_excel(output_file_path, index=True)  # index=True表示保存行索引

得到一个新的表格

2.2处理缺失值

2.2.1检查缺失值

首先需要检查是否存在缺失值。如果存在，需要根据具体情况来处理，比如填补缺失值（使用前一个月的数据，或者使用相邻月份的平均值），或者删除含有缺失值的行。
一般采取填充的方法处理缺失值。时序数据前后相互关联，不建议使用定值填充（作为区分，非时序数据经常用类似-9999这样的值来标识缺失值）。

# 数据处理(检查是否存在缺失值)
import numpy as np
import pandas as pd
import matplotlib.pyplot as pltio = r'D:\\Users\\25705\\Desktop\\乌方每半月兵力损失.xlsx'
df = pd.read_excel(io,sheet_name=0)
df.head()
print(df .isnull().sum())    # 返回每一列中的缺失值（NaN）的数量
print(df .isnull().any(axis=1))  # 返回一个布尔序列，标识每一行是否含有缺失值
# 查看具体某一行的值，可以通过索引或者iloc属性来实现
print(df .iloc[14])

😇这里的十五列指的是数据的第十五列，不包括第一行的文字。

2.2.2邻值填充

邻值填充就是用最近时刻的值来填充，可向前或向后填充。使用fillna()方法，指定method='ffill’进行前向填充，method='bfill’进行后向填充。

df_ffill = df.fillna(method='ffill')  # 前向填充

这是处理之前的数据：

import numpy as np
import pandas as pd
import matplotlib.pyplot as pltio = r'D:\\Users\\25705\\Desktop\\乌方每半月兵力损失.xlsx'
df = pd.read_excel(io,sheet_name=0)
df.head()
df_ffill = df.fillna(method='ffill')  # 前向填充
output_file_path = 'D:\\Users\\25705\\Desktop\\填充后数据.xlsx'  # 替换为你想要保存的文件路径
df_ffill.to_excel(output_file_path, index=True)  # index=True表示保存行索引

执行代码后生成一个新的excel表格：

2.2.3插值填充

针对时序数据的平滑性，使用插值方法填充更加接近真实数据。最常用的是DataFrame.interpolate()线性插值。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import interpolate
file_location = "D:\\Users\\25705\\Desktop\\乌方每半月兵力损失.xlsx"# 使用 pandas 读取 Excel 文件
df = pd.read_excel(file_location, sheet_name=0)  # sheet_name=0 表示读取第一个 sheet# 假设第一列（索引为0）包含日期
# 这里我们假设要插值的列是第二列（索引为1），你需要根据实际情况调整
# 首先，确保日期列是日期类型
df['日期'] = pd.to_datetime(df.iloc[:, 0], errors='coerce')  # 假设第一列是日期# 设置日期列为索引
df.set_index('日期', inplace=True)# 检查数据排序
df = df.sort_index()
df.head()
# 对特定列进行插值（这里以第三列为例）
# 假设列名为 '每半月兵力损失' 或你需要插值的列的实际名称
df['Losses_interpolated'] = df['每半月兵力损失'].interpolate(method='linear')# 如果不需要原始数据，只保存插值后的数据
df_interpolated = df[['Losses_interpolated']]  # 只保留插值后的列# 保存插值后的数据到新的 Excel 文件
output_file_path = 'D:\\Users\\25705\\Desktop\\填充后数据.xlsx'
df_interpolated.to_excel(output_file_path)  # 默认情况下，pandas 会保存索引
# 如果你不想保存索引，可以添加参数 index=False
# df_interpolated.to_excel(output_file_path, index=False)

运行代码前：

运行代码后得到的新excel:

2.3数据类型转换

确保’days’列是日期类型，这对于时间序列分析是必要的。如果’days’列当前不是日期类型，可以使用pandas的to_datetime函数进行转换。

在2.2.3插值填充代码中有涉及

2.4设定索引

对于时间序列分析，需要将时间列（在这个案例中是’days’列）设为索引。

在2.2.3插值填充代码中有涉及

2.5检查数据排序

确认数据是按照时间顺序排列的。如果不是，需要进行排序。

在2.2.3插值填充代码中有涉及

2.6异常值处理

参考：时序数据预处理
时间序列预测中的4大类8种异常值检测方法(从根源上提高预测精度)
异常值的处理一般分为两步：识别异常值和处理异常值。
处理异常值相对来说比较简单：直接删除，或者用处理缺失值的方法填充。DBSCAN将数据分为核心点、边界点和离群点。
核心点：在一定范围的邻域内有一定数量相邻的点。
边界点：一个数据不是核心点，但他在核心点的邻域内。
离群点：以上两者之外的点。

在这里插入代码片

2.6.1Z-score检测

下面部分代码是用Z-score检测出异常值同时用移动平均的思想将异常值进行了替换，从而保证我们数据的正常，不受那些因为异常事件导致的特殊值影响。

# 导入所需模块
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from scipy import stats# 用来正常显示中文标签
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['figure.figsize'] = (10.0, 8.0)  # set default size of plots
plt.rcParams['image.interpolation'] = 'nearest'
plt.rcParams['image.cmap'] = 'gray'# 读取数据
file_path = "D:\\Users\\25705\\Desktop\\乌方每半月兵力损失.xlsx"  # 替换为您的文件路径
data = pd.read_excel(file_path, sheet_name=0)  # sheet_name=0 表示读取第一个 sheet# 计算每半月兵力损失列的Z分数
data['每半月兵力损失_Z_Score'] = stats.zscore(data['每半月兵力损失'])# 定义离群点（这里以Z分数的绝对值大于2为标准）   自己试（后面新产生的表格会有每个数据对应的Z值）
data['每半月兵力损失_Outlier_Z'] = data['每半月兵力损失_Z_Score'].abs() >0.8 #绝对值大于0.8# 创建折线图，并在其中标注离群点
plt.figure(figsize=(12, 6))
plt.style.use('ggplot')# 绘制forecast值的折线图
sns.lineplot(data=data, x=data.index, y='每半月兵力损失', label='乌方每半月兵力损失', color='blue')# 标注离群点
outliers = data[data['每半月兵力损失_Outlier_Z']]
plt.scatter(outliers.index, outliers['每半月兵力损失'], color='red', label='Outliers (Z > 0.8)')# 设置标题和标签
plt.title('每半月兵力损失 Value Over Time with Z-Score > 0.8 Outliers', fontsize=15)
plt.xlabel('Time', fontsize=12)
plt.ylabel('每半月兵力损失', fontsize=12)
plt.legend()# 显示图表
plt.show()# 返回Z分数大于0.8的离群点的数量
outliers_count = outliers['每半月兵力损失'].count()
print(outliers_count)# 计算 moving_avg
data['moving_avg'] = data['每半月兵力损失'].rolling(window=5, min_periods=1).mean()# 计算 Z-Score
data['z_score'] = (data['每半月兵力损失'] - data['每半月兵力损失'].mean()) / data['每半月兵力损失'].std()# 将异常值替换为移动平均值
data.loc[data['z_score'].abs() > 0.8, 'Volu每半月兵力损失e'] = data['moving_avg']# 导出新的excel表格
output_file_path = 'D:\\Users\\25705\\Desktop\\异常值处理数据.xlsx'
data.to_excel(output_file_path)  # 默认情况下，pandas 会保存索引

得到的结果：

第一行的值没有变，不知道什么原因。

Z-Score方法在处理具有高斯（正态）分布特征的数据时特别有效。它可以识别那些与大多数数据显著不同的点。然而，对于不是正态分布的数据集，这种方法可能不太适用，因为它依赖于平均值和标准差。

😎至此，我们完成了数据预处理，现在开始对数据进行分析。

三.数据分析

3.1探索性数据分析（EDA）

首先，我们可以对数据进行初步的分析。可以绘制时间序列图来观察是否有明显的趋势、季节性或者异常点。

import matplotlib.pyplot as plt
import pandas as pd
plt.rcParams['font.sans-serif'] = ['KaiTi', 'SimHei', 'FangSong']  # 汉字字体,优先使用楷体，如果找不到楷体，则使用黑体
plt.rcParams['font.size'] = 12  # 字体大小
plt.rcParams['axes.unicode_minus'] = False  # 正常显示负号#创建一个新的图形，并设置了其宽度10英寸，高6英寸
plt.figure(figsize=(10,6))
file_location ="D:\\Users\\25705\\Desktop\\乌方每半月兵力损失.xlsx"
df = pd.read_excel(file_location, sheet_name=0)
plt.plot(df['每半月兵力损失'])
plt.xlabel('日期')
plt.ylabel('每半月兵力损失')
plt.title('乌方每半月兵力损失')
plt.show()

从图像来看，数据基本上是无法满足平稳性条件的，接下来我们在平稳性检验(ADF单位根检验)上验证这一结论。

3.2ADF单位根检验

#代码运行成功的前提是进行了数据的清洗
# 使用差分消除数据波动
'''
在时间序列中，标签往往具备一定的周期性：例如，标签可能随季节有规律地波动
（比如在夏季标签值高、在冬季标签值较低等），也可能随一周的时间有规律地波动
（比如在周末较高、在工作日较低等）。这种波动可以通过滞后差分来消除，
'''
# 导入必要的库
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.stattools import adfuller
import xlrd# 创建一个函数来检查数据的平稳性
def test_stationarity(timeseries):# 执行Dickey-Fuller测试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] = valueprint(dfoutput)#自己导入数据--不平稳的时间序列
file_location = "D:\\Users\\25705\\Desktop\\乌方每半月兵力损失.xlsx"#避免被转义使用\\
data = xlrd.open_workbook(file_location)#	data是Excel里的数据
sheet = data.sheet_by_index(0)  #根据索引读取,就是说data里有很多sheet	data.sheet_by_index（0）是sheet1，以此类推
days = [sheet.cell_value(r,1) for r in range(1,sheet.nrows)]# 读取该表格中第2列中的数据，（r,n）表示第（n+1）列， range(1,sheet.nrows)从第2行开始
x = list(days)# days为generator，数据转换成list
# 表格从左上角开始是（0，0）和二维数组类似# 把它转换成Pandas的DataFrame格式
df = pd.DataFrame(x, columns=['value'])# 检查原始数据的平稳性
test_stationarity(df['value'])# 进行一阶差分
df['first_difference'] = df['value'] - df['value'].shift(1)# 检查一阶差分后的数据的平稳性
test_stationarity(df['first_difference'].dropna())# 进行二阶差分
df['second_difference'] = df['first_difference'] - df['first_difference'].shift(1)# 检查二阶差分后的数据的平稳性
test_stationarity(df['second_difference'].dropna())# 可视化原始数据和差分后的数据
plt.figure(figsize=(12, 6))
plt.plot(df['value'], label='Original')
plt.plot(df['first_difference'], label='1st Order Difference')
plt.plot(df['second_difference'], label='2nd Order Difference')
plt.legend(loc='best')
plt.title('Original and Differenced Time Series')
plt.show()

查看运行结果（ADF检验）：

Results of Dickey-Fuller Test:
Test Statistic -0.907420
p-value 0.785519
#Lags Used 3.000000
Number of Observations Used 21.000000
Critical Value (1%) -3.788386
Critical Value (5%) -3.013098
Critical Value (10%) -2.646397
dtype: float64
Results of Dickey-Fuller Test:
Test Statistic -6.574510e+00
p-value 7.786070e-09
#Lags Used 1.000000e+00
Number of Observations Used 2.200000e+01
Critical Value (1%) -3.769733e+00
Critical Value (5%) -3.005426e+00
Critical Value (10%) -2.642501e+00
dtype: float64
Results of Dickey-Fuller Test:
Test Statistic -7.482329e+00
p-value 4.736406e-11
#Lags Used 1.000000e+00
Number of Observations Used 2.100000e+01
Critical Value (1%) -3.788386e+00
Critical Value (5%) -3.013098e+00
Critical Value (10%) -2.646397e+00
dtype: float64

原始数据：
Test Statistic -0.907420是adt检验的结果，简称为T值，表示t统计量。
p-value 0.785519简称为p值，表示t统计量对应的概率值。
Lags Used 3.000000表示延迟。
Number of Observations Used 21.000000表示测试的次数。
第五个是配合T值一起看的，是在99%，95%，90%置信区间下的临界的ADF检验的值。
首先，Test Statistic -0.907420大于（若小于，则不存在单位根，不存在单位根即是平稳序列）三个置信区间的临界值，即存在单位根。
其次，p值大于（若小于，则不存在单位根，不存在单 位根即是平稳序列，等于0最好）给定的显著水平（一般是0.05）。本数据中，P值为 0.785519,大于0.05，即存在单位根。
所以综上所述，此序列不是平稳序列。
一阶差分数据和二阶差分数据按如上分析即可。

可以参考下面这篇文章：
时间序列分析

出现报错：
ValueError: could not convert string to float: ‘’
原因是：excel中缺失了值

3.3差分

我们使用pandas的diff()函数进行差分。
通过遍历不同的差分步数，并对每个步数的结果进行Augmented Dickey-Fuller (ADF)检验。函数将打印出每个阶数的ADF检验统计量、p值、标准差和方程，帮助选择最佳的差分步数。代码如下：

# 首先导入了所需要的库和函数
from statsmodels.tsa.stattools import adfuller
import pandas as pdio = r'D:\Users\25705\Desktop\乌方每半月兵力损失.xlsx'
df = pd.read_excel(io,sheet_name=0)
df.head()
# 定义一个名为calculate_diff的函数，接收三个参数：df是要处理的DataFrame，max_diff是最大的差分步数，significance_level是判定平稳性的显著性水平
def calculate_diff(df, max_diff, significance_level=0.05):# 初始化最佳差分阶数和最小p值best_diff = Nonemin_pvalue = 1.0min_variance = float('inf')  # 初始化最小方差min_adf_stat = float('inf')  # 初始化最小ADF统计量# 循环，差分阶数从1到max_difffor i in range(1, max_diff+1):# 对数据进行差分，并去除NA值df_diff = df['每半月兵力损失'].diff(i).dropna()# 对差分后的数据进行ADF单位根检验result = adfuller(df_diff)# 打印出差分阶数，ADF统计量，p值，标准差和方差print(f'{i}步差分')print('ADF统计量: %f' % result[0])print('p值: %.10e' % result[1])print('标准差: %f' % df_diff.std())print('方差: %f' % df_diff.var())# 判断p值是否小于显著性水平，如果小于则认为差分后的数据可能是平稳的if result[1] < significance_level:print('=> 根据这个差分阶数，序列可能是平稳的')# 判断当前的p值是否小于最小p值，如果小于则更新最小p值和最佳差分阶数if result[1] < min_pvalue:min_pvalue = result[1]best_diff = imin_variance = df_diff.var()  # 更新最小方差min_adf_stat = result[0]  # 更新最小ADF统计量else:print('=> 根据这个差分阶数，序列可能是非平稳的')print('--------------------------------')# 如果找到了使数据平稳的差分阶数，打印出最佳差分阶数和其对应的p值if best_diff is not None:print(f'最佳差分阶数是: {best_diff}，p值为: {min_pvalue}，方差为: {min_variance}，ADF统计量为: {min_adf_stat}')# 使用函数对数据进行差分并测试其平稳性
calculate_diff(df, max_diff=24)

这段代码的主要目的是对时间序列进行差分并测试差分后的数据是否平稳。它首先在一个循环中对数据进行不同步数的差分，并对差分后的数据进行ADF（Augmented Dickey-Fuller）单位根检验。然后在所有使数据平稳的差分阶数中，选择p值最小的作为最佳差分步数。部分输出如下：

如果许多差分都能够满足单位根检验，我们通常选择ADF值最小、或差分后数据方差最小的步数。

ADF统计量：ADF统计量越小（越负），说明数据越可能是平稳的。所以在满足平稳性的情况下，可以选择ADF统计量最小的差分阶数。
p值：p值要低于预定的显著性水平（如0.05），才能认定该序列是平稳的。如果多个差分阶数的p值都低于显著性水平，可以选择ADF统计量最小的差分阶数。
标准差和方差：理想的平稳序列，其标准差和方差应该维持在一定范围，不应有大的波动。可以考虑选择使标准差和方差最小的差分阶数。

3.3.1 绘制差分后的图像

我们可以绘制差分后的图像和原始图像作对比，来查看数据的平稳性。
如上：我们选择的最佳差分步数是1步，所以我们就选择1步差分的数据进行绘制（我的数据选的不是很好）

import matplotlib.pyplot as plt
import pandas as pdio = r'D:\Users\25705\Desktop\乌方每半月兵力损失.xlsx'
df = pd.read_excel(io,sheet_name=0)
df.head()
# 设置1步差分
diff_steps = 1
df_diff = df['每半月兵力损失'].diff(diff_steps).dropna()
plt.figure(figsize=(12,6))# 绘制原始数据
plt.plot(df['每半月兵力损失'], label='Original Data')# 绘制12步差分后的数据
plt.plot(df_diff, label=f'Data after {diff_steps}-step differencing', color='orange')# 设置标题和标签
plt.title('Original Data vs Differenced Data')
plt.xlabel('日期')
plt.ylabel('每半月兵力损失')
plt.legend()
plt.show()

在原始数据和差分后的数据中，数值的范围可能有很大的不同。原始数据可能包含趋势，其数值可能会随时间不断增加或减少，导致数据的范围较大。而差分是用来消除这种趋势的一种方法，差分后的数据主要表示相邻观测值之间的变化量，这个变化量可能相对较小，导致差分后的数据范围较小。因此，原始数据和差分后的数据在同一张图中展示时，如果共享同一y轴，可能会使差分后的数据的变化不明显。
差分后的数据(黄线)主要表示相邻观测值之间的变化量，所以我们观察这个变化量的模式。变化量在0附近随机波动，没有明显的趋势和季节性，那么可以认为差分后的数据是平稳的。

所以我们选择1步差分的数据作为后续的建模数据

df['每半月兵力损失_1'] = df['每半月兵力损失'].diff(1)
output_file_path = 'D:\\Users\\25705\\Desktop\\差分后数据.xlsx'
df.to_excel(output_file_path)  # 默认情况下，pandas 会保存索引

[3.4白噪声检验(Ljung-Box检验)]

白噪声是一种具有随机性的信号，其定义为均值为零，方差为常数的时间序列。具有白噪声属性的信号在统计学上是无相关的，意味着它们的时间间隔没有任何关联。白噪声检验方法常用有以下3种方法（自相关图、Box-Pierce检验、Ljung-Box检验），其中Ljung-Box检验相对用的多一些，Ljung-Box检验即LB检验、随机性检验，用来检验序列是否为白噪声，若是白噪声数据，则该数据没有价值提取，即不用继续分析了。
💧白噪声检验的前提是数据是平稳的。

from statsmodels.stats.diagnostic import acorr_ljungbox as lb_test
import xlrdfile_location = "D:\\Users\\25705\\Desktop\\乌方每半月兵力损失.xlsx"#避免被转义使用\\
data = xlrd.open_workbook(file_location)#	data是Excel里的数据
sheet = data.sheet_by_index(0)  #根据索引读取,就是说data里有很多sheet	data.sheet_by_index（0）是sheet1，以此类推
days = [sheet.cell_value(r,1) for r in range(1,sheet.nrows)]# 读取该表格中第2列中的数据，（r,n）表示第（n+1）列， range(1,sheet.nrows)从第2行开始
x = list(days)# days为generator，数据转换成list
re = lb_test(x, lags=20)
print(re)

结果如下：

       lb_stat   lb_pvalue1    0.479715   0.4885522    0.501843   0.7780833    0.606819   0.8948704    0.637297   0.9588305    0.758627   0.9795976    0.781087   0.9925717    1.712527   0.9740288    2.426010   0.9650929    2.993907   0.96453610   3.346418   0.97205811   4.176186   0.96447712   4.179372   0.97997813   5.004039   0.97510414   5.979877   0.96699615   6.138485   0.97732616  13.110057   0.66469417  13.174055   0.72445818  13.174702   0.78109219  13.174787   0.82950420  13.245442   0.866607

我们主要看第二列的P值，lags为检验的延迟数，一般指定是20，或是序列长度，每一个P值都小于0.05或等于0，说明该数据不是白噪声数据，数据有价值，可以继续分析。反之，如果大于0.05，则说明是白噪声序列，是纯随机性序列。参考：时间序列平稳性检验（ADF）和白噪声检验（Ljung-Box）

😥 在时间序列中，白噪声检验除了用于在预测前判断平稳序列是否随机外，还能有哪些用法呢？
检验残差是否为白噪声，判断模型拟合的是否足够好，是否还存在有价值的信息待提取。
残差为白噪声，说明模型拟合的很好，残差部分为无法捕捉的纯随机数据。
残差非白噪声，说明模型哪里出了问题，比如参数没调好，需要继续优化；若如何优化模型也无法使得残差为白噪声，换模型或者集成模型，或者对残差进行二次预测。

四.模型选择

4.1分析ACF及PACF

ACF和PACF图有助于你确定时间序列模型（例如ARIMA模型）的参数。这两种图形分别反映了时间序列的自相关性和偏自相关性。我们通过图像来分析。

import matplotlib.pyplot as plt
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
import pandas as pdio = r'D:\Users\25705\Desktop\差分后数据.xlsx'
df = pd.read_excel(io,sheet_name=0)
df.head()# 绘制ACF图
plt.figure(figsize=(12,8))
plot_acf(df['每半月兵力损失_1'].dropna(), lags=11, ax=plt.gca())# lags根据数据个数设置
plt.title('Autocorrelation for 每半月兵力损失 with 1 lags Difference')
plt.show()# 绘制PACF图
plt.figure(figsize=(12,8))
plot_pacf(df['每半月兵力损失_1'].dropna(), lags=11, ax=plt.gca())
plt.title('Partial Autocorrelation for 每半月兵力损失 with 1 lags Difference')
plt.show()

在这段代码中，我们用到了plot_acf和plot_pacf函数来绘制ACF和PACF图。lags参数用于设定我们想要观察的滞后阶数，ax=plt.gca()用于将图形绘制在当前的Axes对象上。
请注意，因为在计算差分时会引入NaN值，所以在绘制图形之前，我们需要用dropna()函数来去掉这些值。
从ACF和PACF的图像可以看出，两个自相关系数都没有呈现拖尾的情况，因此p和q都不为0，此时我们需要的至少是一个ARIMA(p,d,q)模型。

针对差分后序列根据其自相关系数（ACF）图和偏自相关（PACF）图选择AR、MA、ARMA模型。

💜拖尾：指的是ACF或PACF并不在某阶后均为0，而是呈现出一种衰减的趋势，但并不会完全为0。这通常意味着时间序列数据具有长期记忆性，即过去的数据对未来的数据仍有一定的影响。拖尾的情况在AR模型和MA模型中都有可能出现。

💜截尾：截尾则是指ACF或PACF在某阶后均为0的性质。这意味着时间序列数据在某一阶数后，过去的数据对未来数据的影响可以忽略不计。在AR模型中，PACF通常表现出截尾性，而在MA模型中，ACF则通常表现出截尾性。

💜模型选择：对于ARMA模型的ACF和PACF图，我们可以通过观察其图形特征来判断模型的阶数。如果ACF图呈现出拖尾的特征，而PACF图呈现出截尾的特征，那么可以考虑使用AR模型进行拟合；如果ACF图呈现出截尾的特征，而PACF图呈现出拖尾的特征，那么可以考虑使用MA模型进行拟合。如果ACF和PACF图都呈现出拖尾的特征，那么可能需要考虑使用ARMA模型进行拟合。

需要注意的是，在实际应用中，我们可能需要根据最优信息准则（如AIC、BIC等）来选择最合适的模型阶数，而不仅仅依赖于ACF和PACF图的图形特征。此外，对于复杂的时间序列数据，可能需要结合其他统计方法和模型来进行分析和预测。

4.2使用多阶差分确定d值

当我们确定ARIMA模型的p和q参数都不为0后，我们需要找到一个适合的d值，也就是使序列平稳的最小差分阶数。我们通常会尝试从0到3的阶数，并从中选择一个合适的。我们的目标是找到一个可以消除数据中的趋势和季节性，同时使差分后的序列方差最小、噪声最小的d值。

一种有效的策略是绘制不同阶差分后的自相关函数（ACF）图。如果ACF图中的第一个滞后值是负的，那么我们可能会怀疑数据存在过差分问题，也就是我们对数据做了过多的差分。另外，我们还需要注意的是，我们应该选择能使差分后的序列方差最小的d值。因此，我们将打印出0到3阶差分的序列方差和ACF图，以便进行比较好的选择。

为什么此处不再绘制PACF图了呢？因为我们的数据已经是经过1步差分后的数据，是一个平稳数据了，当尝试对一个已经是平稳的序列再次进行差分操作，特别是当序列中有大量的0时，这会导致计算偏自相关函数（PACF）时出现问题。

from statsmodels.graphics.tsaplots import plot_acf
import matplotlib.pyplot as plt
import numpy as np
import pandas as pdio = r'D:\Users\25705\Desktop\差分后数据.xlsx'
df = pd.read_excel(io,sheet_name=0)
df.head()def diff(data, order):if order == 0:return dataelse:return diff(data.diff().dropna(), order - 1)# 循环不同的d值
for d in range(1, 4):# 创建新的图形plt.figure(figsize=(12, 4))plt.suptitle(f'Order {d} Difference', fontsize=16)# 绘制ACF图plot_acf(diff(df['每半月兵力损失_1'], d), ax=plt.gca(), title='Autocorrelation')# 注意：这里是已经差分后的数据plt.show()

这个函数每次它被调用，它会对数据进行一阶差分，然后将差分后的数据和阶数减一的值作为新的参数再次调用自身，直到阶数降到0为止。然后对每个d值，它都会绘制对应的ACF图。

除了原始数据，其他阶差分都出现了“过差分”的情况（滞后为1的ACF值为负），因此我们可以判断对ARIMA模型而言最佳的d取值为0。

4.3ARIMA模型实现

4.3.1获取不同的p、d、q组合

我们上面已经确定了，d=0，所以我们此处做一下pdq组合，一般p和q不会取很大的值，此处取1~4

from statsmodels.graphics.tsaplots import plot_acf
import matplotlib.pyplot as plt
import numpy as np
import itertools
import pandas as pdio = r'D:\Users\25705\Desktop\差分后数据.xlsx'
df = pd.read_excel(io,sheet_name=0)
df.head()#定义p, d和q参数，取0到4之间的任意值
p = q = range(0, 4)# 生成p, d和q三元组的所有不同组合
pdq = [(x[0], 0, x[1]) for x in list(itertools.product(p, q))]
print(pdq)

4.3.2超参数搜索

我们的目标是找到能够最大化模型性能的p和q值。我们可以通过遍历不同的p和q值，并使用每个值组合来拟合模型，然后比较每个模型的性能，来找到最优的p和q值。模型性能可以使用AIC（Akaike Information Criterion）来评估，AIC越小，模型的性能越好。

以下是对p和q值进行调优的代码，此处我们先选择MSE作为评估标准。

import warnings
warnings.filterwarnings("ignore")
from sklearn.metrics import mean_squared_error
import matplotlib.pyplot as plt
from statsmodels.tsa.arima.model import ARIMA
import itertools
from statsmodels.graphics.tsaplots import plot_acf
import matplotlib.pyplot as plt
import numpy as np
import pandas as pdio = r'D:\Users\25705\Desktop\差分后数据.xlsx'
df = pd.read_excel(io,sheet_name=0)
df.head()#定义p, d和q参数，取0到4之间的任意值
p = q = range(0, 4)# 生成p, d和q三元组的所有不同组合
pdq = [(x[0], 0, x[1]) for x in list(itertools.product(p, q))]
print(pdq)# 将数据分为训练集和测试集
train_data = df['每半月兵力损失_1'].iloc[:-20]
test_data = df['每半月兵力损失_1'].iloc[-20:]# 初始化最佳模型及其参数和MSE值
best_model = None
best_param = None
best_mse = float('inf')# 对每一种参数组合进行迭代
for param in pdq:try:# 实例化ARIMA模型model = ARIMA(train_data, order=param)model_fit = model.fit()# 进行预测predictions = model_fit.predict(start=len(train_data), end=len(train_data) + len(test_data) - 1)# 计算MSEmse = mean_squared_error(test_data, predictions)# 如果当前模型的MSE比最佳MSE小，更新最佳模型、参数和MSEif mse < best_mse:best_model = model_fitbest_param = parambest_mse = mse# 绘制真实值和预测值plt.figure(figsize=(6, 4))plt.plot(df.index, df['每半月兵力损失_1'], label='Actual')plt.plot(df.index[-20:], predictions, color='red', label='Predicted')# 添加标题和标签plt.title('Actual vs Predicted Sales for ARIMA{}'.format(param))plt.xlabel('日期')plt.ylabel('每半月兵力损失_1')plt.legend()# 显示图形plt.show()except Exception as e:print('Error:', e)continue# 打印出最优模型的参数及其MSE
print(f'The best model is ARIMA{best_param}, MSE = {best_mse}')# 打印最优模型的摘要
print(best_model.summary(alpha=0.05))

并打印出最优模型的参数及其MSE。分析见时间序列模型(五)：时间序列案例_实现销售额预测并有手动调优SARIMA模型。