时间序列分析

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from statsmodels.tsa.stattools import adfuller, acf, pacf
from statsmodels.tsa.seasonal import seasonal_decompose
from statsmodels.tsa.holtwinters import ExponentialSmoothing, SimpleExpSmoothing
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
import warnings
warnings.filterwarnings('ignore')

plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
sns.set_theme(style='whitegrid')

# 从列创建日期时间
dates = pd.date_range(start='2023-01-01', periods=365, freq='D')
ts = pd.Series(np.random.randn(365), index=dates)

print(ts.head())
print(f'索引类型: {ts.index.dtype}')

# 字符串转 datetime
df = pd.DataFrame({
    'date': ['2024-01-01', '2024-01-02', '2024-01-03'],
    'value': [10, 20, 15]
})
df['date'] = pd.to_datetime(df['date'])
df = df.set_index('date')

np.random.seed(42)
n = 365 * 2  # 两年数据

t = np.arange(n)

# 趋势：线性增长
trend = 0.05 * t

# 季节：年度周期 + 周周期
seasonal_yearly = 10 * np.sin(2 * np.pi * t / 365)
seasonal_weekly = 3 * np.sin(2 * np.pi * t / 7)

# 噪声
noise = np.random.normal(0, 2, n)

# 合成
values = 50 + trend + seasonal_yearly + seasonal_weekly + noise

dates = pd.date_range(start='2023-01-01', periods=n, freq='D')
df = pd.DataFrame({'value': values}, index=dates)

plt.figure(figsize=(12, 4))
plt.plot(df.index, df['value'], alpha=0.7)
plt.title('合成时间序列（趋势 + 季节 + 噪声）')
plt.xlabel('日期')
plt.ylabel('值')
plt.show()

# 日 -> 月（聚合）
monthly = df.resample('ME').mean()
print('月度均值:')
print(monthly.head())

# 日 -> 周
weekly = df.resample('W').sum()

# 日 -> 季度
quarterly = df.resample('QE').mean()

# 降采样对比
fig, axes = plt.subplots(3, 1, figsize=(12, 8))
df['value'].plot(ax=axes[0], title='日数据 (365 点)')
weekly['value'].plot(ax=axes[1], title='周数据 (52 点)')
monthly['value'].plot(ax=axes[2], title='月数据 (12 点)')
plt.tight_layout()
plt.show()

# 滞后（lag）：今天的值和昨天比
df['yesterday'] = df['value'].shift(1)
df['diff_1'] = df['value'] - df['yesterday']  # 一阶差分

print('滞后和差分:')
print(df.head(10)[['value', 'yesterday', 'diff_1']])

# 多阶滞后
df['lag_7'] = df['value'].shift(7)   # 一周前
df['lag_365'] = df['value'].shift(365) # 一年前

# 日变化率
df['pct_change'] = df['value'].pct_change() * 100

# 7 日移动平均
df['ma_7'] = df['value'].rolling(window=7).mean()
df['ma_30'] = df['value'].rolling(window=30).mean()

# 滚动标准差（波动率）
df['volatility'] = df['value'].rolling(window=30).std()

# 滚动相关系数
df['rolling_corr'] = df['value'].rolling(window=60).corr(df['value'].shift(1))

# 可视化
fig, axes = plt.subplots(2, 1, figsize=(12, 6))

df[['value', 'ma_7', 'ma_30']].plot(ax=axes[0])
axes[0].set_title('原始数据与移动平均')
axes[0].set_ylabel('值')

df['volatility'].plot(ax=axes[1], color='red')
axes[1].set_title('30 日滚动波动率')
axes[1].set_ylabel('标准差')

plt.tight_layout()
plt.show()

# 累计统计量
df['cumsum'] = df['value'].expanding().sum()
df['cummean'] = df['value'].expanding().mean()
df['cummax'] = df['value'].expanding().max()
df['cummin'] = df['value'].expanding().min()

print('截至当前日期的累计统计:')
print(df[['value', 'cummean', 'cummax']].tail())

# 不同粒度的时间视图
fig, axes = plt.subplots(3, 1, figsize=(12, 8))

# 全年概览
df['2023'].plot(ax=axes[0])
axes[0].set_title('2023 年全年')

# 月度对比
df.groupby(df.index.month)['value'].mean().plot(kind='bar', ax=axes[1])
axes[1].set_title('月度均值（跨年聚合）')
axes[1].set_xlabel('月份')

# 周模式
df.groupby(df.index.dayofweek)['value'].mean().plot(kind='bar', ax=axes[2])
axes[2].set_title('周内模式（0=周一）')
axes[2].set_xticklabels(['周一', '周二', '周三', '周四', '周五', '周六', '周日'])
axes[2].set_xlabel('星期')

