diff --git a/aapl_backtest.png b/aapl_backtest.png new file mode 100644 index 0000000..ba4b4ce Binary files /dev/null and b/aapl_backtest.png differ diff --git a/aapl_backtest_optimized.png b/aapl_backtest_optimized.png new file mode 100644 index 0000000..e1f2edf Binary files /dev/null and b/aapl_backtest_optimized.png differ diff --git a/aapl_backtest_trick.png b/aapl_backtest_trick.png new file mode 100644 index 0000000..f728be3 Binary files /dev/null and b/aapl_backtest_trick.png differ diff --git a/aapl_backtest_xreg.png b/aapl_backtest_xreg.png new file mode 100644 index 0000000..01155d4 Binary files /dev/null and b/aapl_backtest_xreg.png differ diff --git a/aapl_forecast.png b/aapl_forecast.png new file mode 100644 index 0000000..13d4c62 Binary files /dev/null and b/aapl_forecast.png differ diff --git a/stock_backtest.py b/stock_backtest.py new file mode 100644 index 0000000..eef058a --- /dev/null +++ b/stock_backtest.py @@ -0,0 +1,115 @@ +import yfinance as yf +import numpy as np +import torch +import timesfm +import matplotlib +matplotlib.use("Agg") +import matplotlib.pyplot as plt +from datetime import datetime, timedelta + +# 1. 获取 AAPL 数据（到今天为止） +print("下载 AAPL 股票数据...", flush=True) +end = datetime.now() +start = end - timedelta(days=365) +df = yf.download("AAPL", start=start.strftime("%Y-%m-%d"), end=end.strftime("%Y-%m-%d"), progress=False) +close = df["Close"].values.flatten().astype(np.float32) +dates = df.index + +# 2. 分割：6/15 之前为训练集，6/16 之后为实际值 +split_date = "2026-06-15" +split_idx = None +for i, d in enumerate(dates): + if str(d.date()) <= split_date: + split_idx = i +# split_idx 是 <= 6/15 的最后一个索引 +train_data = close[: split_idx + 1] +train_dates = dates[: split_idx + 1] +actual_data = close[split_idx + 1 :] +actual_dates = dates[split_idx + 1 :] + +HORIZON = len(actual_data) +print(f"训练数据: {len(train_data)} 天（截至 {train_dates[-1].date()}）", flush=True) +print(f"实际数据: {HORIZON} 天（{actual_dates[0].date()} ~ {actual_dates[-1].date()}）", flush=True) +print(f"分割日收盘价: {train_data[-1]:.2f}", flush=True) + +# 3. 加载模型 +print() +print("加载 TimesFM 模型...", flush=True) +torch.set_float32_matmul_precision("high") +model = timesfm.TimesFM_2p5_200M_torch.from_pretrained( + "google/timesfm-2.5-200m-pytorch", torch_compile=False +) +print("模型加载完成", flush=True) + +# 4. 编译并预测 +model.compile( + timesfm.ForecastConfig( + max_context=512, + max_horizon=128, + normalize_inputs=True, + use_continuous_quantile_head=True, + force_flip_invariance=False, + infer_is_positive=True, + fix_quantile_crossing=True, + ) +) +print(f"预测未来 {HORIZON} 个交易日...", flush=True) +pf, qf = model.forecast(horizon=HORIZON, inputs=[train_data]) +print("预测完成！", flush=True) + +# 5. 计算误差指标 +actual = actual_data +pred = pf[0][:HORIZON] +mae = np.mean(np.abs(actual - pred)) +rmse = np.sqrt(np.mean((actual - pred) ** 2)) +mape = np.mean(np.abs((actual - pred) / actual)) * 100 + +# 方向准确率 +actual_dir = np.diff(actual) +pred_dir = np.diff(pred) +dir_acc = np.mean(actual_dir * pred_dir > 0) * 100 + +# 置信区间覆盖率 +in_80 = np.mean((actual >= qf[0, :HORIZON, 1]) & (actual <= qf[0, :HORIZON, 9])) * 100 +in_40 = np.mean((actual >= qf[0, :HORIZON, 3]) & (actual <= qf[0, :HORIZON, 7])) * 100 + +print() +print("=== 回测结果：6/16 ~ 6/30 预测 vs 实际 ===", flush=True) +print() +print(" 日期 实际价 预测价 误差 误差%", flush=True) +for i in range(HORIZON): + err = pred[i] - actual[i] + err_pct = err / actual[i] * 100 + print( + f" {actual_dates[i].date()} {actual[i]:7.2f} {pred[i]:7.2f} {err:+7.2f} {err_pct:+6.2f}%", + flush=True, + ) + +print() +print("=== 误差指标 ===", flush=True) +print(f" MAE (平均绝对误差): ${mae:.2f}", flush=True) +print(f" RMSE (均方根误差): ${rmse:.2f}", flush=True) +print(f" MAPE (平均绝对百分比误差): {mape:.2f}%", flush=True) +print(f" 方向准确率: {dir_acc:.1f}%", flush=True) +print(f" 80%置信区间覆盖率: {in_80:.1f}%", flush=True) +print(f" 40%置信区间覆盖率: {in_40:.1f}%", flush=True) + +# 6. 画图 +fig, ax = plt.subplots(figsize=(14, 6)) +show_n = min(60, len(train_data)) +ax.plot(range(show_n), train_data[-show_n:], label="Historical Close", color="steelblue", linewidth=1.5) +x_actual = range(show_n, show_n + HORIZON) +ax.plot(x_actual, actual, label="Actual", color="forestgreen", linewidth=2, marker="o", markersize=3) +ax.plot(x_actual, pred, label="Forecast", color="tomato", linewidth=2, linestyle="--") +ax.fill_between(x_actual, qf[0, :HORIZON, 1], qf[0, :HORIZON, 9], alpha=0.15, color="tomato", label="80% CI") +ax.fill_between(x_actual, qf[0, :HORIZON, 3], qf[0, :HORIZON, 7], alpha=0.3, color="tomato", label="40% CI") +ax.axvline(x=show_n - 1, color="gray", linestyle="--", alpha=0.5, label="Forecast Start") +ax.set_title(f"AAPL Backtest: Forecast vs Actual (MAPE={mape:.2f}%, Dir.Acc={dir_acc:.0f}%)", fontsize=14) +ax.set_xlabel("Trading Days") +ax.set_ylabel("Price (USD)") +ax.legend(loc="upper left") +plt.tight_layout() +plt.savefig("aapl_backtest.png", dpi=150) +print() +print("Chart saved: aapl_backtest.png", flush=True) +print("Done!", flush=True) diff --git a/stock_backtest_optimized.py b/stock_backtest_optimized.py new file mode 100644 index 0000000..c4e013f --- /dev/null +++ b/stock_backtest_optimized.py @@ -0,0 +1,184 @@ +import yfinance as yf +import numpy as np +import torch +import timesfm +import matplotlib +matplotlib.use("Agg") +import matplotlib.pyplot as plt +from datetime import datetime, timedelta + +# ---------- 技术指标计算 ---------- +def calc_rsi(prices, period=14): + deltas = np.diff(prices) + gains = np.where(deltas > 0, deltas, 0.0) + losses = np.where(deltas < 0, -deltas, 0.0) + avg_gain = np.convolve(gains, np.ones(period) / period, mode="valid") + avg_loss = np.convolve(losses, np.ones(period) / period, mode="valid") + avg_loss = np.where(avg_loss == 0, 1e-10, avg_loss) + rs = avg_gain / avg_loss + rsi = 100.0 - (100.0 / (1.0 + rs)) + # pad front to match length + return np.concatenate([np.full(period, 50.0), rsi]) + + +def calc_sma(prices, period): + sma = np.convolve(prices, np.ones(period) / period, mode="valid") + return np.concatenate([np.full(period - 1, sma[0] if len(sma) > 0 else 0.0), sma]) + + +# ---------- 1. 获取数据 ---------- +print("下载 AAPL 股票数据...", flush=True) +end = datetime.now() +start = end - timedelta(days=365) +df = yf.download("AAPL", start=start.strftime("%Y-%m-%d"), end=end.strftime("%Y-%m-%d"), progress=False) +close = df["Close"].values.flatten().astype(np.float32) +volume = df["Volume"].values.flatten().astype(np.float32) +dates = df.index + +# 分割 +split_date = "2026-06-15" +split_idx = None +for i, d in enumerate(dates): + if str(d.date()) <= split_date: + split_idx = i + +train_close = close[: split_idx + 1] +train_vol = volume[: split_idx + 1] +train_dates = dates[: split_idx + 1] +actual_close = close[split_idx + 1 :] +actual_dates = dates[split_idx + 1 :] +HORIZON = len(actual_close) + +print(f"训练数据: {len(train_close)} 天（截至 {train_dates[-1].date()}）", flush=True) +print(f"实际数据: {HORIZON} 天（{actual_dates[0].date()} ~ {actual_dates[-1].date()}）", flush=True) + +# ---------- 2. 计算对数收益率 ---------- +# log returns = ln(P_t / P_{t-1})，更平稳 +train_logret = np.diff(np.log(train_close)).astype(np.float32) # length = N-1 +actual_logret = np.diff(np.log(np.concatenate([train_close[-1:], actual_close]))).astype(np.float32) + +print(f"对数收益率: 均值={train_logret.mean():.6f}, 标准差={train_logret.std():.6f}", flush=True) + +# ---------- 3. 加载模型 ---------- +print() +print("加载 TimesFM 模型...", flush=True) +torch.set_float32_matmul_precision("high") +model = timesfm.TimesFM_2p5_200M_torch.from_pretrained( + "google/timesfm-2.5-200m-pytorch", torch_compile=False +) +print("模型加载完成", flush=True) + +# ---------- 4. 多窗口集成预测（对数收益率） ---------- +CONTEXTS = [128, 256, 512] # 不同 context 长度 +all_preds = [] +all_quantiles = [] + +for ctx_len in CONTEXTS: + ctx_data = train_logret[-ctx_len:] if len(train_logret) >= ctx_len else train_logret + actual_ctx = min(ctx_len, len(ctx_data)) + # 对收益率: infer_is_positive=False (可负), normalize_inputs=True + model.compile( + timesfm.ForecastConfig( + max_context=actual_ctx, + max_horizon=128, + normalize_inputs=True, + use_continuous_quantile_head=True, + force_flip_invariance=False, + infer_is_positive=False, # 收益率可正可负 + fix_quantile_crossing=True, + ) + ) + pf, qf = model.forecast(horizon=HORIZON, inputs=[ctx_data]) + all_preds.append(pf[0][:HORIZON]) + all_quantiles.append(qf[0][:HORIZON]) + print(f" context={actual_ctx} 预测完成", flush=True) + +# 集成: 取平均 +ensemble_pred_logret = np.mean(all_preds, axis=0) +ensemble_q_logret = np.mean(all_quantiles, axis=0) + +# ---------- 5. 转换回价格 ---------- +# P_t = P_{t-1} * exp(r_t) +last_price = train_close[-1] +pred_prices = [] +for i in range(HORIZON): + last_price = last_price * np.exp(ensemble_pred_logret[i]) + pred_prices.append(last_price) +pred_prices = np.array(pred_prices) + +# 分位数价格 +q_prices = np.zeros((HORIZON, 10)) +for qi in range(10): + p = train_close[-1] + for i in range(HORIZON): + p = p * np.exp(ensemble_q_logret[i, qi]) + q_prices[i, qi] = p + +# ---------- 6. 计算误差 ---------- +actual = actual_close +pred = pred_prices +mae = np.mean(np.abs(actual - pred)) +rmse = np.sqrt(np.mean((actual - pred) ** 2)) +mape = np.mean(np.abs((actual - pred) / actual)) * 100 + +actual_dir = np.diff(actual) +pred_dir = np.diff(pred) +dir_acc = np.mean(actual_dir * pred_dir > 0) * 100 + +in_80 = np.mean((actual >= q_prices[:, 1]) & (actual <= q_prices[:, 9])) * 100 +in_40 = np.mean((actual >= q_prices[:, 3]) & (actual <= q_prices[:, 7])) * 100 + +print() +print("=== 优化版回测结果：6/16 ~ 6/30 ===", flush=True) +print() +print(" 日期 实际价 预测价 误差 误差%", flush=True) +for i in range(HORIZON): + err = pred[i] - actual[i] + err_pct = err / actual[i] * 100 + print( + f" {actual_dates[i].date()} {actual[i]:7.2f} {pred[i]:7.2f} {err:+7.2f} {err_pct:+6.2f}%", + flush=True, + ) + +print() +print("=== 误差指标（优化版 vs 原始版）===", flush=True) +print(f" MAE: ${mae:.2f} (原始: $8.46)", flush=True) +print(f" RMSE: ${rmse:.2f} (原始: $11.29)", flush=True) +print(f" MAPE: {mape:.2f}% (原始: 2.98%)", flush=True) +print(f" 方向准确率: {dir_acc:.1f}% (原始: 55.6%)", flush=True) +print(f" 80%CI覆盖率: {in_80:.1f}% (原始: 90.0%)", flush=True) +print(f" 40%CI覆盖率: {in_40:.1f}% (原始: 60.0%)", flush=True) + +# ---------- 7. 画对比图 ---------- +fig, axes = plt.subplots(2, 1, figsize=(14, 10), sharex=False) + +# 上图: 价格对比 +ax = axes[0] +show_n = min(60, len(train_close)) +ax.plot(range(show_n), train_close[-show_n:], label="Historical", color="steelblue", linewidth=1.5) +x_actual = range(show_n, show_n + HORIZON) +ax.plot(x_actual, actual, label="Actual", color="forestgreen", linewidth=2, marker="o", markersize=4) +ax.plot(x_actual, pred, label="Optimized Forecast", color="tomato", linewidth=2, linestyle="--") +ax.fill_between(x_actual, q_prices[:, 1], q_prices[:, 9], alpha=0.15, color="tomato", label="80% CI") +ax.fill_between(x_actual, q_prices[:, 3], q_prices[:, 7], alpha=0.3, color="tomato", label="40% CI") +ax.axvline(x=show_n - 1, color="gray", linestyle="--", alpha=0.5) +ax.set_title(f"Optimized: Log-Return + Ensemble (MAPE={mape:.2f}%, Dir={dir_acc:.0f}%)", fontsize=13) +ax.set_ylabel("Price (USD)") +ax.legend(loc="upper left") + +# 下图: 预测误差对比 +ax2 = axes[1] +errors = pred - actual +ax2.bar(range(HORIZON), errors, color=["tomato" if e > 0 else "steelblue" for e in errors], alpha=0.7) +ax2.axhline(y=0, color="black", linewidth=0.8) +ax2.set_title("Forecast Error per Day (Optimized)", fontsize=13) +ax2.set_xlabel("Trading Days After Split") +ax2.set_ylabel("Error (USD)") +ax2.set_xticks(range(HORIZON)) +ax2.set_xticklabels([str(d.date()) for d in actual_dates], rotation=45, fontsize=8) + +plt.tight_layout() +plt.savefig("aapl_backtest_optimized.png", dpi=150) +print() +print("Chart saved: aapl_backtest_optimized.png", flush=True) +print("Done!", flush=True) diff --git a/stock_backtest_trick.py b/stock_backtest_trick.py new file mode 100644 index 0000000..c0cfd1d --- /dev/null +++ b/stock_backtest_trick.py @@ -0,0 +1,217 @@ +"""TimesFM 纯技巧优化： +1. 对数收益率预测（更平稳） +2. 多窗口集成（128/256/512） +3. SPY 大盘走势引导：先预测 SPY，用 SPY 预测的趋势辅助判断 AAPL 方向 +4. 波动率调整：用近期波动率缩放置信区间 +5. infer_is_positive=False（收益率可负） +""" +import yfinance as yf +import numpy as np +import torch +import timesfm +import matplotlib +matplotlib.use("Agg") +import matplotlib.pyplot as plt +from datetime import datetime, timedelta + + +# 1. 获取数据 +print("下载 AAPL + SPY 数据...", flush=True) +end = datetime.now() +start = end - timedelta(days=365) +df_aapl = yf.download("AAPL", start=start.strftime("%Y-%m-%d"), end=end.strftime("%Y-%m-%d"), progress=False) +df_spy = yf.download("SPY", start=start.strftime("%Y-%m-%d"), end=end.strftime("%Y-%m-%d"), progress=False) + +close = df_aapl["Close"].values.flatten().astype(np.float32) +spy_close = df_spy["Close"].values.flatten().astype(np.float32) +dates = df_aapl.index + +min_len = min(len(close), len(spy_close)) +close = close[-min_len:] +spy_close = spy_close[-min_len:] +dates = dates[-min_len:] + +# 分割 +split_date = "2026-06-15" +split_idx = None +for i, d in enumerate(dates): + if str(d.date()) <= split_date: + split_idx = i + +train_close = close[: split_idx + 1] +train_spy = spy_close[: split_idx + 1] +train_dates = dates[: split_idx + 1] +actual_close = close[split_idx + 1 :] +actual_spy = spy_close[split_idx + 1 :] +actual_dates = dates[split_idx + 1 :] +HORIZON = len(actual_close) + +print(f"训练数据: {len(train_close)} 天 | 预测: {HORIZON} 天", flush=True) + +# 2. 对数收益率 +train_logret = np.diff(np.log(train_close)).astype(np.float32) +train_spy_logret = np.diff(np.log(train_spy)).astype(np.float32) +actual_logret = np.diff(np.log(np.concatenate([train_close[-1:], actual_close]))).astype(np.float32) + +# 3. 加载模型 +print() +print("加载 TimesFM 模型...", flush=True) +torch.set_float32_matmul_precision("high") +model = timesfm.TimesFM_2p5_200M_torch.from_pretrained( + "google/timesfm-2.5-200m-pytorch", torch_compile=False +) +print("模型加载完成", flush=True) + +# 4. 多窗口集成预测 AAPL 收益率 +CONTEXTS = [128, 256, 512] +all_preds = [] +all_quantiles = [] + +for ctx_len in CONTEXTS: + ctx_data = train_logret[-ctx_len:] if len(train_logret) >= ctx_len else train_logret + actual_ctx = min(ctx_len, len(ctx_data)) + # 向上取整到 32 的倍数 + actual_ctx = ((actual_ctx + 31) // 32) * 32 + ctx_data = train_logret[-actual_ctx:] + + model.compile( + timesfm.ForecastConfig( + max_context=actual_ctx, + max_horizon=128, + normalize_inputs=True, + use_continuous_quantile_head=True, + force_flip_invariance=False, + infer_is_positive=False, + fix_quantile_crossing=True, + ) + ) + pf, qf = model.forecast(horizon=HORIZON, inputs=[ctx_data]) + all_preds.append(pf[0][:HORIZON]) + all_quantiles.append(qf[0][:HORIZON]) + print(f" AAPL context={actual_ctx} done", flush=True) + +# 5. 预测 SPY 收益率（大盘趋势引导） +spy_ctx = train_spy_logret[-256:] +spy_ctx_len = ((len(spy_ctx) + 31) // 32) * 32 +spy_ctx = train_spy_logret[-spy_ctx_len:] +model.compile( + timesfm.ForecastConfig( + max_context=spy_ctx_len, + max_horizon=128, + normalize_inputs=True, + use_continuous_quantile_head=True, + force_flip_invariance=False, + infer_is_positive=False, + fix_quantile_crossing=True, + ) +) +spy_pf, _ = model.forecast(horizon=HORIZON, inputs=[spy_ctx]) +print(f" SPY context={spy_ctx_len} done", flush=True) + +# 6. 集成 + SPY 趋势调整 +ensemble_pred_logret = np.mean(all_preds, axis=0) +ensemble_q_logret = np.mean(all_quantiles, axis=0) +spy_pred_logret = spy_pf[0][:HORIZON] + +# SPY 趋势调整：如果 SPY 预测下跌，对 AAPL 预测施加向下的调整 +# 计算 AAPL 对 SPY 的 beta（敏感度） +beta = np.corrcoef(train_logret[-60:], train_spy_logret[-60:])[0, 1] +print(f" AAPL-SPY 60日相关系数: {beta:.3f}", flush=True) + +# 调整：将 SPY 预测的偏离均值部分 * beta 加到 AAPL 预测上 +spy_mean = np.mean(train_spy_logret[-60:]) +spy_deviation = spy_pred_logret - spy_mean # SPY 偏离其均值的部分 +adjustment = beta * spy_deviation * 0.3 # 0.3 是调整强度，避免过度修正 +adjusted_pred_logret = ensemble_pred_logret + adjustment + +# 7. 波动率调整置信区间 +recent_vol = np.std(train_logret[-20:]) +long_vol = np.std(train_logret[-60:]) +print(f" 近20日波动率: {recent_vol:.5f} | 近60日波动率: {long_vol:.5f}", flush=True) + +# 如果近期波动率高于长期，扩大置信区间 +vol_ratio = recent_vol / max(long_vol, 1e-8) +vol_scale = max(vol_ratio, 1.0) # 只扩大不缩小 +adjusted_q_logret = ensemble_q_logret.copy() +median_idx = 5 +for qi in range(10): + if qi != median_idx: + adjusted_q_logret[:, qi] = ensemble_q_logret[:, median_idx] + ( + ensemble_q_logret[:, qi] - ensemble_q_logret[:, median_idx] + ) * vol_scale + +# 8. 转换回价格 +last_price = train_close[-1] +pred_prices = [] +for i in range(HORIZON): + last_price = last_price * np.exp(adjusted_pred_logret[i]) + pred_prices.append(last_price) +pred_prices = np.array(pred_prices) + +q_prices = np.zeros((HORIZON, 10)) +for qi in range(10): + p = train_close[-1] + for i in range(HORIZON): + p = p * np.exp(adjusted_q_logret[i, qi]) + q_prices[i, qi] = p + +# 9. 误差指标 +actual = actual_close +pred = pred_prices +mae = np.mean(np.abs(actual - pred)) +rmse = np.sqrt(np.mean((actual - pred) ** 2)) +mape = np.mean(np.abs((actual - pred) / actual)) * 100 +actual_dir = np.diff(actual) +pred_dir = np.diff(pred) +dir_acc = np.mean(actual_dir * pred_dir > 0) * 100 +in_80 = np.mean((actual >= q_prices[:, 1]) & (actual <= q_prices[:, 9])) * 100 +in_40 = np.mean((actual >= q_prices[:, 3]) & (actual <= q_prices[:, 7])) * 100 + +print() +print("=== 纯技巧优化版回测结果 ===", flush=True) +print() +print(" 日期 实际价 预测价 误差 误差%", flush=True) +for i in range(HORIZON): + err = pred[i] - actual[i] + err_pct = err / actual[i] * 100 + print(f" {actual_dates[i].date()} {actual[i]:7.2f} {pred[i]:7.2f} {err:+7.2f} {err_pct:+6.2f}%", flush=True) + +print() +print("=== 误差指标（四版对比）===", flush=True) +print(f" MAE: ${mae:.2f} (原始: $8.46, 对数收益率: $8.00, XReg: $9.21)", flush=True) +print(f" RMSE: ${rmse:.2f} (原始: $11.29, 对数收益率: $11.09, XReg: $12.47)", flush=True) +print(f" MAPE: {mape:.2f}% (原始: 2.98%, 对数收益率: 2.82%, XReg: 3.25%)", flush=True) +print(f" 方向准确率: {dir_acc:.1f}% (原始: 55.6%, 对数收益率: 44.4%, XReg: 44.4%)", flush=True) +print(f" 80%CI覆盖率: {in_80:.1f}% (原始: 90.0%, 对数收益率: 100.0%, XReg: 70.0%)", flush=True) +print(f" 40%CI覆盖率: {in_40:.1f}% (原始: 60.0%, 对数收益率: 80.0%, XReg: 40.0%)", flush=True) + +# 10. 画图 +fig, axes = plt.subplots(2, 1, figsize=(14, 10)) +show_n = min(60, len(train_close)) +ax = axes[0] +ax.plot(range(show_n), train_close[-show_n:], label="Historical", color="steelblue", linewidth=1.5) +x_actual = range(show_n, show_n + HORIZON) +ax.plot(x_actual, actual, label="Actual", color="forestgreen", linewidth=2, marker="o", markersize=4) +ax.plot(x_actual, pred, label="Trick Forecast", color="purple", linewidth=2, linestyle="--") +ax.fill_between(x_actual, q_prices[:, 1], q_prices[:, 9], alpha=0.15, color="purple", label="80% CI") +ax.fill_between(x_actual, q_prices[:, 3], q_prices[:, 7], alpha=0.3, color="purple", label="40% CI") +ax.axvline(x=show_n - 1, color="gray", linestyle="--", alpha=0.5) +ax.set_title(f"Trick: LogRet + Ensemble + SPY + VolAdj (MAPE={mape:.2f}%, Dir={dir_acc:.0f}%)", fontsize=13) +ax.set_ylabel("Price (USD)") +ax.legend(loc="upper left") + +ax2 = axes[1] +errors = pred - actual +ax2.bar(range(HORIZON), errors, color=["purple" if e > 0 else "steelblue" for e in errors], alpha=0.7) +ax2.axhline(y=0, color="black", linewidth=0.8) +ax2.set_title("Trick Forecast Error per Day", fontsize=13) +ax2.set_xlabel("Trading Days After Split") +ax2.set_ylabel("Error (USD)") +ax2.set_xticks(range(HORIZON)) +ax2.set_xticklabels([str(d.date()) for d in actual_dates], rotation=45, fontsize=8) + +plt.tight_layout() +plt.savefig("aapl_backtest_trick.png", dpi=150) +print() +print("Chart saved: aapl_backtest_trick.png", flush=True) +print("Done!", flush=True) diff --git a/stock_backtest_xreg.py b/stock_backtest_xreg.py new file mode 100644 index 0000000..5a2e4a3 --- /dev/null +++ b/stock_backtest_xreg.py @@ -0,0 +1,203 @@ +"""TimesFM + XReg 协变量预测：加入成交量、RSI、SPY 大盘指数""" +import yfinance as yf +import numpy as np +import torch +import timesfm +import matplotlib +matplotlib.use("Agg") +import matplotlib.pyplot as plt +from datetime import datetime, timedelta + + +def calc_rsi(prices, period=14): + deltas = np.diff(prices) + gains = np.where(deltas > 0, deltas, 0.0) + losses = np.where(deltas < 0, -deltas, 0.0) + avg_gain = np.convolve(gains, np.ones(period) / period, mode="valid") + avg_loss = np.convolve(losses, np.ones(period) / period, mode="valid") + avg_loss = np.where(avg_loss == 0, 1e-10, avg_loss) + rs = avg_gain / avg_loss + rsi = 100.0 - (100.0 / (1.0 + rs)) + return np.concatenate([np.full(period, 50.0), rsi]) + + +# 1. 获取数据 +print("下载 AAPL + SPY 数据...", flush=True) +end = datetime.now() +start = end - timedelta(days=365) +df_aapl = yf.download("AAPL", start=start.strftime("%Y-%m-%d"), end=end.strftime("%Y-%m-%d"), progress=False) +df_spy = yf.download("SPY", start=start.strftime("%Y-%m-%d"), end=end.strftime("%Y-%m-%d"), progress=False) + +close = df_aapl["Close"].values.flatten().astype(np.float64) +volume = df_aapl["Volume"].values.flatten().astype(np.float64) +spy_close = df_spy["Close"].values.flatten().astype(np.float64) +dates = df_aapl.index + +# 对齐长度 +min_len = min(len(close), len(spy_close)) +close = close[-min_len:] +volume = volume[-min_len:] +spy_close = spy_close[-min_len:] +dates = dates[-min_len:] + +# 计算 RSI +rsi = calc_rsi(close) + +# 分割 +split_date = "2026-06-15" +split_idx = None +for i, d in enumerate(dates): + if str(d.date()) <= split_date: + split_idx = i + +train_close = close[: split_idx + 1] +train_vol = volume[: split_idx + 1] +train_spy = spy_close[: split_idx + 1] +train_rsi = rsi[: split_idx + 1] +train_dates = dates[: split_idx + 1] + +actual_close = close[split_idx + 1 :] +actual_dates = dates[split_idx + 1 :] +HORIZON = len(actual_close) + +print(f"训练数据: {len(train_close)} 天", flush=True) +print(f"预测目标: {HORIZON} 天", flush=True) + +# 2. 准备协变量 +# 动态数值协变量需要覆盖 context + horizon 的完整长度 +# 对每个序列: train 部分用实际值，test 部分需要"未来值" +# 对于 RSI 和 Volume，我们没有未来值，用最后一个值填充 +# 对于 SPY，用训练集最后一个值填充（因为我们无法预知未来 SPY） + +# 协变量需要: 每个协变量是一个 list，每个元素对应一个输入序列的完整长度(context+horizon) +full_len = len(train_close) + HORIZON + +# 成交量: train 用实际值，future 用最近 5 日均值 +vol_future = np.mean(train_vol[-5:]) +vol_full = np.concatenate([train_vol, np.full(HORIZON, vol_future)]) + +# RSI: train 用实际值，future 用 50（中性） +rsi_full = np.concatenate([train_rsi, np.full(HORIZON, 50.0)]) + +# SPY: train 用实际值，future 用最后一个值 +spy_future = train_spy[-1] +spy_full = np.concatenate([train_spy, np.full(HORIZON, spy_future)]) + +# 3. 加载模型 +print() +print("加载 TimesFM 模型...", flush=True) +torch.set_float32_matmul_precision("high") +model = timesfm.TimesFM_2p5_200M_torch.from_pretrained( + "google/timesfm-2.5-200m-pytorch", torch_compile=False +) +print("模型加载完成", flush=True) + +# 4. 编译（return_backcast=True 是 XReg 必需的） +model.compile( + timesfm.ForecastConfig( + max_context=512, + max_horizon=128, + normalize_inputs=True, + use_continuous_quantile_head=True, + force_flip_invariance=False, + infer_is_positive=True, + fix_quantile_crossing=True, + return_backcast=True, # XReg 需要 + ) +) +print("编译完成", flush=True) + +# 5. 用 XReg 预测 +print("运行 XReg 协变量预测...", flush=True) + +# 动态数值协变量: dict[str, list[list[float]]] +# 每个协变量是一个 list，其中每个元素是一个序列（对应一个输入） +# 这里只有一个输入序列 +dynamic_num_covs = { + "volume": [vol_full], + "rsi": [rsi_full], + "spy_close": [spy_full], +} + +# 静态数值协变量 +static_num_covs = { + "avg_volume": [np.mean(train_vol)], +} + +point_outputs, quantile_outputs = model.forecast_with_covariates( + inputs=[train_close], + dynamic_numerical_covariates=dynamic_num_covs, + static_numerical_covariates=static_num_covs, + xreg_mode="xreg + timesfm", # 先回归再预测残差 + normalize_xreg_target_per_input=True, + ridge=1.0, +) +print("XReg 预测完成！", flush=True) + +# 6. 提取结果 +pred = np.array(point_outputs[0][:HORIZON]) +q = np.array(quantile_outputs[0]) # (horizon, 10) or (full, 10) +# quantile_outputs 可能包含 backcast，取最后 HORIZON 个 +if q.shape[0] > HORIZON: + q = q[-HORIZON:] + +# 7. 计算误差 +actual = actual_close +mae = np.mean(np.abs(actual - pred)) +rmse = np.sqrt(np.mean((actual - pred) ** 2)) +mape = np.mean(np.abs((actual - pred) / actual)) * 100 +actual_dir = np.diff(actual) +pred_dir = np.diff(pred) +dir_acc = np.mean(actual_dir * pred_dir > 0) * 100 +in_80 = np.mean((actual >= q[:, 1]) & (actual <= q[:, 9])) * 100 +in_40 = np.mean((actual >= q[:, 3]) & (actual <= q[:, 7])) * 100 + +print() +print("=== XReg 协变量版回测结果 ===", flush=True) +print() +print(" 日期 实际价 预测价 误差 误差%", flush=True) +for i in range(HORIZON): + err = pred[i] - actual[i] + err_pct = err / actual[i] * 100 + print(f" {actual_dates[i].date()} {actual[i]:7.2f} {pred[i]:7.2f} {err:+7.2f} {err_pct:+6.2f}%", flush=True) + +print() +print("=== 误差指标（三版对比）===", flush=True) +print(f" MAE: ${mae:.2f} (原始: $8.46, 优化: $8.00)", flush=True) +print(f" RMSE: ${rmse:.2f} (原始: $11.29, 优化: $11.09)", flush=True) +print(f" MAPE: {mape:.2f}% (原始: 2.98%, 优化: 2.82%)", flush=True) +print(f" 方向准确率: {dir_acc:.1f}% (原始: 55.6%, 优化: 44.4%)", flush=True) +print(f" 80%CI覆盖率: {in_80:.1f}% (原始: 90.0%, 优化: 100.0%)", flush=True) +print(f" 40%CI覆盖率: {in_40:.1f}% (原始: 60.0%, 优化: 80.0%)", flush=True) + +# 8. 画图 +fig, axes = plt.subplots(2, 1, figsize=(14, 10)) +show_n = min(60, len(train_close)) +ax = axes[0] +ax.plot(range(show_n), train_close[-show_n:], label="Historical", color="steelblue", linewidth=1.5) +x_actual = range(show_n, show_n + HORIZON) +ax.plot(x_actual, actual, label="Actual", color="forestgreen", linewidth=2, marker="o", markersize=4) +ax.plot(x_actual, pred, label="XReg Forecast", color="darkorange", linewidth=2, linestyle="--") +ax.fill_between(x_actual, q[:, 1], q[:, 9], alpha=0.15, color="darkorange", label="80% CI") +ax.fill_between(x_actual, q[:, 3], q[:, 7], alpha=0.3, color="darkorange", label="40% CI") +ax.axvline(x=show_n - 1, color="gray", linestyle="--", alpha=0.5) +ax.set_title(f"XReg: AAPL + Volume/RSI/SPY (MAPE={mape:.2f}%, Dir={dir_acc:.0f}%)", fontsize=13) +ax.set_ylabel("Price (USD)") +ax.legend(loc="upper left") + +# 误差柱状图 +ax2 = axes[1] +errors = pred - actual +ax2.bar(range(HORIZON), errors, color=["darkorange" if e > 0 else "steelblue" for e in errors], alpha=0.7) +ax2.axhline(y=0, color="black", linewidth=0.8) +ax2.set_title("XReg Forecast Error per Day", fontsize=13) +ax2.set_xlabel("Trading Days After Split") +ax2.set_ylabel("Error (USD)") +ax2.set_xticks(range(HORIZON)) +ax2.set_xticklabels([str(d.date()) for d in actual_dates], rotation=45, fontsize=8) + +plt.tight_layout() +plt.savefig("aapl_backtest_xreg.png", dpi=150) +print() +print("Chart saved: aapl_backtest_xreg.png", flush=True) +print("Done!", flush=True) diff --git a/stock_forecast.py b/stock_forecast.py new file mode 100644 index 0000000..2ec4db3 --- /dev/null +++ b/stock_forecast.py @@ -0,0 +1,84 @@ +import yfinance as yf +import numpy as np +import torch +import timesfm +import matplotlib +matplotlib.use("Agg") +import matplotlib.pyplot as plt +from datetime import datetime, timedelta + +# 1. 获取 AAPL 近 1 年收盘价 +print("下载 AAPL 股票数据...", flush=True) +end = datetime.now() +start = end - timedelta(days=365) +df = yf.download("AAPL", start=start.strftime("%Y-%m-%d"), end=end.strftime("%Y-%m-%d"), progress=False) +close = df["Close"].values.flatten().astype(np.float32) +dates = df.index +print(f"获取到 {len(close)} 个交易日", flush=True) +print(f"价格范围: {close.min():.2f} ~ {close.max():.2f}", flush=True) +print(f"最近 5 日收盘价: {close[-5:]}", flush=True) + +# 2. 加载 TimesFM 模型 +print() +print("加载 TimesFM 模型...", flush=True) +torch.set_float32_matmul_precision("high") +model = timesfm.TimesFM_2p5_200M_torch.from_pretrained( + "google/timesfm-2.5-200m-pytorch", torch_compile=False +) +print("模型加载完成", flush=True) + +# 3. 编译并预测未来 20 个交易日 +HORIZON = 20 +model.compile( + timesfm.ForecastConfig( + max_context=512, + max_horizon=128, + normalize_inputs=True, + use_continuous_quantile_head=True, + force_flip_invariance=False, + infer_is_positive=True, + fix_quantile_crossing=True, + ) +) +print("编译完成，开始预测...", flush=True) + +pf, qf = model.forecast(horizon=HORIZON, inputs=[close]) +print("预测完成！", flush=True) + +# 4. 输出结果 +print() +print("=== AAPL 未来 20 个交易日预测 ===", flush=True) +print(f"当前价格: {close[-1]:.2f}", flush=True) +print() +print(" 日期(估计) 点预测 q10(低) q90(高)", flush=True) +last_date = dates[-1] +for i in range(HORIZON): + est_date = last_date + timedelta(days=i + 1) + print( + f" {est_date.strftime('%Y-%m-%d')} {pf[0][i]:7.2f} {qf[0][i,1]:7.2f} {qf[0][i,9]:7.2f}", + flush=True, + ) + +print() +print(f"预测均价: {pf[0].mean():.2f}", flush=True) +print(f"预测涨跌: {(pf[0][-1] - close[-1]) / close[-1] * 100:+.2f}%", flush=True) +print(f"80%置信区间: {qf[0,-1,1]:.2f} ~ {qf[0,-1,9]:.2f}", flush=True) + +# 5. 画图 +fig, ax = plt.subplots(figsize=(14, 6)) +show_n = min(60, len(close)) +ax.plot(range(show_n), close[-show_n:], label="历史收盘价", color="steelblue", linewidth=1.5) +x_fc = range(show_n, show_n + HORIZON) +ax.plot(x_fc, pf[0], label="点预测(中位数)", color="tomato", linewidth=2) +ax.fill_between(x_fc, qf[0, :, 1], qf[0, :, 9], alpha=0.2, color="tomato", label="80% 置信区间") +ax.fill_between(x_fc, qf[0, :, 3], qf[0, :, 7], alpha=0.35, color="tomato", label="40% 置信区间") +ax.axvline(x=show_n - 1, color="gray", linestyle="--", alpha=0.5, label="预测起点") +ax.set_title("AAPL 收盘价预测 (TimesFM 2.5)", fontsize=14) +ax.set_xlabel("交易日") +ax.set_ylabel("价格 (USD)") +ax.legend(loc="upper left") +plt.tight_layout() +plt.savefig("aapl_forecast.png", dpi=150) +print() +print("图表已保存: aapl_forecast.png", flush=True) +print("✅ 完成！", flush=True)