Contents

中文EN

Theory: Factor Models & Styles

1. Single-Factor ModelsCode2. Multi-Factor Models3. Classic Style Factors4. Expanding the Factor Universe

Practice: Factor Production & Management

5. Factor Construction and Neutralization6. Factor Testing FrameworkCode7. Portfolio Optimization and Attribution8. Factor Timing and Crowding

Extension: Macro Quant & Asset Allocation

9. Macro Data and Factors10. Macro Regime Identification11. Asset AllocationCode12. A Decade of Macro Quant NotesCode

6. Factor Testing Framework

Factor testing is a central step in factor research. Its purpose is to evaluate a factor’s predictive power, stability, and practical tradability in a systematic way.

Before a factor is used in portfolio construction or live trading, it should be thoroughly tested: does it have persistent cross-sectional explanatory power, and does it remain implementable under trading-cost constraints? An incomplete testing framework can produce backtests that look strong on the surface but fail to hold up out of sample or in live markets.

Some researchers treat a regression t-statistic above 2 as proof of a valid factor. Strong regression results, however, do not automatically imply stable stock-selection ability—cases exist where a factor looks fine in regression yet group returns are disordered (e.g., group 3 above group 4 in a quantile backtest), suggesting statistical coincidence. Layered (quantile) backtests and IC analysis are therefore essential checks.

This chapter is organized around four core modules: layered backtests, IC/IR analysis, factor-return series analysis, and turnover analysis. These dimensions complement one another and form a basic framework for assessing factor validity. A recommended workflow is: run IC analysis first to discard weak factors quickly, then validate with layered backtests and return-series diagnostics, and finally assess tradability via turnover.

6.1 Layered Backtests

Layered backtests are among the most common tools for initial factor screening. The idea is simple: sort stocks by factor value, assign them to groups, and compare the groups’ subsequent returns.

If the factor is valid, group returns should usually show some monotonicity—e.g., high-factor groups earn more than low-factor groups, with intermediate groups ordered roughly in between. If group returns look unstructured, the factor’s validity is typically weak.

6.1.1 Building Groups

Grouping logic: at each month-end, sort the full cross-section by factor value and split into $N$ equal groups.

In practice, 10 groups (deciles) are common: enough to see monotonicity without making each group too thin. For higher-frequency factors, 5 groups often suffice. The core rule is that groups must be formed strictly from factor values available at the sort date—no look-ahead information.

Concrete steps:

1. Sort: at each month-end, rank all stocks by factor value from high to low
2. Group: split into $N$ equal groups; group 1 has the highest values, group $N$ the lowest
3. Compute returns: within each group, equal-weight or value-weight stocks and compute next-month returns
4. Roll: repeat monthly to obtain a time series of group returns

Note: if the factor distribution is uneven (many extremes), equal-count splits can leave groups with very different sizes. Prefer MAD winsorization before grouping.

Layered-backtest results are also easily contaminated by shared exposures such as size and industry. If one group is mostly small-caps and another mostly large-caps, return gaps may reflect the size effect rather than the factor’s own forecast power. Control relevant risk exposures whenever possible.

You can further run layered backtests on a size-neutralized factor: regress the factor on size, take residuals as the new factor, then form groups. This helps remove size interference and isolate the factor’s independent explanatory power.

If a group ends up with too few names—e.g., fewer than 10—its statistics are unreliable. In practice, aim for at least 20 stocks per group; if that is hard, merge adjacent groups to stabilize the test.

Group labels can be generated as follows:

import pandas as pd
import numpy as np

# Python pseudocode: decile grouping
def factor_quantile(factor_df, n_groups=10):
 """
 factor_df: factor-value matrix with dates as index and tickers as columns
 returns a matrix of group labels
 """
 # Sort and group each month
 groups = factor_df.apply(
 lambda x: pd.qcut(x.rank, q=n_groups, labels=False) + 1,
 axis=1
 )
 return groups

A complete layered-backtest skeleton looks like this:

def layered_backtest(factor_df, return_df, n_groups=10):
 """
 factor_df: factor values, index=date, columns=tickers
 return_df: next-period returns, index=date, columns=tickers
 """
 results = {}
 for date in factor_df.index:
 # Current-period factor values
 factors = factor_df.loc[date].dropna
 # Rank by factor and form groups
 ranked = factors.rank
 group_labels = pd.qcut(ranked, n_groups, labels=False)
 # Next-period return by group
 returns = return_df.loc[date]
 for g in range(n_groups):
 mask = group_labels == g
 group_returns = returns[mask.index[mask]].mean
 results.setdefault(g, []).append(group_returns)
 return pd.DataFrame(results)

  File "<ipython-input-2-ff26dcdeb80a>", line 1
    def分层回测(factor_df, return_df, n_groups=10):
                                              ^
SyntaxError: invalid syntax

6.1.2 Long–Short Portfolio Returns

A long–short portfolio goes long group 1 (highest factor values) and short group $N$ (lowest).

That return answers: if you select stocks purely on this factor, how much excess return can you earn?

Looking at long–short returns helps strip out overall market moves. For example, if the market rises 20%, group 1 rises 25%, and group 5 rises 23%, the factor’s incremental contribution is limited; a 2% long–short return better reflects the factor’s true contribution.

As a rule of thumb, factors with annualized long–short returns below 5% usually do not belong in a multi-factor model—after trading costs, little room remains.

When computing long–short returns, watch two details:

• Rebalance frequency: monthly rebalancing is most common; higher-frequency factors may need weekly rebalancing
• Trading costs: long–short turnover is often high—subtract roughly 0.1%–0.2% one-way cost

# Long–short portfolio return
long_short_ret = group_returns.iloc[:, 0] - group_returns.iloc[:, -1]
# Subtract trading costs (assume 0.15% one-way)
turnover = estimate_turnover(group_assignments)
long_short_ret_net = long_short_ret - turnover * 0.0015 * 2

6.1.3 Monotonicity Tests

For a valid factor, group returns should be monotonic: from group 1 to group $N$, returns should decline (or rise) gradually.

If group 3 earns more than group 2, the factor fails in some regions of the distribution.

Three common checks:

Method	Description	When to use
Return monotonicity	Whether mean returns by group are strictly monotonic	Initial screen; quickly drop obviously invalid factors
Monotonicity statistic	Share of consecutive group-return differences with the expected sign	Formal reports; quantify how monotonic the pattern is
Patton–Timmermann test	Statistical test of monotonic group returns	Academic papers; when rigor is required

Monotonicity is not the same as linearity. Factor returns can be convex (high at both ends) or concave. As long as the direction is consistent, the pattern can still be called monotonic.

Cases exist where the first five groups are monotonic but group 6 jumps up—later traced to a few anomalous names in group 6. Combine monotonicity checks with outlier analysis.

6.1.4 Group Net-Value Curves

Precise numbers still lose to a clear chart. Group net-value curves plot cumulative returns by group.

A good factor’s curves usually show:

• Group 1 on top, group $N$ at the bottom
• Little crossing among groups
• A steadily rising long–short net-value path

Crossing curves suggest periods when the factor failed. A volatile long–short path signals high factor risk.

When plotting net-value curves, log scales can make period-to-period return differences clearer, especially early versus late in the sample.

Practical issues to watch carefully:

• Survivorship bias: use the full sample, including delisted stocks. Surviving-only samples can overstate performance—weak names already exited.
• Look-ahead bias: factor values used for grouping must be available after the rebalance-day close. Future information invalidates the result.
• Industry concentration: if group 1 is dominated by one industry, the factor may be an industry effect rather than stock selection. Prefer industry neutralization.
• Small-cap effect: many factors work better in small caps. Stratify by size to confirm validity across size buckets.
• Size and industry controls: run layered backtests on size-neutralized factors, or control industry exposure when grouping, to isolate independent explanatory power.

Layered backtests are foundational in factor investing. Build groups carefully, net long–short returns of costs, test monotonicity statistically, and show net-value curves clearly. After these steps, you can usually judge whether the factor is worth keeping.

6.2 IC/IR Analysis

Layered backtests focus on group-return performance; IC/IR analysis evaluates predictive power via correlation. IC analysis can discard many invalid factors quickly and is usually the first line of defense in factor investing.

6.2.1 Definition and Computation of IC

IC (Information Coefficient) measures how well factor values predict next-period returns—i.e., do higher-factor stocks perform better next period?

Two mainstream IC definitions:

• Spearman rank correlation: usually preferred in practice. It depends only on ranks, not raw levels, and is less sensitive to extremes. In markets like A-shares that often see extreme moves, rank correlation is more robust than Pearson.
• Pearson correlation: assumes linearity and is sensitive to outliers. Avoid it unless the factor has been carefully normalized.

Steps to compute IC:

1. On each rebalance date, take factor values for all stocks
2. Compute those stocks’ next-period returns (e.g., next day or next week)
3. Compute the correlation between factor values and next-period returns

Implementation:

from scipy.stats import spearmanr

def calc_ic(factor_series, forward_return_series):
 """
 Cross-sectional IC
 factor_series: factor values on one day, index=tickers
 forward_return_series: matching next-period returns
 """
 # Drop missing values
 valid = factor_series.notna & forward_return_series.notna
 f = factor_series[valid]
 r = forward_return_series[valid]

 # Spearman rank correlation
 ic, p_value = spearmanr(f, r)
 return ic, p_value

# Example: factor vs next-period return on one trading day
# ic_value, p_val = calc_ic(factor_data['2024-01-05'],
# forward_ret['2024-01-05'])
# print(f"IC = {ic_value:.4f}, p-value = {p_val:.4f}")

A more compact single-period interface:

def calculate_ic(factor_series, return_series):
 """Single-period IC (Spearman rank correlation)"""
 combined = pd.concat([factor_series, return_series], axis=1).dropna
 if len(combined) < 30:
 return np.nan
 ic, p_value = spearmanr(combined.iloc[:, 0], combined.iloc[:, 1])
 return ic

  File "<ipython-input-4-dd7623c68873>", line 5
    return np.nan
    ^
IndentationError: expected an indented block after 'if' statement on line 4

Key point: IC lies in $[-1, 1]$. Positive IC means the factor covaries positively with future returns; negative IC means the opposite. Larger absolute values imply stronger forecast power.

6.2.2 Statistical Tests on the IC Series

A single-day IC is not persuasive: if today’s IC is 0.05 and tomorrow’s is −0.03, you cannot declare the factor valid. You need the statistical significance of the IC series.

Typical analyses:

• Mean IC: is the long-run average IC significantly different from zero?
• IC t-statistic: test whether the mean is statistically significant
• Positive IC share: what fraction of IC observations are positive? Near 50% is little better than a coin flip

Note: the IC series is autocorrelated. Ordinary t-tests overstate significance. Prefer Newey–West–adjusted standard errors.

from statsmodels.stats.stattools import durbin_watson
from statsmodels.regression.linear_model import OLS

def ic_statistical_test(ic_series):
 """
 Statistical tests on an IC series
 """
 # Basic statistics
 mean_ic = ic_series.mean
 std_ic = ic_series.std
 t_stat = mean_ic / (std_ic / np.sqrt(len(ic_series)))

 # Newey–West adjustment (for autocorrelation)
 X = np.ones((len(ic_series), 1))
 model = OLS(ic_series.values, X).fit(cov_type='HAC', cov_kwds={'maxlags': 5})
 nw_t_stat = model.tvalues[0]

 # Positive share
 positive_ratio = (ic_series > 0).mean

 return {
 'mean_ic': mean_ic,
 'std_ic': std_ic,
 't_stat': t_stat,
 'nw_t_stat': nw_t_stat,
 'positive_ratio': positive_ratio
 }

A common bar is mean IC at least above 0.03 and a t-statistic above 2 (about 95% confidence). A mean IC of only 0.01 may be statistically significant yet still fail to cover trading costs in practice.

Reference ranges:

Metric	Excellent	Good	Fair	Weak
Mean IC	>0.05	0.03–0.05	0.01–0.03	<0.01
IR	>0.5	0.3–0.5	0.1–0.3	<0.1
Positive IC share	>60%	55%–60%	50%–55%	<50%

A common mistake is judging factors by mean IC alone. If mean IC is 0.04 but the standard deviation is 0.15, IR is only 0.27—the signal is unstable. When assessing validity, put more weight on IR (or ICIR) than on mean IC alone.

6.2.3 IC Decay Analysis

IC decay is how quickly predictive power fades with horizon. Causes include changing market structure, more arbitrage, and crowding that dilutes alpha.

Examine IC by holding horizon:

• 1-day IC (daily)
• 5-day IC (weekly)
• 20-day IC (monthly)
• 60-day IC (quarterly)

If the factor works only on day 1 and decays to zero by day 5, it is a short-horizon reversal factor. If IC stays stable out to 20 days, it behaves more like a trend factor.

For example, a momentum factor may show monthly mean IC of 0.05 with t = 3.2—apparently useful—yet decay analysis may show explanatory power only for the next month, with IC turning negative in month 2. That factor suits short holding periods, not long-horizon allocation.

def decay_analysis(factor_df, forward_returns_dict):
 """
 Analyze IC decay
 factor_df: daily factor values, index=date, columns=tickers
 forward_returns_dict: returns by horizon, e.g. {'1d': ..., '5d': ...}
 """
 decay_results = {}
 for horizon, ret_df in forward_returns_dict.items:
 ics = []
 for date in factor_df.index:
 if date in ret_df.index:
 ic, _ = calc_ic(factor_df.loc[date], ret_df.loc[date])
 ics.append(ic)
 decay_results[horizon] = np.mean(ics)

 return decay_results

# Example output
# {'1d': 0.042, '5d': 0.038, '20d': 0.025, '60d': 0.008}

  File "<ipython-input-6-08a0b285071e>", line 9
    ics = []
    ^
IndentationError: expected an indented block after 'for' statement on line 8

Cases exist with daily IC as high as 0.08 but 5-day IC down to 0.01. Such fast-decaying factors usually fit only high-frequency strategies and struggle to cover trading costs at daily or lower frequency.

6.2.4 ICIR: Information Coefficient and Information Ratio

ICIR (Information Coefficient Information Ratio) is mean IC divided by the standard deviation of IC. It measures the stability of predictive power. The text often shortens this to IR.

Formula:

$$ ICIR = mean(IC) / std(IC) $$

Why ICIR matters:

• Factor A: mean IC 0.05, std 0.02 → ICIR = 2.5
• Factor B: mean IC 0.08, std 0.10 → ICIR = 0.8

Factor A has a lower IC but higher stability; B has a higher IC but more volatility. In live trading, A often beats B—you can assign it more weight without sudden monthly failures.

Common ICIR screens:

ICIR range	Rating	Suggestion
> 2.0	Excellent	Usable directly
1.0 – 2.0	Good	Combine with other factors
0.5 – 1.0	Fair	Use cautiously; watch risk
< 0.5	Poor	Prefer to drop

ICIR often matters more than IC itself. A factor with ICIR 2.0 and mean IC 0.03 usually beats one with ICIR 0.5 and mean IC 0.08. Stability is the key metric.

Note: the table above is an ICIR standard based on IC-series stability; earlier IR thresholds (e.g., excellent > 0.5) often refer to monthly-frequency experience. Units and sample frequencies differ—compare metrics at a common frequency in practice.

Common practical issues:

1. Positive IC but losing strategy: the factor may fail in extreme regimes. Inspect the IC time series around bull/bear transitions for sign flips.
2. High ICIR but poor live performance: possible overfitting. In-sample ICIR of 3.0 with out-of-sample 0.5 often signals data mining.
3. Very fast IC decay: the factor may be crowded, or costs may be too high. Check turnover; above 500% annualized, fast decay is unsurprising.

When developing new factors, run IC analysis before layered backtests. IC analysis is the first line of defense in factor investing.

6.3 Factor-Return Series Analysis

Layered backtests and IC analysis mainly evaluate factors in the cross-section; factor-return series analysis adds a time-series lens.