plt.tight_layout()
plt.show()

# 赛季性子图
fig, axes = plt.subplots(4, 3, figsize=(15, 10))

for year in [2023, 2024]:
    for month in range(1, 13):
        ax = axes[year-2023, month-1]
        month_data = df[f'{year}-{month:02d}']
        if len(month_data) > 0:
            ax.plot(month_data.index.day, month_data['value'])
            ax.set_title(f'{year}-{month:02d}')
            ax.set_xlim(0, 31)

plt.tight_layout()
plt.show()

# 生成平稳和非平稳序列
n = 300
t = np.arange(n)

# 平稳：围绕固定均值波动
stationary = np.random.normal(0, 1, n) + 50

# 非平稳：有趋势
non_stationary = 0.1 * t + np.random.normal(0, 1, n)

# 非平稳：方差变化
non_stationary_var = np.random.normal(0, 1 + t/50, n)

fig, axes = plt.subplots(3, 1, figsize=(12, 8))

axes[0].plot(stationary)
axes[0].set_title('平稳序列（均值恒定）')
axes[0].axhline(y=np.mean(stationary), color='r', linestyle='--')

axes[1].plot(non_stationary)
axes[1].set_title('非平稳（趋势）')

axes[2].plot(non_stationary_var)
axes[2].set_title('非平稳（异方差）')

plt.tight_layout()
plt.show()

def check_stationarity(series, name=''):
    result = adfuller(series.dropna())
    print(f'{name}:')
    print(f'  ADF 统计量: {result[0]:.4f}')
    print(f'  p-value: {result[1]:.4f}')
    print(f'  临界值:')
    for key, val in result[4].items():
        print(f'    {key}: {val:.4f}')
    
    if result[1] <= 0.05:
        print(f'  → 结论: 平稳（拒绝 H0，序列无单位根）\n')
    else:
        print(f'  → 结论: 非平稳（无法拒绝 H0，序列有单位根）\n')

check_stationarity(df['value'], '原始序列')

# 一阶差分后的平稳性
df['diff'] = df['value'].diff()
check_stationarity(df['diff'], '一阶差分后')

# 一阶差分
df['diff_1'] = df['value'].diff()

# 季节性差分（年度）
df['diff_seasonal'] = df['value'].diff(365)

# 可视化对比
fig, axes = plt.subplots(3, 1, figsize=(12, 8))

df['value'].plot(ax=axes[0])
axes[0].set_title('原始序列')

df['diff_1'].plot(ax=axes[1])
axes[1].set_title('一阶差分')

df['diff_seasonal'].plot(ax=axes[2])
axes[2].set_title('年度季节差分')

plt.tight_layout()
plt.show()

# 加法分解（假设各成分相加）
decomposition = seasonal_decompose(
    df['value'], 
    model='additive', 
    period=365
)

fig, axes = plt.subplots(4, 1, figsize=(12, 10))

decomposition.observed.plot(ax=axes[0])
axes[0].set_title('原始序列')

decomposition.trend.plot(ax=axes[1])
axes[1].set_title('趋势成分')

decomposition.seasonal.iloc[:365].plot(ax=axes[2])  # 只画一个周期
axes[2].set_title('季节成分（一年周期）')

decomposition.resid.plot(ax=axes[3])
axes[3].set_title('残差（噪声）')

plt.tight_layout()
plt.show()

# 季节成分强度
seasonal_strength = 1 - np.var(decomposition.resid.dropna()) / np.var(decomposition.seasonal.dropna() + decomposition.resid.dropna())
print(f'季节强度: {seasonal_strength:.3f}（越接近 1 季节性越强）')

fig, axes = plt.subplots(2, 2, figsize=(12, 8))

plot_acf(df['value'].dropna(), lags=50, ax=axes[0, 0])
axes[0, 0].set_title('原始序列 ACF')

plot_pacf(df['value'].dropna(), lags=50, method='ywm', ax=axes[0, 1])
axes[0, 1].set_title('原始序列 PACF')

plot_acf(df['diff'].dropna(), lags=50, ax=axes[1, 0])
axes[1, 0].set_title('差分后 ACF')

plot_pacf(df['diff'].dropna(), lags=50, method='ywm', ax=axes[1, 1])
axes[1, 1].set_title('差分后 PACF')

plt.tight_layout()
plt.show()

# 按年分割
train = df['value']['2023']
test = df['value']['2024']

# 朴素预测（用最后一个值预测所有未来值）
naive_forecast = pd.Series(train[-1], index=test.index)

# 移动平均预测
ma_forecast = train.rolling(30).mean().iloc[-1]
ma_forecast = pd.Series(ma_forecast, index=test.index)

