Trong chương này, chúng ta đối mặt với vấn đề kiểm định đa giả thuyết trong định giá tài sản thực nghiệm. Chúng ta tái lập một tập hợp rộng các hiện tượng bất thường (anomaly) đã được công bố trên thị trường Việt Nam, áp dụng các hiệu chỉnh chính thức theo số lượng kiểm định đã thực hiện, ước lượng tỷ lệ các phát hiện sai, và xây dựng một khung phân tích để quyết định nhân tố nào xứng đáng được tin cậy.
Đến năm 2016, tài chính học thuật đã ghi nhận hơn 300 biến được cho là dự báo được mặt cắt ngang của lợi nhuận cổ phiếu. Cochrane (2011) đã gọi sự sinh sôi này bằng một cụm từ đáng nhớ là “vườn thú nhân tố” (zoo of factors). Harvey, Liu, và Zhu (2016) thống kê hơn 300 nhân tố đã công bố và lập luận rằng ngưỡng ý nghĩa thống kê thông thường \(t > 2.0\) là quá thấp so với số lượng giả thuyết đã được kiểm định, dù là công khai hay ngầm định, qua nhiều thập kỷ nghiên cứu. Ngưỡng \(t > 3.0\) (và trong một số trường hợp \(t > 3.78\)) mà họ đề xuất sẽ loại bỏ phần lớn các hiện tượng bất thường đã công bố.
Vấn đề này không chỉ mang tính học thuật. Mỗi nhà nghiên cứu sắp xếp cổ phiếu theo một biến mới và thấy \(t > 2\) đều đang ngầm kiểm định một giả thuyết đã được kiểm định hàng trăm lần với các biến khác. Xác suất để ít nhất một trong 300 kiểm định độc lập cho ra \(t > 2\) do ngẫu nhiên (ngay cả khi không tồn tại hiệu ứng thực) về cơ bản bằng 1,0. Đây chính là vấn đề kiểm định đa giả thuyết, và việc không tính đến nó có nghĩa là một phần đáng kể của “vườn thú nhân tố” là các phát hiện sai.
Đối với nghiên cứu cổ phiếu Việt Nam, vấn đề có lẽ còn nghiêm trọng hơn. Mặt cắt ngang nhỏ hơn (600-800 cổ phiếu so với hơn 4.000 ở Mỹ), chuỗi thời gian ngắn hơn (từ 2008 so với từ 1963), và dữ liệu nhiễu hơn (kém thanh khoản, các ngày không giao dịch, những bất thường về kế toán). Tất cả các yếu tố này làm tăng rủi ro của cả sai lầm loại I (tìm thấy điều không tồn tại) lẫn sai lầm loại II (bỏ sót điều có thật). Chương này cung cấp các công cụ thống kê để xử lý bức tranh đó một cách chặt chẽ.
Harvey, Liu, và Zhu (2016) đếm được hơn 300 nhân tố được công bố trên các tạp chí hàng đầu đến năm 2012. Chen và Zimmermann (2022) xây dựng và công bố công khai hơn 300 hiện tượng bất thường “có ý nghĩa rõ ràng” từ các tài liệu của Mỹ. Hou, Xue, và Zhang (2020) cố gắng tái lập 452 hiện tượng bất thường và thấy rằng 64% (với lợi nhuận theo trọng số vốn hóa, VW) không đạt \(t > 1.96\) trong mẫu của họ. Jensen, Kelly, và Pedersen (2023) có cái nhìn lạc quan hơn, khi thấy rằng phần lớn các nhân tố được tái lập về hướng, dù không phải lúc nào cũng về độ lớn.
Điểm mấu chốt là số nhân tố được kiểm định vượt xa số được công bố. Với mỗi nhân tố được công bố có \(t > 2\), có thể có hàng chục kiểm định không được công bố với \(t < 2\) chưa bao giờ được in ra, đó là vấn đề “ngăn kéo hồ sơ” (file drawer) kinh điển. Tổng số kiểm định mà giới nghiên cứu đã thực hiện là không thể biết, nhưng chắc chắn lên tới hàng nghìn.
Việc đếm các tín hiệu (\(|t| > 1.96\) theo cả hai định nghĩa) cũng đánh giá thấp vấn đề, vì mỗi tín hiệu có thể được đo lường theo nhiều cách hợp lý. Binsbergen, Ma, và Schwert (2026) nhấn mạnh điểm này rất rõ. Họ giữ nguyên các biến dùng để sắp xếp nhưng thay đổi một quy ước tưởng chừng vô hại: thay vì trừ đi lãi suất tín phiếu kho bạc kỳ hạn một tháng, họ trừ đi lợi nhuận của một danh mục trái phiếu chính phủ zero-coupon có kỳ hạn bình quân tương ứng, để một nhân tố được đánh giá so với một kịch bản dòng tiền cố định cùng kỳ hạn thay vì so với tiền mặt. Chạy lại 153 hiện tượng bất thường đã công bố trong các giai đoạn mẫu gốc, họ thấy 66 hiện tượng vẫn vững (\(|t| > 1.96\) theo cả hai định nghĩa) trong khi 20 hiện tượng đổi trạng thái, với khoảng 0,30 phát hiện sai trên mỗi phát hiện đúng. Mở rộng bài tập sang 1.395 cách sắp xếp giả định xây dựng từ các tỷ số Compustat, tỷ lệ giữa phát hiện sai và phát hiện đúng tăng lên khoảng 1,0 và ổn định qua các cửa sổ trượt. Ngay cả nhân tố thị trường cũng bị ảnh hưởng: nhân tố tương ứng đã điều chỉnh theo kỳ hạn có tương quan dưới 0,5 với lợi nhuận vượt mức thị trường thông thường, và phần bù cổ phiếu không còn khác không một cách có ý nghĩa thống kê. Cơ chế là sự sụt giảm dài hạn của lãi suất trong bốn thập kỷ, vốn mang lại một khoản lợi bất ngờ mang tính cơ học cho các vị thế kỳ hạn dài của danh mục mua–bán khống. Kết quả của họ không khẳng định định nghĩa truyền thống là sai, mà cho thấy quá trình phát hiện nhạy cảm với một lựa chọn mà hầu hết các bài báo không bao giờ báo cáo, đó chính là “đa vũ trụ” (multiverse) các bậc tự do của nhà nghiên cứu mà các hiệu chỉnh kiểm định đa giả thuyết không nắm bắt được.
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import statsmodels.api as sm
from scipy import stats
from statsmodels.stats.multitest import multipletests
from itertools import product
import warnings
warnings.filterwarnings('ignore')
plt.rcParams.update({
'figure.figsize': (12, 6),
'figure.dpi': 150,
'font.size': 11,
'axes.spines.top': False,
'axes.spines.right': False
})from datacore import DataCoreClient
client = DataCoreClient()
# Monthly returns (survivorship-bias-free)
monthly = client.get_monthly_returns(
exchanges=['HOSE', 'HNX'],
start_date='2008-01-01',
end_date='2024-12-31',
include_delisted=True,
fields=[
'ticker', 'month_end', 'monthly_return', 'market_cap',
'exchange', 'volume_avg_20d', 'turnover_value_avg_20d',
'n_zero_volume_days'
]
)
# Pre-computed characteristic panel (from factor construction chapter)
characteristics = client.get_characteristics_panel(
exchanges=['HOSE', 'HNX'],
start_date='2008-01-01',
end_date='2024-12-31',
include_delisted=True
)
# Factor returns for spanning tests
factors_ff = client.get_factor_returns(
market='vietnam',
start_date='2008-01-01',
end_date='2024-12-31',
factors=['mkt_excess', 'smb', 'hml', 'rmw', 'cma', 'wml']
)
monthly['month_end'] = pd.to_datetime(monthly['month_end'])
print(f"Monthly returns: {len(monthly):,}")
print(f"Characteristics: {characteristics.shape}")
print(f"Available signals: {[c for c in characteristics.columns if c not in ['ticker', 'month_end']]}")Chúng ta xây dựng một tập hợp toàn diện các biến bất thường bao quát các nhóm chính trong tài liệu nghiên cứu. Mỗi tín hiệu được lấy độ trễ phù hợp để tránh thiên lệch nhìn trước (các tín hiệu kế toán dùng cách căn chỉnh theo thời điểm công bố thực tế; các tín hiệu giá dùng giá trị của tháng trước).
# Merge returns with characteristics
panel = monthly.merge(characteristics, on=['ticker', 'month_end'], how='left')
panel = panel.sort_values(['ticker', 'month_end'])
# ============================================================
# Category 1: VALUE (8 signals)
# ============================================================
# Already in characteristics panel from PIT merge:
# bm, ep, cfp, sp, dp, ev_ebitda
# Add earnings yield variants
value_signals = ['bm', 'ep', 'cfp', 'sp', 'dp']
# ============================================================
# Category 2: SIZE (3 signals)
# ============================================================
panel['log_mcap'] = np.log(panel['market_cap'].clip(lower=1))
panel['log_assets'] = np.log(panel['total_assets'].clip(lower=1))
panel['log_revenue'] = np.log(panel['revenue'].clip(lower=1))
size_signals = ['log_mcap', 'log_assets', 'log_revenue']
# ============================================================
# Category 3: PROFITABILITY (8 signals)
# ============================================================
# GP/A (Novy-Marx 2013), OP (FF 2015), ROE, ROA, etc.
prof_signals = ['gp_at', 'op', 'roe', 'roa', 'profit_margin',
'asset_turnover', 'gross_margin', 'oper_margin']
# ============================================================
# Category 4: INVESTMENT / GROWTH (6 signals)
# ============================================================
inv_signals = ['asset_growth', 'capex_at', 'investment',
'sales_growth', 'equity_issuance', 'debt_issuance']
# ============================================================
# Category 5: MOMENTUM / REVERSAL (6 signals)
# ============================================================
# Compute from returns
for lag_start, lag_end, name in [
(2, 12, 'ret_12_2'), (2, 6, 'ret_6_2'), (7, 12, 'ret_12_7'),
(1, 1, 'ret_1'), (1, 3, 'ret_3_1')
]:
if name not in panel.columns:
panel[name] = (
panel.groupby('ticker')['monthly_return']
.transform(lambda x: x.shift(lag_start).rolling(lag_end - lag_start + 1)
.apply(lambda r: (1 + r).prod() - 1, raw=True))
)
# Earnings momentum (SUE) should already be in characteristics
mom_signals = ['ret_12_2', 'ret_6_2', 'ret_12_7', 'ret_1', 'ret_3_1', 'sue']
# ============================================================
# Category 6: RISK (5 signals)
# ============================================================
risk_signals = ['beta', 'ivol', 'tvol', 'max_ret', 'skewness']
# ============================================================
# Category 7: LIQUIDITY / TRADING (6 signals)
# ============================================================
liq_signals = ['amihud', 'zero_return_pct', 'log_turnover',
'roll_spread', 'cs_spread', 'volume_cv']
# ============================================================
# Category 8: QUALITY / FINANCIAL HEALTH (5 signals)
# ============================================================
qual_signals = ['leverage', 'current_ratio', 'accruals',
'earnings_variability', 'piotroski_f']
# ============================================================
# Category 9: DIVIDEND / PAYOUT (3 signals)
# ============================================================
div_signals = ['div_yield', 'payout_ratio', 'net_repurchase']
# Combine all
all_signals = (value_signals + size_signals + prof_signals +
inv_signals + mom_signals + risk_signals +
liq_signals + qual_signals + div_signals)
# Remove signals with too little coverage
signal_coverage = {}
for sig in all_signals:
if sig in panel.columns:
cov = panel[sig].notna().mean()
signal_coverage[sig] = cov
signal_coverage = pd.Series(signal_coverage).sort_values(ascending=False)
valid_signals = signal_coverage[signal_coverage > 0.30].index.tolist()
print(f"Total signals defined: {len(all_signals)}")
print(f"Signals with >30% coverage: {len(valid_signals)}")
print(f"\nSignal categories:")
for cat_name, sigs in [
('Value', value_signals), ('Size', size_signals),
('Profitability', prof_signals), ('Investment', inv_signals),
('Momentum', mom_signals), ('Risk', risk_signals),
('Liquidity', liq_signals), ('Quality', qual_signals),
('Dividend', div_signals)
]:
n_valid = sum(1 for s in sigs if s in valid_signals)
print(f" {cat_name:<15}: {n_valid}/{len(sigs)}")Chúng ta chạy bộ máy xây dựng nhân tố (từ chương trước) trên tất cả các tín hiệu hợp lệ để tạo ra một “vườn thú” lợi nhuận nhân tố:
# Import the factor engine from previous chapter
# (In practice, this would be imported from a shared module)
def construct_factor_simple(panel_df, signal_col, long_high=True,
n_groups=5, weighting='value',
min_stocks=20):
"""
Simplified factor construction: quintile sort on signal,
long-short = Q5 - Q1 (or reversed if long_high=False).
Returns a Series of monthly long-short returns.
"""
df = panel_df[['ticker', 'month_end', 'monthly_return',
'market_cap', signal_col]].dropna().copy()
df = df.sort_values(['ticker', 'month_end'])
# Lag the signal by one month
df['signal_lag'] = df.groupby('ticker')[signal_col].shift(1)
df = df.dropna(subset=['signal_lag', 'monthly_return'])
results = []
for month, group in df.groupby('month_end'):
if len(group) < min_stocks:
continue
group['quintile'] = pd.qcut(
group['signal_lag'].rank(method='first'),
n_groups, labels=False, duplicates='drop'
)
if weighting == 'value':
def wret(g):
if g['market_cap'].sum() > 0:
return np.average(g['monthly_return'],
weights=g['market_cap'])
return g['monthly_return'].mean()
port_ret = group.groupby('quintile').apply(wret)
else:
port_ret = group.groupby('quintile')['monthly_return'].mean()
if 0 in port_ret.index and (n_groups - 1) in port_ret.index:
if long_high:
ls = port_ret[n_groups - 1] - port_ret[0]
else:
ls = port_ret[0] - port_ret[n_groups - 1]
results.append({'month_end': month, 'ls_return': ls})
if not results:
return None
return pd.DataFrame(results).set_index('month_end')['ls_return']
# Determine expected sign for each signal
# (positive = high signal predicts high returns)
signal_directions = {
# Value: high = higher returns
'bm': True, 'ep': True, 'cfp': True, 'sp': True, 'dp': True,
# Size: small = higher returns (low signal = high returns)
'log_mcap': False, 'log_assets': False, 'log_revenue': False,
# Profitability: high = higher returns
'gp_at': True, 'op': True, 'roe': True, 'roa': True,
'profit_margin': True, 'asset_turnover': True,
'gross_margin': True, 'oper_margin': True,
# Investment: low growth = higher returns
'asset_growth': False, 'capex_at': False, 'investment': False,
'sales_growth': False, 'equity_issuance': False,
'debt_issuance': False,
# Momentum: high = higher returns (except reversal)
'ret_12_2': True, 'ret_6_2': True, 'ret_12_7': True,
'ret_1': False, 'ret_3_1': False, 'sue': True,
# Risk: low risk = higher returns (anomaly direction)
'beta': False, 'ivol': False, 'tvol': False,
'max_ret': False, 'skewness': False,
# Liquidity: illiquid = higher returns
'amihud': True, 'zero_return_pct': True,
'log_turnover': False, 'roll_spread': True,
'cs_spread': True, 'volume_cv': True,
# Quality: high quality = higher returns
'leverage': False, 'current_ratio': True, 'accruals': False,
'earnings_variability': False, 'piotroski_f': True,
# Dividend: high = higher returns
'div_yield': True, 'payout_ratio': True, 'net_repurchase': False,
}
# Construct all factors
zoo_results = {}
for sig in valid_signals:
long_high = signal_directions.get(sig, True)
ls = construct_factor_simple(
panel, sig, long_high=long_high,
n_groups=5, weighting='value'
)
if ls is not None and len(ls) >= 60:
mean_ret = ls.mean()
se = ls.std() / np.sqrt(len(ls))
t = mean_ret / se if se > 0 else 0
zoo_results[sig] = {
'mean_monthly': mean_ret,
'ann_return': mean_ret * 12,
'ann_vol': ls.std() * np.sqrt(12),
'sharpe': mean_ret / ls.std() * np.sqrt(12) if ls.std() > 0 else 0,
't_stat': t,
'abs_t': abs(t),
'p_value': 2 * (1 - stats.t.cdf(abs(t), df=len(ls) - 1)),
'n_months': len(ls),
'direction': 'Long High' if long_high else 'Long Low',
'category': next(
(cat for cat, sigs in [
('Value', value_signals), ('Size', size_signals),
('Profitability', prof_signals),
('Investment', inv_signals),
('Momentum', mom_signals), ('Risk', risk_signals),
('Liquidity', liq_signals),
('Quality', qual_signals),
('Dividend', div_signals)
] if sig in sigs), 'Other'
),
'returns': ls # Store for later use
}
zoo_df = pd.DataFrame({k: {kk: vv for kk, vv in v.items() if kk != 'returns'}
for k, v in zoo_results.items()}).T
zoo_df = zoo_df.sort_values('abs_t', ascending=False)
print(f"Factors successfully constructed: {len(zoo_df)}")
print(f"Factors with |t| > 2.0: {(zoo_df['abs_t'] > 2.0).sum()} "
f"({(zoo_df['abs_t'] > 2.0).mean():.1%})")
print(f"Factors with |t| > 3.0: {(zoo_df['abs_t'] > 3.0).sum()} "
f"({(zoo_df['abs_t'] > 3.0).mean():.1%})")fig, axes = plt.subplots(1, 2, figsize=(16, 6))
# Panel A: Histogram
t_vals = zoo_df['t_stat'].values
axes[0].hist(t_vals, bins=30, color='#2C5F8A', alpha=0.7,
edgecolor='white', density=True)
# Overlay normal distribution under null
x_range = np.linspace(-5, 5, 200)
axes[0].plot(x_range, stats.norm.pdf(x_range), 'k--', linewidth=1.5,
label='N(0,1) null')
# Threshold lines
axes[0].axvline(x=1.96, color='#E67E22', linestyle='--', linewidth=1.5,
label='|t| = 1.96')
axes[0].axvline(x=-1.96, color='#E67E22', linestyle='--', linewidth=1.5)
axes[0].axvline(x=3.0, color='#C0392B', linestyle='--', linewidth=1.5,
label='|t| = 3.0 (HLZ)')
axes[0].axvline(x=-3.0, color='#C0392B', linestyle='--', linewidth=1.5)
axes[0].set_xlabel('t-statistic')
axes[0].set_ylabel('Density')
axes[0].set_title('Panel A: t-Statistic Distribution')
axes[0].legend(fontsize=9)
# Panel B: Ranked bar chart
top_n = min(40, len(zoo_df))
top = zoo_df.head(top_n).copy()
category_colors = {
'Value': '#2C5F8A', 'Size': '#1ABC9C', 'Profitability': '#27AE60',
'Investment': '#E67E22', 'Momentum': '#C0392B', 'Risk': '#8E44AD',
'Liquidity': '#3498DB', 'Quality': '#F1C40F', 'Dividend': '#95A5A6'
}
bar_colors = [category_colors.get(top.iloc[i]['category'], '#BDC3C7')
for i in range(len(top))]
axes[1].barh(range(top_n), top['abs_t'].values,
color=bar_colors, alpha=0.85, edgecolor='white')
axes[1].set_yticks(range(top_n))
axes[1].set_yticklabels(top.index, fontsize=7)
axes[1].axvline(x=1.96, color='#E67E22', linestyle='--', linewidth=1)
axes[1].axvline(x=3.0, color='#C0392B', linestyle='--', linewidth=1)
axes[1].set_xlabel('|t-statistic|')
axes[1].set_title('Panel B: Anomalies Ranked by |t|')
axes[1].invert_yaxis()
# Legend for categories
from matplotlib.patches import Patch
legend_elements = [Patch(facecolor=c, label=cat)
for cat, c in category_colors.items()]
axes[1].legend(handles=legend_elements, fontsize=7,
loc='lower right', ncol=2)
plt.tight_layout()
plt.show()Giả sử bạn kiểm định \(M\) giả thuyết độc lập, mỗi giả thuyết ở mức ý nghĩa \(\alpha = 0.05\). Dưới giả thuyết không toàn cục (không có hiệu ứng thực), số lần bác bỏ sai kỳ vọng là \(M \times \alpha\). Với \(M = 50\) hiện tượng bất thường, bạn kỳ vọng có 2,5 kết quả “có ý nghĩa” hoàn toàn do ngẫu nhiên.
Xác suất có ít nhất một lần bác bỏ sai là:
\[ P(\text{at least one false discovery}) = 1 - (1 - \alpha)^M \tag{21.1}\]
Với \(M = 50\) và \(\alpha = 0.05\), giá trị này là \(1 - 0.95^{50} = 0.923\). Với \(M = 300\), nó về cơ bản bằng 1. Giá trị p của một kiểm định đơn lẻ vô dụng trong việc quyết định nhân tố nào trong số nhiều nhân tố được kiểm định là thật.
Các hiệu chỉnh kiểm định đa giả thuyết kiểm soát một trong hai tỷ lệ sai sót:
Tỷ lệ sai sót theo họ giả thuyết (FWER). Xác suất mắc ít nhất một sai lầm loại I trên toàn bộ các kiểm định. Đây là tiêu chí thận trọng nhất. Nếu FWER = 0,05, bạn có thể tin chắc 95% rằng mọi giả thuyết không đã bị bác bỏ đều là phát hiện đúng. Các phương pháp: hiệu chỉnh Bonferroni (1936), Holm step-down, Romano và Wolf (2005).
Tỷ lệ phát hiện sai (FDR). Tỷ lệ kỳ vọng các lần bác bỏ sai trong số các giả thuyết bị bác bỏ. Tiêu chí này ít thận trọng hơn. Nếu FDR = 0,05, bạn kỳ vọng 5% số phát hiện của mình là sai, còn 95% còn lại là thật. Các phương pháp: Benjamini và Hochberg (1995) (BH), Storey (2003).
Trong định giá tài sản, kiểm soát FDR nhìn chung phù hợp hơn FWER vì chúng ta không cố đảm bảo không có phát hiện sai nào, mà cố xác định một tập các nhân tố trong đó phần lớn là thật. Một nhà nghiên cứu kiểm soát FDR ở mức 5% và phát hiện 10 nhân tố sẽ kỳ vọng khoảng 0,5 nhân tố trong số đó là sai, điều này chấp nhận được khi xây dựng một mô hình nhân tố.
M_range = np.arange(1, 301)
alpha = 0.05
fwer_unadj = 1 - (1 - alpha) ** M_range
fwer_bonf = np.minimum(1, M_range * alpha) # Not actual FWER, just threshold
# Expected false discoveries at alpha = 0.05
expected_false = M_range * alpha
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
axes[0].plot(M_range, fwer_unadj, color='#C0392B', linewidth=2,
label='Unadjusted FWER')
axes[0].axhline(y=0.05, color='gray', linestyle='--', linewidth=1,
label='α = 0.05 target')
axes[0].set_xlabel('Number of Tests (M)')
axes[0].set_ylabel('P(≥1 False Discovery)')
axes[0].set_title('Panel A: Family-Wise Error Rate')
axes[0].legend()
axes[1].plot(M_range, expected_false, color='#2C5F8A', linewidth=2)
axes[1].set_xlabel('Number of Tests (M)')
axes[1].set_ylabel('Expected False Discoveries')
axes[1].set_title('Panel B: Expected False Discoveries at α = 0.05')
axes[1].axhline(y=1, color='gray', linestyle='--', linewidth=0.8)
axes[1].text(250, 1.5, '1 false discovery', color='gray', fontsize=9)
plt.tight_layout()
plt.show()Hiệu chỉnh đơn giản nhất và thận trọng nhất: chỉ bác bỏ giả thuyết \(i\) nếu giá trị p của nó thấp hơn \(\alpha / M\), trong đó \(M\) là tổng số kiểm định:
\[ \text{Reject } H_{0,i} \text{ if } p_i < \frac{\alpha}{M} \tag{21.2}\]
Phương pháp này kiểm soát FWER ở mức \(\alpha\) bất kể cấu trúc phụ thuộc giữa các kiểm định. Cái giá phải trả là mất rất nhiều sức mạnh kiểm định: với \(M = 50\) kiểm định ở \(\alpha = 0.05\), ngưỡng cho mỗi kiểm định trở thành \(0.001\), đòi hỏi \(|t| > 3.29\).
M = len(zoo_df)
alpha = 0.05
# Bonferroni threshold
bonf_threshold = alpha / M
bonf_t_threshold = stats.norm.ppf(1 - bonf_threshold / 2)
zoo_df['bonf_reject'] = zoo_df['p_value'] < bonf_threshold
print(f"Number of tests: {M}")
print(f"Bonferroni p-value threshold: {bonf_threshold:.6f}")
print(f"Equivalent t-statistic threshold: {bonf_t_threshold:.2f}")
print(f"Rejections: {zoo_df['bonf_reject'].sum()} / {M}")
print(f"\nSurviving anomalies:")
survivors = zoo_df[zoo_df['bonf_reject']]
if len(survivors) > 0:
print(survivors[['ann_return', 't_stat', 'category']].to_string())
else:
print(" None survive Bonferroni correction")Thủ tục Holm luôn mạnh hơn Bonferroni mà vẫn kiểm soát FWER. Nó xếp hạng các giá trị p từ nhỏ nhất đến lớn nhất và áp dụng các ngưỡng nới lỏng dần:
p_values = zoo_df['p_value'].values
signal_names = zoo_df.index.values
# Using statsmodels implementation
reject_holm, pvals_corrected_holm, _, _ = multipletests(
p_values, alpha=0.05, method='holm'
)
zoo_df['holm_reject'] = reject_holm
zoo_df['holm_p_corrected'] = pvals_corrected_holm
print(f"Holm rejections: {reject_holm.sum()} / {M}")
holm_survivors = zoo_df[zoo_df['holm_reject']]
if len(holm_survivors) > 0:
print("\nSurviving anomalies (Holm):")
print(holm_survivors[['ann_return', 't_stat', 'holm_p_corrected',
'category']].to_string())Thủ tục Benjamini và Hochberg (1995) (BH) kiểm soát tỷ lệ phát hiện sai (tức tỷ lệ kỳ vọng các lần bác bỏ sai trong số các giả thuyết bị bác bỏ) ở mức \(q\):
Thủ tục BH hợp lệ khi các kiểm định độc lập hoặc phụ thuộc dương (PRDS), một điều kiện nhìn chung được thỏa mãn đối với các kiểm định nhân tố tài chính, vì các nhân tố có xu hướng tương quan dương dưới giả thuyết không.
# BH procedure at q = 0.05 (expect 5% false discoveries)
reject_bh, pvals_corrected_bh, _, _ = multipletests(
p_values, alpha=0.05, method='fdr_bh'
)
zoo_df['bh_reject'] = reject_bh
zoo_df['bh_p_corrected'] = pvals_corrected_bh
# Also at q = 0.10 (more liberal)
reject_bh10, _, _, _ = multipletests(
p_values, alpha=0.10, method='fdr_bh'
)
zoo_df['bh10_reject'] = reject_bh10
print(f"BH rejections (q=0.05): {reject_bh.sum()} / {M}")
print(f"BH rejections (q=0.10): {reject_bh10.sum()} / {M}")
bh_survivors = zoo_df[zoo_df['bh_reject']].sort_values('abs_t', ascending=False)
if len(bh_survivors) > 0:
print(f"\nSurviving anomalies (BH, q=0.05):")
print(bh_survivors[['ann_return', 't_stat', 'bh_p_corrected',
'category']].head(20).to_string())sorted_pvals = np.sort(p_values)
k = np.arange(1, M + 1)
bh_line_05 = k / M * 0.05
bh_line_10 = k / M * 0.10
fig, ax = plt.subplots(figsize=(12, 6))
ax.scatter(k, sorted_pvals, s=15, color='#2C5F8A', alpha=0.7,
label='Sorted p-values', zorder=3)
ax.plot(k, bh_line_05, color='#C0392B', linewidth=2,
label='BH line (q=0.05)')
ax.plot(k, bh_line_10, color='#E67E22', linewidth=2,
linestyle='--', label='BH line (q=0.10)')
# Highlight rejections
n_reject_05 = reject_bh.sum()
if n_reject_05 > 0:
ax.scatter(k[:n_reject_05], sorted_pvals[:n_reject_05],
s=40, color='#C0392B', zorder=4, label=f'Rejected (q=0.05)')
ax.set_xlabel('Rank (k)')
ax.set_ylabel('p-value')
ax.set_title('Benjamini-Hochberg Procedure')
ax.legend()
ax.set_ylim([-0.01, 0.15])
ax.set_xlim([0, min(M, 60)])
plt.tight_layout()
plt.show()Harvey, Liu, và Zhu (2016) tiếp cận theo cách khác: thay vì áp dụng một hiệu chỉnh kiểm định đa giả thuyết chính thức cho một tập kiểm định cụ thể, họ ước lượng thống kê t ngưỡng có tính đến số lượng nhân tố tích lũy đã được toàn giới nghiên cứu kiểm định. Sử dụng khung Bayes được hiệu chỉnh theo 316 nhân tố đã công bố, họ suy ra:
Logic ở đây là ngay cả khi nghiên cứu của bạn chỉ kiểm định một nhân tố, việc hàng trăm nhà nghiên cứu đã kiểm định hàng trăm tín hiệu trước bạn có nghĩa là xác suất tiên nghiệm để tín hiệu cụ thể của bạn là một biến dự báo thực sự là thấp.
hlz_thresholds = {
'Conventional (|t| > 2.0)': 2.0,
'HLZ 100 tests (|t| > 2.8)': 2.8,
'HLZ 316 tests (|t| > 3.0)': 3.0,
'HLZ conservative (|t| > 3.78)': 3.78
}
print("Anomalies Surviving Different t-Thresholds:")
print(f"{'Threshold':<35} {'Surviving':>10} {'% of Zoo':>10}")
print("-" * 55)
for name, thresh in hlz_thresholds.items():
n_survive = (zoo_df['abs_t'] > thresh).sum()
pct = n_survive / len(zoo_df) * 100
print(f"{name:<35} {n_survive:>10} {pct:>10.1f}%")comparison = {
'Unadjusted (|t| > 2)': (zoo_df['abs_t'] > 2.0).sum(),
'Bonferroni': zoo_df['bonf_reject'].sum(),
'Holm': zoo_df['holm_reject'].sum(),
'BH (q=0.05)': zoo_df['bh_reject'].sum(),
'BH (q=0.10)': zoo_df['bh10_reject'].sum(),
'HLZ (|t| > 3.0)': (zoo_df['abs_t'] > 3.0).sum(),
'HLZ (|t| > 3.78)': (zoo_df['abs_t'] > 3.78).sum(),
}
fig, ax = plt.subplots(figsize=(10, 5))
names = list(comparison.keys())
counts = list(comparison.values())
colors = ['#BDC3C7', '#C0392B', '#E74C3C', '#27AE60', '#2ECC71',
'#2C5F8A', '#1A3C5A']
bars = ax.barh(range(len(names)), counts, color=colors,
alpha=0.85, edgecolor='white')
ax.set_yticks(range(len(names)))
ax.set_yticklabels(names)
ax.set_xlabel('Number of Surviving Anomalies')
ax.set_title('Anomalies Surviving Multiple Testing Corrections')
for i, (name, count) in enumerate(zip(names, counts)):
ax.text(count + 0.3, i, str(count), va='center', fontsize=10)
ax.invert_yaxis()
plt.tight_layout()
plt.show()cat_survival = (
zoo_df.groupby('category')
.agg(
n_total=('abs_t', 'count'),
n_survive_bh=('bh_reject', 'sum'),
n_survive_hlz=('abs_t', lambda x: (x > 3.0).sum()),
mean_abs_t=('abs_t', 'mean')
)
.sort_values('mean_abs_t', ascending=False)
)
cat_survival['survival_rate_bh'] = (
cat_survival['n_survive_bh'] / cat_survival['n_total']
)
fig, ax = plt.subplots(figsize=(12, 5))
x = np.arange(len(cat_survival))
ax.bar(x, cat_survival['n_total'], color='#BDC3C7', alpha=0.5,
label='Total tested', edgecolor='white')
ax.bar(x, cat_survival['n_survive_bh'], color='#27AE60', alpha=0.85,
label='Survive BH (q=0.05)', edgecolor='white')
ax.set_xticks(x)
ax.set_xticklabels(cat_survival.index, rotation=30, fontsize=9)
ax.set_ylabel('Number of Anomalies')
ax.set_title('Anomaly Survival by Category')
ax.legend()
# Add survival rate text
for i, (_, row) in enumerate(cat_survival.iterrows()):
if row['n_total'] > 0:
ax.text(i, row['n_total'] + 0.2,
f"{row['survival_rate_bh']:.0%}",
ha='center', fontsize=9, color='#27AE60')
plt.tight_layout()
plt.show()Thay vì kiểm soát FDR ở một mức định trước, chúng ta có thể ước lượng tỷ lệ các giả thuyết được kiểm định thực sự là giả thuyết không (\(\pi_0\)) và tỷ lệ phát hiện sai hàm ý (FDP) tại một ngưỡng bất kỳ. Storey (2003) đề xuất ước lượng \(\pi_0\) từ phân phối của các giá trị p: dưới giả thuyết không, các giá trị p phân phối đều trên \([0, 1]\), nên mật độ của các giá trị p trong vùng “phẳng” (xa số không) cho biết \(\pi_0\).
\[ \hat{\pi}_0(\lambda) = \frac{\#\{p_i > \lambda\}}{M(1 - \lambda)} \tag{21.3}\]
với một tham số điều chỉnh \(\lambda\). Tỷ lệ phát hiện sai ước lượng tại ngưỡng \(t^*\) khi đó là:
\[ \widehat{\text{FDP}}(t^*) = \frac{\hat{\pi}_0 \cdot M \cdot P(|t| > t^* | H_0)}{\#\{|t_i| > t^*\}} \tag{21.4}\]
def estimate_pi0(p_values, lambdas=np.arange(0.05, 0.95, 0.05)):
"""
Estimate pi_0 (proportion of true nulls) using Storey (2003).
Uses bootstrap to select optimal lambda.
"""
M = len(p_values)
pi0_estimates = []
for lam in lambdas:
n_above = (p_values > lam).sum()
pi0 = n_above / (M * (1 - lam))
pi0_estimates.append({'lambda': lam, 'pi0': min(pi0, 1.0)})
pi0_df = pd.DataFrame(pi0_estimates)
# Use the median of estimates for robustness
# (alternatives: bootstrap, spline smoothing)
pi0_hat = pi0_df['pi0'].median()
return pi0_hat, pi0_df
pi0_hat, pi0_df = estimate_pi0(zoo_df['p_value'].values)
print(f"Estimated π₀ (proportion of true nulls): {pi0_hat:.3f}")
print(f"Implied proportion of true factors: {1 - pi0_hat:.3f}")
print(f"Estimated number of true factors: {(1 - pi0_hat) * len(zoo_df):.0f} "
f"out of {len(zoo_df)}")fig, axes = plt.subplots(1, 2, figsize=(14, 5))
# Panel A: p-value histogram
axes[0].hist(zoo_df['p_value'], bins=20, density=True,
color='#2C5F8A', alpha=0.7, edgecolor='white')
axes[0].axhline(y=pi0_hat, color='#C0392B', linewidth=2, linestyle='--',
label=f'π₀ = {pi0_hat:.2f} (null density)')
axes[0].axhline(y=1, color='gray', linewidth=1, linestyle=':',
label='Uniform (all null)')
axes[0].set_xlabel('p-value')
axes[0].set_ylabel('Density')
axes[0].set_title('Panel A: p-Value Distribution')
axes[0].legend()
# Panel B: pi_0 estimates across lambda
axes[1].plot(pi0_df['lambda'], pi0_df['pi0'],
color='#2C5F8A', linewidth=2, marker='o', markersize=4)
axes[1].axhline(y=pi0_hat, color='#C0392B', linestyle='--',
linewidth=1.5, label=f'Median π₀ = {pi0_hat:.2f}')
axes[1].set_xlabel('λ')
axes[1].set_ylabel('π₀(λ)')
axes[1].set_title('Panel B: π₀ Estimates by λ')
axes[1].legend()
axes[1].set_ylim([0, 1.05])
plt.tight_layout()
plt.show()Sử dụng \(\hat{\pi}_0\), chúng ta tính FDP ước lượng tại mọi ngưỡng t có thể:
t_thresholds = np.arange(1.0, 5.01, 0.1)
fdp_curve = []
for t_thresh in t_thresholds:
n_rejected = (zoo_df['abs_t'] > t_thresh).sum()
if n_rejected == 0:
fdp_curve.append({'t_threshold': t_thresh, 'n_rejected': 0,
'fdp': 0})
continue
# Expected false discoveries under null
p_null_reject = 2 * (1 - stats.norm.cdf(t_thresh))
expected_false = pi0_hat * len(zoo_df) * p_null_reject
fdp = min(expected_false / n_rejected, 1.0)
fdp_curve.append({'t_threshold': t_thresh, 'n_rejected': n_rejected,
'fdp': fdp})
fdp_df = pd.DataFrame(fdp_curve)
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
# Panel A: FDP curve
axes[0].plot(fdp_df['t_threshold'], fdp_df['fdp'] * 100,
color='#C0392B', linewidth=2)
axes[0].axhline(y=5, color='gray', linestyle='--', linewidth=1,
label='5% FDP')
axes[0].axhline(y=10, color='gray', linestyle=':', linewidth=1,
label='10% FDP')
axes[0].axvline(x=2.0, color='#E67E22', linestyle='--',
linewidth=1, label='|t| = 2.0')
axes[0].axvline(x=3.0, color='#2C5F8A', linestyle='--',
linewidth=1, label='|t| = 3.0')
axes[0].set_xlabel('|t| Threshold')
axes[0].set_ylabel('Estimated FDP (%)')
axes[0].set_title('Panel A: False Discovery Proportion')
axes[0].legend(fontsize=8)
axes[0].set_ylim([-2, 80])
# Panel B: Number of rejections
axes[1].plot(fdp_df['t_threshold'], fdp_df['n_rejected'],
color='#2C5F8A', linewidth=2)
axes[1].set_xlabel('|t| Threshold')
axes[1].set_ylabel('Number of Discoveries')
axes[1].set_title('Panel B: Discoveries vs Threshold')
plt.tight_layout()
plt.show()Các hiệu chỉnh giải tích giả định các phân phối đã biết. Khi phân phối lợi nhuận không chuẩn (đuôi dày, lệch, tự tương quan chuỗi), như trường hợp ở Việt Nam, các phương pháp bootstrap cho suy luận đáng tin cậy hơn.
White (2000) đề xuất một kiểm tra thực tế bằng bootstrap (reality check) có tính đến việc khai thác dữ liệu (data snooping) trên nhiều chiến lược. Giả thuyết không là chiến lược tốt nhất có lợi nhuận kỳ vọng bằng không. Thủ tục như sau:
def white_reality_check(factor_returns_dict, n_bootstrap=1000,
block_size=6, seed=42):
"""
White (2000) Reality Check for data snooping.
Parameters
----------
factor_returns_dict : dict
{factor_name: pd.Series of monthly returns}
n_bootstrap : int
Number of bootstrap replications
block_size : int
Block length for circular block bootstrap
Returns
-------
DataFrame with original t-stats and bootstrap p-values
"""
rng = np.random.default_rng(seed)
# Align all factor returns to common dates
all_names = list(factor_returns_dict.keys())
common_dates = None
for name in all_names:
dates = set(factor_returns_dict[name].index)
common_dates = dates if common_dates is None else common_dates & dates
common_dates = sorted(common_dates)
T = len(common_dates)
# Build return matrix (T x M)
M = len(all_names)
return_matrix = np.column_stack([
factor_returns_dict[name].reindex(common_dates).values
for name in all_names
])
# Observed t-statistics
means = np.nanmean(return_matrix, axis=0)
ses = np.nanstd(return_matrix, axis=0) / np.sqrt(T)
t_obs = means / np.where(ses > 0, ses, np.nan)
# Block bootstrap
n_blocks = int(np.ceil(T / block_size))
max_t_bootstrap = np.zeros(n_bootstrap)
for b in range(n_bootstrap):
# Circular block bootstrap
block_starts = rng.integers(0, T, size=n_blocks)
indices = np.concatenate([
np.arange(start, start + block_size) % T
for start in block_starts
])[:T]
boot_returns = return_matrix[indices, :]
# Center under null (subtract original mean)
boot_centered = boot_returns - means[np.newaxis, :]
# Compute t-stats under null
boot_means = np.nanmean(boot_centered, axis=0)
boot_ses = np.nanstd(boot_centered, axis=0) / np.sqrt(T)
boot_t = boot_means / np.where(boot_ses > 0, boot_ses, np.nan)
max_t_bootstrap[b] = np.nanmax(np.abs(boot_t))
# Bootstrap p-values (one-sided: is obs_t extreme relative to max?)
boot_pvals = np.array([
np.mean(max_t_bootstrap >= abs(t)) if np.isfinite(t) else 1.0
for t in t_obs
])
results = pd.DataFrame({
'factor': all_names,
't_stat_obs': t_obs,
'boot_pval': boot_pvals
}).set_index('factor')
return results, max_t_bootstrap
# Collect factor return series
factor_returns_dict = {
name: zoo_results[name]['returns']
for name in zoo_results
if 'returns' in zoo_results[name]
}
boot_results, max_t_dist = white_reality_check(
factor_returns_dict, n_bootstrap=1000, block_size=6
)
boot_results = boot_results.sort_values('t_stat_obs', ascending=False)
print("White Reality Check Results (top 20):")
print(boot_results.head(20).round(4).to_string())
print(f"\nFactors surviving (boot p < 0.05): "
f"{(boot_results['boot_pval'] < 0.05).sum()}")fig, ax = plt.subplots(figsize=(12, 5))
ax.hist(max_t_dist, bins=40, density=True, color='#BDC3C7',
alpha=0.7, edgecolor='white', label='Bootstrap null\n(max |t| across all factors)')
# 95th percentile
p95 = np.percentile(max_t_dist, 95)
ax.axvline(x=p95, color='#C0392B', linewidth=2, linestyle='--',
label=f'95th percentile = {p95:.2f}')
# Observed top factor t-statistics
top_5 = boot_results.head(5)
for i, (name, row) in enumerate(top_5.iterrows()):
ax.axvline(x=abs(row['t_stat_obs']),
color=plt.cm.Set1(i), linewidth=1.5, linestyle=':',
label=f'{name}: |t| = {abs(row["t_stat_obs"]):.2f}')
ax.set_xlabel('Maximum |t-statistic|')
ax.set_ylabel('Density')
ax.set_title("White's Reality Check: Bootstrap Null Distribution")
ax.legend(fontsize=8, loc='upper right')
plt.tight_layout()
plt.show()Thủ tục step-down của Romano và Wolf (2005) cải tiến White Reality Check bằng cách lần lượt loại bỏ các giả thuyết đã bị bác bỏ và kiểm định lại các giả thuyết còn lại, nhờ đó tăng sức mạnh ở mỗi bước:
def romano_wolf_stepdown(factor_returns_dict, n_bootstrap=1000,
block_size=6, alpha=0.05, seed=42):
"""
Romano-Wolf (2005) step-down procedure for FWER control.
More powerful than single-step methods because rejected
hypotheses are removed before re-testing.
"""
rng = np.random.default_rng(seed)
all_names = list(factor_returns_dict.keys())
common_dates = None
for name in all_names:
dates = set(factor_returns_dict[name].index)
common_dates = dates if common_dates is None else common_dates & dates
common_dates = sorted(common_dates)
T = len(common_dates)
M = len(all_names)
return_matrix = np.column_stack([
factor_returns_dict[name].reindex(common_dates).values
for name in all_names
])
means = np.nanmean(return_matrix, axis=0)
ses = np.nanstd(return_matrix, axis=0) / np.sqrt(T)
t_obs = np.abs(means / np.where(ses > 0, ses, np.nan))
# Sort by observed t-stat (descending)
order = np.argsort(-t_obs)
t_sorted = t_obs[order]
names_sorted = [all_names[i] for i in order]
# Step-down procedure
rejected = set()
remaining = set(range(M))
for step in range(M):
if not remaining:
break
remaining_idx = sorted(remaining)
# Bootstrap max t-stat among remaining
n_blocks = int(np.ceil(T / block_size))
max_t_boot = np.zeros(n_bootstrap)
for b in range(n_bootstrap):
block_starts = rng.integers(0, T, size=n_blocks)
indices = np.concatenate([
np.arange(start, start + block_size) % T
for start in block_starts
])[:T]
boot_returns = return_matrix[indices][:, remaining_idx]
boot_centered = boot_returns - means[remaining_idx]
boot_means = np.nanmean(boot_centered, axis=0)
boot_ses = np.nanstd(boot_centered, axis=0) / np.sqrt(T)
boot_t = np.abs(boot_means / np.where(boot_ses > 0, boot_ses, np.nan))
max_t_boot[b] = np.nanmax(boot_t)
# Critical value
cv = np.percentile(max_t_boot, (1 - alpha) * 100)
# Test the most significant remaining factor
# Find the factor with highest t among remaining
best_remaining = max(remaining, key=lambda j: t_obs[j])
if t_obs[best_remaining] > cv:
rejected.add(best_remaining)
remaining.remove(best_remaining)
else:
break # Cannot reject any more
results = pd.DataFrame({
'factor': all_names,
't_stat': t_obs,
'rejected': [i in rejected for i in range(M)]
}).sort_values('t_stat', ascending=False)
return results
rw_results = romano_wolf_stepdown(
factor_returns_dict, n_bootstrap=500, block_size=6
)
print(f"Romano-Wolf rejections: {rw_results['rejected'].sum()} / {len(rw_results)}")
print("\nRejected factors:")
print(rw_results[rw_results['rejected']][['factor', 't_stat']].to_string())McLean và Pontiff (2016) thấy rằng lợi nhuận của các hiện tượng bất thường giảm 32% sau khi công bố trên thị trường Mỹ, gợi ý rằng khoảng một phần ba phần bù là do định giá sai được điều chỉnh nhờ giao dịch của các nhà đầu tư có thông tin. Với Việt Nam, chúng ta dùng cách chia theo chuỗi thời gian làm phép xấp xỉ cho kiểm định ngoài mẫu:
# Split at midpoint of available data
all_dates = sorted(panel['month_end'].unique())
mid_date = all_dates[len(all_dates) // 2]
oos_comparison = []
for name in zoo_results:
ls = zoo_results[name]['returns']
in_sample = ls[ls.index <= mid_date]
out_sample = ls[ls.index > mid_date]
if len(in_sample) < 36 or len(out_sample) < 36:
continue
is_mean = in_sample.mean() * 12
is_t = in_sample.mean() / (in_sample.std() / np.sqrt(len(in_sample)))
oos_mean = out_sample.mean() * 12
oos_t = out_sample.mean() / (out_sample.std() / np.sqrt(len(out_sample)))
oos_comparison.append({
'signal': name,
'is_return': is_mean,
'is_t': is_t,
'oos_return': oos_mean,
'oos_t': oos_t,
'decay_pct': (1 - oos_mean / is_mean) * 100 if is_mean != 0 else np.nan,
'same_sign': np.sign(is_mean) == np.sign(oos_mean),
'category': zoo_results[name].get('category', 'Other')
})
oos_df = pd.DataFrame(oos_comparison)
print(f"Split date: {mid_date.strftime('%Y-%m')}")
print(f"Factors with same sign IS and OOS: "
f"{oos_df['same_sign'].mean():.1%}")
print(f"Average OOS decay: {oos_df['decay_pct'].median():.1f}%")fig, axes = plt.subplots(1, 2, figsize=(14, 5))
# Panel A: Scatter
category_colors = {
'Value': '#2C5F8A', 'Size': '#1ABC9C', 'Profitability': '#27AE60',
'Investment': '#E67E22', 'Momentum': '#C0392B', 'Risk': '#8E44AD',
'Liquidity': '#3498DB', 'Quality': '#F1C40F', 'Dividend': '#95A5A6'
}
for cat in oos_df['category'].unique():
subset = oos_df[oos_df['category'] == cat]
axes[0].scatter(subset['is_return'] * 100, subset['oos_return'] * 100,
color=category_colors.get(cat, '#BDC3C7'),
s=50, alpha=0.7, label=cat, edgecolors='white')
lim = max(abs(oos_df['is_return'].max()), abs(oos_df['oos_return'].max())) * 100 + 2
axes[0].plot([-lim, lim], [-lim, lim], 'k--', linewidth=1, alpha=0.5)
axes[0].plot([-lim, lim], [0, 0], color='gray', linewidth=0.5)
axes[0].plot([0, 0], [-lim, lim], color='gray', linewidth=0.5)
axes[0].set_xlabel('In-Sample Return (% ann.)')
axes[0].set_ylabel('Out-of-Sample Return (% ann.)')
axes[0].set_title('Panel A: IS vs OOS Factor Returns')
axes[0].legend(fontsize=7, ncol=2)
# Panel B: Decay by category
cat_decay = (
oos_df.groupby('category')
.agg(
median_decay=('decay_pct', 'median'),
pct_same_sign=('same_sign', 'mean'),
n=('signal', 'count')
)
.sort_values('median_decay')
)
colors_bar = [category_colors.get(cat, '#BDC3C7') for cat in cat_decay.index]
axes[1].barh(range(len(cat_decay)), cat_decay['median_decay'],
color=colors_bar, alpha=0.85, edgecolor='white')
axes[1].set_yticks(range(len(cat_decay)))
axes[1].set_yticklabels(cat_decay.index)
axes[1].set_xlabel('Median OOS Decay (%)')
axes[1].set_title('Panel B: OOS Decay by Category')
axes[1].axvline(x=0, color='gray', linewidth=0.8)
plt.tight_layout()
plt.show()Bây giờ chúng ta kết hợp mọi bằng chứng, gồm các hiệu chỉnh kiểm định đa giả thuyết, hiệu quả ngoài mẫu và tính hợp lý về mặt kinh tế, để xác định những nhân tố xứng đáng được tin cậy trên thị trường Việt Nam.
# Merge all test results
final = zoo_df.copy()
# Add OOS results
oos_lookup = oos_df.set_index('signal')
for col in ['oos_return', 'oos_t', 'same_sign']:
final[col] = final.index.map(lambda x: oos_lookup.loc[x, col]
if x in oos_lookup.index else np.nan)
# Add bootstrap results
boot_lookup = boot_results
final['boot_pval'] = final.index.map(
lambda x: boot_lookup.loc[x, 'boot_pval']
if x in boot_lookup.index else np.nan
)
# Multi-hurdle criteria
final['pass_t2'] = final['abs_t'] > 2.0
final['pass_t3'] = final['abs_t'] > 3.0
final['pass_bh'] = final['bh_reject']
final['pass_boot'] = final['boot_pval'] < 0.05
final['pass_oos_sign'] = final['same_sign'] == True
final['pass_oos_t'] = final['oos_t'].abs() > 1.5 # Relaxed OOS threshold
# Count hurdles passed
hurdle_cols = ['pass_t2', 'pass_t3', 'pass_bh', 'pass_boot',
'pass_oos_sign', 'pass_oos_t']
final['n_hurdles'] = final[hurdle_cols].sum(axis=1)
# "Robust" = passes at least 5 of 6 hurdles
final['robust'] = final['n_hurdles'] >= 5
robust_factors = final[final['robust']].sort_values('abs_t', ascending=False)
print(f"Multi-Hurdle Filter Results:")
print(f" Pass |t| > 2.0: {final['pass_t2'].sum()}")
print(f" Pass |t| > 3.0: {final['pass_t3'].sum()}")
print(f" Pass BH (q=0.05): {final['pass_bh'].sum()}")
print(f" Pass bootstrap: {final['pass_boot'].sum()}")
print(f" Pass OOS sign: {final['pass_oos_sign'].sum()}")
print(f" Pass OOS |t| > 1.5: {final['pass_oos_t'].sum()}")
print(f"\n ROBUST (≥5/6 hurdles): {final['robust'].sum()}")
if len(robust_factors) > 0:
print(f"\nRobust Factors:")
display_cols = ['ann_return', 't_stat', 'bh_p_corrected',
'oos_return', 'oos_t', 'n_hurdles', 'category']
available_cols = [c for c in display_cols if c in robust_factors.columns]
print(robust_factors[available_cols].round(3).to_string())top25 = final.sort_values('abs_t', ascending=False).head(25)
fig, ax = plt.subplots(figsize=(10, 10))
hurdle_matrix = top25[hurdle_cols].astype(float).values
sns.heatmap(hurdle_matrix, ax=ax,
xticklabels=['|t|>2', '|t|>3', 'BH', 'Boot', 'OOS sign', 'OOS |t|'],
yticklabels=top25.index,
cmap=['#E74C3C', '#27AE60'],
cbar=False, linewidths=0.5, linecolor='white',
annot=np.where(hurdle_matrix == 1, '✓', '✗'),
fmt='')
# Add hurdle count column
for i, (_, row) in enumerate(top25.iterrows()):
ax.text(len(hurdle_cols) + 0.3, i + 0.5,
f"{int(row['n_hurdles'])}/6",
va='center', fontsize=9,
fontweight='bold' if row['n_hurdles'] >= 5 else 'normal',
color='#27AE60' if row['n_hurdles'] >= 5 else '#C0392B')
ax.set_title('Multi-Hurdle Test Results (Top 25 Anomalies)')
plt.tight_layout()
plt.show()Phân tích trong chương này đưa ra một số hàm ý thực tiễn:
Ngưỡng thông thường là không đủ. Với \(M \approx 50\) hiện tượng bất thường được kiểm định đồng thời, chỉ yêu cầu \(|t| > 2.0\) gần như chắc chắn dẫn đến các phát hiện sai. Các nhà nghiên cứu Việt Nam nên áp dụng ngưỡng tối thiểu \(|t| > 3.0\) cho từng nhân tố, nhất quán với Harvey, Liu, và Zhu (2016), và nên báo cáo các giá trị p đã điều chỉnh theo BH cho mọi nghiên cứu kiểm định nhiều tín hiệu.
Mặt cắt ngang nhỏ làm khuếch đại nhiễu. Với 600-800 cổ phiếu, các danh mục ngũ phân vị chỉ chứa 120-160 cổ phiếu mỗi danh mục. Điều này tạo ra lợi nhuận danh mục nhiễu hơn và khoảng tin cậy rộng hơn so với các nghiên cứu ở Mỹ với hơn 4.000 cổ phiếu. Những chiến lược có vẻ có ý nghĩa ở Mỹ với \(|t| = 2.5\) có thể chỉ có \(|t| < 1.5\) ở Việt Nam đơn giản vì mẫu nhỏ hơn. Đây không phải là bằng chứng cho thấy nhân tố không tồn tại ở Việt Nam, mà là bằng chứng cho thấy dữ liệu Việt Nam thiếu sức mạnh để phát hiện ra nó.
Kiểm định ngoài mẫu là thiết yếu. Cách chia theo chuỗi thời gian cho thấy sự suy giảm đáng kể ở nhiều hiện tượng bất thường, nhất quán với McLean và Pontiff (2016). Những nhân tố vượt qua cả trong mẫu lẫn ngoài mẫu với dấu và độ lớn nhất quán là hiếm và có giá trị.
Nhóm tín hiệu đóng vai trò quan trọng. Các tín hiệu động lượng và khả năng sinh lời có xu hướng có tỷ lệ tái lập cao nhất giữa các thị trường (Fama và French 2012; Jacobs và Müller 2020). Các tín hiệu giá trị và đầu tư mang tính đặc thù theo quốc gia nhiều hơn (Griffin 2002). Các tín hiệu thanh khoản có khả năng mạnh nhất ở các thị trường mới nổi như Việt Nam, nơi các ma sát giao dịch lớn nhất.
Đường đi của lãi suất là một phần của mẫu. Giai đoạn mẫu của Việt Nam bị chi phối bởi một biến động một chiều của lãi suất: lãi suất chính sách đạt đỉnh gần 14% năm 2011 và giảm xuống khoảng 4,5% vào năm 2023, và lợi suất trái phiếu chính phủ cũng giảm theo. Binsbergen, Ma, và Schwert (2026) chỉ ra rằng ở Mỹ, một sự sụt giảm dài hạn tương tự đủ để phân loại lại một tỷ lệ đáng kể các hiện tượng bất thường khi lợi nhuận vượt mức được đo so với một danh mục trái phiếu có kỳ hạn bình quân tương ứng thay vì lãi suất ngắn hạn. Vì sự sụt giảm ở Việt Nam diễn ra trong một khoảng thời gian ngắn hơn nhiều, bất kỳ danh mục mua–bán khống nào có các vị thế khác nhau về kỳ hạn (giá trị so với tăng trưởng, tỷ suất cổ tức cao so với thấp, đòn bẩy cao so với thấp) có thể đang nắm bắt sự định giá lại đó thay vì một phần bù rủi ro mặt cắt ngang. Tối thiểu, các nhà nghiên cứu nên báo cáo tỷ suất cổ tức hoặc chênh lệch kỳ hạn giữa các vị thế mua và bán, và thận trọng với những nhân tố có phần bù tập trung trong giai đoạn lãi suất sụt giảm 2011-2015.
Báo cáo đầy đủ bộ chỉ tiêu. Mọi nghiên cứu đưa ra một nhân tố mới cho thị trường Việt Nam nên báo cáo: (i) thống kê t thô, (ii) giá trị p đã điều chỉnh theo BH với số kiểm định trong nghiên cứu, (iii) hiệu quả ngoài mẫu trong một giai đoạn được giữ lại, và (iv) độ nhạy với các lựa chọn xây dựng (từ chương trước). Chuẩn báo cáo này sẽ cải thiện đáng kể độ tin cậy của nghiên cứu nhân tố ở Việt Nam.
| Phương pháp | Kiểm soát | Ngưỡng (xấp xỉ) | Số còn lại | Triết lý |
|---|---|---|---|---|
| Không điều chỉnh | Kiểm định đơn | |t| > 2.0 | Nhiều | Không hiệu chỉnh |
| Bonferroni | FWER | |t| > 3.3–4.0 | Ít | Không có dương tính giả |
| Holm | FWER | Theo từng bước | Ít | Ít thận trọng hơn một chút |
| BH | FDR ở mức q | Thích ứng | Trung bình | Chấp nhận một số dương tính giả |
| HLZ | Toàn giới nghiên cứu | |t| > 3.0–3.78 | Trung bình | Phân phối tiên nghiệm trên mọi nhân tố đã kiểm định |
| Bootstrap White/RW | FWER (theo dữ liệu) | Giá trị tới hạn bootstrap | Ít | Tính đến phân phối không chuẩn |
| Storey FDP | Ước lượng FDP | Liên tục | Chẩn đoán | Ước lượng tỷ lệ giả thuyết không thực sự đúng |
| Đa ngưỡng | Kết hợp | Nhiều kiểm định | Vững nhất | Đòi hỏi bằng chứng hội tụ |
Vườn thú nhân tố không phải là hiện tượng riêng của nước Mỹ. Cám dỗ kiểm định nhiều tín hiệu và chỉ báo cáo những tín hiệu cho kết quả tốt nhất tồn tại ở mọi thị trường. Ở Việt Nam, nơi dữ liệu ngắn hơn, nhiễu hơn và cộng đồng học thuật đang phát triển nhanh, rủi ro phát hiện sai là cao. Các công cụ được phát triển trong chương này, gồm điều chỉnh BH, kiểm tra thực tế bằng bootstrap, ước lượng \(\pi_0\) của Storey, cách chia ngoài mẫu và bộ lọc đa ngưỡng, cung cấp một khung nguyên tắc để phân tách các nhân tố thực sự khỏi nhiễu thống kê.
Kết luận trung thực có lẽ sẽ khiêm tốn: trong số hàng chục hiện tượng bất thường được ghi nhận ở các thị trường phát triển, chỉ một số ít vượt qua kiểm định đa giả thuyết nghiêm ngặt trên dữ liệu Việt Nam. Những nhân tố sống sót đó (có thể tập trung ở động lượng, khả năng sinh lời và thanh khoản) tạo nền tảng cho nghiên cứu định giá tài sản đáng tin cậy ở thị trường này.
# Vườn thú nhân tố và kiểm định đa giả thuyết
::: callout-note
Trong chương này, chúng ta đối mặt với vấn đề kiểm định đa giả thuyết trong định giá tài sản thực nghiệm. Chúng ta tái lập một tập hợp rộng các hiện tượng bất thường (anomaly) đã được công bố trên thị trường Việt Nam, áp dụng các hiệu chỉnh chính thức theo số lượng kiểm định đã thực hiện, ước lượng tỷ lệ các phát hiện sai, và xây dựng một khung phân tích để quyết định nhân tố nào xứng đáng được tin cậy.
:::
Đến năm 2016, tài chính học thuật đã ghi nhận hơn 300 biến được cho là dự báo được mặt cắt ngang của lợi nhuận cổ phiếu. @cochrane2011presidential đã gọi sự sinh sôi này bằng một cụm từ đáng nhớ là "vườn thú nhân tố" (zoo of factors). @harvey2016and thống kê hơn 300 nhân tố đã công bố và lập luận rằng ngưỡng ý nghĩa thống kê thông thường $t > 2.0$ là quá thấp so với số lượng giả thuyết đã được kiểm định, dù là công khai hay ngầm định, qua nhiều thập kỷ nghiên cứu. Ngưỡng $t > 3.0$ (và trong một số trường hợp $t > 3.78$) mà họ đề xuất sẽ loại bỏ phần lớn các hiện tượng bất thường đã công bố.
Vấn đề này không chỉ mang tính học thuật. Mỗi nhà nghiên cứu sắp xếp cổ phiếu theo một biến mới và thấy $t > 2$ đều đang ngầm kiểm định một giả thuyết đã được kiểm định hàng trăm lần với các biến khác. Xác suất để *ít nhất một* trong 300 kiểm định độc lập cho ra $t > 2$ do ngẫu nhiên (ngay cả khi không tồn tại hiệu ứng thực) về cơ bản bằng 1,0. Đây chính là vấn đề kiểm định đa giả thuyết, và việc không tính đến nó có nghĩa là một phần đáng kể của "vườn thú nhân tố" là các phát hiện sai.
Đối với nghiên cứu cổ phiếu Việt Nam, vấn đề có lẽ còn nghiêm trọng hơn. Mặt cắt ngang nhỏ hơn (600-800 cổ phiếu so với hơn 4.000 ở Mỹ), chuỗi thời gian ngắn hơn (từ 2008 so với từ 1963), và dữ liệu nhiễu hơn (kém thanh khoản, các ngày không giao dịch, những bất thường về kế toán). Tất cả các yếu tố này làm tăng rủi ro của cả sai lầm loại I (tìm thấy điều không tồn tại) lẫn sai lầm loại II (bỏ sót điều có thật). Chương này cung cấp các công cụ thống kê để xử lý bức tranh đó một cách chặt chẽ.
## Quy mô của vấn đề {#sec-factor-zoo-scale}
### Có bao nhiêu nhân tố?
@harvey2016and đếm được hơn 300 nhân tố được công bố trên các tạp chí hàng đầu đến năm 2012. @chen2022open xây dựng và công bố công khai hơn 300 hiện tượng bất thường "có ý nghĩa rõ ràng" từ các tài liệu của Mỹ. @hou2020replicating cố gắng tái lập 452 hiện tượng bất thường và thấy rằng 64% (với lợi nhuận theo trọng số vốn hóa, VW) không đạt $t > 1.96$ trong mẫu của họ. @jensen2023there có cái nhìn lạc quan hơn, khi thấy rằng phần lớn các nhân tố được tái lập về hướng, dù không phải lúc nào cũng về độ lớn.
Điểm mấu chốt là số nhân tố được *kiểm định* vượt xa số được *công bố*. Với mỗi nhân tố được công bố có $t > 2$, có thể có hàng chục kiểm định không được công bố với $t < 2$ chưa bao giờ được in ra, đó là vấn đề "ngăn kéo hồ sơ" (file drawer) kinh điển. Tổng số kiểm định mà giới nghiên cứu đã thực hiện là không thể biết, nhưng chắc chắn lên tới hàng nghìn.
Việc đếm các *tín hiệu* ($|t| > 1.96$ theo cả hai định nghĩa) cũng đánh giá thấp vấn đề, vì mỗi tín hiệu có thể được đo lường theo nhiều cách hợp lý. @binsbergen2026factor nhấn mạnh điểm này rất rõ. Họ giữ nguyên các biến dùng để sắp xếp nhưng thay đổi một quy ước tưởng chừng vô hại: thay vì trừ đi lãi suất tín phiếu kho bạc kỳ hạn một tháng, họ trừ đi lợi nhuận của một danh mục trái phiếu chính phủ zero-coupon có kỳ hạn bình quân tương ứng, để một nhân tố được đánh giá so với một kịch bản dòng tiền cố định cùng kỳ hạn thay vì so với tiền mặt. Chạy lại 153 hiện tượng bất thường đã công bố trong các giai đoạn mẫu gốc, họ thấy 66 hiện tượng vẫn vững ($|t| > 1.96$ theo cả hai định nghĩa) trong khi 20 hiện tượng đổi trạng thái, với khoảng 0,30 phát hiện sai trên mỗi phát hiện đúng. Mở rộng bài tập sang 1.395 cách sắp xếp giả định xây dựng từ các tỷ số Compustat, tỷ lệ giữa phát hiện sai và phát hiện đúng tăng lên khoảng 1,0 và ổn định qua các cửa sổ trượt. Ngay cả nhân tố thị trường cũng bị ảnh hưởng: nhân tố tương ứng đã điều chỉnh theo kỳ hạn có tương quan dưới 0,5 với lợi nhuận vượt mức thị trường thông thường, và phần bù cổ phiếu không còn khác không một cách có ý nghĩa thống kê. Cơ chế là sự sụt giảm dài hạn của lãi suất trong bốn thập kỷ, vốn mang lại một khoản lợi bất ngờ mang tính cơ học cho các vị thế kỳ hạn dài của danh mục mua–bán khống. Kết quả của họ không khẳng định định nghĩa truyền thống là sai, mà cho thấy quá trình phát hiện nhạy cảm với một lựa chọn mà hầu hết các bài báo không bao giờ báo cáo, đó chính là "đa vũ trụ" (multiverse) các bậc tự do của nhà nghiên cứu mà các hiệu chỉnh kiểm định đa giả thuyết không nắm bắt được.
```{python}
#| label: setup
#| code-summary: "Import libraries and configure environment"
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import statsmodels.api as sm
from scipy import stats
from statsmodels.stats.multitest import multipletests
from itertools import product
import warnings
warnings.filterwarnings('ignore')
plt.rcParams.update({
'figure.figsize': (12, 6),
'figure.dpi': 150,
'font.size': 11,
'axes.spines.top': False,
'axes.spines.right': False
})
```
```{python}
#| label: data-load
#| eval: false
#| code-summary: "Load returns and characteristics panel"
from datacore import DataCoreClient
client = DataCoreClient()
# Monthly returns (survivorship-bias-free)
monthly = client.get_monthly_returns(
exchanges=['HOSE', 'HNX'],
start_date='2008-01-01',
end_date='2024-12-31',
include_delisted=True,
fields=[
'ticker', 'month_end', 'monthly_return', 'market_cap',
'exchange', 'volume_avg_20d', 'turnover_value_avg_20d',
'n_zero_volume_days'
]
)
# Pre-computed characteristic panel (from factor construction chapter)
characteristics = client.get_characteristics_panel(
exchanges=['HOSE', 'HNX'],
start_date='2008-01-01',
end_date='2024-12-31',
include_delisted=True
)
# Factor returns for spanning tests
factors_ff = client.get_factor_returns(
market='vietnam',
start_date='2008-01-01',
end_date='2024-12-31',
factors=['mkt_excess', 'smb', 'hml', 'rmw', 'cma', 'wml']
)
monthly['month_end'] = pd.to_datetime(monthly['month_end'])
print(f"Monthly returns: {len(monthly):,}")
print(f"Characteristics: {characteristics.shape}")
print(f"Available signals: {[c for c in characteristics.columns if c not in ['ticker', 'month_end']]}")
```
## Xây dựng vườn thú các hiện tượng bất thường {#sec-factor-zoo-zoo}
### Định nghĩa tín hiệu
Chúng ta xây dựng một tập hợp toàn diện các biến bất thường bao quát các nhóm chính trong tài liệu nghiên cứu. Mỗi tín hiệu được lấy độ trễ phù hợp để tránh thiên lệch nhìn trước (các tín hiệu kế toán dùng cách căn chỉnh theo thời điểm công bố thực tế; các tín hiệu giá dùng giá trị của tháng trước).
```{python}
#| label: signal-construction
#| eval: false
#| code-summary: "Construct the full set of anomaly signals"
# Merge returns with characteristics
panel = monthly.merge(characteristics, on=['ticker', 'month_end'], how='left')
panel = panel.sort_values(['ticker', 'month_end'])
# ============================================================
# Category 1: VALUE (8 signals)
# ============================================================
# Already in characteristics panel from PIT merge:
# bm, ep, cfp, sp, dp, ev_ebitda
# Add earnings yield variants
value_signals = ['bm', 'ep', 'cfp', 'sp', 'dp']
# ============================================================
# Category 2: SIZE (3 signals)
# ============================================================
panel['log_mcap'] = np.log(panel['market_cap'].clip(lower=1))
panel['log_assets'] = np.log(panel['total_assets'].clip(lower=1))
panel['log_revenue'] = np.log(panel['revenue'].clip(lower=1))
size_signals = ['log_mcap', 'log_assets', 'log_revenue']
# ============================================================
# Category 3: PROFITABILITY (8 signals)
# ============================================================
# GP/A (Novy-Marx 2013), OP (FF 2015), ROE, ROA, etc.
prof_signals = ['gp_at', 'op', 'roe', 'roa', 'profit_margin',
'asset_turnover', 'gross_margin', 'oper_margin']
# ============================================================
# Category 4: INVESTMENT / GROWTH (6 signals)
# ============================================================
inv_signals = ['asset_growth', 'capex_at', 'investment',
'sales_growth', 'equity_issuance', 'debt_issuance']
# ============================================================
# Category 5: MOMENTUM / REVERSAL (6 signals)
# ============================================================
# Compute from returns
for lag_start, lag_end, name in [
(2, 12, 'ret_12_2'), (2, 6, 'ret_6_2'), (7, 12, 'ret_12_7'),
(1, 1, 'ret_1'), (1, 3, 'ret_3_1')
]:
if name not in panel.columns:
panel[name] = (
panel.groupby('ticker')['monthly_return']
.transform(lambda x: x.shift(lag_start).rolling(lag_end - lag_start + 1)
.apply(lambda r: (1 + r).prod() - 1, raw=True))
)
# Earnings momentum (SUE) should already be in characteristics
mom_signals = ['ret_12_2', 'ret_6_2', 'ret_12_7', 'ret_1', 'ret_3_1', 'sue']
# ============================================================
# Category 6: RISK (5 signals)
# ============================================================
risk_signals = ['beta', 'ivol', 'tvol', 'max_ret', 'skewness']
# ============================================================
# Category 7: LIQUIDITY / TRADING (6 signals)
# ============================================================
liq_signals = ['amihud', 'zero_return_pct', 'log_turnover',
'roll_spread', 'cs_spread', 'volume_cv']
# ============================================================
# Category 8: QUALITY / FINANCIAL HEALTH (5 signals)
# ============================================================
qual_signals = ['leverage', 'current_ratio', 'accruals',
'earnings_variability', 'piotroski_f']
# ============================================================
# Category 9: DIVIDEND / PAYOUT (3 signals)
# ============================================================
div_signals = ['div_yield', 'payout_ratio', 'net_repurchase']
# Combine all
all_signals = (value_signals + size_signals + prof_signals +
inv_signals + mom_signals + risk_signals +
liq_signals + qual_signals + div_signals)
# Remove signals with too little coverage
signal_coverage = {}
for sig in all_signals:
if sig in panel.columns:
cov = panel[sig].notna().mean()
signal_coverage[sig] = cov
signal_coverage = pd.Series(signal_coverage).sort_values(ascending=False)
valid_signals = signal_coverage[signal_coverage > 0.30].index.tolist()
print(f"Total signals defined: {len(all_signals)}")
print(f"Signals with >30% coverage: {len(valid_signals)}")
print(f"\nSignal categories:")
for cat_name, sigs in [
('Value', value_signals), ('Size', size_signals),
('Profitability', prof_signals), ('Investment', inv_signals),
('Momentum', mom_signals), ('Risk', risk_signals),
('Liquidity', liq_signals), ('Quality', qual_signals),
('Dividend', div_signals)
]:
n_valid = sum(1 for s in sigs if s in valid_signals)
print(f" {cat_name:<15}: {n_valid}/{len(sigs)}")
```
### Xây dựng nhân tố hàng loạt
Chúng ta chạy bộ máy xây dựng nhân tố (từ chương trước) trên tất cả các tín hiệu hợp lệ để tạo ra một "vườn thú" lợi nhuận nhân tố:
```{python}
#| label: mass-construction
#| eval: false
#| code-summary: "Construct long-short portfolios for every valid signal"
# Import the factor engine from previous chapter
# (In practice, this would be imported from a shared module)
def construct_factor_simple(panel_df, signal_col, long_high=True,
n_groups=5, weighting='value',
min_stocks=20):
"""
Simplified factor construction: quintile sort on signal,
long-short = Q5 - Q1 (or reversed if long_high=False).
Returns a Series of monthly long-short returns.
"""
df = panel_df[['ticker', 'month_end', 'monthly_return',
'market_cap', signal_col]].dropna().copy()
df = df.sort_values(['ticker', 'month_end'])
# Lag the signal by one month
df['signal_lag'] = df.groupby('ticker')[signal_col].shift(1)
df = df.dropna(subset=['signal_lag', 'monthly_return'])
results = []
for month, group in df.groupby('month_end'):
if len(group) < min_stocks:
continue
group['quintile'] = pd.qcut(
group['signal_lag'].rank(method='first'),
n_groups, labels=False, duplicates='drop'
)
if weighting == 'value':
def wret(g):
if g['market_cap'].sum() > 0:
return np.average(g['monthly_return'],
weights=g['market_cap'])
return g['monthly_return'].mean()
port_ret = group.groupby('quintile').apply(wret)
else:
port_ret = group.groupby('quintile')['monthly_return'].mean()
if 0 in port_ret.index and (n_groups - 1) in port_ret.index:
if long_high:
ls = port_ret[n_groups - 1] - port_ret[0]
else:
ls = port_ret[0] - port_ret[n_groups - 1]
results.append({'month_end': month, 'ls_return': ls})
if not results:
return None
return pd.DataFrame(results).set_index('month_end')['ls_return']
# Determine expected sign for each signal
# (positive = high signal predicts high returns)
signal_directions = {
# Value: high = higher returns
'bm': True, 'ep': True, 'cfp': True, 'sp': True, 'dp': True,
# Size: small = higher returns (low signal = high returns)
'log_mcap': False, 'log_assets': False, 'log_revenue': False,
# Profitability: high = higher returns
'gp_at': True, 'op': True, 'roe': True, 'roa': True,
'profit_margin': True, 'asset_turnover': True,
'gross_margin': True, 'oper_margin': True,
# Investment: low growth = higher returns
'asset_growth': False, 'capex_at': False, 'investment': False,
'sales_growth': False, 'equity_issuance': False,
'debt_issuance': False,
# Momentum: high = higher returns (except reversal)
'ret_12_2': True, 'ret_6_2': True, 'ret_12_7': True,
'ret_1': False, 'ret_3_1': False, 'sue': True,
# Risk: low risk = higher returns (anomaly direction)
'beta': False, 'ivol': False, 'tvol': False,
'max_ret': False, 'skewness': False,
# Liquidity: illiquid = higher returns
'amihud': True, 'zero_return_pct': True,
'log_turnover': False, 'roll_spread': True,
'cs_spread': True, 'volume_cv': True,
# Quality: high quality = higher returns
'leverage': False, 'current_ratio': True, 'accruals': False,
'earnings_variability': False, 'piotroski_f': True,
# Dividend: high = higher returns
'div_yield': True, 'payout_ratio': True, 'net_repurchase': False,
}
# Construct all factors
zoo_results = {}
for sig in valid_signals:
long_high = signal_directions.get(sig, True)
ls = construct_factor_simple(
panel, sig, long_high=long_high,
n_groups=5, weighting='value'
)
if ls is not None and len(ls) >= 60:
mean_ret = ls.mean()
se = ls.std() / np.sqrt(len(ls))
t = mean_ret / se if se > 0 else 0
zoo_results[sig] = {
'mean_monthly': mean_ret,
'ann_return': mean_ret * 12,
'ann_vol': ls.std() * np.sqrt(12),
'sharpe': mean_ret / ls.std() * np.sqrt(12) if ls.std() > 0 else 0,
't_stat': t,
'abs_t': abs(t),
'p_value': 2 * (1 - stats.t.cdf(abs(t), df=len(ls) - 1)),
'n_months': len(ls),
'direction': 'Long High' if long_high else 'Long Low',
'category': next(
(cat for cat, sigs in [
('Value', value_signals), ('Size', size_signals),
('Profitability', prof_signals),
('Investment', inv_signals),
('Momentum', mom_signals), ('Risk', risk_signals),
('Liquidity', liq_signals),
('Quality', qual_signals),
('Dividend', div_signals)
] if sig in sigs), 'Other'
),
'returns': ls # Store for later use
}
zoo_df = pd.DataFrame({k: {kk: vv for kk, vv in v.items() if kk != 'returns'}
for k, v in zoo_results.items()}).T
zoo_df = zoo_df.sort_values('abs_t', ascending=False)
print(f"Factors successfully constructed: {len(zoo_df)}")
print(f"Factors with |t| > 2.0: {(zoo_df['abs_t'] > 2.0).sum()} "
f"({(zoo_df['abs_t'] > 2.0).mean():.1%})")
print(f"Factors with |t| > 3.0: {(zoo_df['abs_t'] > 3.0).sum()} "
f"({(zoo_df['abs_t'] > 3.0).mean():.1%})")
```
```{python}
#| label: fig-t-stat-distribution
#| eval: false
#| fig-cap: "Phân phối của thống kê t trên tất cả các hiện tượng bất thường được kiểm định trên thị trường Việt Nam. Panel A hiển thị biểu đồ tần suất cùng các ngưỡng thông thường (|t| > 1.96) và ngưỡng Harvey-Liu-Zhu (|t| > 3.0). Panel B xếp hạng các hiện tượng bất thường theo giá trị tuyệt đối của thống kê t, với màu sắc theo nhóm. Một phần đáng kể các hiện tượng bất thường có |t| > 2 có thể là phát hiện sai do số lượng kiểm định đồng thời."
#| code-summary: "Visualize the distribution of t-statistics across all anomalies"
fig, axes = plt.subplots(1, 2, figsize=(16, 6))
# Panel A: Histogram
t_vals = zoo_df['t_stat'].values
axes[0].hist(t_vals, bins=30, color='#2C5F8A', alpha=0.7,
edgecolor='white', density=True)
# Overlay normal distribution under null
x_range = np.linspace(-5, 5, 200)
axes[0].plot(x_range, stats.norm.pdf(x_range), 'k--', linewidth=1.5,
label='N(0,1) null')
# Threshold lines
axes[0].axvline(x=1.96, color='#E67E22', linestyle='--', linewidth=1.5,
label='|t| = 1.96')
axes[0].axvline(x=-1.96, color='#E67E22', linestyle='--', linewidth=1.5)
axes[0].axvline(x=3.0, color='#C0392B', linestyle='--', linewidth=1.5,
label='|t| = 3.0 (HLZ)')
axes[0].axvline(x=-3.0, color='#C0392B', linestyle='--', linewidth=1.5)
axes[0].set_xlabel('t-statistic')
axes[0].set_ylabel('Density')
axes[0].set_title('Panel A: t-Statistic Distribution')
axes[0].legend(fontsize=9)
# Panel B: Ranked bar chart
top_n = min(40, len(zoo_df))
top = zoo_df.head(top_n).copy()
category_colors = {
'Value': '#2C5F8A', 'Size': '#1ABC9C', 'Profitability': '#27AE60',
'Investment': '#E67E22', 'Momentum': '#C0392B', 'Risk': '#8E44AD',
'Liquidity': '#3498DB', 'Quality': '#F1C40F', 'Dividend': '#95A5A6'
}
bar_colors = [category_colors.get(top.iloc[i]['category'], '#BDC3C7')
for i in range(len(top))]
axes[1].barh(range(top_n), top['abs_t'].values,
color=bar_colors, alpha=0.85, edgecolor='white')
axes[1].set_yticks(range(top_n))
axes[1].set_yticklabels(top.index, fontsize=7)
axes[1].axvline(x=1.96, color='#E67E22', linestyle='--', linewidth=1)
axes[1].axvline(x=3.0, color='#C0392B', linestyle='--', linewidth=1)
axes[1].set_xlabel('|t-statistic|')
axes[1].set_title('Panel B: Anomalies Ranked by |t|')
axes[1].invert_yaxis()
# Legend for categories
from matplotlib.patches import Patch
legend_elements = [Patch(facecolor=c, label=cat)
for cat, c in category_colors.items()]
axes[1].legend(handles=legend_elements, fontsize=7,
loc='lower right', ncol=2)
plt.tight_layout()
plt.show()
```
## Vấn đề kiểm định đa giả thuyết {#sec-factor-zoo-multiple-testing}
### Vì sao giá trị p của một kiểm định đơn lẻ gây hiểu lầm
Giả sử bạn kiểm định $M$ giả thuyết độc lập, mỗi giả thuyết ở mức ý nghĩa $\alpha = 0.05$. Dưới giả thuyết không toàn cục (không có hiệu ứng thực), số lần bác bỏ sai kỳ vọng là $M \times \alpha$. Với $M = 50$ hiện tượng bất thường, bạn kỳ vọng có 2,5 kết quả "có ý nghĩa" hoàn toàn do ngẫu nhiên.
Xác suất có *ít nhất một* lần bác bỏ sai là:
$$
P(\text{at least one false discovery}) = 1 - (1 - \alpha)^M
$$ {#eq-family-error}
Với $M = 50$ và $\alpha = 0.05$, giá trị này là $1 - 0.95^{50} = 0.923$. Với $M = 300$, nó về cơ bản bằng 1. Giá trị p của một kiểm định đơn lẻ vô dụng trong việc quyết định nhân tố nào trong số nhiều nhân tố được kiểm định là thật.
### Hai khái niệm về tỷ lệ sai sót
Các hiệu chỉnh kiểm định đa giả thuyết kiểm soát một trong hai tỷ lệ sai sót:
**Tỷ lệ sai sót theo họ giả thuyết (FWER).** Xác suất mắc *ít nhất một* sai lầm loại I trên toàn bộ các kiểm định. Đây là tiêu chí thận trọng nhất. Nếu FWER = 0,05, bạn có thể tin chắc 95% rằng *mọi* giả thuyết không đã bị bác bỏ đều là phát hiện đúng. Các phương pháp: hiệu chỉnh @bonferroni1936teoria, Holm step-down, @romano2005exact.
**Tỷ lệ phát hiện sai (FDR).** *Tỷ lệ* kỳ vọng các lần bác bỏ sai trong số các giả thuyết bị bác bỏ. Tiêu chí này ít thận trọng hơn. Nếu FDR = 0,05, bạn kỳ vọng 5% số phát hiện của mình là sai, còn 95% còn lại là thật. Các phương pháp: @benjamini1995controlling (BH), @storey2003positive.
Trong định giá tài sản, kiểm soát FDR nhìn chung phù hợp hơn FWER vì chúng ta không cố đảm bảo không có phát hiện sai nào, mà cố xác định một tập các nhân tố trong đó *phần lớn* là thật. Một nhà nghiên cứu kiểm soát FDR ở mức 5% và phát hiện 10 nhân tố sẽ kỳ vọng khoảng 0,5 nhân tố trong số đó là sai, điều này chấp nhận được khi xây dựng một mô hình nhân tố.
```{python}
#| label: fig-error-rates
#| eval: false
#| fig-cap: "Minh họa các khái niệm FWER và FDR. Khi số kiểm định tăng lên, xác suất có ít nhất một phát hiện sai (FWER) tiến nhanh về 1 trừ khi ngưỡng cho mỗi kiểm định được điều chỉnh. FDR, tức tỷ lệ kỳ vọng của các phát hiện sai trong số các giả thuyết bị bác bỏ, là một tiêu chí ít nghiêm ngặt hơn nhưng thực tế hơn để lựa chọn nhân tố."
#| code-summary: "Visualize how error rates scale with number of tests"
M_range = np.arange(1, 301)
alpha = 0.05
fwer_unadj = 1 - (1 - alpha) ** M_range
fwer_bonf = np.minimum(1, M_range * alpha) # Not actual FWER, just threshold
# Expected false discoveries at alpha = 0.05
expected_false = M_range * alpha
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
axes[0].plot(M_range, fwer_unadj, color='#C0392B', linewidth=2,
label='Unadjusted FWER')
axes[0].axhline(y=0.05, color='gray', linestyle='--', linewidth=1,
label='α = 0.05 target')
axes[0].set_xlabel('Number of Tests (M)')
axes[0].set_ylabel('P(≥1 False Discovery)')
axes[0].set_title('Panel A: Family-Wise Error Rate')
axes[0].legend()
axes[1].plot(M_range, expected_false, color='#2C5F8A', linewidth=2)
axes[1].set_xlabel('Number of Tests (M)')
axes[1].set_ylabel('Expected False Discoveries')
axes[1].set_title('Panel B: Expected False Discoveries at α = 0.05')
axes[1].axhline(y=1, color='gray', linestyle='--', linewidth=0.8)
axes[1].text(250, 1.5, '1 false discovery', color='gray', fontsize=9)
plt.tight_layout()
plt.show()
```
## Các phương pháp hiệu chỉnh {#sec-factor-zoo-corrections}
### Hiệu chỉnh Bonferroni (kiểm soát FWER)
Hiệu chỉnh đơn giản nhất và thận trọng nhất: chỉ bác bỏ giả thuyết $i$ nếu giá trị p của nó thấp hơn $\alpha / M$, trong đó $M$ là tổng số kiểm định:
$$
\text{Reject } H_{0,i} \text{ if } p_i < \frac{\alpha}{M}
$$ {#eq-bonferroni}
Phương pháp này kiểm soát FWER ở mức $\alpha$ bất kể cấu trúc phụ thuộc giữa các kiểm định. Cái giá phải trả là mất rất nhiều sức mạnh kiểm định: với $M = 50$ kiểm định ở $\alpha = 0.05$, ngưỡng cho mỗi kiểm định trở thành $0.001$, đòi hỏi $|t| > 3.29$.
```{python}
#| label: bonferroni
#| eval: false
#| code-summary: "Apply Bonferroni correction to the anomaly zoo"
M = len(zoo_df)
alpha = 0.05
# Bonferroni threshold
bonf_threshold = alpha / M
bonf_t_threshold = stats.norm.ppf(1 - bonf_threshold / 2)
zoo_df['bonf_reject'] = zoo_df['p_value'] < bonf_threshold
print(f"Number of tests: {M}")
print(f"Bonferroni p-value threshold: {bonf_threshold:.6f}")
print(f"Equivalent t-statistic threshold: {bonf_t_threshold:.2f}")
print(f"Rejections: {zoo_df['bonf_reject'].sum()} / {M}")
print(f"\nSurviving anomalies:")
survivors = zoo_df[zoo_df['bonf_reject']]
if len(survivors) > 0:
print(survivors[['ann_return', 't_stat', 'category']].to_string())
else:
print(" None survive Bonferroni correction")
```
### Thủ tục Holm step-down (kiểm soát FWER)
Thủ tục Holm luôn mạnh hơn Bonferroni mà vẫn kiểm soát FWER. Nó xếp hạng các giá trị p từ nhỏ nhất đến lớn nhất và áp dụng các ngưỡng nới lỏng dần:
1. Sắp xếp các giá trị p: $p_{(1)} \leq p_{(2)} \leq \ldots \leq p_{(M)}$
2. Với $i = 1, 2, \ldots$, bác bỏ $H_{0,(i)}$ nếu $p_{(i)} < \alpha / (M - i + 1)$
3. Dừng tại $i$ đầu tiên mà $H_{0,(i)}$ không bị bác bỏ
```{python}
#| label: holm
#| eval: false
#| code-summary: "Apply the Holm step-down procedure"
p_values = zoo_df['p_value'].values
signal_names = zoo_df.index.values
# Using statsmodels implementation
reject_holm, pvals_corrected_holm, _, _ = multipletests(
p_values, alpha=0.05, method='holm'
)
zoo_df['holm_reject'] = reject_holm
zoo_df['holm_p_corrected'] = pvals_corrected_holm
print(f"Holm rejections: {reject_holm.sum()} / {M}")
holm_survivors = zoo_df[zoo_df['holm_reject']]
if len(holm_survivors) > 0:
print("\nSurviving anomalies (Holm):")
print(holm_survivors[['ann_return', 't_stat', 'holm_p_corrected',
'category']].to_string())
```
### Thủ tục Benjamini-Hochberg (kiểm soát FDR)
Thủ tục @benjamini1995controlling (BH) kiểm soát tỷ lệ phát hiện sai (tức tỷ lệ kỳ vọng các lần bác bỏ sai trong số các giả thuyết bị bác bỏ) ở mức $q$:
1. Sắp xếp các giá trị p: $p_{(1)} \leq p_{(2)} \leq \ldots \leq p_{(M)}$
2. Tìm $k$ lớn nhất sao cho $p_{(k)} \leq \frac{k}{M} q$
3. Bác bỏ tất cả $H_{0,(i)}$ với $i = 1, \ldots, k$
Thủ tục BH hợp lệ khi các kiểm định độc lập hoặc phụ thuộc dương (PRDS), một điều kiện nhìn chung được thỏa mãn đối với các kiểm định nhân tố tài chính, vì các nhân tố có xu hướng tương quan dương dưới giả thuyết không.
```{python}
#| label: bh-fdr
#| eval: false
#| code-summary: "Apply the Benjamini-Hochberg FDR procedure"
# BH procedure at q = 0.05 (expect 5% false discoveries)
reject_bh, pvals_corrected_bh, _, _ = multipletests(
p_values, alpha=0.05, method='fdr_bh'
)
zoo_df['bh_reject'] = reject_bh
zoo_df['bh_p_corrected'] = pvals_corrected_bh
# Also at q = 0.10 (more liberal)
reject_bh10, _, _, _ = multipletests(
p_values, alpha=0.10, method='fdr_bh'
)
zoo_df['bh10_reject'] = reject_bh10
print(f"BH rejections (q=0.05): {reject_bh.sum()} / {M}")
print(f"BH rejections (q=0.10): {reject_bh10.sum()} / {M}")
bh_survivors = zoo_df[zoo_df['bh_reject']].sort_values('abs_t', ascending=False)
if len(bh_survivors) > 0:
print(f"\nSurviving anomalies (BH, q=0.05):")
print(bh_survivors[['ann_return', 't_stat', 'bh_p_corrected',
'category']].head(20).to_string())
```
```{python}
#| label: fig-bh-procedure
#| eval: false
#| fig-cap: "Biểu diễn trực quan thủ tục Benjamini-Hochberg. Các giá trị p đã sắp xếp được so với đường ngưỡng BH (độ dốc = q/M). Tất cả các giá trị p nằm dưới đường này, cho đến điểm cắt lớn nhất, đều bị bác bỏ. Thủ tục này ít thận trọng hơn Bonferroni vì các ngưỡng về sau nới lỏng dần."
#| code-summary: "Plot the BH procedure graphically"
sorted_pvals = np.sort(p_values)
k = np.arange(1, M + 1)
bh_line_05 = k / M * 0.05
bh_line_10 = k / M * 0.10
fig, ax = plt.subplots(figsize=(12, 6))
ax.scatter(k, sorted_pvals, s=15, color='#2C5F8A', alpha=0.7,
label='Sorted p-values', zorder=3)
ax.plot(k, bh_line_05, color='#C0392B', linewidth=2,
label='BH line (q=0.05)')
ax.plot(k, bh_line_10, color='#E67E22', linewidth=2,
linestyle='--', label='BH line (q=0.10)')
# Highlight rejections
n_reject_05 = reject_bh.sum()
if n_reject_05 > 0:
ax.scatter(k[:n_reject_05], sorted_pvals[:n_reject_05],
s=40, color='#C0392B', zorder=4, label=f'Rejected (q=0.05)')
ax.set_xlabel('Rank (k)')
ax.set_ylabel('p-value')
ax.set_title('Benjamini-Hochberg Procedure')
ax.legend()
ax.set_ylim([-0.01, 0.15])
ax.set_xlim([0, min(M, 60)])
plt.tight_layout()
plt.show()
```
### Ngưỡng Harvey-Liu-Zhu
@harvey2016and tiếp cận theo cách khác: thay vì áp dụng một hiệu chỉnh kiểm định đa giả thuyết chính thức cho một tập kiểm định cụ thể, họ ước lượng thống kê t ngưỡng có tính đến số lượng nhân tố *tích lũy* đã được toàn giới nghiên cứu kiểm định. Sử dụng khung Bayes được hiệu chỉnh theo 316 nhân tố đã công bố, họ suy ra:
- Với một nhân tố mới được kiểm định riêng lẻ: $|t| > 2.0$ (thông thường)
- Tính đến 100 kiểm định trước đó: $|t| > 2.8$
- Tính đến 316 kiểm định trước đó: $|t| > 3.0$
- Tính đến mọi kiểm định có thể hình dung: $|t| > 3.78$
Logic ở đây là ngay cả khi nghiên cứu của *bạn* chỉ kiểm định một nhân tố, việc hàng trăm nhà nghiên cứu đã kiểm định hàng trăm tín hiệu trước bạn có nghĩa là xác suất tiên nghiệm để tín hiệu cụ thể của bạn là một biến dự báo thực sự là thấp.
```{python}
#| label: hlz-threshold
#| eval: false
#| code-summary: "Apply the Harvey-Liu-Zhu t-statistic thresholds"
hlz_thresholds = {
'Conventional (|t| > 2.0)': 2.0,
'HLZ 100 tests (|t| > 2.8)': 2.8,
'HLZ 316 tests (|t| > 3.0)': 3.0,
'HLZ conservative (|t| > 3.78)': 3.78
}
print("Anomalies Surviving Different t-Thresholds:")
print(f"{'Threshold':<35} {'Surviving':>10} {'% of Zoo':>10}")
print("-" * 55)
for name, thresh in hlz_thresholds.items():
n_survive = (zoo_df['abs_t'] > thresh).sum()
pct = n_survive / len(zoo_df) * 100
print(f"{name:<35} {n_survive:>10} {pct:>10.1f}%")
```
## So sánh các phương pháp {#sec-factor-zoo-comparison}
```{python}
#| label: fig-correction-comparison
#| eval: false
#| fig-cap: "Số hiện tượng bất thường vượt qua từng hiệu chỉnh kiểm định đa giả thuyết. Bonferroni và Holm (kiểm soát FWER) thận trọng nhất, trong khi BH (kiểm soát FDR) giữ lại nhiều phát hiện hơn. Ngưỡng HLZ |t| > 3.0 nằm giữa hai nhóm này. Khoảng cách giữa BH và Bonferroni biểu thị các hiện tượng bất thường có khả năng là thật (FDR thấp) nhưng không thể đảm bảo cho từng hiện tượng riêng lẻ (FWER không được kiểm soát)."
#| code-summary: "Compare all correction methods side by side"
comparison = {
'Unadjusted (|t| > 2)': (zoo_df['abs_t'] > 2.0).sum(),
'Bonferroni': zoo_df['bonf_reject'].sum(),
'Holm': zoo_df['holm_reject'].sum(),
'BH (q=0.05)': zoo_df['bh_reject'].sum(),
'BH (q=0.10)': zoo_df['bh10_reject'].sum(),
'HLZ (|t| > 3.0)': (zoo_df['abs_t'] > 3.0).sum(),
'HLZ (|t| > 3.78)': (zoo_df['abs_t'] > 3.78).sum(),
}
fig, ax = plt.subplots(figsize=(10, 5))
names = list(comparison.keys())
counts = list(comparison.values())
colors = ['#BDC3C7', '#C0392B', '#E74C3C', '#27AE60', '#2ECC71',
'#2C5F8A', '#1A3C5A']
bars = ax.barh(range(len(names)), counts, color=colors,
alpha=0.85, edgecolor='white')
ax.set_yticks(range(len(names)))
ax.set_yticklabels(names)
ax.set_xlabel('Number of Surviving Anomalies')
ax.set_title('Anomalies Surviving Multiple Testing Corrections')
for i, (name, count) in enumerate(zip(names, counts)):
ax.text(count + 0.3, i, str(count), va='center', fontsize=10)
ax.invert_yaxis()
plt.tight_layout()
plt.show()
```
```{python}
#| label: fig-survival-by-category
#| eval: false
#| fig-cap: "Tỷ lệ sống sót của các hiện tượng bất thường theo nhóm dưới thủ tục BH (q=0.05). Các nhóm có tỷ lệ sống sót cao (chẳng hạn động lượng, khả năng sinh lời) chứa nhiều tín hiệu vững hơn. Các nhóm có tỷ lệ sống sót thấp có thể chủ yếu chứa các mối liên hệ giả. Tổng chiều cao của mỗi cột cho biết số hiện tượng bất thường được kiểm định trong nhóm đó."
#| code-summary: "Plot survival rates by anomaly category"
cat_survival = (
zoo_df.groupby('category')
.agg(
n_total=('abs_t', 'count'),
n_survive_bh=('bh_reject', 'sum'),
n_survive_hlz=('abs_t', lambda x: (x > 3.0).sum()),
mean_abs_t=('abs_t', 'mean')
)
.sort_values('mean_abs_t', ascending=False)
)
cat_survival['survival_rate_bh'] = (
cat_survival['n_survive_bh'] / cat_survival['n_total']
)
fig, ax = plt.subplots(figsize=(12, 5))
x = np.arange(len(cat_survival))
ax.bar(x, cat_survival['n_total'], color='#BDC3C7', alpha=0.5,
label='Total tested', edgecolor='white')
ax.bar(x, cat_survival['n_survive_bh'], color='#27AE60', alpha=0.85,
label='Survive BH (q=0.05)', edgecolor='white')
ax.set_xticks(x)
ax.set_xticklabels(cat_survival.index, rotation=30, fontsize=9)
ax.set_ylabel('Number of Anomalies')
ax.set_title('Anomaly Survival by Category')
ax.legend()
# Add survival rate text
for i, (_, row) in enumerate(cat_survival.iterrows()):
if row['n_total'] > 0:
ax.text(i, row['n_total'] + 0.2,
f"{row['survival_rate_bh']:.0%}",
ha='center', fontsize=9, color='#27AE60')
plt.tight_layout()
plt.show()
```
## Ước lượng tỷ lệ phát hiện sai {#sec-factor-zoo-fdp}
### Cách tiếp cận của Storey
Thay vì kiểm soát FDR ở một mức định trước, chúng ta có thể *ước lượng* tỷ lệ các giả thuyết được kiểm định thực sự là giả thuyết không ($\pi_0$) và tỷ lệ phát hiện sai hàm ý (FDP) tại một ngưỡng bất kỳ. @storey2003positive đề xuất ước lượng $\pi_0$ từ phân phối của các giá trị p: dưới giả thuyết không, các giá trị p phân phối đều trên $[0, 1]$, nên mật độ của các giá trị p trong vùng "phẳng" (xa số không) cho biết $\pi_0$.
$$
\hat{\pi}_0(\lambda) = \frac{\#\{p_i > \lambda\}}{M(1 - \lambda)}
$$ {#eq-pi0}
với một tham số điều chỉnh $\lambda$. Tỷ lệ phát hiện sai ước lượng tại ngưỡng $t^*$ khi đó là:
$$
\widehat{\text{FDP}}(t^*) = \frac{\hat{\pi}_0 \cdot M \cdot P(|t| > t^* | H_0)}{\#\{|t_i| > t^*\}}
$$ {#eq-fdp}
```{python}
#| label: storey-pi0
#| eval: false
#| code-summary: "Estimate the proportion of true nulls using Storey's method"
def estimate_pi0(p_values, lambdas=np.arange(0.05, 0.95, 0.05)):
"""
Estimate pi_0 (proportion of true nulls) using Storey (2003).
Uses bootstrap to select optimal lambda.
"""
M = len(p_values)
pi0_estimates = []
for lam in lambdas:
n_above = (p_values > lam).sum()
pi0 = n_above / (M * (1 - lam))
pi0_estimates.append({'lambda': lam, 'pi0': min(pi0, 1.0)})
pi0_df = pd.DataFrame(pi0_estimates)
# Use the median of estimates for robustness
# (alternatives: bootstrap, spline smoothing)
pi0_hat = pi0_df['pi0'].median()
return pi0_hat, pi0_df
pi0_hat, pi0_df = estimate_pi0(zoo_df['p_value'].values)
print(f"Estimated π₀ (proportion of true nulls): {pi0_hat:.3f}")
print(f"Implied proportion of true factors: {1 - pi0_hat:.3f}")
print(f"Estimated number of true factors: {(1 - pi0_hat) * len(zoo_df):.0f} "
f"out of {len(zoo_df)}")
```
```{python}
#| label: fig-pvalue-distribution
#| eval: false
#| fig-cap: "Phân phối của các giá trị p trên toàn bộ vườn thú các hiện tượng bất thường. Panel A hiển thị biểu đồ tần suất cùng tỷ lệ giả thuyết không ước lượng (π₀ = đường ngang). Phần mật độ dôi ra gần số không biểu thị các phát hiện thật. Panel B hiển thị π₀ ước lượng như một hàm của tham số điều chỉnh λ. Vùng tương đối phẳng ở các giá trị λ cao hơn cho ước lượng ổn định nhất."
#| code-summary: "Plot p-value distribution and pi_0 estimation"
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
# Panel A: p-value histogram
axes[0].hist(zoo_df['p_value'], bins=20, density=True,
color='#2C5F8A', alpha=0.7, edgecolor='white')
axes[0].axhline(y=pi0_hat, color='#C0392B', linewidth=2, linestyle='--',
label=f'π₀ = {pi0_hat:.2f} (null density)')
axes[0].axhline(y=1, color='gray', linewidth=1, linestyle=':',
label='Uniform (all null)')
axes[0].set_xlabel('p-value')
axes[0].set_ylabel('Density')
axes[0].set_title('Panel A: p-Value Distribution')
axes[0].legend()
# Panel B: pi_0 estimates across lambda
axes[1].plot(pi0_df['lambda'], pi0_df['pi0'],
color='#2C5F8A', linewidth=2, marker='o', markersize=4)
axes[1].axhline(y=pi0_hat, color='#C0392B', linestyle='--',
linewidth=1.5, label=f'Median π₀ = {pi0_hat:.2f}')
axes[1].set_xlabel('λ')
axes[1].set_ylabel('π₀(λ)')
axes[1].set_title('Panel B: π₀ Estimates by λ')
axes[1].legend()
axes[1].set_ylim([0, 1.05])
plt.tight_layout()
plt.show()
```
### Đường cong FDP
Sử dụng $\hat{\pi}_0$, chúng ta tính FDP ước lượng tại mọi ngưỡng t có thể:
```{python}
#| label: fig-fdp-curve
#| eval: false
#| fig-cap: "Tỷ lệ phát hiện sai ước lượng như một hàm của ngưỡng thống kê t. Tại ngưỡng thông thường |t| > 2.0, một phần đáng kể các phát hiện có thể là sai. FDP giảm xuống dưới 5% quanh mức |t| > 3.0, gần khớp với khuyến nghị của Harvey-Liu-Zhu."
#| code-summary: "Compute and plot the FDP curve"
t_thresholds = np.arange(1.0, 5.01, 0.1)
fdp_curve = []
for t_thresh in t_thresholds:
n_rejected = (zoo_df['abs_t'] > t_thresh).sum()
if n_rejected == 0:
fdp_curve.append({'t_threshold': t_thresh, 'n_rejected': 0,
'fdp': 0})
continue
# Expected false discoveries under null
p_null_reject = 2 * (1 - stats.norm.cdf(t_thresh))
expected_false = pi0_hat * len(zoo_df) * p_null_reject
fdp = min(expected_false / n_rejected, 1.0)
fdp_curve.append({'t_threshold': t_thresh, 'n_rejected': n_rejected,
'fdp': fdp})
fdp_df = pd.DataFrame(fdp_curve)
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
# Panel A: FDP curve
axes[0].plot(fdp_df['t_threshold'], fdp_df['fdp'] * 100,
color='#C0392B', linewidth=2)
axes[0].axhline(y=5, color='gray', linestyle='--', linewidth=1,
label='5% FDP')
axes[0].axhline(y=10, color='gray', linestyle=':', linewidth=1,
label='10% FDP')
axes[0].axvline(x=2.0, color='#E67E22', linestyle='--',
linewidth=1, label='|t| = 2.0')
axes[0].axvline(x=3.0, color='#2C5F8A', linestyle='--',
linewidth=1, label='|t| = 3.0')
axes[0].set_xlabel('|t| Threshold')
axes[0].set_ylabel('Estimated FDP (%)')
axes[0].set_title('Panel A: False Discovery Proportion')
axes[0].legend(fontsize=8)
axes[0].set_ylim([-2, 80])
# Panel B: Number of rejections
axes[1].plot(fdp_df['t_threshold'], fdp_df['n_rejected'],
color='#2C5F8A', linewidth=2)
axes[1].set_xlabel('|t| Threshold')
axes[1].set_ylabel('Number of Discoveries')
axes[1].set_title('Panel B: Discoveries vs Threshold')
plt.tight_layout()
plt.show()
```
## Kiểm định đa giả thuyết dựa trên bootstrap {#sec-factor-zoo-bootstrap}
### White Reality Check
Các hiệu chỉnh giải tích giả định các phân phối đã biết. Khi phân phối lợi nhuận không chuẩn (đuôi dày, lệch, tự tương quan chuỗi), như trường hợp ở Việt Nam, các phương pháp bootstrap cho suy luận đáng tin cậy hơn.
@white2000reality đề xuất một kiểm tra thực tế bằng bootstrap (reality check) có tính đến việc khai thác dữ liệu (data snooping) trên nhiều chiến lược. Giả thuyết không là chiến lược tốt nhất có lợi nhuận kỳ vọng bằng không. Thủ tục như sau:
1. Tạo $B$ mẫu bootstrap của chuỗi thời gian (bootstrap theo khối để bảo toàn sự phụ thuộc chuỗi).
2. Với mỗi mẫu bootstrap, tính thống kê t cho mọi nhân tố.
3. Giá trị p của nhân tố $i$ là tỷ lệ các mẫu bootstrap trong đó thống kê t lớn nhất trên tất cả các nhân tố vượt quá thống kê t quan sát được của nhân tố $i$.
```{python}
#| label: bootstrap-rc
#| eval: false
#| code-summary: "Implement White's Reality Check via block bootstrap"
def white_reality_check(factor_returns_dict, n_bootstrap=1000,
block_size=6, seed=42):
"""
White (2000) Reality Check for data snooping.
Parameters
----------
factor_returns_dict : dict
{factor_name: pd.Series of monthly returns}
n_bootstrap : int
Number of bootstrap replications
block_size : int
Block length for circular block bootstrap
Returns
-------
DataFrame with original t-stats and bootstrap p-values
"""
rng = np.random.default_rng(seed)
# Align all factor returns to common dates
all_names = list(factor_returns_dict.keys())
common_dates = None
for name in all_names:
dates = set(factor_returns_dict[name].index)
common_dates = dates if common_dates is None else common_dates & dates
common_dates = sorted(common_dates)
T = len(common_dates)
# Build return matrix (T x M)
M = len(all_names)
return_matrix = np.column_stack([
factor_returns_dict[name].reindex(common_dates).values
for name in all_names
])
# Observed t-statistics
means = np.nanmean(return_matrix, axis=0)
ses = np.nanstd(return_matrix, axis=0) / np.sqrt(T)
t_obs = means / np.where(ses > 0, ses, np.nan)
# Block bootstrap
n_blocks = int(np.ceil(T / block_size))
max_t_bootstrap = np.zeros(n_bootstrap)
for b in range(n_bootstrap):
# Circular block bootstrap
block_starts = rng.integers(0, T, size=n_blocks)
indices = np.concatenate([
np.arange(start, start + block_size) % T
for start in block_starts
])[:T]
boot_returns = return_matrix[indices, :]
# Center under null (subtract original mean)
boot_centered = boot_returns - means[np.newaxis, :]
# Compute t-stats under null
boot_means = np.nanmean(boot_centered, axis=0)
boot_ses = np.nanstd(boot_centered, axis=0) / np.sqrt(T)
boot_t = boot_means / np.where(boot_ses > 0, boot_ses, np.nan)
max_t_bootstrap[b] = np.nanmax(np.abs(boot_t))
# Bootstrap p-values (one-sided: is obs_t extreme relative to max?)
boot_pvals = np.array([
np.mean(max_t_bootstrap >= abs(t)) if np.isfinite(t) else 1.0
for t in t_obs
])
results = pd.DataFrame({
'factor': all_names,
't_stat_obs': t_obs,
'boot_pval': boot_pvals
}).set_index('factor')
return results, max_t_bootstrap
# Collect factor return series
factor_returns_dict = {
name: zoo_results[name]['returns']
for name in zoo_results
if 'returns' in zoo_results[name]
}
boot_results, max_t_dist = white_reality_check(
factor_returns_dict, n_bootstrap=1000, block_size=6
)
boot_results = boot_results.sort_values('t_stat_obs', ascending=False)
print("White Reality Check Results (top 20):")
print(boot_results.head(20).round(4).to_string())
print(f"\nFactors surviving (boot p < 0.05): "
f"{(boot_results['boot_pval'] < 0.05).sum()}")
```
```{python}
#| label: fig-bootstrap-distribution
#| eval: false
#| fig-cap: "Phân phối bootstrap của thống kê t lớn nhất trên tất cả các nhân tố dưới giả thuyết không là không có hiệu ứng thực. Các đường thẳng đứng đánh dấu các thống kê t quan sát được của các hiện tượng bất thường hàng đầu. Chỉ những hiện tượng bất thường có t quan sát được vượt phân vị thứ 95 của phân phối này mới vượt qua hiệu chỉnh khai thác dữ liệu."
#| code-summary: "Plot the bootstrap null distribution"
fig, ax = plt.subplots(figsize=(12, 5))
ax.hist(max_t_dist, bins=40, density=True, color='#BDC3C7',
alpha=0.7, edgecolor='white', label='Bootstrap null\n(max |t| across all factors)')
# 95th percentile
p95 = np.percentile(max_t_dist, 95)
ax.axvline(x=p95, color='#C0392B', linewidth=2, linestyle='--',
label=f'95th percentile = {p95:.2f}')
# Observed top factor t-statistics
top_5 = boot_results.head(5)
for i, (name, row) in enumerate(top_5.iterrows()):
ax.axvline(x=abs(row['t_stat_obs']),
color=plt.cm.Set1(i), linewidth=1.5, linestyle=':',
label=f'{name}: |t| = {abs(row["t_stat_obs"]):.2f}')
ax.set_xlabel('Maximum |t-statistic|')
ax.set_ylabel('Density')
ax.set_title("White's Reality Check: Bootstrap Null Distribution")
ax.legend(fontsize=8, loc='upper right')
plt.tight_layout()
plt.show()
```
### Bootstrap step-down Romano-Wolf
Thủ tục step-down của @romano2005exact cải tiến White Reality Check bằng cách lần lượt loại bỏ các giả thuyết đã bị bác bỏ và kiểm định lại các giả thuyết còn lại, nhờ đó tăng sức mạnh ở mỗi bước:
```{python}
#| label: romano-wolf
#| eval: false
#| code-summary: "Implement the Romano-Wolf step-down bootstrap procedure"
def romano_wolf_stepdown(factor_returns_dict, n_bootstrap=1000,
block_size=6, alpha=0.05, seed=42):
"""
Romano-Wolf (2005) step-down procedure for FWER control.
More powerful than single-step methods because rejected
hypotheses are removed before re-testing.
"""
rng = np.random.default_rng(seed)
all_names = list(factor_returns_dict.keys())
common_dates = None
for name in all_names:
dates = set(factor_returns_dict[name].index)
common_dates = dates if common_dates is None else common_dates & dates
common_dates = sorted(common_dates)
T = len(common_dates)
M = len(all_names)
return_matrix = np.column_stack([
factor_returns_dict[name].reindex(common_dates).values
for name in all_names
])
means = np.nanmean(return_matrix, axis=0)
ses = np.nanstd(return_matrix, axis=0) / np.sqrt(T)
t_obs = np.abs(means / np.where(ses > 0, ses, np.nan))
# Sort by observed t-stat (descending)
order = np.argsort(-t_obs)
t_sorted = t_obs[order]
names_sorted = [all_names[i] for i in order]
# Step-down procedure
rejected = set()
remaining = set(range(M))
for step in range(M):
if not remaining:
break
remaining_idx = sorted(remaining)
# Bootstrap max t-stat among remaining
n_blocks = int(np.ceil(T / block_size))
max_t_boot = np.zeros(n_bootstrap)
for b in range(n_bootstrap):
block_starts = rng.integers(0, T, size=n_blocks)
indices = np.concatenate([
np.arange(start, start + block_size) % T
for start in block_starts
])[:T]
boot_returns = return_matrix[indices][:, remaining_idx]
boot_centered = boot_returns - means[remaining_idx]
boot_means = np.nanmean(boot_centered, axis=0)
boot_ses = np.nanstd(boot_centered, axis=0) / np.sqrt(T)
boot_t = np.abs(boot_means / np.where(boot_ses > 0, boot_ses, np.nan))
max_t_boot[b] = np.nanmax(boot_t)
# Critical value
cv = np.percentile(max_t_boot, (1 - alpha) * 100)
# Test the most significant remaining factor
# Find the factor with highest t among remaining
best_remaining = max(remaining, key=lambda j: t_obs[j])
if t_obs[best_remaining] > cv:
rejected.add(best_remaining)
remaining.remove(best_remaining)
else:
break # Cannot reject any more
results = pd.DataFrame({
'factor': all_names,
't_stat': t_obs,
'rejected': [i in rejected for i in range(M)]
}).sort_values('t_stat', ascending=False)
return results
rw_results = romano_wolf_stepdown(
factor_returns_dict, n_bootstrap=500, block_size=6
)
print(f"Romano-Wolf rejections: {rw_results['rejected'].sum()} / {len(rw_results)}")
print("\nRejected factors:")
print(rw_results[rw_results['rejected']][['factor', 't_stat']].to_string())
```
## Kiểm định ngoài mẫu {#sec-factor-zoo-oos}
### Suy giảm trước và sau khi công bố
@mclean2016does thấy rằng lợi nhuận của các hiện tượng bất thường giảm 32% sau khi công bố trên thị trường Mỹ, gợi ý rằng khoảng một phần ba phần bù là do định giá sai được điều chỉnh nhờ giao dịch của các nhà đầu tư có thông tin. Với Việt Nam, chúng ta dùng cách chia theo chuỗi thời gian làm phép xấp xỉ cho kiểm định ngoài mẫu:
```{python}
#| label: oos-split
#| eval: false
#| code-summary: "Compare factor premia in-sample vs out-of-sample"
# Split at midpoint of available data
all_dates = sorted(panel['month_end'].unique())
mid_date = all_dates[len(all_dates) // 2]
oos_comparison = []
for name in zoo_results:
ls = zoo_results[name]['returns']
in_sample = ls[ls.index <= mid_date]
out_sample = ls[ls.index > mid_date]
if len(in_sample) < 36 or len(out_sample) < 36:
continue
is_mean = in_sample.mean() * 12
is_t = in_sample.mean() / (in_sample.std() / np.sqrt(len(in_sample)))
oos_mean = out_sample.mean() * 12
oos_t = out_sample.mean() / (out_sample.std() / np.sqrt(len(out_sample)))
oos_comparison.append({
'signal': name,
'is_return': is_mean,
'is_t': is_t,
'oos_return': oos_mean,
'oos_t': oos_t,
'decay_pct': (1 - oos_mean / is_mean) * 100 if is_mean != 0 else np.nan,
'same_sign': np.sign(is_mean) == np.sign(oos_mean),
'category': zoo_results[name].get('category', 'Other')
})
oos_df = pd.DataFrame(oos_comparison)
print(f"Split date: {mid_date.strftime('%Y-%m')}")
print(f"Factors with same sign IS and OOS: "
f"{oos_df['same_sign'].mean():.1%}")
print(f"Average OOS decay: {oos_df['decay_pct'].median():.1f}%")
```
```{python}
#| label: fig-is-oos
#| eval: false
#| fig-cap: "Lợi nhuận nhân tố trong mẫu so với ngoài mẫu. Panel A là biểu đồ phân tán trong đó mỗi điểm là một hiện tượng bất thường; các điểm nằm trên đường 45 độ mạnh lên ngoài mẫu, còn các điểm nằm dưới thì yếu đi. Panel B hiển thị tốc độ suy giảm trung bình theo nhóm. Các nhân tố được tái lập ngoài mẫu với độ lớn tương tự là đáng tin cậy nhất. Một mẫu hình suy giảm chung phù hợp với McLean-Pontiff (2016) là bằng chứng cho việc định giá sai trong quá khứ đang được điều chỉnh."
#| code-summary: "Scatter in-sample vs out-of-sample returns"
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
# Panel A: Scatter
category_colors = {
'Value': '#2C5F8A', 'Size': '#1ABC9C', 'Profitability': '#27AE60',
'Investment': '#E67E22', 'Momentum': '#C0392B', 'Risk': '#8E44AD',
'Liquidity': '#3498DB', 'Quality': '#F1C40F', 'Dividend': '#95A5A6'
}
for cat in oos_df['category'].unique():
subset = oos_df[oos_df['category'] == cat]
axes[0].scatter(subset['is_return'] * 100, subset['oos_return'] * 100,
color=category_colors.get(cat, '#BDC3C7'),
s=50, alpha=0.7, label=cat, edgecolors='white')
lim = max(abs(oos_df['is_return'].max()), abs(oos_df['oos_return'].max())) * 100 + 2
axes[0].plot([-lim, lim], [-lim, lim], 'k--', linewidth=1, alpha=0.5)
axes[0].plot([-lim, lim], [0, 0], color='gray', linewidth=0.5)
axes[0].plot([0, 0], [-lim, lim], color='gray', linewidth=0.5)
axes[0].set_xlabel('In-Sample Return (% ann.)')
axes[0].set_ylabel('Out-of-Sample Return (% ann.)')
axes[0].set_title('Panel A: IS vs OOS Factor Returns')
axes[0].legend(fontsize=7, ncol=2)
# Panel B: Decay by category
cat_decay = (
oos_df.groupby('category')
.agg(
median_decay=('decay_pct', 'median'),
pct_same_sign=('same_sign', 'mean'),
n=('signal', 'count')
)
.sort_values('median_decay')
)
colors_bar = [category_colors.get(cat, '#BDC3C7') for cat in cat_decay.index]
axes[1].barh(range(len(cat_decay)), cat_decay['median_decay'],
color=colors_bar, alpha=0.85, edgecolor='white')
axes[1].set_yticks(range(len(cat_decay)))
axes[1].set_yticklabels(cat_decay.index)
axes[1].set_xlabel('Median OOS Decay (%)')
axes[1].set_title('Panel B: OOS Decay by Category')
axes[1].axvline(x=0, color='gray', linewidth=0.8)
plt.tight_layout()
plt.show()
```
## Những nhân tố nào còn tồn tại? {#sec-factor-zoo-survivors}
### Bộ lọc đa ngưỡng
Bây giờ chúng ta kết hợp mọi bằng chứng, gồm các hiệu chỉnh kiểm định đa giả thuyết, hiệu quả ngoài mẫu và tính hợp lý về mặt kinh tế, để xác định những nhân tố xứng đáng được tin cậy trên thị trường Việt Nam.
```{python}
#| label: multi-hurdle
#| eval: false
#| code-summary: "Apply a multi-hurdle filter to identify robust factors"
# Merge all test results
final = zoo_df.copy()
# Add OOS results
oos_lookup = oos_df.set_index('signal')
for col in ['oos_return', 'oos_t', 'same_sign']:
final[col] = final.index.map(lambda x: oos_lookup.loc[x, col]
if x in oos_lookup.index else np.nan)
# Add bootstrap results
boot_lookup = boot_results
final['boot_pval'] = final.index.map(
lambda x: boot_lookup.loc[x, 'boot_pval']
if x in boot_lookup.index else np.nan
)
# Multi-hurdle criteria
final['pass_t2'] = final['abs_t'] > 2.0
final['pass_t3'] = final['abs_t'] > 3.0
final['pass_bh'] = final['bh_reject']
final['pass_boot'] = final['boot_pval'] < 0.05
final['pass_oos_sign'] = final['same_sign'] == True
final['pass_oos_t'] = final['oos_t'].abs() > 1.5 # Relaxed OOS threshold
# Count hurdles passed
hurdle_cols = ['pass_t2', 'pass_t3', 'pass_bh', 'pass_boot',
'pass_oos_sign', 'pass_oos_t']
final['n_hurdles'] = final[hurdle_cols].sum(axis=1)
# "Robust" = passes at least 5 of 6 hurdles
final['robust'] = final['n_hurdles'] >= 5
robust_factors = final[final['robust']].sort_values('abs_t', ascending=False)
print(f"Multi-Hurdle Filter Results:")
print(f" Pass |t| > 2.0: {final['pass_t2'].sum()}")
print(f" Pass |t| > 3.0: {final['pass_t3'].sum()}")
print(f" Pass BH (q=0.05): {final['pass_bh'].sum()}")
print(f" Pass bootstrap: {final['pass_boot'].sum()}")
print(f" Pass OOS sign: {final['pass_oos_sign'].sum()}")
print(f" Pass OOS |t| > 1.5: {final['pass_oos_t'].sum()}")
print(f"\n ROBUST (≥5/6 hurdles): {final['robust'].sum()}")
if len(robust_factors) > 0:
print(f"\nRobust Factors:")
display_cols = ['ann_return', 't_stat', 'bh_p_corrected',
'oos_return', 'oos_t', 'n_hurdles', 'category']
available_cols = [c for c in display_cols if c in robust_factors.columns]
print(robust_factors[available_cols].round(3).to_string())
```
```{python}
#| label: fig-hurdle-heatmap
#| eval: false
#| fig-cap: "Kết quả kiểm định đa ngưỡng cho 25 hiện tượng bất thường hàng đầu, xếp hạng theo thống kê t trong mẫu. Mỗi cột đại diện cho một ngưỡng thống kê; ô màu xanh lá cho biết vượt qua, ô màu đỏ cho biết không đạt. Chỉ những hiện tượng bất thường vượt qua nhiều kiểm định độc lập mới đáng được xem xét nghiêm túc. Cột ngoài cùng bên phải tổng hợp số ngưỡng đã vượt qua."
#| code-summary: "Heatmap of multi-hurdle results"
top25 = final.sort_values('abs_t', ascending=False).head(25)
fig, ax = plt.subplots(figsize=(10, 10))
hurdle_matrix = top25[hurdle_cols].astype(float).values
sns.heatmap(hurdle_matrix, ax=ax,
xticklabels=['|t|>2', '|t|>3', 'BH', 'Boot', 'OOS sign', 'OOS |t|'],
yticklabels=top25.index,
cmap=['#E74C3C', '#27AE60'],
cbar=False, linewidths=0.5, linecolor='white',
annot=np.where(hurdle_matrix == 1, '✓', '✗'),
fmt='')
# Add hurdle count column
for i, (_, row) in enumerate(top25.iterrows()):
ax.text(len(hurdle_cols) + 0.3, i + 0.5,
f"{int(row['n_hurdles'])}/6",
va='center', fontsize=9,
fontweight='bold' if row['n_hurdles'] >= 5 else 'normal',
color='#27AE60' if row['n_hurdles'] >= 5 else '#C0392B')
ax.set_title('Multi-Hurdle Test Results (Top 25 Anomalies)')
plt.tight_layout()
plt.show()
```
## Hàm ý cho nghiên cứu nhân tố ở Việt Nam {#sec-factor-zoo-implications}
Phân tích trong chương này đưa ra một số hàm ý thực tiễn:
**Ngưỡng thông thường là không đủ.** Với $M \approx 50$ hiện tượng bất thường được kiểm định đồng thời, chỉ yêu cầu $|t| > 2.0$ gần như chắc chắn dẫn đến các phát hiện sai. Các nhà nghiên cứu Việt Nam nên áp dụng ngưỡng tối thiểu $|t| > 3.0$ cho từng nhân tố, nhất quán với @harvey2016and, và nên báo cáo các giá trị p đã điều chỉnh theo BH cho mọi nghiên cứu kiểm định nhiều tín hiệu.
**Mặt cắt ngang nhỏ làm khuếch đại nhiễu.** Với 600-800 cổ phiếu, các danh mục ngũ phân vị chỉ chứa 120-160 cổ phiếu mỗi danh mục. Điều này tạo ra lợi nhuận danh mục nhiễu hơn và khoảng tin cậy rộng hơn so với các nghiên cứu ở Mỹ với hơn 4.000 cổ phiếu. Những chiến lược có vẻ có ý nghĩa ở Mỹ với $|t| = 2.5$ có thể chỉ có $|t| < 1.5$ ở Việt Nam đơn giản vì mẫu nhỏ hơn. Đây không phải là bằng chứng cho thấy nhân tố không tồn tại ở Việt Nam, mà là bằng chứng cho thấy dữ liệu Việt Nam thiếu sức mạnh để phát hiện ra nó.
**Kiểm định ngoài mẫu là thiết yếu.** Cách chia theo chuỗi thời gian cho thấy sự suy giảm đáng kể ở nhiều hiện tượng bất thường, nhất quán với @mclean2016does. Những nhân tố vượt qua cả trong mẫu lẫn ngoài mẫu với dấu và độ lớn nhất quán là hiếm và có giá trị.
**Nhóm tín hiệu đóng vai trò quan trọng.** Các tín hiệu động lượng và khả năng sinh lời có xu hướng có tỷ lệ tái lập cao nhất giữa các thị trường [@fama2012size; @jacobs2020anomalies]. Các tín hiệu giá trị và đầu tư mang tính đặc thù theo quốc gia nhiều hơn [@griffin2002fama]. Các tín hiệu thanh khoản có khả năng mạnh nhất ở các thị trường mới nổi như Việt Nam, nơi các ma sát giao dịch lớn nhất.
**Đường đi của lãi suất là một phần của mẫu.** Giai đoạn mẫu của Việt Nam bị chi phối bởi một biến động một chiều của lãi suất: lãi suất chính sách đạt đỉnh gần 14% năm 2011 và giảm xuống khoảng 4,5% vào năm 2023, và lợi suất trái phiếu chính phủ cũng giảm theo. @binsbergen2026factor chỉ ra rằng ở Mỹ, một sự sụt giảm dài hạn tương tự đủ để phân loại lại một tỷ lệ đáng kể các hiện tượng bất thường khi lợi nhuận vượt mức được đo so với một danh mục trái phiếu có kỳ hạn bình quân tương ứng thay vì lãi suất ngắn hạn. Vì sự sụt giảm ở Việt Nam diễn ra trong một khoảng thời gian ngắn hơn nhiều, bất kỳ danh mục mua–bán khống nào có các vị thế khác nhau về kỳ hạn (giá trị so với tăng trưởng, tỷ suất cổ tức cao so với thấp, đòn bẩy cao so với thấp) có thể đang nắm bắt sự định giá lại đó thay vì một phần bù rủi ro mặt cắt ngang. Tối thiểu, các nhà nghiên cứu nên báo cáo tỷ suất cổ tức hoặc chênh lệch kỳ hạn giữa các vị thế mua và bán, và thận trọng với những nhân tố có phần bù tập trung trong giai đoạn lãi suất sụt giảm 2011-2015.
**Báo cáo đầy đủ bộ chỉ tiêu.** Mọi nghiên cứu đưa ra một nhân tố mới cho thị trường Việt Nam nên báo cáo: (i) thống kê t thô, (ii) giá trị p đã điều chỉnh theo BH với số kiểm định trong nghiên cứu, (iii) hiệu quả ngoài mẫu trong một giai đoạn được giữ lại, và (iv) độ nhạy với các lựa chọn xây dựng (từ chương trước). Chuẩn báo cáo này sẽ cải thiện đáng kể độ tin cậy của nghiên cứu nhân tố ở Việt Nam.
## Tóm tắt {#sec-factor-zoo-summary}
| Phương pháp | Kiểm soát | Ngưỡng (xấp xỉ) | Số còn lại | Triết lý |
|---------------|---------------|---------------|---------------|---------------|
| Không điều chỉnh | Kiểm định đơn | \|t\| \> 2.0 | Nhiều | Không hiệu chỉnh |
| Bonferroni | FWER | \|t\| \> 3.3--4.0 | Ít | Không có dương tính giả |
| Holm | FWER | Theo từng bước | Ít | Ít thận trọng hơn một chút |
| BH | FDR ở mức q | Thích ứng | Trung bình | Chấp nhận một số dương tính giả |
| HLZ | Toàn giới nghiên cứu | \|t\| \> 3.0--3.78 | Trung bình | Phân phối tiên nghiệm trên mọi nhân tố đã kiểm định |
| Bootstrap White/RW | FWER (theo dữ liệu) | Giá trị tới hạn bootstrap | Ít | Tính đến phân phối không chuẩn |
| Storey FDP | Ước lượng FDP | Liên tục | Chẩn đoán | Ước lượng tỷ lệ giả thuyết không thực sự đúng |
| Đa ngưỡng | Kết hợp | Nhiều kiểm định | Vững nhất | Đòi hỏi bằng chứng hội tụ |
: Tóm tắt các phương pháp kiểm định đa giả thuyết và đặc tính của chúng. {#tbl-summary}
Vườn thú nhân tố không phải là hiện tượng riêng của nước Mỹ. Cám dỗ kiểm định nhiều tín hiệu và chỉ báo cáo những tín hiệu cho kết quả tốt nhất tồn tại ở mọi thị trường. Ở Việt Nam, nơi dữ liệu ngắn hơn, nhiễu hơn và cộng đồng học thuật đang phát triển nhanh, rủi ro phát hiện sai là cao. Các công cụ được phát triển trong chương này, gồm điều chỉnh BH, kiểm tra thực tế bằng bootstrap, ước lượng $\pi_0$ của Storey, cách chia ngoài mẫu và bộ lọc đa ngưỡng, cung cấp một khung nguyên tắc để phân tách các nhân tố thực sự khỏi nhiễu thống kê.
Kết luận trung thực có lẽ sẽ khiêm tốn: trong số hàng chục hiện tượng bất thường được ghi nhận ở các thị trường phát triển, chỉ một số ít vượt qua kiểm định đa giả thuyết nghiêm ngặt trên dữ liệu Việt Nam. Những nhân tố sống sót đó (có thể tập trung ở động lượng, khả năng sinh lời và thanh khoản) tạo nền tảng cho nghiên cứu định giá tài sản đáng tin cậy ở thị trường này.
<!-- ## Bài tập {#sec-factor-zoo-exercises}
1. **Cách tiếp cận Bayes.** Triển khai khung Bayes của @harvey2021lucky để phân biệt các nhân tố may mắn với các nhân tố thực sự. Gán một xác suất tiên nghiệm $\pi$ rằng một nhân tố bất kỳ là biến dự báo thực sự, và tính xác suất hậu nghiệm để nó là nhân tố thật, dựa trên thống kê t và số kiểm định. Việc thay đổi $\pi$ ảnh hưởng thế nào đến việc nhân tố nào còn tồn tại?
2. **BH điều chỉnh theo phụ thuộc.** Thủ tục BH tiêu chuẩn giả định các thống kê kiểm định độc lập (hoặc PRDS). Hãy ước lượng cấu trúc tương quan của lợi nhuận nhân tố và triển khai thủ tục Benjamini-Yekutieli (BY), vốn kiểm soát FDR dưới sự phụ thuộc bất kỳ. So với BH, số nhân tố còn tồn tại ít hơn bao nhiêu?
3. **Mô phỏng thiên lệch công bố.** Mô phỏng một nhà nghiên cứu kiểm định 100 tín hiệu ngẫu nhiên trên dữ liệu Việt Nam, chỉ công bố những tín hiệu có $|t| > 2$, và báo cáo thống kê t trung bình của các tín hiệu được công bố. Lặp lại 1.000 lần. Phân phối của thống kê t "được công bố" như thế nào? So sánh với độ lớn hiệu ứng thực (bằng không theo cách xây dựng) ra sao?
4. **Tái lập giữa các sàn giao dịch.** Kiểm định xem các nhân tố có ý nghĩa trên HOSE có được tái lập trên HNX hay không. Nếu cùng một nhân tố có ý nghĩa một cách độc lập trên cả hai sàn, xác suất phát hiện sai thấp hơn đáng kể (vì hai mẫu độc lập một phần).
5. **Ngưỡng chi phí giao dịch.** Với mỗi nhân tố vượt qua bộ lọc đa ngưỡng, hãy tính chi phí giao dịch khứ hồi tối thiểu sẽ triệt tiêu phần bù sau chi phí. Những nhân tố có chi phí hòa vốn thấp hơn chi phí giao dịch thực tế ở Việt Nam là không thể triển khai, bất kể ý nghĩa thống kê của chúng.
6. **Tầm quan trọng của đặc trưng trong học máy.** Sử dụng mô hình rừng ngẫu nhiên (random forest) hoặc cây tăng cường độ dốc (gradient-boosted tree) để dự báo lợi nhuận tháng tiếp theo từ tất cả các tín hiệu có sẵn đồng thời. Trích xuất các thước đo tầm quan trọng của đặc trưng. So sánh xếp hạng tín hiệu theo học máy với xếp hạng theo thống kê t đơn biến. Chúng có nhất quán không?
7. **Tái lập quốc tế.** Chọn 10 nhân tố hàng đầu từ vườn thú nhân tố của Việt Nam. Kiểm định xem cùng các nhân tố đó có ý nghĩa ở Thái Lan (SET), Indonesia (IDX) hoặc Philippines (PSE) hay không, sử dụng dữ liệu tương đương. Những nhân tố được tái lập qua nhiều thị trường ASEAN ít có khả năng là kết quả giả tạo do dữ liệu đặc thù theo quốc gia.
8. **Giám sát FDP theo cửa sổ trượt.** Triển khai phiên bản cửa sổ trượt của ước lượng $\pi_0$ của Storey. Tính $\hat{\pi}_0$ cho mỗi cửa sổ 60 tháng. Tỷ lệ các nhân tố thực có thay đổi theo thời gian không? Một $\pi_0$ tăng dần (nhiều giả thuyết không thực sự đúng hơn) có thể cho thấy các hiện tượng bất thường đang bị kinh doanh chênh lệch giá triệt tiêu dần khi thị trường Việt Nam trưởng thành. -->