Concretely, form a factor-mimicking portfolio: each month long the highest-factor portfolio and short the lowest, then study that long–short return path (construction matches Section 6.1.2; here the focus is time-series statistics and risk features).

Use the series to ask:

• Does the cumulative-return path keep rising?
• How large is max drawdown, and when does it occur?
• Is return volatility acceptable?
• Are there clear style exposures or regime failures?

For example, a value long–short book may earn 8% annualized with a 25% max drawdown concentrated in late-2015 bull conditions and the 2017 blue-chip rally—evidence of style-dependent failure. Acceptable for broad-market selection; riskier if the strategy covers only one style sleeve.

Key statistics include annualized return, annualized volatility, Sharpe ratio, max drawdown, Calmar ratio (annualized return / max drawdown), and monthly hit rate.

6.4 Turnover Analysis

Turnover analysis is easy to underweight. Even a high-return factor can be hollowed out by trading costs if turnover is extreme—sometimes to the point of impracticality.

Turnover is defined as:

$$ \text{Turnover} = (\text{buys} + \text{sells}) / \text{portfolio market value} $$

For long–short books, measure long- and short-side turnover separately. Focus on:

• One-way turnover: long or short turnover per period
• Two-way turnover: long plus short combined
• Annualized turnover: one-way turnover × rebalance frequency
• Return-per-turnover: annualized return / annualized turnover—return earned per unit turnover

Evaluate turnover against real costs. In A-shares, two-way costs (commission, stamp tax, impact) are often about 0.1%–0.3%. High turnover can erase even a strong IC.

In one case, a high-frequency reversal factor had ~800% raw turnover and 12% gross annualized return; after costs, net return fell to 3%. Adding screens cut turnover to 200%: gross return fell to 9%, but net rose to 6%. Lowering turnover can raise realized return.

In practice, set a turnover cap. For monthly factors, annualized turnover above 300% warrants caution; for weekly factors, above 1000% is usually infeasible. Exact thresholds depend on costs and capacity.

Code example for turnover:

def calculate_turnover(weights_prev, weights_curr):
 """
 Single-period turnover
 weights_prev: prior holdings weights
 weights_curr: current holdings weights
 """
 # Buys: weight increases
 buy = (weights_curr - weights_prev).clip(lower=0).sum
 # Sells: weight decreases
 sell = (weights_prev - weights_curr).clip(lower=0).sum
 # Turnover as average of buys and sells
 turnover = (buy + sell) / 2
 return turnover

Finally, a brief summary of the testing framework. It has four modules with different emphases:

• Layered backtests: check monotonicity of group returns and stable cross-sectional explanatory power (including long–short returns, monotonicity tests, net-value curves, and bias controls);
• IC/IR analysis: assess predictive power and its stability (including statistical tests, decay analysis, and ICIR);
• Factor-return series: examine realized return paths and risk features;
• Turnover analysis: measure trading costs and implementability.

Recommended order: IC first to discard weak factors, then layered backtests and return-series checks, then turnover for tradability. Only after all four dimensions pass should a factor be treated as provisionally valid and moved into multi-factor combination.

Factor testing is not one-and-done. As markets change, validity can drift. Revisit the framework periodically. Widespread out-of-sample failure is a normal research risk.

6.5 Formal Portfolio Sorts: Univariate and Bivariate

Layered backtests and IC analysis are intuitive, but academic replication and strict tests also require clarity on breakpoints, controls for another dimension (e.g., size), and how long–short returns are defined. Formal portfolio sorts are the standard workflow in empirical asset pricing.

6.5.1 Univariate Portfolio Sorts

Basic steps

1. At each rebalance date (e.g., month-end), sort stocks by the factor;
2. Split into $Q$ groups using quantile breakpoints (often 5 or 10);
3. Compute holding-period returns within groups, equal- or value-weighted;
4. Form a high-minus-low (or low-minus-high) long–short portfolio to get a factor-return series;
5. Test the long–short mean with t-statistics (often Newey–West–adjusted).

This matches earlier layered backtests in spirit (“group—hold—compare”); academic write-ups emphasize breakpoint samples (whether breakpoints use only a subset of exchanges), lag alignment (avoid look-ahead), and the sign convention of the long–short portfolio.

Implementation notes

A complete sort usually includes breakpoint computation, portfolio assignment, within-group return aggregation, and high-minus-low construction. Modularize these steps in your own code.

For A-shares, panel fields typically include ticker, date, return, market cap, and factor value, with filters for ST, suspensions, limit-up/down days, and days since listing.

# Demo: minimal runnable univariate-sort skeleton (synthetic data)
import numpy as np
import pandas as pd

rng = np.random.default_rng(42)
dates = pd.date_range("2020-01-31", periods=24, freq="ME")
rows = []
for dt in dates:
 n = 300
 factor = rng.normal(size=n)
 ret = 0.02 * factor + rng.normal(0, 0.05, size=n)
 mktcap = np.exp(rng.normal(22, 1, size=n))
 rows.append(pd.DataFrame({
 "date": dt, "id": np.arange(n), "factor": factor, "ret": ret, "mktcap": mktcap,
 }))
panel = pd.concat(rows, ignore_index=True)

def univariate_sort(df, q=5, weight="equal"):
 out = []
 for dt, g in df.groupby("date"):
 g = g.copy
 g["port"] = pd.qcut(g["factor"], q=q, labels=False, duplicates="drop") + 1
 if weight == "equal":
 port_ret = g.groupby("port")["ret"].mean
 else:
 def wavg(x):
 w = x["mktcap"] / x["mktcap"].sum
 return np.sum(w * x["ret"])
 port_ret = g.groupby("port").apply(wavg, include_groups=False)
 long_short = port_ret.loc[q] - port_ret.loc[1]
 out.append({"date": dt, "long_short": long_short})
 return pd.DataFrame(out).set_index("date")

res = univariate_sort(panel, q=5)
print(res[["long_short"]].mean)
print(
 "t-stat (iid approx):",
 res["long_short"].mean / (res["long_short"].std(ddof=1) / np.sqrt(len(res))),
)

  File "<ipython-input-8-a86b468f3358>", line 21
    g = g.copy()
    ^
IndentationError: expected an indented block after 'for' statement on line 20

6.5.2 Bivariate Portfolio Sorts

Univariate sorts cannot answer whether factor returns merely proxy for size, industry, or another style. A common bivariate approach:

1. First group by the control variable (e.g., size);
2. Then sort by the target factor within each group—or form a 2D grid via independent sorts;
3. Average across the control dimension to obtain “controlled” factor-portfolio returns.

Fama–French Size–BM 2×3 or 5×5 portfolios are classic examples.

Practical implication: if the target factor’s long–short return remains significant after size-neutral bivariate sorts, its incremental information is more likely independent of the size effect.

6.6 Fama–MacBeth Regressions

Fama–MacBeth (1973) regressions estimate risk premia on characteristics/exposures and complement portfolio sorts: sorts speak to economic magnitude and monotonicity; FM speaks to statistical pricing and multivariate robustness.

Step 1 (cross-section): for each period $t$, estimate

$$ R_{i,t}= \gamma_{0,t}+\gamma_{1,t}X_{i,t-1}+\varepsilon_{i,t} $$

where $X_{i,t-1}$ is a lagged factor exposure or firm characteristic.

Step 2 (time series): average $\{\gamma_{1,t}\}$ to estimate the risk premium, and use time-series standard errors (often Newey–West) to test whether it differs from 0.

Method	Main question	Strengths	Limits
Portfolio sorts	Do high-exposure portfolios earn more?	Intuitive; can show nonlinearity	Hard to control many variables at once
Fama–MacBeth	After controls, is the exposure priced?	Multivariate controls	Relies on linearity; sensitive to outliers

# Demo: handwritten two-step Fama–MacBeth (reuses panel from previous cell)
import statsmodels.api as sm

def fama_macbeth(panel, y="ret", x_cols=("factor",), date_col="date"):
 gammas = []
 for dt, g in panel.groupby(date_col):
 if len(g) < 30:
 continue
 X = sm.add_constant(g.loc[:, list(x_cols)])
 model = sm.OLS(g[y], X).fit
 gammas.append({"date": dt, **model.params.to_dict})
 gdf = pd.DataFrame(gammas).set_index("date").sort_index
 summary = pd.DataFrame({
 "risk_premium": gdf.mean,
 "std_error_iid": gdf.std(ddof=1) / np.sqrt(len(gdf)),
 })
 summary["t_iid"] = summary["risk_premium"] / summary["std_error_iid"]
 return gdf, summary

gdf, summary = fama_macbeth(panel, y="ret", x_cols=("factor",))
print(summary)

  File "<ipython-input-9-9068e354d5e3>", line 7
    if len(g) < 30:
    ^
IndentationError: expected an indented block after 'for' statement on line 6

Multivariate Controls and Standard Errors

In practice, put several characteristics in the same cross-sectional regression—e.g., market beta, book-to-market, and log market cap—to isolate the target factor’s incremental pricing. For the second-step coefficient series, prefer Newey–West or other HAC standard errors.

Multivariate sketch:

# Cross-section at t: ret_excess ~ factor + beta + bm + log_mktcap
# Then time-series means of coefficients with Newey–West t-tests

Implementation Checklist

1. Lags: regressors must be strictly predetermined relative to returns—avoid look-ahead;
2. Standardization / winsorization: cross-sectional regressions are outlier-sensitive—clean and standardize first (Chapter 5);
3. Industry or size controls: add controls in the regression, or neutralize first then regress (Chapter 5);
4. Standard errors: coefficient series are often autocorrelated—prefer Newey–West t-stats;
5. Relation to IC: Rank IC measures order correlation; FM estimates conditional risk premia—they are not substitutes.

6.7 Clustered Standard Errors in Panel Regressions (Essentials)

When stacking many cross-sections into a panel, or estimating stock-level time-series regressions, residuals often correlate within stocks or within dates. Ordinary OLS standard errors then overstate significance.

In panel or pooled regressions, residual dependence shapes standard errors. Common fixes:

• Cluster by firm: allow within-stock residual correlation over time;
• Cluster by time: allow within-date residual correlation across stocks;
• Two-way clustering: allow both (common in empirical asset pricing).

Relative to FM: classic FM already partly handles cross-sectional dependence by inferring from the time series of cross-sectional coefficients; a one-shot panel regression needs clustering made explicit. In Python, libraries such as linearmodels support fixed effects and clustered SEs.

Course suggestion: prefer FM + Newey–West for factor-pricing tests; for panel characteristic regressions report at least one-way clustered SEs, and check two-way clustering for important conclusions. Multiple testing and p-hacking in sort designs appear in Chapter 4 (expanding the factor map).

6.8 A-Share Replication: Univariate Portfolio Sorts

This section fully replicates the tidyfinance chapter Univariate Portfolio Sorts, replacing the data with an A-share monthly sample. Portfolio sorts are among the most common tools in empirical asset pricing: sort stocks by an observable characteristic into portfolios and compare subsequent returns to test cross-sectional predictability.

A univariate sort uses a single sorting variable $x_{i,t-1}$ (the $t-1$ subscript means the signal is known at the start of period $t$). Here we use the market beta estimated in Section 2.6 as the sorting variable and study its link to next-period excess returns $r_{i,t}$.

import warnings
from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import statsmodels.api as sm
from IPython.display import display
from matplotlib import font_manager

warnings.filterwarnings("ignore", category=UserWarning)

# CJK-capable fonts (same plotting setup as Chapter 2)
_cn_candidates = ["Microsoft YaHei", "SimHei", "PingFang SC", "Noto Sans CJK SC", "Arial Unicode MS"]
_available = {f.name for f in font_manager.fontManager.ttflist}
for _name in _cn_candidates:
    if _name in _available:
        plt.rcParams["font.sans-serif"] = [_name]
        break
plt.rcParams["axes.unicode_minus"] = False

DATA_DIR = Path(r"D:/A_Topics/202607_03_TidyFinanceAShare/202608_02_传统量化/Data")
CACHE_DIR = DATA_DIR / "_cache_beta"

print("DATA_DIR =", DATA_DIR)
print("beta exists =", (DATA_DIR / "beta_ashare.parquet").exists)
print("monthly cache exists =", (CACHE_DIR / "crsp_monthly_ashare.parquet").exists)

DATA_DIR = D:\A_Topics\202607_03_TidyFinanceAShare\202608_02_传统量化\Data
beta exists = True
monthly cache exists = True

6.8.1 Data Preparation

Corresponds to Data Preparation in the original. The asset universe uses A-share monthly stock excess returns and lagged free-float market cap; market excess returns are used for CAPM alphas. beta comes from the monthly rolling estimates in Section 2.6.

crsp_monthly = (
    pd.read_parquet(CACHE_DIR / "crsp_monthly_ashare.parquet")
    [["permno", "date", "ret_excess", "mkt_excess", "mktcap_lag"]]
    .copy
)
crsp_monthly["date"] = pd.to_datetime(crsp_monthly["date"])
crsp_monthly["permno"] = crsp_monthly["permno"].astype(str).str.zfill(6)

# Market factor: one mkt_excess per month (same role as factors_ff3_monthly in the original)
factors_mkt_monthly = (
    crsp_monthly[["date", "mkt_excess"]]
    .drop_duplicates("date")
    .sort_values("date")
    .reset_index(drop=True)
)

beta = (
    pd.read_parquet(DATA_DIR / "beta_ashare.parquet")
    .query("return_type == 'monthly'")
    [["permno", "date", "beta"]]
    .copy
)
beta["date"] = pd.to_datetime(beta["date"])
beta["permno"] = beta["permno"].astype(str).str.zfill(6)

print("crsp_monthly:", crsp_monthly.shape)
print("beta:", beta.shape)
print("months:", factors_mkt_monthly["date"].nunique)
display(crsp_monthly.head(3))
display(beta.head(3))

crsp_monthly: (890537, 5)
beta: (623256, 3)
months: 319

	permno	date	ret_excess	mkt_excess	mktcap_lag
0	000001	2000-01-01	0.060035	0.158982	NaN
1	000001	2000-02-01	-0.013189	0.120168	19843822.88
2	000001	2000-03-01	0.000873	0.054070	19618933.36

	permno	date	beta
0	000001	2003-12-01	0.945177
1	000001	2004-01-01	0.956954
2	000001	2004-02-01	0.978056

6.8.2 Sorting by Market Beta

Corresponds to Sorting by Market Beta. Use a one-period-lagged beta as the sorting variable so the information is known at portfolio formation. Shift beta dates forward by one month and merge to returns on (permno, date), rather than a naive groupby.shift(1)—the latter misaligns when months are implicitly missing.

beta_lag = beta.copy
beta_lag["date"] = beta_lag["date"] + pd.DateOffset(months=1)
beta_lag = beta_lag.rename(columns={"beta": "beta_lag"}).dropna(subset=["beta_lag"])

data_for_sorts = crsp_monthly.merge(
    beta_lag, on=["permno", "date"], how="inner"
).dropna(subset=["ret_excess", "mktcap_lag", "beta_lag"])

print("data_for_sorts:", data_for_sorts.shape)
print(
    "date range:",
    data_for_sorts["date"].min.date,
    "->",
    data_for_sorts["date"].max.date,
)
display(data_for_sorts.head(3))

data_for_sorts: (614523, 6)
date range: 2004-01-01 -> 2026-07-01

	permno	date	ret_excess	mkt_excess	mktcap_lag	beta_lag
0	000001	2004-01-01	0.088847	0.073093	11993670.32	0.945177
1	000001	2004-02-01	0.114744	0.069838	13078879.04	0.956954
2	000001	2004-03-01	0.027323	0.030699	14600989.96	0.978056

The first step in a portfolio sort is computing breakpoints. Split stocks into low / high using the median of lagged beta, then compute each group’s value-weighted excess return (weights = mktcap_lag).

def value_weighted_ret(g: pd.DataFrame) -> float:
    w = g["mktcap_lag"]
    return float((g["ret_excess"] * w).sum / w.sum)