# 可视化
plt.figure(figsize=(12, 5))
plt.plot(train.index, train, label='训练集')
plt.plot(test.index, test, label='测试集')
plt.plot(test.index, naive_forecast, '--', label='朴素预测')
plt.plot(test.index, ma_forecast, '--', label='30日移动平均')
plt.legend()
plt.title('简单预测方法对比')
plt.show()

# 评估
from sklearn.metrics import mean_absolute_error, mean_squared_error

print('朴素预测:')
print(f'  MAE: {mean_absolute_error(test, naive_forecast):.2f}')
print(f'  RMSE: {np.sqrt(mean_squared_error(test, naive_forecast)):.2f}')

print('移动平均预测:')
print(f'  MAE: {mean_absolute_error(test, ma_forecast):.2f}')
print(f'  RMSE: {np.sqrt(mean_squared_error(test, ma_forecast)):.2f}')

# 简单指数平滑（SES）
ses = SimpleExpSmoothing(train).fit(smoothing_level=0.3, optimized=False)
ses_forecast = ses.forecast(len(test))

# Holt 线性趋势
holt = ExponentialSmoothing(
    train,
    trend='add',
    seasonal=None
).fit()
holt_forecast = holt.forecast(len(test))

# Holt-Winters 季节性
hw = ExponentialSmoothing(
    train,
    trend='add',
    seasonal='add',
    seasonal_periods=365
).fit()
hw_forecast = hw.forecast(len(test))

# 对比
fig, axes = plt.subplots(2, 1, figsize=(12, 8))

axes[0].plot(train.index, train, label='训练集')
axes[0].plot(test.index, test, label='测试集', alpha=0.7)
axes[0].plot(test.index, ses_forecast, '--', label='SES (α=0.3)')
axes[0].plot(test.index, holt_forecast, '--', label='Holt')
axes[0].plot(test.index, hw_forecast, '--', label='Holt-Winters')
axes[0].legend()
axes[0].set_title('指数平滑方法对比')

# 残差
axes[1].plot(test.index, test - hw_forecast, label='Holt-Winters 残差')
axes[1].axhline(y=0, color='r', linestyle='--')
axes[1].set_title('预测残差')
axes[1].legend()

plt.tight_layout()
plt.show()

# 评估
models = {
    'SES (α=0.3)': ses_forecast,
    'Holt': holt_forecast,
    'Holt-Winters': hw_forecast
}

for name, fcast in models.items():
    mae = mean_absolute_error(test, fcast)
    rmse = np.sqrt(mean_squared_error(test, fcast))
    print(f'{name}: MAE={mae:.2f}, RMSE={rmse:.2f}')

# 手动选择参数
# p = PACF 显著滞后的数量（约 2）
# d = 差分阶数（1）
# q = ACF 显著滞后的数量（约 2）

model = ARIMA(train, order=(2, 1, 2))
arima_result = model.fit()

print(arima_result.summary())

arima_forecast = arima_result.forecast(len(test))

# 自动选择参数（使用信息准则）
from statsmodels.tsa.arima.model import ARIMA
import itertools

# 网格搜索最优参数
p = d = q = range(0, 4)
pdq = list(itertools.product(p, d, q))

best_aic = np.inf
best_pdq = None

for param in pdq:
    try:
        model = ARIMA(train, order=param)
        result = model.fit()
        if result.aic < best_aic:
            best_aic = result.aic
            best_pdq = param
    except:
        continue

print(f'最优参数 (p,d,q): {best_pdq}, AIC: {best_aic:.2f}')

# 用最优参数重新拟合
best_model = ARIMA(train, order=best_pdq)
best_result = best_model.fit()
best_forecast = best_result.forecast(len(test))

# 评估
print(f'ARIMA{best_pdq} 预测:')
print(f'  MAE: {mean_absolute_error(test, best_forecast):.2f}')
print(f'  RMSE: {np.sqrt(mean_squared_error(test, best_forecast)):.2f}')

# 使用 yfinance 获取股票数据（需先安装: pip install yfinance）
import yfinance as yf

# 下载中国平安股票数据
stock = yf.download('601318.SS', start='2023-01-01', end='2024-12-31')
stock.columns = ['_'.join(col).strip() for col in stock.columns.values]

print(stock.head())

close = stock['Close_601318.SS']

# 日收益率
returns = close.pct_change().dropna()

fig, axes = plt.subplots(2, 2, figsize=(12, 8))

# 收盘价
close.plot(ax=axes[0, 0])
axes[0, 0].set_title('中国平安 收盘价')

# 收益率
returns.plot(ax=axes[0, 1])
axes[0, 1].set_title('日收益率')
axes[0, 1].axhline(y=0, color='r', linestyle='--')

