Python打卡训练营Day56

DAY 56 时序数据的检验

知识点回顾：

1. 假设检验基础知识
   1. 原假设与备择假设
   2. P值、统计量、显著水平、置信区间
2. 白噪声
   1. 白噪声的定义
   2. 自相关性检验：ACF检验和Ljung-Box 检验
   3. 偏自相关性检验：PACF检验
3. 平稳性
   1. 平稳性的定义
   2. 单位根检验
4. 季节性检验
   1. ACF检验
   2. 序列分解：趋势+季节性+残差

记忆口诀：p越小，落在置信区间外，越拒绝原假设。

时序部分需要铺垫的知识非常多，相信这次应该说清楚了假设检验相关的基础知识。

import numpy as np
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings("ignore")
# 中文显示设置
plt.rcParams['font.sans-serif'] = ['SimHei']  # 设置中文字体
plt.rcParams['axes.unicode_minus'] = False  # 解决负号显示为方块的问题# --- 1. 生成随机序列数据 ---# 为了让每次运行的结果都一样，设置一个随机种子（可选）
np.random.seed(42)# 定义序列的长度
num_points = 200# 生成一个包含 200 个点的随机序列
# np.random.randn() 从标准正态分布（均值为0，方差为1）中抽取随机样本
random_sequence = np.random.randn(num_points)print("生成的前10个数据点:")
print(random_sequence[:10])# --- 2. 可视化序列 ---# 设置图形大小
plt.figure(figsize=(12, 6))# 绘制线图
plt.plot(random_sequence, label='Random Sequence (White Noise)')# 添加标题和标签
plt.title('Visualization of a Randomly Generated Sequence', fontsize=16)
plt.xlabel('Time Step (时间步)', fontsize=12)
plt.ylabel('Value (值)', fontsize=12)# 添加一条水平线，表示序列的均值（接近于0）
plt.axhline(y=0, color='r', linestyle='--', label='Mean (均值 ≈ 0)')# 显示网格和图例
plt.grid(True, linestyle='--', alpha=0.6)
plt.legend()# 显示图形
plt.show()from statsmodels.graphics.tsaplots import plot_acf 
print("--- 开始检验白噪声属性 ---")# 检验 1: 均值是否接近 0
mean = np.mean(random_sequence)
print(f"1. 序列的均值: {mean:.4f}")
if -0.1 < mean < 0.1:print("   (结论: 均值非常接近0，满足条件。)\n")
else:print("   (结论: 均值偏离0较远。)\n")# 检验 2: 方差是否恒定（且接近理论值1）
# 对于我们生成的数据，方差恒定是与生俱来的。我们主要检查其值。
variance = np.var(random_sequence)
print(f"2. 序列的方差: {variance:.4f}")
if 0.8 < variance < 1.2:print("   (结论: 方差接近于1，满足条件。np.random.randn理论方差为1)\n")
else:print("   (结论: 方差偏离1较远。)\n")# 检验 3: 自相关性是否为 0
# 这是最核心的检验。我们通过绘制ACF图来完成。
print("3. 检验自相关性 (使用ACF图):")
print("   - ACF图展示了序列与它过去值之间的相关性。")
print("   - 对于白噪声，只有lag=0时相关性为1，其他所有lag的相关性都应在蓝色置信区间内（统计上不显著）。")# 创建一个新的图形来绘制ACF图
fig, ax = plt.subplots(figsize=(12, 5))
plot_acf(random_sequence, lags=30, ax=ax) # 我们查看前30个滞后的相关性
ax.set_title('序列的自相关函数图 (ACF Plot)')
ax.set_xlabel('Lag (滞后阶数)')
ax.set_ylabel('Autocorrelation (自相关系数)')
plt.grid(True, linestyle='--', alpha=0.6)
plt.show()from statsmodels.graphics.tsaplots import plot_pacf # 引入PACF图
# --- 绘制PACF图 ---
fig, ax = plt.subplots(figsize=(12, 5))
plot_pacf(random_sequence, lags=30, ax=ax)
ax.set_title('序列的偏自相关函数图 (PACF Plot)')
ax.set_xlabel('Lag (滞后阶数)')
ax.set_ylabel('Partial Autocorrelation')
plt.grid(True, linestyle='--', alpha=0.6)
plt.show()# --- 新增：使用 Ljung-Box 检验进行严格的白噪声检验 ---
# 引入Ljung-Box检验的函数
from statsmodels.stats.diagnostic import acorr_ljungbox 
print("\n" + "="*50)
print("4. 进行严格的白噪声检验 (Ljung-Box Test)")
print("="*50)
print("   - 原假设(H₀): 序列是白噪声。")
print("   - 判断标准: 如果 p-value > 0.05，则接受原假设，认为序列是白噪声。")# 执行Ljung-Box检验
# 我们通常会检查一系列的滞后项，比如前10、20、30个
# 函数返回一个包含统计量和p值的DataFrame
ljung_box_result = acorr_ljungbox(random_sequence, lags=[10, 20, 30], return_df=True)print("\nLjung-Box检验结果:")
print(ljung_box_result)# --- 结论解释 ---
print("\n--- 检验结论 ---")
# 我们可以检查最后一个（最严格的）p值
# .iloc[-1] 获取最后一行, .loc['lb_pvalue'] 获取p值
last_p_value = ljung_box_result.iloc[-1]['lb_pvalue']if last_p_value < 0.05:print(f"在滞后30阶时，p-value ({last_p_value:.4f}) 小于 0.05。")print("结论：我们拒绝原假设，该序列不是白噪声。")
else:print(f"在滞后30阶时，p-value ({last_p_value:.4f}) 大于 0.05。")print("结论：我们无法拒绝原假设，该序列是白噪声。")# 引入ADF检验的函数
from statsmodels.tsa.stattools import adfuller # --- 新增：使用ADF检验来判断平稳性 ---print("开始进行ADF平稳性检验...")# 执行ADF检验
# adfuller()函数会返回一个包含多个结果的元组
adf_result = adfuller(random_sequence)# 提取并展示主要结果
adf_statistic = adf_result[0]
p_value = adf_result[1]
critical_values = adf_result[4]print(f"ADF统计量 (ADF Statistic): {adf_statistic:.4f}")
print(f"p值 (p-value): {p_value:.4f}")
print("临界值 (Critical Values):")
for key, value in critical_values.items():print(f'    {key}: {value:.4f}')print("\n--- 检验结论 ---")
# 根据p值进行判断
if p_value < 0.05:print(f"p-value ({p_value:.4f}) 小于 0.05，我们强烈拒绝原假设(H₀)。")print("结论：该序列是平稳的 (Stationary)。")
else:print(f"p-value ({p_value:.4f}) 大于或等于 0.05，我们无法拒绝原假设(H₀)。")print("结论：该序列是非平稳的 (Non-stationary)。")# 也可以通过比较ADF统计量和临界值来判断，结论是一致的
if adf_statistic < critical_values['5%']:print("\n补充判断：ADF统计量小于5%的临界值，同样表明序列是平稳的。")import numpy as np
import matplotlib.pyplot as plt
from statsmodels.tsa.seasonal import seasonal_decompose
from statsmodels.graphics.tsaplots import plot_acf
# 显示中文
plt.rcParams['font.sans-serif'] = ['SimHei']  # 设置中文字体
plt.rcParams['axes.unicode_minus'] = False# --- 1. 创建一个带季节性的序列 ---
# 我们模拟一个为期5年的月度数据（60个点）
num_points = 60
time = np.arange(num_points)# a. 创建一个线性趋势
trend = 0.5 * time# b. 创建一个季节性成分（周期为12个月）
# 使用sin函数来模拟年度周期性波动
seasonal_component = 15 * np.sin(2 * np.pi * time / 12)# c. 创建一些随机噪声
np.random.seed(10)
noise = np.random.randn(num_points) * 2# d. 合成最终序列
seasonal_data = trend + seasonal_component + noise# --- 2. 开始检验季节性 ---# 方法一：肉眼观察
print("--- 方法一：肉眼观察 ---")
plt.figure(figsize=(14, 6))
plt.plot(seasonal_data)
plt.title('带趋势和季节性的时间序列图', fontsize=16)
plt.xlabel('时间步 (月)')
plt.ylabel('值')
plt.grid(True, linestyle='--', alpha=0.6)
plt.show()
# 观察：我们可以清晰地看到一个整体上升的趋势，以及每年重复的波峰和波谷。# 方法二：ACF图
print("\n--- 方法二：ACF图 ---")
fig, ax = plt.subplots(figsize=(14, 6))
plot_acf(seasonal_data, lags=30, ax=ax)
ax.set_title('季节性序列的ACF图')
plt.grid(True, linestyle='--', alpha=0.6)
plt.show()
# 观察：ACF图不仅整体缓慢下降（表明有趋势），更重要的是在lag=12和24的位置出现了明显的峰值！# 方法三：序列分解
print("\n--- 方法三：序列分解 ---")
# 使用statsmodels进行分解（假设为加法模型）
decomposition = seasonal_decompose(seasonal_data, model='additive', period=12)# 绘制分解图
fig = decomposition.plot()
fig.set_size_inches(14, 8)
plt.suptitle('时间序列分解图', y=1.02, fontsize=16)
plt.show()
# 观察：分解图清晰地将数据拆分成了趋势、季节性和残差。季节性部分呈现完美的年度周期，而残差看起来像随机噪声。