rows = []
for dt, g in data_for_sorts.groupby("date", sort=True):
    g = g.copy
    # Median breakpoint: two groups
    g["portfolio"] = pd.qcut(
        g["beta_lag"].rank(method="first"),
        q=2,
        labels=["low", "high"],
    )
    for port, pg in g.groupby("portfolio", observed=True):
        if pg["mktcap_lag"].sum <= 0:
            continue
        rows.append(
            {
                "date": dt,
                "portfolio": str(port),
                "ret": value_weighted_ret(pg),
            }
        )

beta_portfolios_2 = pd.DataFrame(rows).sort_values(["date", "portfolio"])
print(beta_portfolios_2.head)
print("n months:", beta_portfolios_2["date"].nunique)

        date portfolio       ret
1 2004-01-01      high  0.099191
0 2004-01-01       low  0.056309
3 2004-02-01      high  0.087797
2 2004-02-01       low  0.068211
5 2004-03-01      high  0.035140
n months: 271

6.8.3 Performance Evaluation (Two-Way Long–Short)

Corresponds to Performance Evaluation. Form a long–short book: long high beta, short low beta. Ignoring frictions, net market exposure of the long–short can be near zero (here we first look at raw long–short returns).

beta_longshort_2 = (
    beta_portfolios_2.pivot(index="date", columns="portfolio", values="ret")
    .sort_index
    .dropna(subset=["low", "high"])
)
beta_longshort_2["long_short"] = beta_longshort_2["high"] - beta_longshort_2["low"]
beta_longshort_2 = beta_longshort_2.reset_index

display(beta_longshort_2.head)
print(
    "mean long_short =",
    round(beta_longshort_2["long_short"].mean, 6),
    "std =",
    round(beta_longshort_2["long_short"].std, 6),
)

portfolio	date	high	low	long_short
0	2004-01-01	0.099191	0.056309	0.042882
1	2004-02-01	0.087797	0.068211	0.019586
2	2004-03-01	0.035140	0.025616	0.009524
3	2004-04-01	-0.103090	-0.114164	0.011074
4	2004-05-01	-0.014376	-0.028354	0.013978

mean long_short = -0.001431 std = 0.041392

Test whether mean long–short excess returns differ from zero. Asset-pricing papers often use Newey–West (HAC) $t$-statistics for autocorrelation; the original defaults to a 6-month lag. Here we use statsmodels with cov_type="HAC".

y = beta_longshort_2["long_short"]
X = np.ones((len(y), 1))
model_fit = sm.OLS(y, X, missing="drop").fit(
    cov_type="HAC", cov_kwds={"maxlags": 6}
)
print("Long-short ~ 1  (Newey-West, lags=6)")
print(model_fit.summary)

Long-short ~ 1  (Newey-West, lags=6)
                            OLS Regression Results                            
==============================================================================
Dep. Variable:             long_short   R-squared:                       0.000
Model:                            OLS   Adj. R-squared:                  0.000
Method:                 Least Squares   F-statistic:                       nan
Date:                Sat, 08 Aug 2026   Prob (F-statistic):                nan
Time:                        21:15:11   Log-Likelihood:                 479.01
No. Observations:                 271   AIC:                            -956.0
Df Residuals:                     270   BIC:                            -952.4
Df Model:                           0                                         
Covariance Type:                  HAC                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
const         -0.0014      0.002     -0.722      0.470      -0.005       0.002
==============================================================================
Omnibus:                       37.393   Durbin-Watson:                   1.942
Prob(Omnibus):                  0.000   Jarque-Bera (JB):              187.790
Skew:                          -0.363   Prob(JB):                     1.67e-41
Kurtosis:                       7.013   Cond. No.                         1.00
==============================================================================

Notes:
[1] Standard Errors are heteroscedasticity and autocorrelation robust (HAC) using 6 lags and without small sample correction

If CAPM holds, high-beta portfolios should have higher expected returns, so “long high / short low” mean excess returns should be significantly positive. Insignificant or negative results conflict with that CAPM intuition—this is the low-beta / betting-against-beta thread developed below.

6.8.4 Functional Programming for Portfolio Sorts

Corresponds to Functional Programming for Portfolio Sorts. Wrap the sort logic in functions so stocks can be assigned to any $N$ groups. Be careful when clustering is severe, breakpoints repeat, or early samples have few constituents (empty portfolios, extreme weights).

def assign_portfolio(series: pd.Series, n_portfolios: int) -> pd.Series:
    # Map the sorting variable to labels 1..n_portfolios by quantile
    return pd.qcut(
        series.rank(method="first"),
        q=n_portfolios,
        labels=[str(i) for i in range(1, n_portfolios + 1)],
        duplicates="drop",
    )

def form_value_weighted_portfolios(
    data: pd.DataFrame,
    sorting_variable: str = "beta_lag",
    n_portfolios: int = 10,
) -> pd.DataFrame:
    # Each month, split into n groups by sorting_variable; value-weighted excess returns
    rows = []
    for dt, g in data.groupby("date", sort=True):
        g = g.copy
        g["portfolio"] = assign_portfolio(g[sorting_variable], n_portfolios)
        for port, pg in g.groupby("portfolio", observed=True):
            if pg["mktcap_lag"].sum <= 0:
                continue
            rows.append(
                {
                    "date": dt,
                    "portfolio": str(port),
                    "ret": value_weighted_ret(pg),
                }
            )
    out = pd.DataFrame(rows)
    # Ordered categories for plotting in 1..10 order
    cats = [str(i) for i in range(1, n_portfolios + 1)]
    out["portfolio"] = pd.Categorical(out["portfolio"], categories=cats, ordered=True)
    out = out.merge(factors_mkt_monthly, on="date", how="left")
    return out.sort_values(["date", "portfolio"]).reset_index(drop=True)

beta_portfolios = form_value_weighted_portfolios(
    data_for_sorts, sorting_variable="beta_lag", n_portfolios=10
)
print(beta_portfolios.head)
print("portfolios:", sorted(beta_portfolios["portfolio"].dropna.unique.tolist))

        date portfolio       ret  mkt_excess
0 2004-01-01         1  0.024797    0.073093
1 2004-01-01         2  0.066696    0.073093
2 2004-01-01         3  0.054822    0.073093
3 2004-01-01         4  0.070948    0.073093
4 2004-01-01         5  0.082065    0.073093
portfolios: ['1', '10', '2', '3', '4', '5', '6', '7', '8', '9']

6.8.5 Decile Performance: Alpha, Portfolio Beta, and Average Returns

Corresponds to More Performance Evaluation. For each beta portfolio, estimate CAPM:

$$ r_{p,t}=\alpha_p+\beta_p\,r_{m,t}+\varepsilon_{p,t}, $$

and summarize $\hat\alpha_p$, $\hat\beta_p$, and average excess returns.

summary_rows = []
for port, g in beta_portfolios.groupby("portfolio", observed=True):
    g = g.dropna(subset=["ret", "mkt_excess"])
    if len(g) < 24:
        continue
    X = sm.add_constant(g["mkt_excess"])
    m = sm.OLS(g["ret"], X).fit
    summary_rows.append(
        {
            "portfolio": str(port),
            "alpha": float(m.params["const"]),
            "beta": float(m.params["mkt_excess"]),
            "ret": float(g["ret"].mean),
            "n": int(len(g)),
        }
    )

beta_portfolios_summary = pd.DataFrame(summary_rows)
cats = [str(i) for i in range(1, 11)]
beta_portfolios_summary["portfolio"] = pd.Categorical(
    beta_portfolios_summary["portfolio"], categories=cats, ordered=True
)
beta_portfolios_summary = beta_portfolios_summary.sort_values("portfolio").reset_index(drop=True)
display(beta_portfolios_summary)

	portfolio	alpha	beta	ret	n
0	1	0.003172	0.715856	0.009123	271
1	2	0.001743	0.922710	0.009414	271
2	3	0.000964	1.012170	0.009378	271
3	4	0.000680	1.041719	0.009340	271
4	5	0.000177	1.070959	0.009080	271
5	6	0.000239	1.095577	0.009346	271
6	7	-0.000868	1.146563	0.008663	271
7	8	-0.000416	1.182178	0.009411	271
8	9	-0.001407	1.215726	0.008699	271
9	10	-0.006544	1.278707	0.004086	271

The figure below shows CAPM alphas of the ten beta portfolios. A pattern of “low-beta alphas high, high-beta alphas low” contradicts CAPM (risk-adjusted alphas near 0 and returns rising with beta).

fig, ax = plt.subplots(figsize=(9, 4.5))
ax.bar(
    beta_portfolios_summary["portfolio"].astype(str),
    beta_portfolios_summary["alpha"],
    color="#4C78A8",
    edgecolor="white",
)
ax.axhline(0, color="black", linewidth=0.8)
ax.set_xlabel("Portfolio (1=low beta … 10=high beta)")
ax.set_ylabel("CAPM alpha (monthly)")
ax.set_title("CAPM alphas of beta-sorted portfolios (A-share)")
ax.yaxis.set_major_formatter(plt.FuncFormatter(lambda x, _: f"{100*x:.2f}%"))
plt.tight_layout
plt.show

6.8.6 The Security Market Line (SML) and Beta Portfolios

Corresponds to The Security Market Line and Beta Portfolios. CAPM predicts portfolios lie on the SML with slope equal to the market risk premium. The figure shows:

1. The theoretical SML through the origin with slope equal to average market excess return (solid);
2. An empirical fit through the ten portfolios’ (portfolio beta, average return) points (dashed).

# Theoretical SML: through origin, slope = sample mean market excess return
mkt_premium = float(factors_mkt_monthly["mkt_excess"].mean)

# Empirical line: ret ~ a + b * beta
X = sm.add_constant(beta_portfolios_summary["beta"])
sml_fit = sm.OLS(beta_portfolios_summary["ret"], X).fit
intercept_hat = float(sml_fit.params["const"])
slope_hat = float(sml_fit.params["beta"])

fig, ax = plt.subplots(figsize=(7.5, 5.5))
ax.scatter(
    beta_portfolios_summary["beta"],
    beta_portfolios_summary["ret"],
    c=np.arange(len(beta_portfolios_summary)),
    cmap="viridis",
    s=60,
    zorder=3,
)
for _, r in beta_portfolios_summary.iterrows:
    ax.annotate(
        str(r["portfolio"]),
        (r["beta"], r["ret"]),
        textcoords="offset points",
        xytext=(5, 4),
        fontsize=9,
    )

x_line = np.linspace(0, max(2.0, beta_portfolios_summary["beta"].max * 1.05), 100)
ax.plot(x_line, mkt_premium * x_line, color="black", linewidth=1.5, label="Theoretical SML")
ax.plot(
    x_line,
    intercept_hat + slope_hat * x_line,
    color="crimson",
    linestyle="--",
    linewidth=1.5,
    label="Fitted line on portfolios",
)
ax.set_xlim(0, max(2.0, beta_portfolios_summary["beta"].max * 1.05))
ymax = max(mkt_premium * 2, beta_portfolios_summary["ret"].max * 1.2, 0.02)
ax.set_ylim(min(0, beta_portfolios_summary["ret"].min * 1.2), ymax)
ax.set_xlabel("Beta")
ax.set_ylabel("Average excess return")
ax.set_title("Average portfolio excess returns and beta estimates (A-share)")
ax.yaxis.set_major_formatter(plt.FuncFormatter(lambda x, _: f"{100*x:.2f}%"))
ax.legend(frameon=False)
ax.grid(True, alpha=0.3)
plt.tight_layout
plt.show

print(f"market premium (mean mkt_excess) = {mkt_premium:.4%}")
print(f"fitted intercept = {intercept_hat:.4%}, slope = {slope_hat:.4f}")

market premium (mean mkt_excess) = 0.6856%
fitted intercept = 1.4026%, slope = -0.0050

Next form the extreme long–short within deciles: long the highest-beta portfolio, short the lowest, and run Newey–West mean tests plus CAPM alpha tests.

wide = (
    beta_portfolios.pivot(index="date", columns="portfolio", values="ret")
    .sort_index
)
# Portfolio labels '1'..'10'
low_col, high_col = "1", "10"
beta_longshort = wide[[low_col, high_col]].dropna.copy
beta_longshort = beta_longshort.rename(columns={low_col: "low", high_col: "high"})
beta_longshort["long_short"] = beta_longshort["high"] - beta_longshort["low"]
beta_longshort = beta_longshort.reset_index.merge(
    factors_mkt_monthly, on="date", how="left"
)

print("long-short months:", len(beta_longshort))
display(beta_longshort.head)

# Is mean return zero?
y = beta_longshort["long_short"]
model_mu = sm.OLS(y, np.ones((len(y), 1))).fit(
    cov_type="HAC", cov_kwds={"maxlags": 6}
)
print("\nLong-short ~ 1  (NW lags=6)")
print(model_mu.summary)

# CAPM alpha
X = sm.add_constant(beta_longshort["mkt_excess"])
model_capm = sm.OLS(beta_longshort["long_short"], X, missing="drop").fit(
    cov_type="HAC", cov_kwds={"maxlags": 6}
)
print("\nLong-short ~ 1 + mkt_excess  (NW lags=6)")
print(model_capm.summary)

long-short months: 271

	date	low	high	long_short	mkt_excess
0	2004-01-01	0.024797	0.138344	0.113547	0.073093
1	2004-02-01	0.041176	0.122426	0.081249	0.069838
2	2004-03-01	0.008092	0.012647	0.004555	0.030699
3	2004-04-01	-0.144966	-0.124091	0.020875	-0.100944
4	2004-05-01	-0.036202	-0.007386	0.028816	-0.024572

Long-short ~ 1  (NW lags=6)
                            OLS Regression Results                            
==============================================================================
Dep. Variable:             long_short   R-squared:                       0.000
Model:                            OLS   Adj. R-squared:                  0.000
Method:                 Least Squares   F-statistic:                       nan
Date:                Sat, 08 Aug 2026   Prob (F-statistic):                nan
Time:                        21:15:13   Log-Likelihood:                 324.30
No. Observations:                 271   AIC:                            -646.6
Df Residuals:                     270   BIC:                            -643.0
Df Model:                           0                                         
Covariance Type:                  HAC                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
const         -0.0050      0.004     -1.419      0.156      -0.012       0.002
==============================================================================
Omnibus:                       29.855   Durbin-Watson:                   1.980
Prob(Omnibus):                  0.000   Jarque-Bera (JB):              152.214
Skew:                          -0.128   Prob(JB):                     8.85e-34
Kurtosis:                       6.663   Cond. No.                         1.00
==============================================================================

Notes:
[1] Standard Errors are heteroscedasticity and autocorrelation robust (HAC) using 6 lags and without small sample correction

Long-short ~ 1 + mkt_excess  (NW lags=6)
                            OLS Regression Results                            
==============================================================================
Dep. Variable:             long_short   R-squared:                       0.346
Model:                            OLS   Adj. R-squared:                  0.343
Method:                 Least Squares   F-statistic:                     63.87
Date:                Sat, 08 Aug 2026   Prob (F-statistic):           3.93e-14
Time:                        21:15:13   Log-Likelihood:                 381.75
No. Observations:                 271   AIC:                            -759.5
Df Residuals:                     269   BIC:                            -752.3
Df Model:                           1                                         
Covariance Type:                  HAC                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
const         -0.0097      0.004     -2.618      0.009      -0.017      -0.002
mkt_excess     0.5629      0.070      7.992      0.000       0.425       0.701
==============================================================================
Omnibus:                       60.286   Durbin-Watson:                   1.841
Prob(Omnibus):                  0.000   Jarque-Bera (JB):              262.919
Skew:                          -0.826   Prob(JB):                     8.09e-58
Kurtosis:                       7.534   Cond. No.                         13.1
==============================================================================

Notes:
[1] Standard Errors are heteroscedasticity and autocorrelation robust (HAC) using 6 lags and without small sample correction

