flowchart LR
prior["Previous beta β(i,t)"] --> fit["Fit today's cross-section"]
fit --> target["Use observed market return g(t)"]
target --> update["EWLS update with decay λ"]
update --> next["Freeze beta β(i,t+1)"]
next --> fit
Walkthrough
A ten-year point-in-time model run and evaluation
Model configuration
The evaluated run uses local marcap observations, a daily top-1,000 universe selected from the previous close, and four factors:
- MKT: recursive market beta against an observed market proxy.
- SIZE: lagged log market cap.
- HML: a 12–60 month long-term-reversal value proxy.
- MOM: a 3–9 month trend, measured three months before today.
Sequential WLS removes them in that order after 1% exposure winsorization, weighting each name by the square root of its lagged market cap. Style exposures are weighted and standardized before an orthogonalization step; MKT remains a beta exposure. The resulting ordered basis is jointly orthogonal, so the final residual is orthogonal to the complete factor span.
Factor definitions
Let \(r_{i,t}\) be a stock’s decimal return and \(C_{i,t}\) its market cap. All inputs are lagged to the start of the interval, so today’s exposure cannot use today’s close.
| Factor | Exposure formula | Meaning |
|---|---|---|
| MKT | \(\beta_{i,t}\) | Recursive beta against the observed lagged-cap-weighted market return. |
| SIZE | \(\log C_{i,t-1}\) | Larger companies have larger exposure. |
| HML | \(-R_{i,t}(1008,252)\) | Long-term reversal: the negative return from 12–60 months ago. |
| MOM | \(R_{i,t}(189,63)\) | Compounded return over roughly 3–9 months ago. |
Here \(R_{i,t}(L,S)\) is the compounded return from the \(L\) sessions ending \(S\) sessions before \(t\) (implemented with log1p, a shift by S, and a rolling window of L). HML negates this older return, making it a long-term reversal exposure.
Before each cross-sectional fit, exposures are winsorized and centered. SIZE, HML, and MOM are then standardized; MKT stays in beta units. Sequential output columns are the weighted orthogonal basis, not the correlated raw columns.
For the regression below, \(w_{i,t}\) is the positive weight for stock \(i\) on date \(t\). The main run uses the square root of lagged market cap. \(x_{j,i,t}\) is the normalized exposure, \(q_{j,i,t}\) is its weighted orthogonalized version, \(f_{j,t}\) is the return attributed to factor \(j\), and \(e^{(j)}_{i,t}\) is the residual after removing the first \(j\) factors. With \(j\) running in the configured order, the regression applies \[ q_{j}=x_{j}-\sum_{k<j} \frac{\langle x_j,q_k\rangle_w}{\langle q_k,q_k\rangle_w}q_k, \qquad e^{(0)}_{i,t}=r_{i,t}-\frac{\sum_i w_{i,t}r_{i,t}}{\sum_i w_{i,t}}, \] \[ f_{j,t}=\frac{\langle q_j,e^{(j-1)}\rangle_w} {\langle q_j,q_j\rangle_w},\qquad e^{(j)}_{i,t}=e^{(j-1)}_{i,t}-q_{j,i,t}f_{j,t}. \]
The factors are removed in the order MKT, SIZE, HML, MOM, so each one explains what remains after the earlier factors, while the orthogonalized basis makes the final residual orthogonal to all four factors. Changing that order changes the basis and the attribution, even though the fitted factor span is the same.
Recursive market beta
At the first usable date, each beta is bootstrapped from a historical covariance with the market return. Thereafter the previous beta is frozen for the current cross-section. Once the observed lagged-cap-weighted market return \(g_t\) is known, an exponentially weighted update produces the beta used on the next date. The fitted cross-sectional MKT slope is kept as an output, but it is not fed back as the time-series shock:
\[ D_{i,t}=\lambda D_{i,t-1}+g_t^2,\qquad N_{i,t}=\lambda N_{i,t-1}+g_t y_{i,t},\qquad \beta_{i,t+1}=\operatorname{clip}(N_{i,t}/D_{i,t}). \]
Here \(y_{i,t}=r_{i,t}-a_t+c_tg_t\), where \(a_t\) is the fitted intercept and \(c_t\) is the beta center. Only the intercept is stripped. This avoids feeding estimated, beta-correlated style returns back into the beta state, and the update applies to every name with a finite return and observed market shock, even when a long-window style excluded that name from the cross-section.
The loop is stateful across dates, but it uses the same cross-sectional fit as the batch pipeline. The flow is:
Show code
from pathlib import Path
from time import perf_counter
import numpy as np
import pandas as pd
import plotly.graph_objects as go
import resid
pd.set_option("display.max_columns", 12)
data_dir = Path("../../marcap/data")
source = resid.MarcapDataSource(data_dir)
base_factor_model = resid.CharacteristicFactorModel(
(
resid.size_factor(name="SIZE"),
resid.long_term_reversal_factor(
lookback_days=1008,
skip_days=252,
min_periods=756,
name="HML",
),
resid.momentum_factor(
lookback_days=189,
skip_days=63,
min_periods=126,
name="MOM",
),
)
)
factor_model = resid.RecursiveMarketBetaModel(
base=base_factor_model,
lookback_days=252,
min_periods=126,
decay=0.97,
name="MKT",
)
window = resid.analysis_window(source, years=10, end=None)
universe_builder = resid.DailyTopMarketCapUniverse(size=1000)
return_calculator = resid.PercentageReturns()
factor_order = ("MKT", "SIZE", "HML", "MOM")
residualizer = resid.SequentialWLSResidualizer(
factor_order=factor_order,
winsor_quantile=0.01,
unscaled_factors=("MKT",),
)
regression_weight_model = resid.SquareRootMarketCapWeights()
validations = (
resid.RegressionCoverageValidation(minimum_date_coverage=0.99),
resid.FiniteOutputValidation(),
resid.ReturnReconstructionValidation(absolute_tolerance=1e-12),
resid.ResidualOrthogonalityValidation(absolute_tolerance=1e-12),
resid.SequentialOrthogonalityValidation(absolute_tolerance=1e-12),
)
started = perf_counter()
result = resid.run_pipeline(
window=window,
source=source,
universe_builder=universe_builder,
return_calculator=return_calculator,
factor_model=factor_model,
residualizer=residualizer,
regression_weight_model=regression_weight_model,
validations=validations,
)
elapsed_seconds = perf_counter() - startedShow code
raw = result.returns.set_index(["date", "ticker"])["return"]
residual = result.specific_returns.set_index(["date", "ticker"])[
"specific_return"
]
diagnostics = result.diagnostics.set_index("date")
regression_weights = result.regression_weights.set_index(["date", "ticker"])[
"regression_weight"
].reindex(raw.index)
weight_sums = regression_weights.groupby(level="date").sum()
raw_dispersion = raw.groupby(level="date").std(ddof=0)
residual_dispersion = residual.groupby(level="date").std(ddof=0)
members = {
date: frozenset(group.index.get_level_values("ticker"))
for date, group in result.universe.loc[result.universe].groupby(level="date")
}
membership_dates = sorted(members)
turnover = pd.Series(
[
np.nan
if index == 0
else len(members[date] - members[membership_dates[index - 1]])
/ len(members[date])
for index, date in enumerate(membership_dates)
],
index=pd.DatetimeIndex(membership_dates),
name="turnover",
)
residual_wide = residual.unstack("ticker")
consecutive_pairs = residual_wide.notna() & residual_wide.shift(1).notna()
ticker_autocorrelation = residual_wide.corrwith(residual_wide.shift(1))
ticker_autocorrelation = ticker_autocorrelation.loc[
consecutive_pairs.sum() >= 60
].dropna()
r_squared = diagnostics["r_squared"]
factor_columns = [name for name in result.model_columns if name != "INTERCEPT"]
all_factor_returns = result.factor_returns.set_index("date")
display_names = {
"INTERCEPT": "Market level",
"MKT": "Beta premium",
"SIZE": "SIZE",
"HML": "HML",
"MOM": "MOM",
}
summary = pd.DataFrame(
{
"Result": [
f"{result.factor_returns['date'].min().date()} to "
f"{result.factor_returns['date'].max().date()}",
f"{len(result.factor_returns):,}",
f"{len(result.specific_returns):,}",
f"{result.specific_returns['ticker'].nunique():,}",
f"{diagnostics['n_observations'].mean():.1f}",
f"{diagnostics['n_observations'].min():.0f}",
f"{diagnostics['rank'].eq(len(result.model_columns)).mean():.1%}",
f"{r_squared.mean():.2%} / {r_squared.median():.2%}",
f"{r_squared.quantile(0.05):.2%} / {r_squared.quantile(0.95):.2%}",
f"{raw_dispersion.mean():.2%} / {residual_dispersion.mean():.2%}",
f"{residual_dispersion.mean() / raw_dispersion.mean():.2%}",
f"{ticker_autocorrelation.mean():.2%}",
f"{residual.abs().quantile(0.95):.2%} / "
f"{residual.abs().quantile(0.99):.2%}",
f"{residual.abs().max():.2%}",
f"{turnover.mean():.2%} / {turnover.quantile(0.95):.2%}",
f"{elapsed_seconds:.1f} seconds",
]
},
index=[
"Window",
"Regression dates",
"Residual-return rows",
"Unique tickers",
"Mean observations per regression",
"Minimum observations",
"Full-rank dates",
"Mean / median daily R²",
"5th / 95th percentile daily R²",
"Mean raw / residual dispersion",
"Residual / raw dispersion",
"Mean ticker lag-one residual autocorrelation",
"95th / 99th percentile absolute residual",
"Maximum absolute residual",
"Mean / 95th percentile universe turnover",
"Runtime on this machine",
],
)
summary
latest_beta_date = result.exposures["date"].max()
latest_betas = result.exposures.loc[
result.exposures["date"] == latest_beta_date, "MKT"
]
latest_beta_quantiles = latest_betas.quantile([0.01, 0.25, 0.50, 0.75, 0.99])
latest_beta_quantile_text = ", ".join(
f"{value:.2f}" for value in latest_beta_quantiles
)The lag-one statistic includes only tickers with at least 60 consecutive observation pairs. The observation count reflects the target universe after missing or non-finite returns, weights, and factor histories are removed.
For each date, R² compares weighted residual sum of squares with the weighted cross-sectional dispersion of raw returns around their weighted mean. The WLS run uses square-root-cap weights; the OLS comparison uses equal weights. It is not a time-series R² against the market proxy.
Universe choice
Show code
latest_caps = (
source.load(window.end, window.end)["market_cap"]
.dropna()
.loc[lambda values: values > 0]
.sort_values(ascending=False)
)
total_market_cap = latest_caps.sum()
coverage_rows = []
for size in (200, 500, 1000):
selected = latest_caps.iloc[:size]
coverage_rows.append(
{
"Names": size,
"Latest market-cap coverage": selected.sum() / total_market_cap,
"Smallest included market cap": selected.iloc[-1],
}
)
coverage = pd.DataFrame(coverage_rows).set_index("Names")
# Derive the figures quoted below the table so prose and table cannot diverge.
target_size = universe_builder.size
listed_count = f"{len(latest_caps):,}"
target_coverage = f"{coverage.loc[target_size, 'Latest market-cap coverage']:.2%}"
target_floor = (
f"₩{coverage.loc[target_size, 'Smallest included market cap'] / 1e9:,.0f} billion"
)
mid_coverage = f"{coverage.loc[500, 'Latest market-cap coverage']:.2%}"
excluded_count = f"{len(latest_caps) - target_size:,}"
excluded_share = (
f"{1 - coverage.loc[target_size, 'Latest market-cap coverage']:.2%}"
)
coverage.style.format(
{
"Latest market-cap coverage": "{:.2%}",
"Smallest included market cap": "₩{:,.0f}",
}
)| Latest market-cap coverage | Smallest included market cap | |
|---|---|---|
| Names | ||
| 200 | 90.38% | ₩2,243,502,841,470 |
| 500 | 95.60% | ₩473,039,949,450 |
| 1000 | 98.05% | ₩167,118,555,474 |
There are 2,872 active securities in the latest snapshot:
- Top 1,000 covers 98.05% of market cap and currently ends near ₩167 billion.
- Top 500 covers 95.60%.
- The excluded 1,872 names contribute only 1.95% in aggregate, and are the most exposed to sparse trading and stale closes.
Top 1,000 broadens the output while excluding the least liquid tail. The main run uses square-root-cap weights, giving larger names more influence. Equal weights make the fit describe the typical included stock and are the OLS alternative.
Ordered factor removal
For each date, the intercept first removes the mean return. Exposures are winsorized, centered, scaled where configured, and orthogonalized in the explicit factor order before the residual loop:
flowchart LR
raw["Raw return r"] --> center["Remove intercept"]
center --> market["1 · MKT"]
market --> size["2 · SIZE"]
size --> value["3 · HML"]
value --> momentum["4 · MOM"]
momentum --> epsilon["Residual return ε"]
MKT gets first claim on shared variation, followed by SIZE, HML, and MOM. Changing that order changes the basis and factor attribution, while the fitted span and final residual remain the same when the columns are full rank.
Show code
ordered_x = result.exposures.set_index(["date", "ticker"])
ordered_f = result.factor_returns.set_index("date")
aligned_ordered_f = ordered_f.reindex(
ordered_x.index.get_level_values("date")
).set_axis(ordered_x.index)
weights = regression_weights.reindex(ordered_x.index)
daily_weight_sums = weights.groupby(level="date").sum()
weighted_mean = (raw * weights).groupby(level="date").sum() / daily_weight_sums
aligned_weighted_mean = weighted_mean.reindex(
ordered_x.index.get_level_values("date")
).set_axis(ordered_x.index)
remaining = raw - aligned_ordered_f["INTERCEPT"]
total_ss = (
weights * (raw - aligned_weighted_mean).pow(2)
).groupby(level="date").sum()
stages = []
for name in factor_order:
exposure = ordered_x[name]
remaining = remaining - exposure * aligned_ordered_f[name]
stage_moment = (
(weights * remaining * exposure).groupby(level="date").sum()
/ daily_weight_sums
).abs()
residual_ss = (weights * remaining.pow(2)).groupby(level="date").sum()
stages.append(
{
"factor removed": name,
"mean cumulative R²": (1 - residual_ss / total_ss).mean(),
"maximum |exposure × stage residual|": stage_moment.max(),
}
)
stage_summary = pd.DataFrame(stages).set_index("factor removed")
stage_summary| mean cumulative R² | maximum |exposure × stage residual| | |
|---|---|---|
| factor removed | ||
| MKT | 0.033894 | 1.045680e-17 |
| SIZE | 0.046815 | 1.314098e-17 |
| HML | 0.055757 | 8.967887e-18 |
| MOM | 0.068537 | 1.045489e-17 |
Each stage moment is measured immediately after removal and should be near zero. Because the basis is orthogonalized, the same calculation against final ε should also be near zero for every factor:
Show code
final_moments = (
ordered_x.loc[:, list(factor_order)]
.mul(residual * weights, axis=0)
.groupby(level="date")
.sum()
.div(daily_weight_sums, axis=0)
.abs()
)
pd.DataFrame(
{
"mean |exposure × ε|": final_moments.mean(),
"maximum |exposure × ε|": final_moments.max(),
}
)| mean |exposure × ε| | maximum |exposure × ε| | |
|---|---|---|
| MKT | 6.040057e-19 | 1.223153e-17 |
| SIZE | 7.864892e-19 | 1.344659e-17 |
| HML | 6.112839e-19 | 8.757221e-18 |
| MOM | 6.961096e-19 | 9.755138e-18 |
Factor sensitivity
The table compares two universe sizes and two estimators. Returns are annualized mean factor returns; the R² columns are mean daily cross-sectional R².
| Universe and estimator | MKT, per unit beta | SIZE, per standardized unit | MKT stage R² | Final R² |
|---|---|---|---|---|
| Top 200, sequential OLS | 12.68% | 1.55% | 2.95% | 9.26% |
| Top 200, sequential WLS | 15.65% | 3.25% | 2.50% | 8.63% |
| Top 1,000, sequential OLS | 15.04% | 3.78% | 2.28% | 4.73% |
| Top 1,000, sequential WLS | 15.62% | 4.92% | 2.44% | 6.56% |
MKT stage R² is measured after the intercept and before SIZE, HML, and MOM. Final R² includes all four factors. The top-200 runs have higher R² because that group is more homogeneous. The top-1,000 runs cover more names but also contain more noisy, idiosyncratic returns. WLS gives larger companies more influence, which is a better fit for that broad universe than giving every name equal weight; it is the estimator used for the main run. The comparison rows are a snapshot; the main run and diagnostics below are recomputed when the page is rendered.
The chart compounds the market-level intercept and cumulatively sums the four factor slopes, which are contributions rather than portfolio returns.
Diagnostics
Show code
plot_columns = ["INTERCEPT", *factor_columns]
colors = ["#212529", "#2780e3", "#3fb618", "#9954bb", "#ff7518"]
cumulative_series = all_factor_returns[plot_columns].cumsum()
cumulative_series["INTERCEPT"] = (
1 + all_factor_returns["INTERCEPT"]
).cumprod() - 1
legend_names = {**display_names, "SIZE": "SIZE"}
def show_chart(figure: go.Figure) -> None:
figure.update_layout(
template="plotly_white",
hovermode="x",
autosize=True,
height=360,
margin={"l": 45, "r": 20, "t": 50, "b": 55},
legend={"orientation": "h", "y": -0.2},
)
figure.show(config={"displaylogo": False, "responsive": True})Show code
figure = go.Figure()
figure.add_trace(
go.Scatter(
x=r_squared.index,
y=r_squared,
mode="lines",
name="Daily R²",
line={"color": "#8ca8d8", "width": 1},
opacity=0.45,
)
)
figure.add_trace(
go.Scatter(
x=r_squared.index,
y=r_squared.rolling(63).mean(),
mode="lines",
name="63-day R²",
line={"color": "#2780e3", "width": 2},
)
)
figure.update_yaxes(tickformat=".0%")
show_chart(figure)
figure = go.Figure()
for name, color in zip(plot_columns, colors, strict=True):
figure.add_trace(
go.Scatter(
x=cumulative_series.index,
y=cumulative_series[name],
mode="lines",
name=legend_names[name],
line={"color": color, "width": 2 if name == "INTERCEPT" else 1.5},
)
)
figure.update_yaxes(tickformat=".0%")
show_chart(figure)
figure = go.Figure()
for name, values, color in (
("Raw dispersion", raw_dispersion, "#ff7518"),
("Residual dispersion", residual_dispersion, "#9954bb"),
):
figure.add_trace(
go.Scatter(
x=values.index,
y=values.rolling(21).mean(),
mode="lines",
name=name,
line={"color": color, "width": 1.5},
)
)
figure.update_yaxes(tickformat=".1%")
show_chart(figure)
figure = go.Figure(
go.Scatter(
x=turnover.index,
y=turnover.rolling(21).mean(),
mode="lines",
name="Universe turnover",
line={"color": "#613d7c", "width": 1.5},
)
)
figure.update_yaxes(tickformat=".1%")
show_chart(figure)Show code
factor_summary = pd.DataFrame(
{
"Annualized mean": all_factor_returns[plot_columns].mean() * 252,
"Annualized volatility": all_factor_returns[plot_columns].std(ddof=0)
* np.sqrt(252),
"Minimum daily return": all_factor_returns[plot_columns].min(),
"Maximum daily return": all_factor_returns[plot_columns].max(),
}
).rename(index=display_names)
factor_summary.style.format("{:.2%}")| Annualized mean | Annualized volatility | Minimum daily return | Maximum daily return | |
|---|---|---|---|---|
| Market level | 8.00% | 21.03% | -12.59% | 9.83% |
| Beta premium | 16.02% | 22.54% | -12.09% | 22.88% |
| SIZE | 5.11% | 5.60% | -1.84% | 2.37% |
| HML | 3.19% | 4.47% | -1.55% | 1.22% |
| MOM | 0.89% | 5.82% | -2.77% | 2.02% |
Show code
validation_rows = []
for validation in result.validation_results:
for metric, value in validation.metrics.items():
validation_rows.append(
{
"Validation": validation.name,
"Passed": validation.passed,
"Metric": metric,
"Value": value,
}
)
pd.DataFrame(validation_rows).style.format({"Value": "{:.6g}"})| Validation | Passed | Metric | Value | |
|---|---|---|---|---|
| 0 | regression_coverage | True | completed_dates | 2452 |
| 1 | regression_coverage | True | expected_dates | 2452 |
| 2 | regression_coverage | True | date_coverage | 1 |
| 3 | finite_outputs | True | nonfinite_outputs | 0 |
| 4 | return_reconstruction | True | maximum_reconstruction_error | 1.11022e-16 |
| 5 | residual_orthogonality | True | maximum_absolute_factor_moment | 5.24577e-16 |
| 6 | sequential_orthogonality | True | maximum_absolute_stage_moment | 2.38586e-17 |
Customize the analysis
Window and universe
window = resid.analysis_window(source, years=5, end="2025-12-31")
universe_builder = resid.DailyTopMarketCapUniverse(size=500)end=None selects the latest local observation. DailyTopMarketCapUniverse reselects from the previous close; FixedTopMarketCapUniverse selects from the close immediately before the first output date and holds that membership. Another point-in-time rule only needs to implement UniverseBuilder.build.
Factors and ordering
def volatility(returns: pd.DataFrame, _: pd.DataFrame) -> pd.DataFrame:
return returns.shift(1).rolling(63, min_periods=42).std()
styles = resid.CharacteristicFactorModel(
(
*base_factor_model.factors,
resid.Factor(
"VOLATILITY",
volatility,
history_business_days=64,
),
)
)
factor_model = resid.RecursiveMarketBetaModel(
base=styles,
lookback_days=252,
min_periods=126,
decay=0.97,
name="MKT",
)
residualizer = resid.SequentialOLSResidualizer(
factor_order=("MKT", "SIZE", "HML", "VOLATILITY", "MOM"),
winsor_quantile=0.01,
unscaled_factors=("MKT",),
)Add or remove characteristics in CharacteristicFactorModel, then make the same change to factor_order. Lookback, skip, minimum history, beta decay, and factor order are model choices. The recursive wrapper is optional when market beta is not part of the model.
Regression and weighting
regression_weight_model = resid.SquareRootMarketCapWeights()
residualizer = resid.SequentialWLSResidualizer(
factor_order=factor_order,
winsor_quantile=0.01,
unscaled_factors=("MKT",),
)The evaluated configuration pairs SquareRootMarketCapWeights with SequentialWLSResidualizer. EqualRegressionWeights with SequentialOLSResidualizer is the equal-influence alternative. For a standard multivariate fit, use OLSResidualizer with either equal or square-root-cap weights. Weight construction and residualization are separate replaceable boundaries.
Validation policy
validations = (
resid.RegressionCoverageValidation(minimum_date_coverage=0.995),
resid.FiniteOutputValidation(),
resid.ReturnReconstructionValidation(absolute_tolerance=1e-12),
resid.ResidualOrthogonalityValidation(absolute_tolerance=1e-12),
resid.SequentialOrthogonalityValidation(absolute_tolerance=1e-12),
)Coverage and numerical tolerances can vary by use case. Joint OLS uses ResidualOrthogonalityValidation; the sequential residualizers use the stage-wise check shown above. Additional policies implement RegressionValidation.
Output destination
result.specific_returns
resid.CsvArtifactWriter().write(result, Path("output/experiment-02"))The result can be consumed in memory or written to a different experiment directory. Another format can write directly from ResidualizationResult.
Interpretation
The mechanical checks pass:
- all expected dates have full-rank fits and finite outputs;
- stored returns reconstruct as \(r = Xf + \varepsilon\); and
- each factor is orthogonal to the residual immediately after its removal and in the final residual.
Those are regression identities, not proof of a complete risk model. Most single-name dispersion remains, and the largest residuals still reflect events and data tails. The small mean lag-one autocorrelation is more encouraging: there is no obvious broad one-day pattern left in ε.
The corrected recursive beta starts from a pre-window bootstrap or the 1.0 fallback, then uses the observed lagged-cap-weighted market return after each final fit. In the latest WLS cross-section (datetime.date(2026, 8, 7)), the beta quantiles from 1st through 99th percentile are -0.65, -0.39, -0.16, 0.14, 0.81. The MKT stage explains 3.39% of daily cross-sectional variance on average; all four factors explain 6.85% on average, with a 4.87% median.
Verdict
The output is a useful baseline for exploratory signal research with tail controls. It is not yet a production risk-model residual because the factor set is sparse and the recursive beta state is noisy.