作业：自行构造数据集，来检查是否符合这个要求。

记忆口诀：p越小，落在置信区间外，越拒绝原假设。

什么叫做白噪声呢？他需要满足以下条件：

1. 均值为0

2. 方差恒定

3. 自相关性为0（即过去的值对未来的值没有影响）

# 构造数据集 24h动态变化周期性的函数 模拟5day中每个24h人口流动的数据集 
num_points = 120
time = np.arange(num_points)# a. 创建一个线性趋势
trend = 1.0 * time# b. 创建一个季节性成分（周期为24h）
# 使用sin函数来模拟周期性波动
day_component = 15 * np.sin(2 * np.pi * time / 24)# c. 创建一些随机噪声
np.random.seed(10)
noise = np.random.randn(num_points) * 2# d. 合成最终序列
daily_data = trend + day_component + noise#step0 可视化
# 设置图形大小
plt.figure(figsize=(12, 6))# 绘制线图
plt.plot(daily_data, label='Daily Data')# 添加标题和标签
plt.title('Visualization of a 24h Dynamic Periodic Time Series', fontsize=16)
plt.xlabel('Time Step (时间步)', fontsize=12)
plt.ylabel('Value (值)', fontsize=12)# 添加一条水平线，表示序列的均值（接近于0）
plt.axhline(y=0, color='r', linestyle='--', label='Mean (均值 ≈ 0)')# 显示网格和图例
plt.grid(True, linestyle='--', alpha=0.6)
plt.legend()# 显示图形
plt.show()#step1 可预测性检验 使用 Ljung-Box 检验进行严格的白噪声检验 
print("\n" + "="*50)
print("4. 进行严格的白噪声检验 (Ljung-Box Test)")
print("="*50)
print("   - 原假设(H₀): 序列是白噪声。")
print("   - 判断标准: 如果 p-value > 0.05，则接受原假设，认为序列是白噪声。")# 执行Ljung-Box检验
# 我们通常会检查一系列的滞后项，比如前10、20、30个
# 函数返回一个包含统计量和p值的DataFrame
ljung_box_result = acorr_ljungbox(daily_data, lags=[10, 20, 30], return_df=True)print("\nLjung-Box检验结果:")
print(ljung_box_result)# --- 结论解释 ---
print("\n--- 检验结论 ---")
# 我们可以检查最后一个（最严格的）p值
# .iloc[-1] 获取最后一行, .loc['lb_pvalue'] 获取p值
last_p_value = ljung_box_result.iloc[-1]['lb_pvalue']if last_p_value < 0.05:print(f"在滞后30阶时，p-value ({last_p_value:.4f}) 小于 0.05。")print("结论：我们拒绝原假设，该序列不是白噪声。")
else:print(f"在滞后30阶时，p-value ({last_p_value:.4f}) 大于 0.05。")print("结论：我们无法拒绝原假设，该序列是白噪声。")#step2 平稳性检验
# --- 新增：使用ADF检验来判断平稳性 ---print("开始进行ADF平稳性检验...")# 执行ADF检验
# adfuller()函数会返回一个包含多个结果的元组
adf_result = adfuller(daily_data)# 提取并展示主要结果
adf_statistic = adf_result[0]
p_value = adf_result[1]
critical_values = adf_result[4]print(f"ADF统计量 (ADF Statistic): {adf_statistic:.4f}")
print(f"p值 (p-value): {p_value:.4f}")
print("临界值 (Critical Values):")
for key, value in critical_values.items():print(f'    {key}: {value:.4f}')print("\n--- 检验结论 ---")
# 根据p值进行判断
if p_value < 0.05:print(f"p-value ({p_value:.4f}) 小于 0.05，我们强烈拒绝原假设(H₀)。")print("结论：该序列是平稳的 (Stationary)。")
else:print(f"p-value ({p_value:.4f}) 大于或等于 0.05，我们无法拒绝原假设(H₀)。")print("结论：该序列是非平稳的 (Non-stationary)。")# 也可以通过比较ADF统计量和临界值来判断，结论是一致的
if adf_statistic < critical_values['5%']:print("\n补充判断：ADF统计量小于5%的临界值，同样表明序列是平稳的。")#step3 结构识别（序列分解，ACF/PACF）
print("\n--- 方法二：ACF图 ---")
fig, ax = plt.subplots(figsize=(14, 6))
plot_acf(daily_data, lags=30, ax=ax)
ax.set_title('24h序列的ACF图')
plt.grid(True, linestyle='--', alpha=0.6)
plt.show()
# 观察：ACF图不仅整体缓慢下降（表明有趋势）# 方法三：序列分解
print("\n--- 方法三：序列分解 ---")
# 使用statsmodels进行分解（假设为加法模型）
decomposition = seasonal_decompose(daily_data, model='additive', period=24)# 绘制分解图
fig = decomposition.plot()
fig.set_size_inches(14, 8)
plt.suptitle('时间序列分解图', y=1.02, fontsize=16)
plt.show()
# 观察：分解图清晰地将数据拆分成了趋势、季节性和残差。季节性部分呈现完美的年度周期，而残差看起来像随机噪声。#stpe4 数据清洗 箱线图和Z分数识别异常值
# 绘制箱线图
plt.figure(figsize=(12, 6))
plt.boxplot(daily_data, vert=False, widths=0.5, patch_artist=True, boxprops=dict(facecolor='blue', color='blue'))
plt.title('24h动态周期性时间序列的箱线图', fontsize=16)
plt.xlabel('值', fontsize=12)
plt.yticks([])  # 隐藏y轴刻度
plt.grid(True, linestyle='--', alpha=0.6)
plt.show()
# 观察：箱线图显示了数据的四分位数、中位数和异常值。这里没有明显的异常值，数据分布相对比较均匀。# 计算Z分数并识别异常值
# 计算Z分数
z_scores = np.abs((daily_data - daily_data.mean()) / daily_data.std())# 定义Z分数阈值（例如：超过3个标准差）
z_threshold = 3# 识别异常值
outlier_indices = np.where(z_scores > z_threshold)[0]# 打印异常值的索引和值
if outlier_indices.size > 0:print(f"发现 {outlier_indices.size} 个异常值（超过Z分数阈值）：")for idx in outlier_indices:print(f"索引 {idx}: 值 {daily_data[idx]:.4f}, Z分数 {z_scores[idx]:.4f}")
else:print("没有发现异常值。")

==================================================

进行严格的白噪声检验 (Ljung-Box Test) ================================================== -

原假设(H₀): 序列是白噪声。

- 判断标准: 如果 p-value > 0.05，则接受原假设，认为序列是白噪声。

Ljung-Box检验结果: lb_stat lb_pvalue 10 869.534899 2.289377e-180 20 1307.450690 7.507512e-265 30 1529.131369 2.447689e-303

--- 检验结论 --- 在滞后30阶时，p-value (0.0000) 小于 0.05。

结论：我们拒绝原假设，该序列不是白噪声。

开始进行ADF平稳性检验

... ADF统计量 (ADF Statistic): -0.2350

p值 (p-value): 0.9342

临界值 (Critical Values): 1%: -3.4936    5%: -2.8892   10%: -2.5815

--- 检验结论 --- p-value (0.9342) 大于或等于 0.05，我们无法拒绝原假设(H₀)。

结论：该序列是非平稳的 (Non-stationary)。

Z-scores 没有发现异常值。