python打卡day56

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

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

import numpy as np
import matplotlib.pyplot as plt
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from statsmodels.tsa.stattools import adfuller
from statsmodels.tsa.seasonal import seasonal_decompose
import pandas as pd
from scipy import stats
import warnings
warnings.filterwarnings('ignore')

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

# 设置随机种子以便结果可重现
np.random.seed(42)

# 时间设置
num_points = 200
time = np.arange(num_points)

# ==================== 1. 构造不同类型的数据集 ====================

# 1.1 白噪声数据
white_noise = np.random.normal(0, 1, num_points)

# 1.2 非平稳数据（带趋势）
trend_data = 0.05 * time + np.random.normal(0, 1, num_points)

# 1.3 季节性数据
seasonal_data = 3 * np.sin(2 * np.pi * time / 12) + np.random.normal(0, 0.5, num_points)

# 1.4 复合数据（趋势+季节性+噪声）
complex_data = 0.02 * time + 2 * np.sin(2 * np.pi * time / 24) + np.random.normal(0, 0.8, num_points)

# 1.5 随机游走（非平稳）
random_walk = np.cumsum(np.random.normal(0, 1, num_points))

# ==================== 2. 检验函数定义 ====================

def check_white_noise(data, name):
    """检验白噪声属性"""
    print(f"\n=== {name} 白噪声检验 ===")
    
    # 1. 均值检验
    mean_val = np.mean(data)
    print(f"均值: {mean_val:.4f} (接近0为好)")
    
    # 2. 方差稳定性检验
    std_val = np.std(data)
    print(f"标准差: {std_val:.4f}")
    
    # 3. 正态性检验 (Shapiro-Wilk)
    stat, p_value = stats.shapiro(data[:50])  # 取前50个点检验
    print(f"正态性检验 p值: {p_value:.4f} (>0.05为正态分布)")
    
    # 4. 自相关检验 - 绘制ACF图
    fig, ax = plt.subplots(figsize=(10, 4))
    plot_acf(data, lags=30, ax=ax, alpha=0.05)
    ax.set_title(f'{name} - 自相关函数图 (ACF)')
    plt.tight_layout()
    plt.show()
    
    # 5. Ljung-Box检验（检验是否为白噪声）
    from statsmodels.stats.diagnostic import acorr_ljungbox
    lb_stat = acorr_ljungbox(data, lags=10, return_df=True)
    print(f"Ljung-Box检验 p值: {lb_stat['lb_pvalue'].iloc[-1]:.4f} (>0.05为白噪声)")
    
    return mean_val, std_val, p_value

def check_stationarity(data, name):
    """检验平稳性"""
    print(f"\n=== {name} 平稳性检验 ===")
    
    # ADF检验
    adf_result = adfuller(data)
    print(f"ADF统计量: {adf_result[0]:.4f}")
    print(f"p值: {adf_result[1]:.4f} (<0.05拒绝原假设，即序列平稳)")
    print(f"临界值:")
    for key, value in adf_result[4].items():
        print(f"\t{key}: {value:.4f}")
    
    if adf_result[1] < 0.05:
        print("结论: 序列是平稳的")
    else:
        print("结论: 序列是非平稳的")
    
    # 绘制时间序列图
    plt.figure(figsize=(12, 4))
    plt.plot(time, data)
    plt.title(f'{name} - 时间序列图')
    plt.xlabel('时间')
    plt.ylabel('数值')
    plt.grid(True, alpha=0.3)
    plt.show()
    
    return adf_result[1] < 0.05

def check_seasonality(data, name, period=12):
    """检验季节性"""
    print(f"\n=== {name} 季节性检验 ===")
    
    # 创建时间索引
    dates = pd.date_range('2020-01-01', periods=len(data), freq='M')
    ts = pd.Series(data, index=dates)
    
    # 季节性分解
    try:
        decomposition = seasonal_decompose(ts, model='additive', period=period)
        
        # 绘制分解图
        fig, axes = plt.subplots(4, 1, figsize=(12, 10))
        
        decomposition.observed.plot(ax=axes[0], title=f'{name} - 原始序列')
        decomposition.trend.plot(ax=axes[1], title='趋势')
        decomposition.seasonal.plot(ax=axes[2], title='季节性')
        decomposition.resid.plot(ax=axes[3], title='残差')
        
        plt.tight_layout()
        plt.show()
        
        # 计算季节性强度
        seasonal_strength = np.var(decomposition.seasonal.dropna()) / np.var(decomposition.observed.dropna())
        print(f"季节性强度: {seasonal_strength:.4f} (>0.1表示有明显季节性)")
        
        return seasonal_strength > 0.1
        
    except Exception as e:
        print(f"季节性分解失败: {e}")
        return False

def comprehensive_analysis(data, name):
    """综合分析"""
    print(f"\n{'='*50}")
    print(f"正在分析: {name}")
    print(f"{'='*50}")
    
    # 白噪声检验
    mean_val, std_val, normality_p = check_white_noise(data, name)
    
    # 平稳性检验
    is_stationary = check_stationarity(data, name)
    
    # 季节性检验
    has_seasonality = check_seasonality(data, name)
    
    # 综合结论
    print(f"\n=== {name} 综合结论 ===")
    print(f"是否为白噪声: {'是' if abs(mean_val) < 0.1 and normality_p > 0.05 else '否'}")
    print(f"是否平稳: {'是' if is_stationary else '否'}")
    print(f"是否有季节性: {'是' if has_seasonality else '否'}")
    
    return {
        'name': name,
        'is_white_noise': abs(mean_val) < 0.1 and normality_p > 0.05,
        'is_stationary': is_stationary,
        'has_seasonality': has_seasonality
    }

# ==================== 3. 执行分析 ====================

# 分析所有数据集
datasets = {
    '白噪声': white_noise,
    '趋势数据': trend_data,
    '季节性数据': seasonal_data,
    '复合数据': complex_data,
    '随机游走': random_walk
}

results = []
for name, data in datasets.items():
    result = comprehensive_analysis(data, name)
    results.append(result)

# ==================== 4. 结果汇总 ====================

print("\n" + "="*60)
print("所有数据集分析结果汇总")
print("="*60)
print(f"{'数据集':<10} {'白噪声':<8} {'平稳性':<8} {'季节性':<8}")
print("-"*40)
for result in results:
    print(f"{result['name']:<10} {'✓' if result['is_white_noise'] else '✗':<8} {'✓' if result['is_stationary'] else '✗':<8} {'✓' if result['has_seasonality'] else '✗':<8}")

print("\n说明:")
print("✓ = 是/有")
print("✗ = 否/无")

==================================================
正在分析: 白噪声
==================================================

=== 白噪声 白噪声检验 ===
均值: -0.0408 (接近0为好)
标准差: 0.9287
正态性检验 p值: 0.6722 (>0.05为正态分布)