# -*- coding: utf-8 -*-
"""雙因子橫斷面選股完整檢驗:月營收年增率 × 自由現金流殖利率。

對應文章:
  https://finlab.finance/blog/python-investing-benefits-tutorial

執行方式(安裝 finlab 後直接跑,首次執行會自動引導登入):
  pip install finlab
  python strategy.py

研究問題:
  「營收成長 × 自由現金流為正」這個聽起來合理的雙因子組合,
  在 2018-01 ~ 2026-06 的台股能不能贏過 0050 買進持有?

四個設計變體(A→D 是文章揭露的迭代過程,不是暗中調參):
  A. 兩因子各做橫斷面 z-score(1%/99% 縮尾),等權合成,取前 20 名
  B. 改用橫斷面百分位排名合成(對極端值穩健),取前 20 名
  C. B 再加上兩因子皆為正的門檻,取前 20 名
  D. C 的持股檔數放寬到 40 檔(分散個股風險)

回測結果(資料截止 2026-06-09,finlab sim() 台股預設交易成本):
  A:CAGR 7.08%、日夏普 0.42、MDD -47.17%
  B:CAGR 11.11%、日夏普 0.58、MDD -43.05%
  C:CAGR 9.88%、日夏普 0.53、MDD -42.52%
  D:CAGR 17.21%、日夏普 0.83、MDD -36.69%
  0050 含息買進持有:CAGR 25.05%、日夏普 1.22、MDD -33.96%
  → 四個變體在本窗口全部輸給 0050,結論詳見文章。

投資警語:本程式僅供量化研究與教學用途,過去績效不代表未來表現,
不構成任何投資建議;實際交易前請自行評估風險、滑價與交易容量。
"""
import pandas as pd
from finlab import data
from finlab.backtest import sim

START = "2018-01-01"
END = "2026-06-09"  # 文章發佈時的資料截止日;重跑最新資料可移除此限制


def cap(df):
    """把資料截到 END,確保結果可重現。"""
    return df[df.index <= END]


# ---------- 載入資料 ----------
close = cap(data.get("price:收盤價"))
vol = cap(data.get("price:成交股數"))
market_value = cap(data.get("etl:market_value"))

# 月營收年增率:finlab 的索引已是「公布截止日」,直接用就不會偷看未來
rev_yoy = cap(data.get("monthly_revenue:去年同月增減(%)"))

# 自由現金流(季資料):index_str_to_date() 把季度索引換成財報公布截止日
fcf_quarterly = cap(data.get("fundamental_features:自由現金流量").index_str_to_date())

# ---------- 股票池:60 日均成交金額在全市場前 50% ----------
amount = (close * vol).rolling(60).mean()
pool = (amount.rank(axis=1, pct=True) > 0.5) & close.notna()

# ---------- 因子計算 ----------
# 因子一:月營收年增率,向前填值對齊到每個交易日
rev_daily = rev_yoy.reindex(close.index, method="ffill")

# 因子二:自由現金流殖利率 = 近四季 FCF 合計 / 市值
fcf_ttm = fcf_quarterly.rolling(4).sum()
fcf_yield = fcf_ttm.reindex(close.index, method="ffill") / market_value

# 兩個因子都要有值才納入評分
valid = pool & rev_daily.notna() & fcf_yield.notna()


# ---------- 打分方法 ----------
def zscore(factor, mask):
    """橫斷面 z 分數:先以 1%/99% 分位縮尾,再標準化。"""
    d = factor.where(mask)
    lower = d.quantile(0.01, axis=1)
    upper = d.quantile(0.99, axis=1)
    d = d.clip(lower, upper, axis=0)
    return d.sub(d.mean(axis=1), axis=0).div(d.std(axis=1), axis=0)


def rank_pct(factor, mask):
    """橫斷面百分位排名(0~1),對極端值不敏感。"""
    return factor.where(mask).rank(axis=1, pct=True)


# 變體 A:z-score 合成,前 20 名
score_a = (zscore(rev_daily, valid) + zscore(fcf_yield, valid)) / 2
position_a = score_a.where(valid).is_largest(20)

# 變體 B:百分位排名合成,前 20 名
score_b = (rank_pct(rev_daily, valid) + rank_pct(fcf_yield, valid)) / 2
position_b = score_b.where(valid).is_largest(20)

# 變體 C / D:加上「營收有成長、自由現金流為正」門檻
mask_positive = valid & (rev_daily > 0) & (fcf_yield > 0)
score_c = (rank_pct(rev_daily, mask_positive) + rank_pct(fcf_yield, mask_positive)) / 2
position_c = score_c.where(mask_positive).is_largest(20)
position_d = score_c.where(mask_positive).is_largest(40)

# ---------- 回測(每月再平衡,finlab 預設台股交易成本) ----------
variants = {
    "A z-score 合成前 20": position_a,
    "B 百分位排名前 20": position_b,
    "C 排名+正值門檻前 20": position_c,
    "D 排名+正值門檻前 40": position_d,
}

for name, position in variants.items():
    position = position[position.index >= START]
    report = sim(position, resample="M", upload=False)
    # sim() 的 creturn 會延伸到執行當日;統計前先雙端截斷到 START~END,
    # 並用純算術公式計算,與 0050 基準同一口徑
    creturn = report.creturn
    creturn = creturn[(creturn.index >= START) & (creturn.index <= END)]
    daily_return = creturn.pct_change().dropna()
    total = creturn.iloc[-1] / creturn.iloc[0] - 1
    n_years = (creturn.index[-1] - creturn.index[0]).days / 365.25
    print(
        name,
        "| CAGR:", round(((1 + total) ** (1 / n_years) - 1) * 100, 2), "%",
        "| 日夏普:", round(daily_return.mean() / daily_return.std() * 252 ** 0.5, 2),
        "| MDD:", round((creturn / creturn.cummax() - 1).min() * 100, 2), "%",
    )

# ---------- 0050 含息基準(還原價純指數算術,不含交易成本) ----------
adj_close = cap(data.get("etl:adj_close"))
benchmark = adj_close["0050"]
benchmark = benchmark[benchmark.index >= START].dropna()
total_return = benchmark.iloc[-1] / benchmark.iloc[0] - 1
years = (benchmark.index[-1] - benchmark.index[0]).days / 365.25
cagr = (1 + total_return) ** (1 / years) - 1
print("0050 含息買進持有 | CAGR:", round(cagr * 100, 2), "%")