""" 經典學術異常在台股：10 個因子的三項驗證（2014-06 起至最新資料） 三項驗證（「贏」＝日夏普高於同期 0050 含息買進持有）： 1. 參數高原：每個因子 2 個主要參數各取 3 個值，9 組裡至少 7 組要贏。 2. 前後兩段：回測期間對半切，兩段各自要贏。 3. 交易成本：sim() 預設手續費 0.1425%（買賣各一次）與證交稅 0.3%（賣出）。 共同設定：上市櫃普通股（含已下市），排除金融股與存託憑證；20 日平均成交金額至少 2,000 萬元； 每月最後一個交易日依因子排名換股，等權持有。財報依各公司實際公布日才可用。 執行：安裝 finlab 後直接跑，資料需要時會自動引導登入。 輸出：metrics.json、data.csv（各因子代表參數的月底淨值與 0050）、report_strategy.html。 """ import datetime import itertools import json import numpy as np import pandas as pd import finlab from finlab import data from finlab.backtest import sim START = "2014-06-01" # 財報從 2013 年第 1 季起，這天起 10 個因子都算得出來 LONG_START = "2008-01-01" # 只用股價的因子，另從這天起跑一次較長的期間 MIN_TURNOVER_VALUE = 2e7 # 20 日平均成交金額至少 2,000 萬元 TRADING_DAYS = 252 # --------------------------------------------------------------------------- # 1. 股票池、價格與 0050 含息基準 # --------------------------------------------------------------------------- def plain(key): df = pd.DataFrame(data.get(key)) df.index = pd.to_datetime(df.index) return df adj_close = plain("etl:adj_close") bench = adj_close["0050"].dropna() categories = pd.DataFrame(data.get("security_categories")).set_index("stock_id")["category"] EXCLUDED_CATEGORIES = {"金融保險業", "金融業", "存託憑證"} def is_common_stock(sid): """4 位數普通股；已下市的公司沒有產業別，用代號排除金融（28xx）與存託憑證（91xx）。""" return ( len(sid) == 4 and sid.isdigit() and not sid.startswith(("0", "28", "91")) and categories.get(sid) not in EXCLUDED_CATEGORIES ) stocks = [s for s in adj_close.columns if is_common_stock(s)] adj_close = adj_close[stocks] days = adj_close.loc[LONG_START:].index daily_return = adj_close.pct_change(fill_method=None).reindex(days) market_return = plain("benchmark_return:發行量加權股價報酬指數").iloc[:, 0].pct_change().reindex(days) avg_turnover_value = plain("price:成交金額").reindex(columns=stocks).rolling(20).mean() liquid = avg_turnover_value.reindex(days) > MIN_TURNOVER_VALUE # 每個月最後一個交易日（換股日） month_id = days.to_period("M") month_ends = days[np.append(month_id[1:] != month_id[:-1], True)] def on_days(df): """把低頻資料展開到每個交易日：每天只看得到當天以前已公布的最新一期。""" df = df.reindex(columns=stocks) return df.reindex(df.index.union(days)).ffill().reindex(days) def quarterly_on_days(df): """季資料換成各公司實際公布日，再展開到每個交易日。""" dated = pd.DataFrame(df.index_str_to_date()) dated.index = pd.to_datetime(dated.index) return on_days(dated) # --------------------------------------------------------------------------- # 2. 十個因子：參數 → 每日分數 # --------------------------------------------------------------------------- # 季財報為單季值，index 為 2013-Q1 這類季別。 gross_profit = data.get("financial_statement:營業毛利") total_assets = data.get("financial_statement:資產總額") net_income = data.get("financial_statement:歸屬母公司淨利_損") + data.get("financial_statement:歸屬非控制權益淨利_損").fillna(0) operating_cash_flow = data.get("financial_statement:營業活動之淨現金流入_流出") def gross_profitability(quarters): """Novy-Marx (2013)：近 quarters 季毛利（年化）÷ 最新總資產。""" annual_gp = gross_profit.rolling(quarters).sum() * 4 / quarters return quarterly_on_days(annual_gp / total_assets) def asset_growth(quarters): """Cooper, Gulen & Schill (2008)：總資產比 quarters 季前的成長率。""" return quarterly_on_days(total_assets / total_assets.shift(quarters) - 1) def accruals(quarters): """Sloan (1996) 的現金流量版：近 quarters 季（淨利 − 營業現金流），年化後 ÷ 平均總資產。""" accrual = (net_income - operating_cash_flow).rolling(quarters).sum() * 4 / quarters average_assets = (total_assets + total_assets.shift(quarters)) / 2 return quarterly_on_days(accrual / average_assets) # 淨股票發行：發行股數（月資料，index 為可取得的月底）扣掉無償配股 issued = plain("internal_equity_changes:發行股數").reindex(columns=stocks) dividends = pd.DataFrame(data.get("dividend_announcement")) dividends = dividends[dividends["除權交易日"].notna() & dividends["stock_id"].astype(str).isin(stocks)].copy() par = np.where(dividends["普通股每股面額"].astype(str).str.contains("5.0000"), 5.0, 10.0) stock_dividend = (dividends["盈餘轉增資配股(元/股)"].fillna(0) + dividends["法定盈餘公積、資本公積轉增資配股(元/股)"].fillna(0)) / par dividends = dividends.assign(ratio=1 + stock_dividend, sid=dividends["stock_id"].astype(str)) dividends = dividends[dividends["ratio"] > 1] # 某月月底的發行股數，反映上個月以前基準日的配股；基準日落在最近兩個月的，那個月不算（股數更新時點不確定） row_month = issued.index.to_period("M") record_month = dividends["權利分派基準日"].dt.to_period("M") split_factor = pd.DataFrame(1.0, index=issued.index, columns=stocks) near_record_date = pd.DataFrame(False, index=issued.index, columns=stocks) for sid, month, ratio in zip(dividends["sid"], record_month, dividends["ratio"]): split_factor.loc[row_month > month, sid] *= ratio near_record_date.loc[(row_month >= month) & (row_month <= month + 1), sid] = True # 彌補虧損減資與面額變更（股票分割）只改變每股單位，不是發行或買回；回看期間內有這類事件的公司不列入排名 unit_change = pd.DataFrame(False, index=issued.index, columns=stocks) events = [] for key in ("capital_reduction_tse:減資原因", "capital_reduction_otc:減資原因"): reason = plain(key).reindex(columns=stocks).stack() events += list(reason[reason.astype(str).str.contains("彌補虧損")].index) for key in ("par_value_change_tse:恢復買賣日期", "par_value_change_otc:恢復買賣日期"): events += list(plain(key).reindex(columns=stocks).stack().dropna().index) for day, sid in events: unit_change.loc[row_month == pd.Period(day, "M"), sid] = True adjusted_shares = np.log(issued / split_factor) def net_issuance(months): """Pontiff & Woodgate (2008)：months 個月來、扣掉配股後的發行股數 log 變化。""" change = adjusted_shares - adjusted_shares.shift(months) # 股數在事件後 1～2 個月才更新，所以事件之後多排除 2 個月 excluded = ( near_record_date | near_record_date.shift(months, fill_value=False) | (unit_change.astype(float).rolling(months + 3, min_periods=1).max() > 0) ) # 展開到每日之後才排除；先遮再展開，ffill 會補回上一期的分數 return on_days(change).mask(on_days(excluded.astype(float)) > 0) def high_52_week(window): """George & Hwang (2004)：還原收盤價 ÷ 近 window 日最高價。""" price = adj_close.reindex(days) return price / price.rolling(window, min_periods=window).max() def _market_moments(window): mean_m = market_return.rolling(window).mean() var_m = market_return.rolling(window).var(ddof=0) cov = daily_return.mul(market_return, axis=0).rolling(window).mean() - daily_return.rolling(window).mean().mul(mean_m, axis=0) return cov, var_m def idiosyncratic_volatility(window): """Ang, Hodrick, Xing & Zhang (2006)：對加權報酬指數迴歸後的殘差波動（年化）。""" cov, var_m = _market_moments(window) residual_var = daily_return.rolling(window).var(ddof=0) - cov.pow(2).div(var_m, axis=0) return np.sqrt(residual_var.clip(lower=0) * TRADING_DAYS) def max_daily_return(window): """Bali, Cakici & Whitelaw (2011)：近 window 日最大單日報酬。""" return daily_return.rolling(window).max() def market_beta(window): """Frazzini & Pedersen (2014) 只做多的版本：近 window 日對加權報酬指數的 beta。""" cov, var_m = _market_moments(window) return cov.div(var_m, axis=0) monthly_price = adj_close.reindex(month_ends) monthly_return = monthly_price.pct_change(fill_method=None) def return_seasonality(years): """Heston & Sadka (2008)：過去 years 年「下個月同一個月份」的平均報酬。""" same_month = [monthly_return.shift(12 * j - 1) for j in range(1, years + 1)] return on_days(sum(same_month) / years) market_value = plain("etl:market_value").reindex(columns=stocks) industry = categories.reindex(stocks) def industry_momentum(months, industries): """Moskowitz & Grinblatt (1999)：近 months 個月市值加權報酬最高的 industries 個產業，各買市值最大的 10 檔。""" past_return = monthly_price / monthly_price.shift(months) - 1 cap = market_value.reindex(days).ffill().reindex(month_ends) past_cap = cap.shift(months) eligible = liquid.reindex(month_ends) & past_return.notna() & past_cap.notna() position = pd.DataFrame(False, index=month_ends, columns=stocks) for day in month_ends: members = eligible.loc[day] members = members[members].index if len(members) == 0: continue frame = pd.DataFrame({ "industry": industry.reindex(members), "ret": past_return.loc[day, members], "weight": past_cap.loc[day, members], "cap": cap.loc[day, members], }).dropna() sizes = frame.groupby("industry")["ret"].size() frame = frame[frame["industry"].isin(sizes[sizes >= 5].index)] weighted = (frame["ret"] * frame["weight"]).groupby(frame["industry"]).sum() / frame.groupby("industry")["weight"].sum() winners = weighted.nlargest(industries).index picks = frame[frame["industry"].isin(winners)].groupby("industry")["cap"].nlargest(10).index.get_level_values(1) position.loc[day, picks] = True return position # --------------------------------------------------------------------------- # 3. 排名選股：分數 → 每月持股 # --------------------------------------------------------------------------- avg_turnover_value = avg_turnover_value.reindex(month_ends) def top_n(score, n, ascending): """分數排名前 n 檔。同分（例如都在 52 週高點、股數都沒變）時，成交金額大的優先。""" score = score.reindex(index=month_ends, columns=stocks).where(liquid.reindex(month_ends)) rank = score.rank(axis=1, ascending=ascending, method="average") tie_break = avg_turnover_value.where(score.notna()).rank(axis=1, ascending=False, pct=True) * 0.5 return (rank + tie_break).rank(axis=1, method="first") <= n def ranked(factor, ascending): def build(param, n): return top_n(factor(param), n, ascending) return build HOLDINGS = [20, 30, 50] FACTORS = { # 名稱: (建立持股, 參數網格, 原始參數) "gross_profitability": (ranked(gross_profitability, False), {"quarters": [2, 4, 8], "n": HOLDINGS}, {"quarters": 4, "n": 30}), "asset_growth": (ranked(asset_growth, True), {"quarters": [2, 4, 8], "n": HOLDINGS}, {"quarters": 4, "n": 30}), "accruals": (ranked(accruals, True), {"quarters": [2, 4, 8], "n": HOLDINGS}, {"quarters": 4, "n": 30}), "net_issuance": (ranked(net_issuance, True), {"months": [6, 12, 24], "n": HOLDINGS}, {"months": 12, "n": 30}), "high_52_week": (ranked(high_52_week, False), {"window": [126, 252, 504], "n": HOLDINGS}, {"window": 252, "n": 30}), "idiosyncratic_vol": (ranked(idiosyncratic_volatility, True), {"window": [10, 21, 42], "n": HOLDINGS}, {"window": 21, "n": 30}), "max_effect": (ranked(max_daily_return, True), {"window": [10, 21, 42], "n": HOLDINGS}, {"window": 21, "n": 30}), "low_beta": (ranked(market_beta, True), {"window": [126, 252, 504], "n": HOLDINGS}, {"window": 252, "n": 30}), "seasonality": (ranked(return_seasonality, False), {"years": [3, 5, 7], "n": HOLDINGS}, {"years": 5, "n": 30}), "industry_momentum": (industry_momentum, {"months": [3, 6, 12], "industries": [2, 3, 5]}, {"months": 6, "industries": 3}), } PRICE_ONLY = ["high_52_week", "idiosyncratic_vol", "max_effect", "low_beta"] # --------------------------------------------------------------------------- # 4. 績效與三項驗證 # --------------------------------------------------------------------------- def perf(equity): equity = equity.dropna() ret = equity.pct_change().dropna() years = (equity.index[-1] - equity.index[0]).days / 365.25 return { "cagr": float((equity.iloc[-1] / equity.iloc[0]) ** (1 / years) - 1), "sharpe": float(ret.mean() / ret.std() * np.sqrt(TRADING_DAYS)), "mdd": float((equity / equity.cummax() - 1).min()), "start": str(equity.index[0].date()), "end": str(equity.index[-1].date()), } def versus_0050(equity): """全期與前後兩段，策略與 0050 在完全相同的日期上比較。""" b = bench.reindex(equity.index).ffill() mid = equity.index[0] + (equity.index[-1] - equity.index[0]) / 2 out = {} for name, part in [("full", slice(None)), ("first_half", slice(None, mid)), ("second_half", slice(mid, None))]: s, k = perf(equity.loc[part]), perf(b.loc[part]) out[name] = {"strategy": s, "bench": k, "win": s["sharpe"] > k["sharpe"]} return out def monthly_turnover(position): """每月換股時，平均有多少比例的部位被換掉（單邊）。""" weights = position.div(position.sum(axis=1).replace(0, np.nan), axis=0).fillna(0) return float(weights.diff().abs().sum(axis=1).iloc[1:].mean() / 2) def backtest(position, start, **cost): position = position.loc[start:] report = sim(position, resample="M", upload=False, **cost) first_day = position.index[position.any(axis=1)][0] return report, report.creturn.loc[first_day:] results, curves, reports, positions = {}, {}, {}, {} for name, (build, grid, center) in FACTORS.items(): rows = [] for values in itertools.product(*grid.values()): params = dict(zip(grid, values)) position = build(*values) report, equity = backtest(position, START) held = position.loc[START:] row = { "params": params, "avg_holdings": float(held[held.any(axis=1)].sum(axis=1).mean()), "turnover": monthly_turnover(held), **versus_0050(equity), } if params == center: curves[name] = equity reports[name] = report positions[name] = held _, zero_cost = backtest(position, START, fee_ratio=0, tax_ratio=0) row["zero_cost"] = versus_0050(zero_cost)["full"] if name in PRICE_ONLY: _, long_equity = backtest(position, LONG_START) row["since_2008"] = versus_0050(long_equity) rows.append(row) print(name, params, f"sharpe {row['full']['strategy']['sharpe']:.2f} vs 0050 {row['full']['bench']['sharpe']:.2f}", flush=True) base = next(r for r in rows if r["params"] == center) wins = sum(r["full"]["win"] for r in rows) checks = { "neighbors": wins / len(rows) >= 0.7, "halves": base["first_half"]["win"] and base["second_half"]["win"], "cost": base["full"]["win"], } checks["all"] = all(checks.values()) results[name] = {"grid": grid, "center": center, "wins": wins, "checks": checks, "combos": rows} print(name, "wins", wins, "/", len(rows), checks, flush=True) # --------------------------------------------------------------------------- # 5. 因子本身：五分位報酬與 Rank IC（原始參數，流動性股票池） # --------------------------------------------------------------------------- SCORES = { "gross_profitability": (gross_profitability(4), False), "asset_growth": (asset_growth(4), True), "accruals": (accruals(4), True), "net_issuance": (net_issuance(12), True), "high_52_week": (high_52_week(252), False), "idiosyncratic_vol": (idiosyncratic_volatility(21), True), "max_effect": (max_daily_return(21), True), "low_beta": (market_beta(252), True), "seasonality": (return_seasonality(5), False), } next_month_return = monthly_return.shift(-1) def annualize(monthly_returns): r = pd.Series(monthly_returns).fillna(0) return float((1 + r).prod() ** (12 / len(r)) - 1) def quintiles(score, ascending): """每月把流動性股票依分數分 5 組（第 5 組＝因子看好的一端），算下個月等權報酬。""" score = score.reindex(index=month_ends, columns=stocks).where(liquid.reindex(month_ends)) if ascending: score = -score groups, medians, ics = [], [], [] for day in month_ends[(month_ends >= START)][:-1]: x = score.loc[day].dropna() y = next_month_return.loc[day].reindex(x.index).dropna() x = x.reindex(y.index) if len(x) < 50: continue q = pd.qcut(x.rank(method="first"), 5, labels=False) groups.append(y.groupby(q).mean()) medians.append(y.groupby(q).median()) ics.append(x.corr(y, method="spearman")) table, median_table = pd.DataFrame(groups), pd.DataFrame(medians) spread = table[4] - table[0] ics = pd.Series(ics) return { "annual_by_quintile": [annualize(table[q]) for q in range(5)], "annual_median_stock_by_quintile": [annualize(median_table[q]) for q in range(5)], "spread_monthly_mean": float(spread.mean()), "spread_t": float(spread.mean() / spread.std() * np.sqrt(len(spread))), "rank_ic_mean": float(ics.mean()), "rank_ic_t": float(ics.mean() / ics.std() * np.sqrt(len(ics))), "months": len(table), } signals = {name: quintiles(score, asc) for name, (score, asc) in SCORES.items()} # 對照組：整個流動性股票池等權持有、每月再平衡 _, universe_equity = backtest(liquid.reindex(month_ends), START) signals["universe_equal_weight"] = versus_0050(universe_equity) curves["universe_equal_weight"] = universe_equity # 各因子原始參數平均持股的產業分布（前 3 大） industry_mix = {} for name, held in positions.items(): share = held.astype(float).T.groupby(industry.reindex(held.columns).fillna("已下市／無分類")).sum().T share = share.div(held.sum(axis=1).replace(0, np.nan), axis=0).mean().nlargest(3) industry_mix[name] = {k: round(float(v), 4) for k, v in share.items()} signals["industry_mix"] = industry_mix # 台股特有的背景數字：漲停板、增資 rebalance_days = month_ends[month_ends >= START] eligible = liquid.reindex(rebalance_days) max_21 = max_daily_return(21).reindex(rebalance_days).where(eligible) issuance_12 = net_issuance(12).reindex(rebalance_days).where(eligible) signals["share_with_limit_up_day_21d"] = float((max_21 >= 0.095).sum(axis=1).div(max_21.notna().sum(axis=1)).mean()) signals["share_issuance_12m"] = { "above_5pct": float((issuance_12 > np.log(1.05)).sum(axis=1).div(issuance_12.notna().sum(axis=1)).mean()), "negative": float((issuance_12 < 0).sum(axis=1).div(issuance_12.notna().sum(axis=1)).mean()), } # --------------------------------------------------------------------------- # 6. 輸出 # --------------------------------------------------------------------------- best = max(reports, key=lambda name: perf(curves[name])["sharpe"]) reports[best].to_html("report_strategy.html") meta = { "run_date": str(datetime.date.today()), "finlab_version": finlab.__version__, "data_end": str(days[-1].date()), "start": START, "long_start": LONG_START, "universe_size_latest": int(liquid.iloc[-1].sum()), "liquidity_filter_twd": MIN_TURNOVER_VALUE, "report_strategy": best, "benchmark": "0050 含息買進持有（etl:adj_close，不含交易成本）", "sharpe": "日報酬平均 / 標準差 × √252，未扣無風險利率", } json.dump({"meta": meta, "factors": results, "signals": signals}, open("metrics.json", "w"), ensure_ascii=False, indent=1, default=str) table = pd.DataFrame(curves) table = table / table.bfill().iloc[0] table["0050"] = bench.reindex(table.index).ffill() table["0050"] = table["0050"] / table["0050"].dropna().iloc[0] table.reindex(month_ends.intersection(table.index)).round(4).to_csv("data.csv", index_label="date") print("done", meta)