Time Series Analysis by aj-geddes/useful-ai-prompts
npx skills add https://github.com/aj-geddes/useful-ai-prompts --skill 'Time Series Analysis'时间序列分析通过检查随时间收集的数据点来识别模式、趋势和季节性,以进行预测并理解时间动态。
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from statsmodels.tsa.seasonal import seasonal_decompose
from statsmodels.tsa.stattools import adfuller, acf, pacf
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.tsa.holtwinters import ExponentialSmoothing
# 创建示例时间序列数据
dates = pd.date_range('2020-01-01', periods=365, freq='D')
values = 100 + np.sin(np.arange(365) * 2*np.pi / 365) * 20 + np.random.normal(0, 5, 365)
ts = pd.Series(values, index=dates)
# 可视化时间序列
fig, axes = plt.subplots(2, 2, figsize=(14, 8))
axes[0, 0].plot(ts)
axes[0, 0].set_title('原始时间序列')
axes[0, 0].set_ylabel('值')
# 分解
decomposition = seasonal_decompose(ts, model='additive', period=30)
axes[0, 1].plot(decomposition.trend)
axes[0, 1].set_title('趋势分量')
axes[1, 0].plot(decomposition.seasonal)
axes[1, 0].set_title('季节性分量')
axes[1, 1].plot(decomposition.resid)
axes[1, 1].set_title('残差分分量')
plt.tight_layout()
plt.show()
# 平稳性检验(增强迪基-福勒检验)
result = adfuller(ts)
print(f"ADF 检验统计量: {result[0]:.6f}")
print(f"P 值: {result[1]:.6f}")
print(f"临界值: {result[4]}")
if result[1] <= 0.05:
print("时间序列是平稳的")
else:
print("时间序列是非平稳的 - 需要进行差分")
# 一阶差分以实现平稳性
ts_diff = ts.diff().dropna()
result_diff = adfuller(ts_diff)
print(f"\n差分后 - ADF p 值: {result_diff[1]:.6f}")
# 自相关和偏自相关
fig, axes = plt.subplots(1, 2, figsize=(12, 4))
plot_acf(ts_diff, lags=40, ax=axes[0])
axes[0].set_title('ACF')
plot_pacf(ts_diff, lags=40, ax=axes[1])
axes[1].set_title('PACF')
plt.tight_layout()
plt.show()
# ARIMA 模型
arima_model = ARIMA(ts, order=(1, 1, 1))
arima_result = arima_model.fit()
print(arima_result.summary())
# 预测
forecast_steps = 30
forecast = arima_result.get_forecast(steps=forecast_steps)
forecast_df = forecast.conf_int()
forecast_mean = forecast.predicted_mean
# 绘制预测图
fig, ax = plt.subplots(figsize=(12, 5))
ax.plot(ts.index[-90:], ts[-90:], label='历史数据')
ax.plot(forecast_df.index, forecast_mean, label='预测', color='red')
ax.fill_between(
forecast_df.index,
forecast_df.iloc[:, 0],
forecast_df.iloc[:, 1],
color='red', alpha=0.2
)
ax.set_title('带置信区间的 ARIMA 预测')
ax.legend()
ax.grid(True, alpha=0.3)
plt.show()
# 指数平滑
exp_smooth = ExponentialSmoothing(
ts, seasonal_periods=30, trend='add', seasonal='add', initialization_method='estimated'
)
exp_result = exp_smooth.fit()
# 模型诊断
fig = exp_result.plot_diagnostics(figsize=(12, 8))
plt.tight_layout()
plt.show()
# 自定义移动平均分析
window_sizes = [7, 30, 90]
fig, ax = plt.subplots(figsize=(12, 5))
ax.plot(ts.index, ts.values, label='原始数据', alpha=0.7)
for window in window_sizes:
ma = ts.rolling(window=window).mean()
ax.plot(ma.index, ma.values, label=f'MA({window})')
ax.set_title('移动平均')
ax.legend()
ax.grid(True, alpha=0.3)
plt.show()
# 季节性子序列图
fig, axes = plt.subplots(2, 2, figsize=(12, 8))
for i, month in enumerate(range(1, 5)):
month_data = ts[ts.index.month == month]
axes[i // 2, i % 2].plot(month_data.values)
axes[i // 2, i % 2].set_title(f'月份 {month} 模式')
plt.tight_layout()
plt.show()
# 预测准确性指标
def calculate_forecast_metrics(actual, predicted):
mae = np.mean(np.abs(actual - predicted))
rmse = np.sqrt(np.mean((actual - predicted) ** 2))
mape = np.mean(np.abs((actual - predicted) / actual)) * 100
return {'MAE': mae, 'RMSE': rmse, 'MAPE': mape}
metrics = calculate_forecast_metrics(ts[-30:], forecast_mean[:30])
print(f"\n预测指标:\n{metrics}")
# 附加分析技术
# 步骤 10: 季节性子序列图
fig, axes = plt.subplots(2, 2, figsize=(12, 8))
for i, season in enumerate([1, 2, 3, 4]):
seasonal_ts = ts[ts.index.month % 4 == season % 4]
axes[i // 2, i % 2].plot(seasonal_ts.values)
axes[i // 2, i % 2].set_title(f'季节 {season}')
plt.tight_layout()
plt.show()
# 步骤 11: 格兰杰因果检验(用于多序列)
from statsmodels.tsa.stattools import grangercausalitytests
# 创建另一个序列用于测试
ts2 = ts.shift(1).fillna(method='bfill')
try:
print("\n格兰杰因果检验:")
print(f"检验 ts2 是否格兰杰导致 ts:")
gc_result = grangercausalitytests(np.column_stack([ts.values, ts2.values]), maxlag=3)
except Exception as e:
print(f"未执行格兰杰因果检验: {str(e)[:50]}")
# 步骤 12: 自相关和偏自相关分析
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
acf_values = acf(ts.dropna(), nlags=20)
pacf_values = pacf(ts.dropna(), nlags=20)
# 步骤 13: 季节性强度
def seasonal_strength(series, seasonal_period=30):
seasonal = seasonal_decompose(series, model='additive', period=seasonal_period)
var_residual = np.var(seasonal.resid.dropna())
var_seasonal = np.var(seasonal.seasonal)
return 1 - (var_residual / (var_residual + var_seasonal)) if (var_residual + var_seasonal) > 0 else 0
ss = seasonal_strength(ts)
print(f"\n季节性强度: {ss:.3f}")
# 步骤 14: 带不确定性的预测
fig, ax = plt.subplots(figsize=(12, 5))
ax.plot(ts.index[-60:], ts.values[-60:], label='历史数据', linewidth=2)
# 多步长预测
for steps_ahead in [10, 20, 30]:
try:
fc = arima_result.get_forecast(steps=steps_ahead)
fc_mean = fc.predicted_mean
ax.plot(pd.date_range(ts.index[-1], periods=steps_ahead+1)[1:],
fc_mean.values, marker='o', label=f'预测 (+{steps_ahead})')
except:
pass
ax.set_title('多步长预测')
ax.set_xlabel('日期')
ax.set_ylabel('值')
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
# 步骤 15: 模型比较总结
print("\n时间序列分析完成!")
print(f"原始序列长度: {len(ts)}")
print(f"趋势强度: {1 - np.var(decomposition.resid.dropna()) / np.var((ts - ts.mean()).dropna()):.3f}")
print(f"季节性强度: {ss:.3f}")
广告位招租
在这里展示您的产品或服务
触达数万 AI 开发者,精准高效
每周安装量
0
代码仓库
GitHub 星标数
121
首次出现时间
1970年1月1日
安全审计
Time series analysis examines data points collected over time to identify patterns, trends, and seasonality for forecasting and understanding temporal dynamics.
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from statsmodels.tsa.seasonal import seasonal_decompose
from statsmodels.tsa.stattools import adfuller, acf, pacf
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.tsa.holtwinters import ExponentialSmoothing
# Create sample time series data
dates = pd.date_range('2020-01-01', periods=365, freq='D')
values = 100 + np.sin(np.arange(365) * 2*np.pi / 365) * 20 + np.random.normal(0, 5, 365)
ts = pd.Series(values, index=dates)
# Visualize time series
fig, axes = plt.subplots(2, 2, figsize=(14, 8))
axes[0, 0].plot(ts)
axes[0, 0].set_title('Original Time Series')
axes[0, 0].set_ylabel('Value')
# Decomposition
decomposition = seasonal_decompose(ts, model='additive', period=30)
axes[0, 1].plot(decomposition.trend)
axes[0, 1].set_title('Trend Component')
axes[1, 0].plot(decomposition.seasonal)
axes[1, 0].set_title('Seasonal Component')
axes[1, 1].plot(decomposition.resid)
axes[1, 1].set_title('Residual Component')
plt.tight_layout()
plt.show()
# Test for stationarity (Augmented Dickey-Fuller)
result = adfuller(ts)
print(f"ADF Test Statistic: {result[0]:.6f}")
print(f"P-value: {result[1]:.6f}")
print(f"Critical Values: {result[4]}")
if result[1] <= 0.05:
print("Time series is stationary")
else:
print("Time series is non-stationary - differencing needed")
# First differencing for stationarity
ts_diff = ts.diff().dropna()
result_diff = adfuller(ts_diff)
print(f"\nAfter differencing - ADF p-value: {result_diff[1]:.6f}")
# Autocorrelation and Partial Autocorrelation
fig, axes = plt.subplots(1, 2, figsize=(12, 4))
plot_acf(ts_diff, lags=40, ax=axes[0])
axes[0].set_title('ACF')
plot_pacf(ts_diff, lags=40, ax=axes[1])
axes[1].set_title('PACF')
plt.tight_layout()
plt.show()
# ARIMA Model
arima_model = ARIMA(ts, order=(1, 1, 1))
arima_result = arima_model.fit()
print(arima_result.summary())
# Forecast
forecast_steps = 30
forecast = arima_result.get_forecast(steps=forecast_steps)
forecast_df = forecast.conf_int()
forecast_mean = forecast.predicted_mean
# Plot forecast
fig, ax = plt.subplots(figsize=(12, 5))
ax.plot(ts.index[-90:], ts[-90:], label='Historical')
ax.plot(forecast_df.index, forecast_mean, label='Forecast', color='red')
ax.fill_between(
forecast_df.index,
forecast_df.iloc[:, 0],
forecast_df.iloc[:, 1],
color='red', alpha=0.2
)
ax.set_title('ARIMA Forecast with Confidence Interval')
ax.legend()
ax.grid(True, alpha=0.3)
plt.show()
# Exponential Smoothing
exp_smooth = ExponentialSmoothing(
ts, seasonal_periods=30, trend='add', seasonal='add', initialization_method='estimated'
)
exp_result = exp_smooth.fit()
# Model diagnostics
fig = exp_result.plot_diagnostics(figsize=(12, 8))
plt.tight_layout()
plt.show()
# Custom moving average analysis
window_sizes = [7, 30, 90]
fig, ax = plt.subplots(figsize=(12, 5))
ax.plot(ts.index, ts.values, label='Original', alpha=0.7)
for window in window_sizes:
ma = ts.rolling(window=window).mean()
ax.plot(ma.index, ma.values, label=f'MA({window})')
ax.set_title('Moving Averages')
ax.legend()
ax.grid(True, alpha=0.3)
plt.show()
# Seasonal subseries plot
fig, axes = plt.subplots(2, 2, figsize=(12, 8))
for i, month in enumerate(range(1, 5)):
month_data = ts[ts.index.month == month]
axes[i // 2, i % 2].plot(month_data.values)
axes[i // 2, i % 2].set_title(f'Month {month} Pattern')
plt.tight_layout()
plt.show()
# Forecast accuracy metrics
def calculate_forecast_metrics(actual, predicted):
mae = np.mean(np.abs(actual - predicted))
rmse = np.sqrt(np.mean((actual - predicted) ** 2))
mape = np.mean(np.abs((actual - predicted) / actual)) * 100
return {'MAE': mae, 'RMSE': rmse, 'MAPE': mape}
metrics = calculate_forecast_metrics(ts[-30:], forecast_mean[:30])
print(f"\nForecast Metrics:\n{metrics}")
# Additional analysis techniques
# Step 10: Seasonal subseries plots
fig, axes = plt.subplots(2, 2, figsize=(12, 8))
for i, season in enumerate([1, 2, 3, 4]):
seasonal_ts = ts[ts.index.month % 4 == season % 4]
axes[i // 2, i % 2].plot(seasonal_ts.values)
axes[i // 2, i % 2].set_title(f'Season {season}')
plt.tight_layout()
plt.show()
# Step 11: Granger causality (for multiple series)
from statsmodels.tsa.stattools import grangercausalitytests
# Create another series for testing
ts2 = ts.shift(1).fillna(method='bfill')
try:
print("\nGranger Causality Test:")
print(f"Test whether ts2 Granger-causes ts:")
gc_result = grangercausalitytests(np.column_stack([ts.values, ts2.values]), maxlag=3)
except Exception as e:
print(f"Granger causality not performed: {str(e)[:50]}")
# Step 12: Autocorrelation and partial autocorrelation analysis
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
acf_values = acf(ts.dropna(), nlags=20)
pacf_values = pacf(ts.dropna(), nlags=20)
# Step 13: Seasonal strength
def seasonal_strength(series, seasonal_period=30):
seasonal = seasonal_decompose(series, model='additive', period=seasonal_period)
var_residual = np.var(seasonal.resid.dropna())
var_seasonal = np.var(seasonal.seasonal)
return 1 - (var_residual / (var_residual + var_seasonal)) if (var_residual + var_seasonal) > 0 else 0
ss = seasonal_strength(ts)
print(f"\nSeasonal Strength: {ss:.3f}")
# Step 14: Forecasting with uncertainty
fig, ax = plt.subplots(figsize=(12, 5))
ax.plot(ts.index[-60:], ts.values[-60:], label='Historical', linewidth=2)
# Multiple horizon forecasts
for steps_ahead in [10, 20, 30]:
try:
fc = arima_result.get_forecast(steps=steps_ahead)
fc_mean = fc.predicted_mean
ax.plot(pd.date_range(ts.index[-1], periods=steps_ahead+1)[1:],
fc_mean.values, marker='o', label=f'Forecast (+{steps_ahead})')
except:
pass
ax.set_title('Multi-step Ahead Forecasts')
ax.set_xlabel('Date')
ax.set_ylabel('Value')
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
# Step 15: Model comparison summary
print("\nTime Series Analysis Complete!")
print(f"Original series length: {len(ts)}")
print(f"Trend strength: {1 - np.var(decomposition.resid.dropna()) / np.var((ts - ts.mean()).dropna()):.3f}")
print(f"Seasonal strength: {ss:.3f}")
Weekly Installs
0
Repository
GitHub Stars
121
First Seen
Jan 1, 1970
Security Audits
专业SEO审计工具:全面网站诊断、技术SEO优化与页面分析指南
62,600 周安装