If the long–short mean return is insignificant but CAPM alpha is significantly negative after controlling for the market, the result conflicts with CAPM and aligns with betting-against-beta (long low beta, short high beta) in the literature [@Frazzini2014]. Note that the Frazzini–Pedersen factor construction is not a standard value-weighted decile sort; whether its anomaly returns are robust remains debated [@NovyMarx2022].

6.8.7 Annual Return Decomposition

Corresponds to the annual bar charts in the original: inspect cumulative yearly returns of low beta, high beta, and the long–short book to see whether a few years dominate the full-sample conclusion.

tmp = beta_longshort.copy
tmp["year"] = tmp["date"].dt.year

def ann_comp(s: pd.Series) -> float:
    return float(np.prod(1.0 + s.dropna) - 1.0)

annual = (
    tmp.groupby("year")
    .agg(
        low=("low", ann_comp),
        high=("high", ann_comp),
        long_short=("long_short", ann_comp),
    )
    .reset_index
)

plot_df = annual.melt(
    id_vars="year",
    value_vars=["low", "high", "long_short"],
    var_name="name",
    value_name="value",
)

fig, axes = plt.subplots(3, 1, figsize=(10, 8), sharex=True)
for ax, name in zip(axes, ["low", "high", "long_short"]):
    sub = plot_df[plot_df["name"] == name]
    ax.bar(sub["year"], sub["value"], color="#4C78A8", width=0.8)
    ax.axhline(0, color="black", linewidth=0.8)
    ax.set_ylabel(name)
    ax.yaxis.set_major_formatter(plt.FuncFormatter(lambda x, _: f"{100*x:.0f}%"))
    ax.grid(True, axis="y", alpha=0.3)
axes[0].set_title("Annual returns of beta portfolios (A-share)")
axes[-1].set_xlabel("")
plt.tight_layout
plt.show

display(annual.tail(10))

	year	low	high	long_short
13	2017	0.287904	-0.157718	-0.353551
14	2018	-0.169761	-0.321624	-0.194517
15	2019	0.242932	0.247311	0.012042
16	2020	0.097370	0.157694	0.056586
17	2021	0.018111	0.079606	0.055313
18	2022	-0.081621	-0.204390	-0.125512
19	2023	0.121660	-0.132788	-0.236273
20	2024	0.281498	0.026509	-0.206064
21	2025	0.093440	0.269349	0.149724
22	2026	0.010071	-0.081584	-0.173771

6.8.8 Key Takeaways

• Univariate portfolio sorts test whether a single characteristic (e.g., lagged market beta) forecasts future excess returns.
• Form groups each month with quantile breakpoints and compute value-weighted portfolio returns—closer to investable strategies.
• Beta-sorted long–short books often fail to deliver stably positive mean excess returns, conflicting with the CAPM intuition that high beta should earn higher expected returns.
• Low-beta portfolios often show relatively higher CAPM alphas (and high-beta the opposite)—evidence consistent with betting-against-beta / the low-beta anomaly.
• Once assign_portfolio and value-weighted returns are functions, the same pipeline migrates to any sorting variable (value, momentum, quality, etc.).

Exercises (adapted from the original for A-shares)

1. Compare median two-way long–short versus extreme-decile long–short: mean return, volatility, Sharpe, and Newey–West $t$—are they materially different?
2. Test CAPM alpha significance for each of the ten portfolios (HAC is fine). Which alphas are significant?
3. Replace the sorting variable with daily beta from Section 2.6 (return_type=="daily" in beta_ashare.parquet), repeat the analysis, and note differences versus monthly beta.
4. Based on this section, design a long–short direction that may earn positive alpha (hint: opposite of “long high / short low”), plot cumulative excess returns, and compare with market excess returns.

6.9 A-Share Replication: Size Sorts and p-Hacking

This section replicates the tidyfinance chapter Size Sorts and p-Hacking: within the univariate-sort framework, switch the sorting variable to firm size (free-float market cap) to recover the classic size premium (long small, short big); then show how research-design choices—breakpoint sample, weighting, number of groups, sample period—change premium estimates. That is the non-standard-errors / p-hacking risk in sort settings.

import warnings
from itertools import product
from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from IPython.display import display
from matplotlib import font_manager

warnings.filterwarnings("ignore", category=UserWarning)

_cn_candidates = ["Microsoft YaHei", "SimHei", "PingFang SC", "Noto Sans CJK SC", "Arial Unicode MS"]
_available = {f.name for f in font_manager.fontManager.ttflist}
for _name in _cn_candidates:
    if _name in _available:
        plt.rcParams["font.sans-serif"] = [_name]
        break
plt.rcParams["axes.unicode_minus"] = False

DATA_DIR = Path(r"D:/A_Topics/202607_03_TidyFinanceAShare/202608_02_传统量化/Data")
CACHE_DIR = DATA_DIR / "_cache_beta"

print("monthly cache exists =", (CACHE_DIR / "crsp_monthly_ashare.parquet").exists)

monthly cache exists = True

6.9.1 Data Preparation

Corresponds to Data Preparation. Size uses lagged free-float market cap mktcap_lag (known at formation); returns use ret_excess. Boards are mapped from permno rules as the A-share analogue of U.S. exchange.

def map_board(permno: str) -> str:
    # Infer listing board from security code (simplified rules)
    code = str(permno).zfill(6)
    if code.startswith("688"):
        return "STAR"
    if code.startswith(("300", "301")):
        return "ChiNext"
    if code.startswith("6"):
        return "SSE Main"
    if code.startswith(("000", "001", "002", "003")):
        return "SZSE Main"
    return "Other"

crsp_monthly = pd.read_parquet(CACHE_DIR / "crsp_monthly_ashare.parquet").copy
crsp_monthly["date"] = pd.to_datetime(crsp_monthly["date"])
crsp_monthly["permno"] = crsp_monthly["permno"].astype(str).str.zfill(6)
crsp_monthly["board"] = crsp_monthly["permno"].map(map_board)

# Sorting/weighting need lagged cap; missing => exclude that month
crsp_monthly = crsp_monthly.dropna(subset=["ret_excess", "mktcap", "mktcap_lag"]).copy
crsp_monthly = crsp_monthly[crsp_monthly["mktcap_lag"] > 0].copy

print(crsp_monthly.shape)
print(crsp_monthly["board"].value_counts)
display(crsp_monthly.head(3))

(884560, 11)
board
上交所主板    344321
深交所主板    318516
创业板      143454
其他        44021
科创板       34248
Name: count, dtype: int64

	permno	month	ret	mktcap	ret_m	rf	date	ret_excess	mkt_excess	mktcap_lag	board
1	000001	2000-02	-0.011333	19618933.36	0.122024	0.001856	2000-02-01	-0.013189	0.120168	19843822.88	深交所主板
2	000001	2000-03	0.002729	19672478.48	0.055926	0.001856	2000-03-01	0.000873	0.054070	19618933.36	深交所主板
3	000001	2000-04	0.037017	20400692.17	0.013014	0.001856	2000-04-01	0.035161	0.011158	19672478.48	深交所主板

6.9.2 Size Distribution

Corresponds to Size Distribution. First check concentration: share of free-float market cap held by the largest 1%/5%/10%/25% of firms. High head concentration means value-weighted portfolios are dominated by mega-caps.

def top_cap_share(g: pd.DataFrame, q: float) -> float:
    thr = g["mktcap"].quantile(q)
    total = g["mktcap"].sum
    if total <= 0:
        return np.nan
    return float(g.loc[g["mktcap"] >= thr, "mktcap"].sum / total)

rows = []
for dt, g in crsp_monthly.groupby("date", sort=True):
    rows.append(
        {
            "date": dt,
            "Largest 1%": top_cap_share(g, 0.99),
            "Largest 5%": top_cap_share(g, 0.95),
            "Largest 10%": top_cap_share(g, 0.90),
            "Largest 25%": top_cap_share(g, 0.75),
        }
    )
conc = pd.DataFrame(rows).sort_values("date")
plot_df = conc.melt(id_vars="date", var_name="name", value_name="value")

fig, ax = plt.subplots(figsize=(10, 4.5))
for name, sub in plot_df.groupby("name"):
    ax.plot(sub["date"], sub["value"], label=name, linewidth=1.6)
ax.set_title("Share of free-float market cap held by largest firms (A-share)")
ax.set_ylabel("Share")
ax.yaxis.set_major_formatter(plt.FuncFormatter(lambda x, _: f"{100*x:.0f}%"))
ax.legend(frameon=False, ncol=2)
ax.grid(True, alpha=0.3)
plt.tight_layout
plt.show

display(conc.tail(3))

	date	Largest 1%	Largest 5%	Largest 10%	Largest 25%
315	2026-05-01	0.300931	0.529262	0.648327	0.809577
316	2026-06-01	0.300993	0.539631	0.662786	0.823773
317	2026-07-01	0.320657	0.552181	0.668810	0.822761

Next plot each listing board’s market-cap share over time (analogue of the original NYSE/NASDAQ/AMEX stacked area chart).

share = (
    crsp_monthly.groupby(["date", "board"], as_index=False)["mktcap"]
    .sum
    .rename(columns={"mktcap": "board_cap"})
)
share["total"] = share.groupby("date")["board_cap"].transform("sum")
share["share"] = share["board_cap"] / share["total"]

pivot = share.pivot(index="date", columns="board", values="share").fillna(0).sort_index
# Fixed stacking order
cols = [c for c in ["SSE Main", "SZSE Main", "ChiNext", "STAR", "Other"] if c in pivot.columns]
pivot = pivot[cols]

fig, ax = plt.subplots(figsize=(10, 4.5))
ax.stackplot(pivot.index, [pivot[c].values for c in cols], labels=cols, alpha=0.85)
ax.set_title("Listing-board share of free-float market cap (A-share)")
ax.set_ylabel("Share")
ax.set_ylim(0, 1)
ax.yaxis.set_major_formatter(plt.FuncFormatter(lambda x, _: f"{100*x:.0f}%"))
ax.legend(loc="upper left", ncol=3, frameon=False, fontsize=9)
ax.grid(True, axis="y", alpha=0.3)
plt.tight_layout
plt.show

Descriptive statistics of market cap by board in the latest month: main-board firms are usually larger on average; ChiNext/STAR have many names but more dispersed size. That is why “main-board-only breakpoints” affect how thick the small-cap portfolio is.

latest = crsp_monthly["date"].max
snap = crsp_monthly[crsp_monthly["date"] == latest].copy

def board_summary(df: pd.DataFrame) -> pd.DataFrame:
    def one(g):
        s = g["mktcap"]
        return pd.Series(
            {
                "count": len(s),
                "mean": s.mean,
                "std": s.std,
                "min": s.min,
                "p5": s.quantile(0.05),
                "p50": s.quantile(0.50),
                "p95": s.quantile(0.95),
                "max": s.max,
            }
        )

    by = df.groupby("board").apply(one, include_groups=False)
    overall = one(df)
    overall.name = "Overall"
    out = pd.concat([by, overall.to_frame.T])
    return out

print("snapshot month:", latest.date)
summary = board_summary(snap)
display(summary.round(0))

snapshot month: 2026-07-01

	count	mean	std	min	p5	p50	p95	max
上交所主板	1699.0	29093131.0	119091033.0	722864.0	1825229.0	6953583.0	104098370.0	2.202785e+09
其他	396.0	1315870.0	2552979.0	16781.0	78636.0	716231.0	4035889.0	3.429679e+07
创业板	1397.0	10177802.0	56703424.0	24043.0	1193492.0	3487125.0	26520174.0	1.684010e+09
深交所主板	1494.0	14654251.0	35785984.0	37467.0	1543550.0	5414599.0	56215101.0	6.028877e+08
科创板	608.0	16589284.0	46424082.0	209392.0	1581280.0	5953805.0	53572655.0	6.948920e+08
Overall	5594.0	17187780.0	75890924.0	16781.0	861935.0	4741815.0	56425655.0	2.202785e+09

6.9.3 Univariate Size Portfolios with Flexible Breakpoints

Corresponds to Univariate Size Portfolios with Flexible Breakpoints. Key extension: breakpoints can be computed on a subset (e.g., SSE Main only) and then applied to all stocks that month—analogue of U.S. “NYSE breakpoints.”

def assign_portfolio(
    data: pd.DataFrame,
    boards: list[str],
    sorting_variable: str,
    n_portfolios: int,
) -> pd.Series:
    # Compute quantile breakpoints on the boards subset, then map to the full sample
    subset = data[data["board"].isin(boards)]
    if len(subset) < n_portfolios:
        # Degenerate: everyone in one group
        return pd.Series(1, index=data.index, dtype=int)

    quantiles = np.linspace(0, 1, n_portfolios + 1)
    breakpoints = subset[sorting_variable].quantile(quantiles).to_numpy(dtype=float)
    breakpoints = np.unique(breakpoints)
    if len(breakpoints) < 2:
        return pd.Series(1, index=data.index, dtype=int)

    # cut needs interior breakpoints; extend ends to -inf/inf
    inner = breakpoints[1:-1]
    labels = list(range(1, len(inner) + 2))
    try:
        port = pd.cut(
            data[sorting_variable],
            bins=[-np.inf, *inner, np.inf],
            labels=labels,
            include_lowest=True,
        )
    except ValueError:
        return pd.Series(1, index=data.index, dtype=int)
    return port.astype(int)

# Example: one month, SSE-Main breakpoints vs all-board breakpoints — count of names in the small group
demo_date = crsp_monthly["date"].max
g = crsp_monthly[crsp_monthly["date"] == demo_date].copy
g["port_all"] = assign_portfolio(g, g["board"].unique.tolist, "mktcap_lag", 2)
g["port_sh"] = assign_portfolio(g, ["SSE Main"], "mktcap_lag", 2)
print("month:", demo_date.date)
print("small-group count | all-board breakpoints:", int((g["port_all"] == 1).sum))
print("small-group count | SSE-main breakpoints:", int((g["port_sh"] == 1).sum))

month: 2026-07-01
small-group count | all-board breakpoints: 2797
small-group count | SH-main breakpoints: 3408

6.9.4 Portfolio Weighting: Value-Weighted vs Equal-Weighted

Corresponds to Weighting Schemes. Value-weighting is closer to a passive investable book; equal-weighting requires monthly rebalancing, is costlier in practice, and often amplifies the small-cap premium. The size premium is defined as: smallest-size-group return − largest-size-group return (small minus big).

def compute_portfolio_returns(
    n_portfolios: int = 10,
    boards: list[str] | None = None,
    value_weighted: bool = True,
    data: pd.DataFrame | None = None,
) -> float:
    # Full-sample average size premium (time-series mean of monthly small-big)
    if data is None:
        data = crsp_monthly
    if boards is None:
        boards = sorted(data["board"].dropna.unique.tolist)

    premia = []
    for dt, g in data.groupby("date", sort=True):
        g = g.copy
        g["portfolio"] = assign_portfolio(g, boards, "mktcap_lag", n_portfolios)
        # Group returns
        rets = {}
        for port, pg in g.groupby("portfolio"):
            if value_weighted:
                w = pg["mktcap_lag"]
                if w.sum <= 0:
                    continue
                rets[int(port)] = float((pg["ret_excess"] * w).sum / w.sum)
            else:
                rets[int(port)] = float(pg["ret_excess"].mean)
        if len(rets) < 2:
            continue
        pmin, pmax = min(rets), max(rets)
        premia.append(rets[pmin] - rets[pmax])

    if not premia:
        return np.nan
    return float(np.mean(premia))

boards_all = ["SSE Main", "SZSE Main", "ChiNext", "STAR"]
ret_all = compute_portfolio_returns(
    n_portfolios=2, boards=boards_all, value_weighted=True, data=crsp_monthly
)
ret_sh = compute_portfolio_returns(
    n_portfolios=2, boards=["SSE Main"], value_weighted=True, data=crsp_monthly
)

cmp = pd.DataFrame(
    {
        "breakpoint_sample": ["All major boards", "SSE Main only"],
        "avg_size_premium_m": [ret_all, ret_sh],
        "ann_approx_pct": [ret_all * 12 * 100, ret_sh * 12 * 100],
    }
)
display(cmp.round(4))
print("Note: with SSE-Main-only breakpoints, the market-wide small group is usually thicker, and premium estimates often change materially.")

	断点样本	平均规模溢价(月)	折年约%
