Python打卡訓練營Day56

DAY 56 時序數據的檢驗

知識點回顧：

假設檢驗基礎知識
1. 原假設與備擇假設
2. P值、統計量、顯著水平、置信區間
白噪聲
1. 白噪聲的定義
2. 自相關性檢驗：ACF檢驗和Ljung-Box 檢驗
3. 偏自相關性檢驗：PACF檢驗
平穩性
1. 平穩性的定義
2. 單位根檢驗
季節性檢驗
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("沒有發現異常值。")