# 收益率分布
returns.hist(ax=axes[1, 0], bins=50)
axes[1, 0].set_title('收益率分布')

# Q-Q 图
from scipy import stats
stats.probplot(returns.dropna(), dist='norm', plot=axes[1, 1])
axes[1, 1].set_title('正态 Q-Q 图')

plt.tight_layout()
plt.show()

# 检查平稳性
check_stationarity(close, '收盘价')
check_stationarity(returns, '收益率（对数差分）')

# 用前 80% 的数据训练，后 20% 测试
split_idx = int(len(close) * 0.8)
train_stock = close.iloc[:split_idx]
test_stock = close.iloc[split_idx:]

# Holt-Winters 预测
hw_stock = ExponentialSmoothing(
    train_stock,
    trend='add',
    seasonal='add',
    seasonal_periods=5  # 周周期（5 个交易日）
).fit()

hw_stock_forecast = hw_stock.forecast(len(test_stock))

# 绘图
plt.figure(figsize=(12, 5))
plt.plot(train_stock.index, train_stock, label='训练集')
plt.plot(test_stock.index, test_stock, label='测试集')
plt.plot(test_stock.index, hw_stock_forecast, '--', label='Holt-Winters 预测')
plt.legend()
plt.title('股票价格预测')
plt.show()

print(f'预测 MAE: {mean_absolute_error(test_stock, hw_stock_forecast):.2f}')

# 使用 Seaborn 的航班数据集做月度分析
flights = sns.load_dataset('flights')
flights['year_month'] = pd.to_datetime(flights['year'].astype(str) + '-' + flights['month'] + '-01')
flights = flights.set_index('year_month')
flights = flights.sort_index()

passengers = flights['passengers']

# 分解
decomp = seasonal_decompose(passengers, model='additive', period=12)

fig, axes = plt.subplots(4, 1, figsize=(12, 8))
decomp.observed.plot(ax=axes[0], title='原始')
decomp.trend.plot(ax=axes[1], title='趋势')
decomp.seasonal.plot(ax=axes[2], title='季节')
decomp.resid.plot(ax=axes[3], title='残差')
plt.tight_layout()
plt.show()

# 预测
train_pass = passengers.iloc[:120]  # 前 10 年
test_pass = passengers.iloc[120:]  # 后 2 年

hw_pass = ExponentialSmoothing(
    train_pass,
    trend='add',
    seasonal='add',
    seasonal_periods=12
).fit()

hw_pass_forecast = hw_pass.forecast(len(test_pass))

plt.figure(figsize=(10, 4))
plt.plot(train_pass.index, train_pass, label='训练集')
plt.plot(test_pass.index, test_pass, label='测试集')
plt.plot(test_pass.index, hw_pass_forecast, '--', label='预测')
plt.legend()
plt.title('航班乘客数预测（Holt-Winters）')
plt.show()

print(f'Holt-Winters MAE: {mean_absolute_error(test_pass, hw_pass_forecast):.2f}')

def forecast_accuracy(forecast, actual):
    mae = mean_absolute_error(actual, forecast)
    mse = mean_squared_error(actual, forecast)
    rmse = np.sqrt(mse)
    mape = np.mean(np.abs((actual - forecast) / actual)) * 100
    
    # 方向准确率
    direction_actual = np.sign(np.diff(actual))
    direction_forecast = np.sign(np.diff(forecast.values))
    direction_acc = np.mean(direction_actual == direction_forecast) * 100
    
    return {
        'MAE': mae,
        'MSE': mse,
        'RMSE': rmse,
        'MAPE(%)': mape,
        '方向准确率(%)': direction_acc
    }

# 评估所有模型
all_forecasts = {
    '朴素': naive_forecast,
    '移动平均': ma_forecast,
    'SES': ses_forecast,
    'Holt': holt_forecast,
    'Holt-Winters': hw_forecast,
    'ARIMA': best_forecast,
}

results = []
for name, fcast in all_forecasts.items():
    metrics = forecast_accuracy(fcast, test.values)
    metrics['模型'] = name
    results.append(metrics)

df_results = pd.DataFrame(results)
df_results = df_results.set_index('模型')
print(df_results.round(2))

陷阱 = '''
1. 过拟合：不要用复杂的模型去拟合噪声
2. 泄漏：预测时必须只用已知信息，不能"看到未来"
3. 伪相关：两个不相关的时间序列也可能出现高相关
4. 非平稳：忘记处理趋势和季节，模型会失效
5. 评估方式错误：时间序列不能随机分割，要按时间顺序
6. 忽视外部因素：模型只看到历史模式，不知道新冠、政策变化
7. 预测区间：点预测不够，要给出置信区间
'''
print(陷阱)