0	全部主要板块	0.0081	9.6901
1	仅上交所主板	0.0068	8.1625

说明：仅用上交所主板断点时，全市场小市值组通常更‘厚’，溢价估计往往变化明显。

6.9.5 p-Hacking and Non-Standard Errors

Corresponds to P-Hacking and Non-Standard Errors. A sort at least requires choosing: number of groups, breakpoint sample, equal vs value weights, and whether to truncate the sample period or drop boards. Each choice has precedents in top journals; none is necessarily “wrong,” but they create non-standard errors (researcher-choice variation). Reporting only the most significant specification slides into p-hacking.

Below we scan a grid of specifications and inspect the distribution of size premia. To keep runtime manageable the grid is slightly tighter than the original, but it still covers the key decision nodes.

n_portfolios_grid = [2, 5, 10]
boards_grid = [
    ["SSE Main"],
    ["SSE Main", "SZSE Main", "ChiNext", "STAR"],
]
value_weighted_grid = [True, False]

# Data slices: full sample / main boards only / pre-2010 / 2010+
data_main = crsp_monthly[crsp_monthly["board"].isin(["SSE Main", "SZSE Main"])].copy
data_pre = crsp_monthly[crsp_monthly["date"] < "2010-01-01"].copy
data_post = crsp_monthly[crsp_monthly["date"] >= "2010-01-01"].copy
data_grid = [
    ("full_sample", crsp_monthly),
    ("main_boards", data_main),
    ("pre_2010", data_pre),
    ("from_2010", data_post),
]

setup = list(product(n_portfolios_grid, boards_grid, value_weighted_grid, data_grid))
print("specifications:", len(setup))

records = []
for n_port, boards, vw, (tag, df) in setup:
    if len(df) < 1000 or df["date"].nunique < 24:
        prem = np.nan
    else:
        prem = compute_portfolio_returns(
            n_portfolios=n_port,
            boards=boards,
            value_weighted=vw,
            data=df,
        )
    records.append(
        {
            "n_portfolios": n_port,
            "breakpoint_boards": "+".join(boards) if len(boards) == 1 else "multi_board",
            "value_weighted": vw,
            "sample": tag,
            "size_premium": prem,
        }
    )

p_hacking_results = pd.DataFrame(records).dropna(subset=["size_premium"])
p_hacking_results = p_hacking_results.sort_values("size_premium", ascending=False)
print(p_hacking_results.head(10))
print("...")
print(p_hacking_results.tail(5))
print(
    "mean=",
    round(p_hacking_results["size_premium"].mean, 5),
    "std=",
    round(p_hacking_results["size_premium"].std, 5),
    "min=",
    round(p_hacking_results["size_premium"].min, 5),
    "max=",
    round(p_hacking_results["size_premium"].max, 5),
)

specifications: 48
    n_portfolios breakpoint_boards  value_weighted sample  size_premium
38            10             上交所主板           False  2010前      0.021159
46            10               多板块           False  2010前      0.020778
34            10             上交所主板            True  2010前      0.017018
45            10               多板块           False    仅主板      0.016913
42            10               多板块            True  2010前      0.016659
44            10               多板块           False    全样本      0.016365
36            10             上交所主板           False    全样本      0.015902
30             5               多板块           False  2010前      0.015687
37            10             上交所主板           False    仅主板      0.015563
22             5             上交所主板           False  2010前      0.014876
...
    n_portfolios breakpoint_boards  value_weighted   sample  size_premium
10             2               多板块            True    2010前      0.007132
3              2             上交所主板            True  2010及以后      0.006864
0              2             上交所主板            True      全样本      0.006802
2              2             上交所主板            True    2010前      0.006699
1              2             上交所主板            True      仅主板      0.006532
mean= 0.01224 std= 0.00363 min= 0.00653 max= 0.02116

6.9.6 Size-Premium Variation across Specs

Corresponds to Size-Premium Variation. The histogram shows monthly average size premia under different specs; the dashed vertical line is a “benchmark” premium (deciles, SSE-Main breakpoints, value-weighted, full sample), playing the role of the FF-SMB mean reference line in the original.

# Benchmark spec (common literature-style choice; not uniquely correct)
benchmark_premium = compute_portfolio_returns(
    n_portfolios=10,
    boards=["SSE Main"],
    value_weighted=True,
    data=crsp_monthly,
)

fig, ax = plt.subplots(figsize=(9, 4.5))
vals = p_hacking_results["size_premium"].values
ax.hist(vals, bins=min(40, max(10, len(vals) // 2)), color="#4C78A8", edgecolor="white")
ax.axvline(benchmark_premium, color="crimson", linestyle="--", linewidth=1.8, label="Benchmark spec")
ax.set_title("Distribution of size premia across sort specs (A-share)")
ax.set_xlabel("Average monthly size premium (small - big)")
ax.set_ylabel("Number of specs")
ax.xaxis.set_major_formatter(plt.FuncFormatter(lambda x, _: f"{100*x:.2f}%"))
ax.legend(frameon=False)
ax.grid(True, axis="y", alpha=0.3)
plt.tight_layout
plt.show

print(f"benchmark monthly size premium = {benchmark_premium:.4%}  (ann. ~ {benchmark_premium*12:.2%})")

# Which decision node matters most? dispersion of premia across levels
impact = []
for col in ["n_portfolios", "breakpoint_boards", "value_weighted", "sample"]:
    grp = p_hacking_results.groupby(col)["size_premium"].mean
    impact.append({"choice": col, "range_across_levels": float(grp.max - grp.min)})
impact_df = pd.DataFrame(impact).sort_values("range_across_levels", ascending=False)
print("\nDecision nodes: range of average premia across levels (larger => more influential)")
display(impact_df)

benchmark monthly size premium = 1.4274%  (ann. ~ 17.13%)

各决策节点：不同水平下平均溢价的极差（越大说明该选择越能左右结果）

	choice	range_across_levels
0	n_portfolios	0.007701
3	sample	0.002164
2	value_weighted	0.001620
1	breakpoint_boards	0.000731

6.9.7 Key Takeaways

• Size is an important sorting variable; small-minus-big returns are the size premium and one empirical source of SMB in the FF three-factor model.
• Portfolio returns are sensitive to design choices such as breakpoint sample, weighting, number of groups, and sample-period / board filters.
• Variation from reasonable alternative pipelines is non-standard errors; reporting only the most significant spec approaches p-hacking.
• A robust practice: pre-register a primary spec and report a result range under several reasonable alternatives (see Section 4.1.4).

Exercises (adapted from the original for A-shares)

1. Plot average market cap over time for SSE Main, SZSE Main, and ChiNext; interpret structural change.
2. Plot the time series of decile size-sort breakpoints (or extreme-group breakpoint gaps) and check for abnormal jumps.
3. Modify compute_portfolio_returns to output CAPM alpha of the long–short book (vs mkt_excess) and compare conclusions with mean excess returns.
4. In the specification grid, use regression or grouped statistics to find the single decision variable with the largest impact on the premium.
5. Which specification is closest to a common SMB convention—“main-board breakpoints + value weights + small-minus-big”? Where does it sit in the distribution?

6.10 A-Share Replication: Value and Bivariate Sorts

This section replicates the tidyfinance chapter Value and Bivariate Sorts: on top of univariate sorts, introduce book-to-market (BM) and run Size × BM independent and dependent (conditional) bivariate sorts; finally compare stock counts and market-cap shares in the 5×5 grids under both protocols.

import warnings
from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from IPython.display import display
from matplotlib import font_manager

warnings.filterwarnings("ignore", category=UserWarning)

_cn_candidates = ["Microsoft YaHei", "SimHei", "PingFang SC", "Noto Sans CJK SC", "Arial Unicode MS"]
_available = {f.name for f in font_manager.fontManager.ttflist}
for _name in _cn_candidates:
    if _name in _available:
        plt.rcParams["font.sans-serif"] = [_name]
        break
plt.rcParams["axes.unicode_minus"] = False

DATA_DIR = Path(r"D:/A_Topics/202607_03_TidyFinanceAShare/202608_02_传统量化/Data")
CACHE_DIR = DATA_DIR / "_cache_beta"
BE_CACHE = CACHE_DIR / "book_equity_ashare.parquet"

print("crsp cache =", (CACHE_DIR / "crsp_monthly_ashare.parquet").exists)
print("BE cache   =", BE_CACHE.exists)

crsp cache = True
BE cache   = True

6.10.1 Data Preparation

Corresponds to Data Preparation. Monthly returns/caps come from the existing panel; if book equity is not cached locally, pull annual reports via akshare.stock_zcfz_em at year-end (YYYY1231) and write a parquet (slow the first time; thereafter read the cache).

def map_board(permno: str) -> str:
    code = str(permno).zfill(6)
    if code.startswith("688"):
        return "STAR"
    if code.startswith(("300", "301")):
        return "ChiNext"
    if code.startswith("6"):
        return "SSE Main"
    if code.startswith(("000", "001", "002", "003")):
        return "SZSE Main"
    return "Other"

def load_or_download_book_equity(
    cache_path: Path,
    years: range | list[int] = range(1999, 2026),
) -> pd.DataFrame:
    """Load cached annual book equity; download via akshare if missing."""
    if cache_path.exists:
        be = pd.read_parquet(cache_path)
        be["permno"] = be["permno"].astype(str).str.zfill(6)
        be["datadate"] = pd.to_datetime(be["datadate"])
        return be

    import time
    import akshare as ak

    print("Cache missing; downloading annual shareholders' equity (first run ~10+ minutes)…")
    frames: list[pd.DataFrame] = []
    partial = cache_path.with_suffix(".partial.parquet")
    done: set[int] = set
    if partial.exists:
        old = pd.read_parquet(partial)
        frames.append(old)
        done = set(pd.to_datetime(old["datadate"]).dt.year.unique.tolist)
        print("resume years:", sorted(done))

    for y in years:
        if y in done:
            continue
        date = f"{y}1231"
        print("fetch", date)
        try:
            raw = ak.stock_zcfz_em(date=date)
        except Exception as exc:  # noqa: BLE001
            print("  FAIL", date, exc)
            time.sleep(2)
            continue
        cols = list(raw.columns)
        # Match Eastmoney Chinese column names from akshare
        code_col = [c for c in cols if "代码" in str(c)][0]
        be_col = [c for c in cols if "股东权益" in str(c)][0]
        out = pd.DataFrame(
            {
                "permno": raw[code_col].astype(str).str.zfill(6),
                "datadate": pd.to_datetime(date),
                "be": pd.to_numeric(raw[be_col], errors="coerce"),
            }
        )
        out = out.dropna(subset=["be"])
        out = out[out["be"] > 0].copy
        # Eastmoney unit: yuan; CSMAR Msmvosd: thousand yuan
        out["be"] = out["be"] / 1000.0
        frames.append(out)
        done.add(y)
        pd.concat(frames, ignore_index=True).to_parquet(partial, index=False)
        time.sleep(0.5)

    be = pd.concat(frames, ignore_index=True)
    be = be.drop_duplicates(["permno", "datadate"], keep="last")
    cache_path.parent.mkdir(parents=True, exist_ok=True)
    be.to_parquet(cache_path, index=False)
    if partial.exists:
        partial.unlink
    print("saved", cache_path, be.shape)
    return be

crsp_monthly = pd.read_parquet(CACHE_DIR / "crsp_monthly_ashare.parquet").copy
crsp_monthly["date"] = pd.to_datetime(crsp_monthly["date"])
crsp_monthly["permno"] = crsp_monthly["permno"].astype(str).str.zfill(6)
crsp_monthly["board"] = crsp_monthly["permno"].map(map_board)
crsp_monthly = crsp_monthly.dropna(subset=["ret_excess", "mktcap", "mktcap_lag"]).copy
crsp_monthly = crsp_monthly[crsp_monthly["mktcap_lag"] > 0].copy

book_equity = load_or_download_book_equity(BE_CACHE)
book_equity = book_equity[book_equity["be"] > 0].dropna(subset=["be"]).copy
# Same as original: align accounting date to month (annual reports as month-start of December)
book_equity["date"] = book_equity["datadate"].dt.to_period("M").dt.to_timestamp

print("crsp:", crsp_monthly.shape, "BE rows:", len(book_equity),
      "years:", book_equity["datadate"].dt.year.min, "-", book_equity["datadate"].dt.year.max)
display(book_equity.head(3))

crsp: (884560, 11) BE rows: 88343 years: 1999 - 2025

	permno	datadate	be	date
0	601963	1999-12-31	4.083935e+05	1999-12-01
1	601229	1999-12-31	4.182000e+06	1999-12-01
2	601187	1999-12-31	3.140239e+05	1999-12-01

6.10.2 Book-to-Market Ratio

Corresponds to Book-to-Market Ratio. The core is avoiding look-ahead bias:

1. Size is lagged by 1 month (sorting_date = date + 1M);
2. In the filing month, match be to that month’s mktcap to get BM, then lag the whole series by 6 months;
3. Forward-fill BM by stock and drop observations whose accounting information is older than 12 months.

High BM = value; low BM = growth.

size = (
    crsp_monthly.assign(sorting_date=lambda d: d["date"] + pd.DateOffset(months=1))
    .rename(columns={"mktcap": "size"})
    [["permno", "sorting_date", "size"]]
)

bm = (
    book_equity.merge(
        crsp_monthly[["permno", "date", "mktcap"]],
        on=["permno", "date"],
        how="inner",
    )
    .assign(
        bm=lambda d: d["be"] / d["mktcap"],
        sorting_date=lambda d: d["date"] + pd.DateOffset(months=6),
        accounting_date=lambda d: d["date"] + pd.DateOffset(months=6),
    )
    [["permno", "sorting_date", "accounting_date", "bm"]]
)
bm = bm.replace([np.inf, -np.inf], np.nan).dropna(subset=["bm"])
bm = bm[bm["bm"] > 0].copy

data_for_sorts = (
    crsp_monthly.merge(
        bm,
        left_on=["permno", "date"],
        right_on=["permno", "sorting_date"],
        how="left",
        suffixes=("", "_bm"),
    )
    .drop(columns=["sorting_date"], errors="ignore")
    .merge(
        size,
        left_on=["permno", "date"],
        right_on=["permno", "sorting_date"],
        how="left",
    )
    .drop(columns=["sorting_date"], errors="ignore")
)

data_for_sorts = data_for_sorts.sort_values(["permno", "date"]).copy
data_for_sorts["bm"] = data_for_sorts.groupby("permno")["bm"].ffill
data_for_sorts["accounting_date"] = data_for_sorts.groupby("permno")["accounting_date"].ffill
data_for_sorts["threshold_date"] = data_for_sorts["date"] - pd.DateOffset(months=12)
data_for_sorts = data_for_sorts[
    data_for_sorts["accounting_date"] > data_for_sorts["threshold_date"]
].copy
data_for_sorts = data_for_sorts.dropna(
    subset=["ret_excess", "mktcap_lag", "size", "bm", "board"]
).copy
data_for_sorts = data_for_sorts.drop(columns=["threshold_date", "accounting_date"])

print(data_for_sorts.shape)
print(
    "date range:",
    data_for_sorts["date"].min.date,
    "→",
    data_for_sorts["date"].max.date,
)
display(data_for_sorts[["permno", "date", "size", "bm", "ret_excess", "board"]].head(3))

(703738, 13)
date range: 2001-06-01 → 2026-07-01

	permno	date	size	bm	ret_excess	board
16	000001	2001-06-01	22568621.18	0.173894	-0.056794	深交所主板
17	000001	2001-07-01	21328740.14	0.173894	-0.098525	深交所主板
18	000001	2001-08-01	19266915.49	0.173894	-0.082839	深交所主板

Breakpoint function is the same as in 6.9: compute quantiles on a board subset, then map to the full cross-section that month.

def assign_portfolio(
    data: pd.DataFrame,
    boards: list[str],
    sorting_variable: str,
    n_portfolios: int,
) -> pd.Series:
    subset = data[data["board"].isin(boards)]
    if len(subset) < n_portfolios:
        return pd.Series(1, index=data.index, dtype=int)

    quantiles = np.linspace(0, 1, n_portfolios + 1)
    breakpoints = subset[sorting_variable].quantile(quantiles).to_numpy(dtype=float)
    breakpoints = np.unique(breakpoints)
    if len(breakpoints) < 2:
        return pd.Series(1, index=data.index, dtype=int)

    inner = breakpoints[1:-1]
    labels = list(range(1, len(inner) + 2))
    try:
        port = pd.cut(
            data[sorting_variable],
            bins=[-np.inf, *inner, np.inf],
            labels=labels,
            include_lowest=True,
        )
    except ValueError:
        return pd.Series(1, index=data.index, dtype=int)
    return port.astype(int)

BREAKPOINT_BOARDS = ["SSE Main"]
print("default breakpoint boards:", BREAKPOINT_BOARDS)

default breakpoint boards: ['上交所主板']

6.10.3 Independent Bivariate Sorts

Corresponds to Independent Sorts. Each month, independently split Size and BM into 5 groups using main-board breakpoints and cross them into 25 portfolios; within portfolios, value-weight by mktcap_lag. Value premium = equal-weighted average of the five high-BM portfolios − equal-weighted average of the five low-BM portfolios (first collapsing across size).

def monthly_independent_returns(df: pd.DataFrame) -> pd.DataFrame:
    rows = []
    for dt, g in df.groupby("date", sort=True):
        g = g.copy
        g["portfolio_bm"] = assign_portfolio(g, BREAKPOINT_BOARDS, "bm", 5)
        g["portfolio_size"] = assign_portfolio(g, BREAKPOINT_BOARDS, "size", 5)
        for (p_bm, p_sz), pg in g.groupby(["portfolio_bm", "portfolio_size"]):
            w = pg["mktcap_lag"]
            if w.sum <= 0:
                continue
            rows.append(
                {
                    "date": dt,
                    "portfolio_bm": int(p_bm),
                    "portfolio_size": int(p_sz),
                    "ret": float((pg["ret_excess"] * w).sum / w.sum),
                }
            )
    return pd.DataFrame(rows)

value_portfolios_ind = monthly_independent_returns(data_for_sorts)

# Each month: equal-weight the 5 size portfolios within each BM group, then high-low
vp = (
    value_portfolios_ind.groupby(["date", "portfolio_bm"], as_index=False)["ret"]
    .mean
)
prem_ind = []
for dt, g in vp.groupby("date"):
    hi = g.loc[g["portfolio_bm"] == g["portfolio_bm"].max, "ret"].mean
    lo = g.loc[g["portfolio_bm"] == g["portfolio_bm"].min, "ret"].mean
    prem_ind.append(hi - lo)
value_premium_ind = float(np.mean(prem_ind))

print(
    f"independent value premium (monthly) = {value_premium_ind:.4%}  "
    f"(ann. ~ {value_premium_ind * 12:.2%})"
)
display(value_portfolios_ind.head(3))

independent value premium (monthly) = 0.1419%  (ann. ~ 1.70%)

	date	portfolio_bm	portfolio_size	ret
0	2001-06-01	1	1	0.022497
1	2001-06-01	1	2	-0.014334
2	2001-06-01	1	3	-0.009245

6.10.4 Dependent Bivariate Sorts

Corresponds to Dependent Sorts. First split Size into 5 groups, then within each size group split BM into 5—BM breakpoints vary by size bucket. Dependent sorts often produce more even stock counts when breakpoints use the full market; if breakpoints still use main boards only, within-group counts may remain uneven, but relative evenness along BM is usually better than under independent sorts.

def monthly_dependent_returns(df: pd.DataFrame) -> pd.DataFrame:
    rows = []
    for dt, g in df.groupby("date", sort=True):
        g = g.copy
        g["portfolio_size"] = assign_portfolio(g, BREAKPOINT_BOARDS, "size", 5)
        parts = []
        for _, sg in g.groupby("portfolio_size"):
            sg = sg.copy
            sg["portfolio_bm"] = assign_portfolio(sg, BREAKPOINT_BOARDS, "bm", 5)
            parts.append(sg)
        g2 = pd.concat(parts, ignore_index=False)
        for (p_bm, p_sz), pg in g2.groupby(["portfolio_bm", "portfolio_size"]):
            w = pg["mktcap_lag"]
            if w.sum <= 0:
                continue
            rows.append(
                {
                    "date": dt,
                    "portfolio_bm": int(p_bm),
                    "portfolio_size": int(p_sz),
                    "ret": float((pg["ret_excess"] * w).sum / w.sum),
                }
            )
    return pd.DataFrame(rows)

value_portfolios_dep = monthly_dependent_returns(data_for_sorts)
vp_d = (
    value_portfolios_dep.groupby(["date", "portfolio_bm"], as_index=False)["ret"]
    .mean
)
prem_dep = []
for dt, g in vp_d.groupby("date"):
    hi = g.loc[g["portfolio_bm"] == g["portfolio_bm"].max, "ret"].mean
    lo = g.loc[g["portfolio_bm"] == g["portfolio_bm"].min, "ret"].mean
    prem_dep.append(hi - lo)
value_premium_dep = float(np.mean(prem_dep))

cmp = pd.DataFrame(
    {
        "sort_method": ["Independent", "Dependent"],
        "avg_value_premium_m": [value_premium_ind, value_premium_dep],
        "ann_approx_pct": [value_premium_ind * 12 * 100, value_premium_dep * 12 * 100],
    }
)
display(cmp.round(4))

	排序方式	月均价值溢价	折年约%
0	独立排序	0.0014	1.7027
1	依赖排序	0.0013	1.5334

6.10.5 Portfolio Composition: Stock Counts and Cap Shares

Corresponds to Portfolio Composition. Beyond returns, inspect how stocks and market cap concentrate in the 5×5 grid—independent vs dependent differs mainly in whether BM piles into corners.

def assign_independent(df: pd.DataFrame) -> pd.DataFrame:
    parts = []
    for dt, g in df.groupby("date", sort=True):
        g = g.copy
        g["portfolio_size"] = assign_portfolio(g, BREAKPOINT_BOARDS, "size", 5)
        g["portfolio_bm"] = assign_portfolio(g, BREAKPOINT_BOARDS, "bm", 5)
        g["sorting_method"] = "Independent"
        parts.append(g)
    return pd.concat(parts, ignore_index=True)

def assign_dependent(df: pd.DataFrame) -> pd.DataFrame:
    parts = []
    for dt, g in df.groupby("date", sort=True):
        g = g.copy
        g["portfolio_size"] = assign_portfolio(g, BREAKPOINT_BOARDS, "size", 5)
        sub = []
        for _, sg in g.groupby("portfolio_size"):
            sg = sg.copy
            sg["portfolio_bm"] = assign_portfolio(sg, BREAKPOINT_BOARDS, "bm", 5)
            sub.append(sg)
        g2 = pd.concat(sub, ignore_index=False)
        g2["sorting_method"] = "Dependent"
        parts.append(g2)
    return pd.concat(parts, ignore_index=True)

assignments = pd.concat(
    [assign_independent(data_for_sorts), assign_dependent(data_for_sorts)],
    ignore_index=True,
)

monthly_char = (
    assignments.groupby(
        ["sorting_method", "date", "portfolio_size", "portfolio_bm"], as_index=False
    )
    .agg(n_stocks=("permno", "size"), mktcap=("mktcap_lag", "sum"))
)

portfolio_characteristics = (
    monthly_char.groupby(
        ["sorting_method", "portfolio_size", "portfolio_bm"], as_index=False
    )
    .agg(n_stocks=("n_stocks", "mean"), mktcap=("mktcap", "mean"))
)
portfolio_characteristics["mktcap_share"] = portfolio_characteristics[
    "mktcap"
] / portfolio_characteristics.groupby("sorting_method")["mktcap"].transform("sum")

display(portfolio_characteristics.head(6))

	sorting_method	portfolio_size	portfolio_bm	n_stocks	mktcap	mktcap_share
0	Dependent	1	1	106.258278	1.883282e+08	0.006357
1	Dependent	1	2	129.678808	2.398458e+08	0.008097
2	Dependent	1	3	132.099338	2.389889e+08	0.008068
3	Dependent	1	4	154.298013	2.474387e+08	0.008353
4	Dependent	1	5	144.781457	2.081802e+08	0.007028
5	Dependent	2	1	108.443709	3.549448e+08	0.011982

def plot_heatmap(df: pd.DataFrame, value_col: str, title: str, fmt) -> None:
    methods = ["Independent", "Dependent"]
    fig, axes = plt.subplots(1, 2, figsize=(11, 4.5), sharey=True)
    for ax, method in zip(axes, methods):
        sub = df[df["sorting_method"] == method]
        mat = sub.pivot(
            index="portfolio_bm", columns="portfolio_size", values=value_col
        ).sort_index(ascending=True)
        # Display: size 1→5 left to right; keep matrix rows = BM
        im = ax.imshow(mat.values, aspect="auto", origin="lower", cmap="Blues")
        ax.set_xticks(range(mat.shape[1]))
        ax.set_yticks(range(mat.shape[0]))
        ax.set_xticklabels(mat.columns.tolist)
        ax.set_yticklabels(mat.index.tolist)
        ax.set_xlabel("Size portfolio")
        ax.set_ylabel("Book-to-market portfolio")
        ax.set_title(method)
        for i in range(mat.shape[0]):
            for j in range(mat.shape[1]):
                val = mat.values[i, j]
                if np.isnan(val):
                    continue
                ax.text(j, i, fmt(val), ha="center", va="center", color="black", fontsize=8)
        fig.colorbar(im, ax=ax, fraction=0.046, pad=0.04)
    fig.suptitle(title, y=1.02)
    plt.tight_layout
    plt.show

plot_heatmap(
    portfolio_characteristics,
    "n_stocks",
    "Average number of stocks per portfolio (A-share)",
    lambda v: f"{v:.0f}",
)
plot_heatmap(
    portfolio_characteristics,
    "mktcap_share",
    "Market-cap share per portfolio (A-share)",
    lambda v: f"{100 * v:.1f}%",
)

Interpretation (same structure as the original, in an A-share setting):

• With SSE-Main breakpoints, small-size groups often absorb many SZSE / ChiNext / STAR names, so stock counts tilt toward the small-size columns.
• Independent sorts reuse one BM breakpoint schedule across all size groups, so BM-group counts within a size column can be uneven; dependent sorts recompute BM breakpoints within size groups, evening out the BM dimension.
• Market-cap share usually concentrates in the large-size × low-BM (growth) corner: a few mega-caps dominate value-weighted returns—a reminder to look at composition alongside returns.

6.10.6 Key Takeaways

• Bivariate sorts build portfolios from two characteristics (here Size and BM)—the empirical basis of Fama–French-style value/size factors.
• Independent sorts assign groups separately then cross; dependent sorts sort the second dimension inside first-dimension groups.
• Accounting data must be lagged and aligned; a 6-month BM lag plus a 12-month validity window is a common conservative look-ahead control.
• Value premia are sensitive to independent/dependent design, breakpoint samples, and weighting—continuing the non-standard-errors narrative from 6.9.
• Composition heatmaps show stock count ≠ market-cap share; value-weighting is dominated by large-cap growth names.

Exercises (adapted from the original for A-shares)

1. Compute breakpoints with “all major boards” vs “SSE Main only” and compare average stock counts in each Size×BM cell.
2. Following 6.9, plot the cross-sectional distribution of book equity and BM (or the time series of decile breakpoints).
3. Write a function that reports within-cell means of Size and BM for all 25 independent/dependent portfolios (characteristic profiles, not just returns).
4. Mimic a common HML convention (e.g., 2×3 sorts) in a small specification grid and see how the value premium moves with design choices (echoing 6.9’s p-hacking).

6.11 A-Share Replication: Replicating Fama–French Factors

This section replicates the tidyfinance chapter Replicating Fama-French Factors: following the official FF protocol (annual June sorts, hold from July for one year; Size×BM 2×3 sorts) construct SMB / HML, then evaluate replication quality against CSMAR’s published monthly three-factor series. If profitability and investment variables are cached locally, extend to five-factor RMW / CMA.

import warnings
from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import statsmodels.api as sm
from IPython.display import display
from matplotlib import font_manager

warnings.filterwarnings("ignore", category=UserWarning)

_cn_candidates = ["Microsoft YaHei", "SimHei", "PingFang SC", "Noto Sans CJK SC", "Arial Unicode MS"]
_available = {f.name for f in font_manager.fontManager.ttflist}
for _name in _cn_candidates:
    if _name in _available:
        plt.rcParams["font.sans-serif"] = [_name]
        break
plt.rcParams["axes.unicode_minus"] = False

DATA_DIR = Path(r"D:/A_Topics/202607_03_TidyFinanceAShare/202608_02_传统量化/Data")
CACHE_DIR = DATA_DIR / "_cache_beta"
BE_CACHE = CACHE_DIR / "book_equity_ashare.parquet"
ANN_CACHE = CACHE_DIR / "compustat_annual_ashare.parquet"  # be/at/op/inv (for FF5)
FF3_XLSX = DATA_DIR / "Fama_French_Factor_Monthly.xlsx"

# CSMAR: SH+SZ A-shares + ChiNext + STAR; free-float value-weighted
CSMAR_MARKET = "P9714"

print("crsp =", (CACHE_DIR / "crsp_monthly_ashare.parquet").exists)
print("BE   =", BE_CACHE.exists)
print("FF3  =", FF3_XLSX.exists)
print("ANN  =", ANN_CACHE.exists, "(optional for FF5)")

crsp = True
BE   = True
FF3  = True
ANN  = True (optional for FF5)

6.11.1 Data Preparation

Corresponds to Data Preparation. Load the monthly panel, annual book equity, and CSMAR three factors; board mapping is used for main-board breakpoints.

def map_board(permno: str) -> str:
    code = str(permno).zfill(6)
    if code.startswith("688"):
        return "STAR"
    if code.startswith(("300", "301")):
        return "ChiNext"
    if code.startswith("6"):
        return "SSE Main"
    if code.startswith(("000", "001", "002", "003")):
        return "SZSE Main"
    return "Other"

def load_csmar_ff3(path: Path, market: str = CSMAR_MARKET) -> pd.DataFrame:
    """Load CSMAR monthly FF3; keep float-weighted factors for one market."""
    raw = pd.read_excel(path)
    # First two rows: Chinese names / units
    df = raw.iloc[2:].copy
    df.columns = ["MarkettypeID", "TradingMonth", "mkt_excess", "mkt_excess_total",
                  "smb", "smb_total", "hml", "hml_total"]
    df = df[df["MarkettypeID"].astype(str) == market].copy
    df["date"] = pd.to_datetime(df["TradingMonth"].astype(str) + "-01")
    for c in ["mkt_excess", "smb", "hml"]:
        df[c] = pd.to_numeric(df[c], errors="coerce")
    out = df[["date", "mkt_excess", "smb", "hml"]].dropna(subset=["smb", "hml"]).sort_values("date")
    return out.reset_index(drop=True)

crsp_monthly = pd.read_parquet(CACHE_DIR / "crsp_monthly_ashare.parquet").copy
crsp_monthly["date"] = pd.to_datetime(crsp_monthly["date"])
crsp_monthly["permno"] = crsp_monthly["permno"].astype(str).str.zfill(6)
crsp_monthly["board"] = crsp_monthly["permno"].map(map_board)
crsp_monthly = crsp_monthly.dropna(subset=["ret_excess", "mktcap", "mktcap_lag"]).copy
crsp_monthly = crsp_monthly[crsp_monthly["mktcap_lag"] > 0].copy

book_equity = pd.read_parquet(BE_CACHE).copy
book_equity["permno"] = book_equity["permno"].astype(str).str.zfill(6)
book_equity["datadate"] = pd.to_datetime(book_equity["datadate"])
book_equity = book_equity[book_equity["be"] > 0].copy

factors_ff3_monthly = load_csmar_ff3(FF3_XLSX, CSMAR_MARKET)

print("crsp:", crsp_monthly.shape)
print("BE years:", book_equity["datadate"].dt.year.min, "-", book_equity["datadate"].dt.year.max)
print("CSMAR FF3:", factors_ff3_monthly["date"].min.date, "→", factors_ff3_monthly["date"].max.date,
      "| market=", CSMAR_MARKET)
display(factors_ff3_monthly.tail(3))

crsp: (884560, 11)
BE years: 1999 - 2025
CSMAR FF3: 1991-07-01 → 2026-07-01 | market= P9714

	date	mkt_excess	smb	hml
416	2026-05-01	-0.003817	-0.009226	-0.022075
417	2026-06-01	0.002287	-0.008796	-0.035981
418	2026-07-01	-0.082799	0.004447	0.204287

Unlike 6.10, the official FF timeline is:

1. Size: use June market cap in year $t$, hold from July through next June;
2. BM: use year $t-1$ book equity ÷ December year $t-1$ market cap, also put into use in July of year $t$ (accounting information can lag up to about 18 months).

Below, sorting_date (each July) aligns size and BM onto one annual grouping key.

# June market cap → July sorting_date
size = (
    crsp_monthly.loc[crsp_monthly["date"].dt.month == 6, ["permno", "board", "date", "mktcap"]]
    .assign(sorting_date=lambda d: d["date"] + pd.DateOffset(months=1))
    .rename(columns={"mktcap": "size"})
    [["permno", "board", "sorting_date", "size"]]
)

# December market cap → +7 months = next July (align with annual BE)
market_equity = (
    crsp_monthly.loc[crsp_monthly["date"].dt.month == 12, ["permno", "date", "mktcap"]]
    .assign(sorting_date=lambda d: d["date"] + pd.DateOffset(months=7))
    .rename(columns={"mktcap": "me"})
    [["permno", "sorting_date", "me"]]
)

book_to_market = book_equity.copy
book_to_market["sorting_date"] = pd.to_datetime(
    dict(
        year=book_to_market["datadate"].dt.year + 1,
        month=7,
        day=1,
    )
)
book_to_market = book_to_market.merge(market_equity, on=["permno", "sorting_date"], how="inner")
book_to_market["bm"] = book_to_market["be"] / book_to_market["me"]
book_to_market = book_to_market.replace([np.inf, -np.inf], np.nan)
book_to_market = book_to_market.dropna(subset=["bm"])
book_to_market = book_to_market[book_to_market["bm"] > 0]

sorting_variables = (
    size.merge(
        book_to_market[["permno", "sorting_date", "me", "bm"]],
        on=["permno", "sorting_date"],
        how="inner",
    )
    .dropna
    .drop_duplicates(["permno", "sorting_date"], keep="first")
)

print(sorting_variables.shape)
print(
    "sorting years:",
    sorting_variables["sorting_date"].dt.year.min,
    "-",
    sorting_variables["sorting_date"].dt.year.max,
)
display(sorting_variables.head(3))

(63650, 6)
sorting years: 2001 - 2026

	permno	board	sorting_date	size	me	bm
0	000001	深交所主板	2001-07-01	21328740.14	20228171.57	0.173894
1	000001	深交所主板	2002-07-01	21140429.48	17264684.07	0.210121
2	000001	深交所主板	2003-07-01	15629824.19	14784207.01	0.266529

6.11.2 Portfolio Sorts (2×3)

Corresponds to Portfolio Sorts. Size splits into Small / Big at the main-board median; BM into Low / Mid / High at main-board 30% / 70% percentiles. Crossing yields six portfolios.

def assign_portfolio(
    data: pd.DataFrame,
    sorting_variable: str,
    percentiles: list[float],
    boards: list[str] | None = None,
) -> pd.Series:
    if boards is None:
        boards = ["SSE Main"]
    subset = data[data["board"].isin(boards)]
    if len(subset) < max(2, len(percentiles) - 1):
        return pd.Series(1, index=data.index, dtype=int)

    breakpoints = subset[sorting_variable].quantile(percentiles).to_numpy(dtype=float)
    breakpoints = np.unique(breakpoints)
    if len(breakpoints) < 2:
        return pd.Series(1, index=data.index, dtype=int)

    inner = breakpoints[1:-1]
    labels = list(range(1, len(inner) + 2))
    try:
        port = pd.cut(
            data[sorting_variable],
            bins=[-np.inf, *inner, np.inf],
            labels=labels,
            include_lowest=True,
        )
    except ValueError:
        return pd.Series(1, index=data.index, dtype=int)
    return port.astype(int)

parts = []
for dt, g in sorting_variables.groupby("sorting_date", sort=True):
    g = g.copy
    g["portfolio_size"] = assign_portfolio(g, "size", [0, 0.5, 1])
    g["portfolio_bm"] = assign_portfolio(g, "bm", [0, 0.3, 0.7, 1])
    parts.append(g[["permno", "sorting_date", "portfolio_size", "portfolio_bm"]])
portfolios_annual = pd.concat(parts, ignore_index=True)

# Monthly returns: Jan–Jun use prior July groups; Jul–Dec use current July groups
crsp = crsp_monthly.copy
y = crsp["date"].dt.year
m = crsp["date"].dt.month
crsp["sorting_date"] = pd.to_datetime(
    dict(year=np.where(m <= 6, y - 1, y), month=7, day=1)
)
portfolios = crsp.merge(portfolios_annual, on=["permno", "sorting_date"], how="inner")
print("monthly portfolio-months:", portfolios.shape)
print(portfolios.groupby(["portfolio_size", "portfolio_bm"]).size.unstack(fill_value=0))

monthly portfolio-months: (689085, 14)
portfolio_bm         1       2       3
portfolio_size                        
1               122372  179872  109645
2               106402  101274   69520

6.11.3 Three Factors: SMB and HML

Corresponds to Fama-French Three-Factor Model. Within each 2×3 cell, value-weight by mktcap_lag;
SMB = equal-weighted average of the three small portfolios − equal-weighted average of the three big portfolios;
HML = equal-weighted average of the two high-BM portfolios − equal-weighted average of the two low-BM portfolios.

def vw_ret(g: pd.DataFrame) -> float:
    w = g["mktcap_lag"]
    return float((g["ret_excess"] * w).sum / w.sum) if w.sum > 0 else np.nan

rows = []
for keys, g in portfolios.groupby(["date", "portfolio_size", "portfolio_bm"], sort=True):
    rows.append(
        {
            "date": keys[0],
            "portfolio_size": int(keys[1]),
            "portfolio_bm": int(keys[2]),
            "ret": vw_ret(g),
        }
    )
port_rets = pd.DataFrame(rows).dropna(subset=["ret"])

factors_replicated = []
for dt, g in port_rets.groupby("date"):
    small = g.loc[g["portfolio_size"] == 1, "ret"].mean
    big = g.loc[g["portfolio_size"] == 2, "ret"].mean
    high = g.loc[g["portfolio_bm"] == g["portfolio_bm"].max, "ret"].mean
    low = g.loc[g["portfolio_bm"] == g["portfolio_bm"].min, "ret"].mean
    factors_replicated.append(
        {"date": dt, "smb_replicated": small - big, "hml_replicated": high - low}
    )
factors_replicated = pd.DataFrame(factors_replicated).sort_values("date")

print(
    factors_replicated[["smb_replicated", "hml_replicated"]]
    .agg(["mean", "std", "count"])
    .round(4)
)
display(factors_replicated.tail(3))

       smb_replicated  hml_replicated
mean          -0.0027         -0.0018
std            0.0434          0.0308
count        318.0000        318.0000

	date	smb_replicated	hml_replicated
315	2026-05-01	0.005685	0.000339
316	2026-06-01	0.023147	-0.039969
317	2026-07-01	-0.140000	0.111628

6.11.4 Replication Evaluation (vs CSMAR)

Corresponds to Replication Evaluation: time-series regress home-built factors on CSMAR official factors
replicated ~ α + β · official. Ideally α≈0, β≈1, and adjusted $R^2$ near 1.

def replication_regression(y: pd.Series, x: pd.Series, name: str) -> dict:
    df = pd.concat([y, x], axis=1).dropna
    df.columns = ["y", "x"]
    if len(df) < 24:
        return {"factor": name, "n": len(df), "alpha": np.nan, "beta": np.nan, "adj_r2": np.nan}
    model = sm.OLS(df["y"], sm.add_constant(df["x"])).fit
    return {
        "factor": name,
        "n": int(model.nobs),
        "alpha": float(model.params["const"]),
        "alpha_t": float(model.tvalues["const"]),
        "beta": float(model.params["x"]),
        "beta_t": float(model.tvalues["x"]),
        "adj_r2": float(model.rsquared_adj),
        "corr": float(df["y"].corr(df["x"])),
    }

test = factors_replicated.merge(factors_ff3_monthly, on="date", how="inner")
# Round to 4 decimals to reduce float noise (same as original)
test["smb_replicated"] = test["smb_replicated"].round(4)
test["hml_replicated"] = test["hml_replicated"].round(4)

rows = [
    replication_regression(test["smb_replicated"], test["smb"], "SMB"),
    replication_regression(test["hml_replicated"], test["hml"], "HML"),
]
eval_df = pd.DataFrame(rows)
display(eval_df.round(4))

fig, axes = plt.subplots(1, 2, figsize=(10, 4.2))
for ax, lab, ycol, xcol in [
    (axes[0], "SMB", "smb_replicated", "smb"),
    (axes[1], "HML", "hml_replicated", "hml"),
]:
    ax.scatter(test[xcol], test[ycol], s=12, alpha=0.55, color="#4C78A8")
    lims = [
        np.nanmin([test[xcol].min, test[ycol].min]),
        np.nanmax([test[xcol].max, test[ycol].max]),
    ]
    ax.plot(lims, lims, "r--", linewidth=1, label="45°")
    ax.set_xlabel(f"CSMAR {lab}")
    ax.set_ylabel(f"Replicated {lab}")
    ax.set_title(lab)
    ax.grid(True, alpha=0.3)
    ax.legend(frameon=False)
plt.suptitle(f"A-share FF3 replication vs CSMAR ({CSMAR_MARKET})", y=1.02)
plt.tight_layout
plt.show

print(
    f"overlap months = {len(test)}; "
    f"SMB corr = {test['smb_replicated'].corr(test['smb']):.3f}; "
    f"HML corr = {test['hml_replicated'].corr(test['hml']):.3f}"
)

	factor	n	alpha	alpha_t	beta	beta_t	adj_r2	corr
0	SMB	318	-0.0080	-9.7005	0.8133	50.0896	0.8878	0.9424
1	HML	318	-0.0036	-3.6164	0.7315	25.6705	0.6749	0.8221

overlap months = 318; SMB corr = 0.942; HML corr = 0.822

If β is far from 1 or $R^2$ is low, common causes include: stock universe mismatch with CSMAR MarkettypeID, ST/financial-industry filter differences, market-cap definition (free-float vs total), and BE definition / filing-date differences. Try codes such as P9706 (main boards only) for robustness.

6.11.5 Five-Factor Extension: RMW and CMA

Corresponds to Fama-French Five-Factor Model. Within Size groups, further split operating profitability (OP) and investment (INV) at 30/70 percentiles to form:

• RMW = high profitability − low profitability;
• CMA = low investment − high investment;
• Five-factor SMB often averages small-minus-big across Size×BM / Size×OP / Size×INV grids.

Requires annual variables op and inv (cache compustat_annual_ashare.parquet). If the cache is missing, the cell below skips; you can also fill via akshare income/balance sheets. This CSMAR file has only three factors, so five-factor results are not matched to an official series.

def load_annual_fundamentals(path: Path) -> pd.DataFrame | None:
    if not path.exists:
        print("Not found:", path, "→ skip five-factor construction.")
        print("Run the Data download script to build be/at/op/inv, then re-run this subsection.")
        return None
    ann = pd.read_parquet(path).copy
    ann["permno"] = ann["permno"].astype(str).str.zfill(6)
    ann["datadate"] = pd.to_datetime(ann["datadate"])
    need = {"be", "op", "inv"}
    if not need.issubset(ann.columns):
        print("Cache missing columns", need - set(ann.columns), "→ skip five-factor.")
        return None
    return ann

ann = load_annual_fundamentals(ANN_CACHE)
if ann is not None:
    other = ann.copy
    other["sorting_date"] = pd.to_datetime(
        dict(year=other["datadate"].dt.year + 1, month=7, day=1)
    )
    other = other.merge(market_equity, on=["permno", "sorting_date"], how="inner")
    other["bm"] = other["be"] / other["me"]
    other = other.replace([np.inf, -np.inf], np.nan)

    sorting5 = (
        size.merge(
            other[["permno", "sorting_date", "me", "be", "bm", "op", "inv"]],
            on=["permno", "sorting_date"],
            how="inner",
        )
        .dropna(subset=["size", "bm", "op", "inv"])
        .drop_duplicates(["permno", "sorting_date"], keep="first")
    )
    # Dependent sorts: Size first, then BM/OP/INV within size
    parts = []
    for dt, g in sorting5.groupby("sorting_date", sort=True):
        g = g.copy
        g["portfolio_size"] = assign_portfolio(g, "size", [0, 0.5, 1])
        sub = []
        for _, sg in g.groupby("portfolio_size"):
            sg = sg.copy
            sg["portfolio_bm"] = assign_portfolio(sg, "bm", [0, 0.3, 0.7, 1])
            sg["portfolio_op"] = assign_portfolio(sg, "op", [0, 0.3, 0.7, 1])
            sg["portfolio_inv"] = assign_portfolio(sg, "inv", [0, 0.3, 0.7, 1])
            sub.append(sg)
        parts.append(pd.concat(sub, ignore_index=False))
    port_ann5 = pd.concat(parts, ignore_index=True)[
        ["permno", "sorting_date", "portfolio_size", "portfolio_bm", "portfolio_op", "portfolio_inv"]
    ]
    port5 = crsp.merge(port_ann5, on=["permno", "sorting_date"], how="inner")

    def factor_from_grid(df, char_col, high_minus_low=True):
        rows = []
        for keys, g in df.groupby(["date", "portfolio_size", char_col], sort=True):
            rows.append(
                {
                    "date": keys[0],
                    "portfolio_size": int(keys[1]),
                    char_col: int(keys[2]),
                    "ret": vw_ret(g),
                }
            )
        rets = pd.DataFrame(rows).dropna(subset=["ret"])
        out = []
        for dt, g in rets.groupby("date"):
            hi = g.loc[g[char_col] == g[char_col].max, "ret"].mean
            lo = g.loc[g[char_col] == g[char_col].min, "ret"].mean
            out.append({"date": dt, "ret": (hi - lo) if high_minus_low else (lo - hi)})
        return pd.DataFrame(out)

    hml5 = factor_from_grid(port5, "portfolio_bm", True).rename(columns={"ret": "hml_replicated"})
    rmw = factor_from_grid(port5, "portfolio_op", True).rename(columns={"ret": "rmw_replicated"})
    cma = factor_from_grid(port5, "portfolio_inv", False).rename(columns={"ret": "cma_replicated"})

    # SMB: average small-minus-big across three size grids
    smb_parts = []
    for col in ["portfolio_bm", "portfolio_op", "portfolio_inv"]:
        rows = []
        for keys, g in port5.groupby(["date", "portfolio_size", col], sort=True):
            rows.append(
                {
                    "date": keys[0],
                    "portfolio_size": int(keys[1]),
                    "ret": vw_ret(g),
                }
            )
        rets = pd.DataFrame(rows).dropna(subset=["ret"])
        for dt, g in rets.groupby("date"):
            s = g.loc[g["portfolio_size"] == 1, "ret"].mean
            b = g.loc[g["portfolio_size"] == 2, "ret"].mean
            smb_parts.append({"date": dt, "ret": s - b})
    smb5 = (
        pd.DataFrame(smb_parts).groupby("date", as_index=False)["ret"].mean
        .rename(columns={"ret": "smb_replicated"})
    )

    factors5 = (
        smb5.merge(hml5, on="date", how="outer")
        .merge(rmw, on="date", how="outer")
        .merge(cma, on="date", how="outer")
        .sort_values("date")
    )
    print("FF5 replicated — monthly means:")
    display(factors5.filter(like="_replicated").agg(["mean", "std", "count"]).T.round(4))
    display(factors5.tail(3))

FF5 replicated — monthly means:

	mean	std	count
smb_replicated	-0.0021	0.0441	318.0
hml_replicated	-0.0019	0.0311	318.0
rmw_replicated	0.0022	0.0241	318.0
cma_replicated	-0.0029	0.0163	318.0

	date	smb_replicated	hml_replicated	rmw_replicated	cma_replicated
315	2026-05-01	0.002664	-0.003790	0.032230	-0.026152
316	2026-06-01	0.015789	-0.045325	-0.012644	-0.055543
317	2026-07-01	-0.119145	0.112453	0.058547	0.047182

6.11.6 Key Takeaways

• Beyond the CAPM market factor, FF3 adds SMB (size) and HML (value); FF5 further adds RMW (profitability) and CMA (investment).
• Official construction emphasizes: June/July annual rebalancing, specific accounting lags, and (U.S.) NYSE / (A-share) main-board breakpoints with 2×3 sorts.
• Replication quality is judged by time-series regressions on official series: focus on α, β, and adjusted $R^2$.
• Sections 6.9/6.10 show sort degrees of freedom; this section shows how locking one FF protocol yields standard factors that can be matched.

Exercises (adapted from the original for A-shares)

1. Replicate RiskPremium1 with a home-built value-weighted market excess return and evaluate against CSMAR mkt_excess with the same regression.
2. Compare P9706 (no ChiNext/STAR) vs P9714 (includes them) for SMB/HML β and $R^2$.
3. Switch breakpoints from “SSE Main only” to full-market quantiles and see how replication correlations change.
4. After the cache is complete, finish five factors and match to CSMAR five-factor data if available.

6.12 A-Share Replication: Fama–MacBeth Regressions

This section replicates the tidyfinance chapter Fama-MacBeth Regressions: treating individual stocks as test assets, estimate risk premia on characteristics linked to the Fama–French three factors (market beta, log market cap, book-to-market). Conceptual background is in Section 6.6; here we give a full A-share pipeline isomorphic to the original.

Fama–MacBeth is a two-step procedure:

1. Cross-sectional regression: each month $t$, explain next-month excess returns $r_{i,t+1}$ with characteristics known at $t$;
2. Time-series aggregation: average the monthly $\hat\lambda^f_t$ and run $t$-tests (often with Newey–West).

Linear sketch (characteristics as proxies for risk exposures):

$$ r_{i,t+1}=\alpha_t+\lambda^{M}_t\beta^{M}_{i,t}+\lambda^{\mathrm{Size}}_t\log(\mathrm{ME})_{i,t}+\lambda^{\mathrm{BM}}_t\mathrm{BM}_{i,t}+\epsilon_{i,t+1}. $$

import warnings
from pathlib import Path

import numpy as np
import pandas as pd
import statsmodels.api as sm
from IPython.display import display

warnings.filterwarnings("ignore", category=UserWarning)

DATA_DIR = Path(r"D:/A_Topics/202607_03_TidyFinanceAShare/202608_02_传统量化/Data")
CACHE_DIR = DATA_DIR / "_cache_beta"

print("crsp =", (CACHE_DIR / "crsp_monthly_ashare.parquet").exists)
print("BE   =", (CACHE_DIR / "book_equity_ashare.parquet").exists)
print("beta =", (DATA_DIR / "beta_ashare.parquet").exists)

crsp = True
BE   = True
beta = True

6.12.1 Data Preparation

Corresponds to Data Preparation. Load monthly returns and market cap, positive book equity, and monthly CAPM beta.

crsp_monthly = pd.read_parquet(
    CACHE_DIR / "crsp_monthly_ashare.parquet",
    columns=["permno", "date", "ret_excess", "mktcap"],
).copy
crsp_monthly["date"] = pd.to_datetime(crsp_monthly["date"])
crsp_monthly["permno"] = crsp_monthly["permno"].astype(str).str.zfill(6)
crsp_monthly = crsp_monthly.dropna(subset=["ret_excess", "mktcap"])
crsp_monthly = crsp_monthly[crsp_monthly["mktcap"] > 0].copy

compustat_annual = pd.read_parquet(CACHE_DIR / "book_equity_ashare.parquet").copy
compustat_annual["permno"] = compustat_annual["permno"].astype(str).str.zfill(6)
compustat_annual["datadate"] = pd.to_datetime(compustat_annual["datadate"])
compustat_annual = compustat_annual[compustat_annual["be"] > 0].copy

beta = pd.read_parquet(DATA_DIR / "beta_ashare.parquet").copy
beta["permno"] = beta["permno"].astype(str).str.zfill(6)
beta["date"] = pd.to_datetime(beta["date"])
beta = beta.loc[beta["return_type"] == "monthly", ["permno", "date", "beta"]].copy

print(crsp_monthly.shape, compustat_annual.shape, beta.shape)
display(crsp_monthly.head(2))

(890537, 4) (88343, 3) (623256, 3)

	permno	date	ret_excess	mktcap
0	000001	2000-01-01	0.060035	19843822.88
1	000001	2000-02-01	-0.013189	19618933.36

Construct BM and $\log(\mathrm{ME})$ from annual BE and same-month market cap, merge contemporaneous beta; lag characteristics by 6 months into the return panel and forward-fill by stock. Finally lead returns by one period to get ret_excess_lead (period-$t$ characteristics forecast period-$t+1$ returns).

# Align annual report dates to month (December)
compustat_annual = compustat_annual.assign(
    date=lambda d: d["datadate"].dt.to_period("M").dt.to_timestamp
)

characteristics = (
    compustat_annual.merge(
        crsp_monthly[["permno", "date", "mktcap"]],
        on=["permno", "date"],
        how="left",
    )
    .merge(beta, on=["permno", "date"], how="left")
    .assign(
        bm=lambda d: d["be"] / d["mktcap"],
        log_mktcap=lambda d: np.log(d["mktcap"]),
        sorting_date=lambda d: d["date"] + pd.DateOffset(months=6),
    )
    [["permno", "bm", "log_mktcap", "beta", "sorting_date"]]
)
characteristics = characteristics.replace([np.inf, -np.inf], np.nan)

data_fama_macbeth = (
    crsp_monthly.merge(
        characteristics,
        left_on=["permno", "date"],
        right_on=["permno", "sorting_date"],
        how="left",
    )
    .drop(columns=["sorting_date"], errors="ignore")
    .sort_values(["permno", "date"])
)
for col in ["beta", "bm", "log_mktcap"]:
    data_fama_macbeth[col] = data_fama_macbeth.groupby("permno")[col].ffill

# Lead returns: shift date back 1 month so current characteristics match next-month returns
lead = data_fama_macbeth[["permno", "date", "ret_excess"]].copy
lead["date"] = lead["date"] - pd.DateOffset(months=1)
lead = lead.rename(columns={"ret_excess": "ret_excess_lead"})

data_fama_macbeth = (
    data_fama_macbeth.merge(lead, on=["permno", "date"], how="left")
    [["permno", "date", "ret_excess_lead", "beta", "log_mktcap", "bm"]]
    .dropna
    .copy
)

print(data_fama_macbeth.shape)
print(
    "date range:",
    data_fama_macbeth["date"].min.date,
    "→",
    data_fama_macbeth["date"].max.date,
)
display(data_fama_macbeth.head(3))

(502694, 6)
date range: 2004-06-01 → 2026-06-01

	permno	date	ret_excess_lead	beta	log_mktcap	bm
53	000001	2004-06-01	-0.058612	0.945177	16.29989	0.366434
54	000001	2004-07-01	0.004530	0.945177	16.29989	0.366434
55	000001	2004-08-01	-0.001635	0.945177	16.29989	0.366434

6.12.2 Cross-Sectional Regression

Corresponds to Cross-Sectional Regression. Each month estimate on the full cross-section

ret_excess_lead ~ beta + log_mktcap + bm

to obtain that month’s risk-premium estimates $\hat\lambda_t$.

def estimate_cross_section(group: pd.DataFrame) -> pd.Series | None:
    # Skip if sample too small or collinear
    if len(group) < 50:
        return None
    y = group["ret_excess_lead"]
    X = sm.add_constant(group[["beta", "log_mktcap", "bm"]], has_constant="add")
    try:
        model = sm.OLS(y, X).fit
    except Exception:  # noqa: BLE001
        return None
    out = model.params.rename({"const": "Intercept"})
    out["date"] = group["date"].iloc[0]
    return out

rows = []
for _, g in data_fama_macbeth.groupby("date", sort=True):
    res = estimate_cross_section(g)
    if res is not None:
        rows.append(res)

risk_premiums = pd.DataFrame(rows).sort_values("date").reset_index(drop=True)
risk_premiums = risk_premiums[["date", "Intercept", "beta", "log_mktcap", "bm"]]
print(risk_premiums.shape)
display(risk_premiums.tail(3))

(265, 5)

	date	Intercept	beta	log_mktcap	bm
262	2026-04-01	-0.173627	-0.010647	0.009697	-0.006927
263	2026-05-01	-0.325958	0.029290	0.015854	-0.015826
264	2026-06-01	-0.066273	-0.111225	0.007892	0.031225

6.12.3 Time-Series Aggregation

Corresponds to Time-Series Aggregation. Average $\hat\lambda_t$ over time as the risk premium and compute ordinary $t$-statistics:

[ t=\frac{\overline{\hat\lambda}}{\mathrm{sd}(\hat\lambda)/\sqrt{T}}. ]

long = risk_premiums.melt(
    id_vars="date", var_name="factor", value_name="estimate"
)

price_of_risk = (
    long.groupby("factor", as_index=False)
    .agg(
        risk_premium=("estimate", "mean"),
        t_statistic=(
            "estimate",
            lambda s: s.mean / s.std(ddof=1) * np.sqrt(len(s)),
        ),
    )
    .sort_values("factor")
)
display(price_of_risk.round(4))

	factor	risk_premium	t_statistic
0	Intercept	0.0565	2.9072
1	beta	-0.0018	-1.1494
2	bm	-0.0003	-0.5250
3	log_mktcap	-0.0027	-2.4319

When reporting risk premia, practitioners often apply Newey–West (HAC) to the $\hat\lambda_t$ series for autocorrelation. Below we regress on a constant with HAC standard errors (lag 6) to get NW $t$-stats and tabulate them beside ordinary $t$—matching the original’s pyfixest / tidyfinance approach.

def estimate_newey_west(group: pd.DataFrame, lags: int = 6) -> pd.Series:
    g = group.sort_values("date")
    y = g["estimate"].to_numpy(dtype=float)
    X = np.ones((len(y), 1))
    fit = sm.OLS(y, X).fit(cov_type="HAC", cov_kwds={"maxlags": lags})
    return pd.Series(
        {
            "factor": g["factor"].iloc[0],
            "t_statistic_newey_west": float(y.mean / fit.bse[0]),
        }
    )

nw_rows = [estimate_newey_west(g) for _, g in long.groupby("factor")]
price_of_risk_newey_west = pd.DataFrame(nw_rows)

fm_table = (
    price_of_risk.merge(price_of_risk_newey_west, on="factor")
    .assign(
        risk_premium=lambda d: d["risk_premium"].round(3),
        t_statistic=lambda d: d["t_statistic"].round(3),
        t_statistic_newey_west=lambda d: d["t_statistic_newey_west"].round(3),
    )
)
display(fm_table)

	factor	risk_premium	t_statistic	t_statistic_newey_west
0	Intercept	0.057	2.907	3.300
1	beta	-0.002	-1.149	-1.146
2	bm	-0.000	-0.525	-0.672
3	log_mktcap	-0.003	-2.432	-2.824

6.12.4 Interpreting Results and a One-Call Wrapper

Interpretation matches the original (signs depend on the sample):

• BM coefficient positive → value premium (high BM, higher expected returns);
• log_mktcap negative → size premium (small caps, higher expected returns);
• beta often insignificant or weak → consistent with earlier univariate beta sorts and a flat SML.

Below is a helper with the same role as tidyfinance.estimate_fama_macbeth (default reports Newey–West $t$).

def estimate_fama_macbeth(
    data: pd.DataFrame,
    model: str = "ret_excess_lead ~ beta + bm + log_mktcap",
    vcov_lags: int = 6,
) -> pd.DataFrame:
    """Two-step Fama-MacBeth with HAC t-stats (A-share helper)."""
    # Parse a minimal formula: y ~ x1 + x2 + ...
    left, right = [s.strip for s in model.split("~")]
    x_cols = [c.strip for c in right.split("+")]
    y_col = left

    prem_rows = []
    for dt, g in data.groupby("date", sort=True):
        if len(g) < 50:
            continue
        y = g[y_col]
        X = sm.add_constant(g[x_cols], has_constant="add")
        try:
            fit = sm.OLS(y, X).fit
        except Exception:  # noqa: BLE001
            continue
        params = fit.params.rename({"const": "Intercept"})
        params["date"] = dt
        prem_rows.append(params)
    prem = pd.DataFrame(prem_rows)
    if prem.empty:
        return prem

    long_ = prem.melt(id_vars="date", var_name="factor", value_name="estimate")
    out_rows = []
    for fac, g in long_.groupby("factor"):
        y = g.sort_values("date")["estimate"].to_numpy(dtype=float)
        fit = sm.OLS(y, np.ones((len(y), 1))).fit(
            cov_type="HAC", cov_kwds={"maxlags": vcov_lags}
        )
        se = float(fit.bse[0])
        mu = float(y.mean)
        out_rows.append(
            {
                "factor": fac if fac != "Intercept" else "intercept",
                "risk_premium": mu,
                "n": len(y),
                "standard_error": se,
                "t_statistic": mu / se if se > 0 else np.nan,
            }
        )
    return pd.DataFrame(out_rows).sort_values("factor").reset_index(drop=True)

display(
    estimate_fama_macbeth(
        data_fama_macbeth,
        model="ret_excess_lead ~ beta + bm + log_mktcap",
        vcov_lags=6,
    ).round(4)
)

	factor	risk_premium	n	standard_error	t_statistic
0	beta	-0.0018	265	0.0016	-1.1460
1	bm	-0.0003	265	0.0005	-0.6723
2	intercept	0.0565	265	0.0171	3.3003
3	log_mktcap	-0.0027	265	0.0010	-2.8240

6.12.5 Key Takeaways

• Fama–MacBeth estimates characteristic risk premia by time-averaging cross-sectional regression coefficients, complementing portfolio sorts.
• Commonly used to test whether size, value, market beta, etc. are priced in the return cross-section.
• Pre-estimated betas have measurement error; bias shrinks in large samples. In practice, at least use Newey–West / HAC SEs for $\lambda_t$; stock-level regressions may also use value-weighted WLS (Hou et al., 2020).
• This section implements the two-step procedure with statsmodels and wraps it in estimate_fama_macbeth, matching the original’s manual workflow and tidyfinance one-call interface.

Exercises (adapted from the original for A-shares)

1. Estimate stock-level size/value loadings with rolling windows against home-built or CSMAR SMB/HML, replace characteristics log_mktcap/bm, re-run FM, and compare coefficients.
2. If industry-portfolio returns are available, estimate rolling beta/size/value exposures on industry portfolios instead of stocks, then run FM and check robustness.
3. Following Hou et al. (2020), run value-weighted WLS cross-sectional regressions; also run OLS after dropping the smallest 20% of firms each month. Compare risk premia across the three specs.