Factor Returns Computation Validation#

In this tutorial we are going to compute factor returns from a standard Bayesline factor model and compare the results to those we manually compute with statsmodels. The specific steps we will follow are:

Create a basic risk model and extract factor returns.
Manually compute factor returns with statsmodels.
Compare results of extracted factor returns and statsmodels regression.

Throughout this notebook we work with a randomly generated dataset. The results should generalize to real data, but we do not show any real data on our public API. Bayesline clients can run this notebook on real data.

0. Imports & Setup#

For this tutorial notebook, you will need to import the following packages.

import numpy as np
import numpy.testing as tst
import pandas as pd
from statsmodels.api import WLS
from statsmodels.stats.weightstats import DescrStatsW

from bayesline.api.equity import (
    CategoricalExposureGroupSettings,
    ContinuousExposureGroupSettings,
    ExposureSettings,
    FactorRiskModelSettings,
    HierarchyLevel,
    ModelConstructionSettings,
    UniverseSettings,
)
from bayesline.apiclient import BayeslineApiClient

We will also need to have a Bayesline API client configured.

bln = BayeslineApiClient.new_client(
    endpoint="https://[ENDPOINT]",
    api_key="[API-KEY]",
)

1. Create a basic risk model and extract factor returns#

Set up the risk model settings#

We will set up a risk model to use InvIdioVar weights and otherwise default settings.

factorriskmodel_settings = FactorRiskModelSettings(
    universe=UniverseSettings(dataset="Bayesline-US-All-1y"),
    exposures=ExposureSettings(
        exposures=[
            ContinuousExposureGroupSettings(hierarchy="market"),
            CategoricalExposureGroupSettings(hierarchy="trbc"),
            ContinuousExposureGroupSettings(
                hierarchy="style", standardize_method="equal_weighted"
            ),
        ],
    ),
    modelconstruction=ModelConstructionSettings(
        weights="InvIdioVar",
        return_clip_bounds=(None, None),
        zero_sum_constraints={"trbc": "mcap_weighted"},
    ),
)

Let’s verify the risk model settings we configured above.

print(factorriskmodel_settings.model_dump_json(indent=2))

{
  "universe": [
    {
      "dataset": "Bayesline-US-All-1y",
      "id_type": "bayesid",
      "calendar": {
        "dataset": "Bayesline-US-All-1y",
        "filters": [
          [
            "XNYS"
          ]
        ]
      },
      "categorical_filters": [],
      "portfolio_filter": null,
      "mcap_filter": {
        "lower": 0.0,
        "upper": 1e20
      }
    }
  ],
  "exposures": [
    {
      "exposures": [
        {
          "exposure_type": "continuous",
          "hierarchy": {
            "hierarchy_type": "level",
            "name": "market",
            "level": 1
          },
          "factor_group": "market",
          "include": "All",
          "exclude": [],
          "standardize_method": "none"
        },
        {
          "exposure_type": "categorical",
          "hierarchy": {
            "hierarchy_type": "level",
            "name": "trbc",
            "level": 1
          },
          "factor_group": "trbc",
          "include": "All",
          "exclude": []
        },
        {
          "exposure_type": "continuous",
          "hierarchy": {
            "hierarchy_type": "level",
            "name": "style",
            "level": 1
          },
          "factor_group": "style",
          "include": "All",
          "exclude": [],
          "standardize_method": "equal_weighted"
        }
      ]
    }
  ],
  "modelconstruction": [
    {
      "currency": "USD",
      "weights": "InvIdioVar",
      "return_clip_bounds": [
        null,
        null
      ],
      "thin_category_shrinkage": {},
      "thin_category_shrinkage_overrides": {},
      "zero_sum_constraints": {
        "trbc": "mcap_weighted"
      },
      "known_factor_map": {},
      "fx_convert_returns": true
    }
  ],
  "halflife_idio_vra": null
}

Construct risk model with settings#

risk_model = bln.equity.riskmodels.load(factorriskmodel_settings).get_model()

Extract factor returns from risk model#

factor_returns = (
    (risk_model.fret())
    .to_pandas()
    .set_index("date")
)
factor_returns.columns = pd.MultiIndex.from_tuples(
    factor_returns.columns.str.split(".").to_list(), names=["factor_group", "factor"]
)

Let’s take a peek at the factor returns that our risk model computed. We have several groups of factors (market, industry, region, and style), and one or two factors in each group.

factor_returns.head()

factor_group	market	trbc													style
factor	Market	Energy	Basic Materials	Industrials	Consumer Cyclicals	Consumer Non-Cyclicals	Financials	Healthcare	Technology	Utilities	...	Institutions, Associations & Organizations	Government Activity	Academic & Educational Services	Size	Value	Growth	Volatility	Momentum	Dividend	Leverage
date
2024-08-27	-0.001146	-0.006938	0.000545	-0.001435	-0.001540	0.001061	0.001554	-0.000899	0.000959	-0.005658	...	0.0	-0.021948	-0.008925	-0.000081	7.483370e-05	-4.414787e-07	-0.000658	-9.355579e-07	-0.000432	0.000089
2024-08-28	-0.004645	-0.002795	0.001472	0.001878	-0.000945	-0.000169	0.002164	0.000441	-0.001346	0.000131	...	0.0	-0.051230	0.002912	0.000242	8.580666e-05	-7.595508e-08	-0.005433	1.385296e-05	-0.002449	0.000024
2024-08-29	0.004473	0.007536	-0.001839	0.002499	-0.004165	-0.003382	-0.003082	0.002344	0.001913	-0.001034	...	0.0	0.048979	-0.001993	-0.000221	9.156034e-07	8.478744e-07	0.000643	-6.729405e-06	0.001615	-0.000178
2024-08-30	0.006078	-0.003329	-0.001179	0.003683	-0.002518	0.004085	-0.002175	0.001650	-0.000175	0.004160	...	0.0	-0.072346	-0.003713	-0.000396	-9.593043e-05	3.463018e-07	0.004755	-1.841490e-05	0.001826	-0.000166
2024-09-03	-0.019542	-0.011158	0.002842	-0.010038	0.006229	0.002907	0.010235	-0.003587	-0.003931	0.000904	...	0.0	-0.119801	-0.006943	0.001270	2.319981e-05	-3.115371e-06	-0.019810	4.590908e-05	-0.005244	0.000612

5 rows × 21 columns

2. Manually compute factor returns with `statsmodels`#

Now, we will manually compute factor returns for numerical comparison to our results.

Write manual regression code#

This is our implementation of a basic factor returns regression with statsmodels which we will use for comparison. This implementation uses weighted constrained linear regression to model factor returns. Specifically, it solves the following optimization problem to compute factor returns \(f_t\)

(1)#\[\begin{align} \min_{f_t\in\mathbb{R}^k}{\sum_{i=1}^{n_t}{W_{i,t}\left(r_{i,t}-X_{i,t}^\top f_t\right)^2}}, \quad\mathrm{subject\ to} \quad\sum_{j\in\mathrm{ind}}{\left(\sum_{i=1}^{n_t}\mathrm{MCap}_{i,t}X_{i,j,t}\right)f_{j,t}}=0, \end{align}\]

where \(j\in\mathrm{ind}\) are all industry factors and \(W_{i,t}\) are the weights of each stock. The constraint effectively says that the (market-cap weighted) industry factor returns sum to zero.

We use weighted linear regression, and from an econometric perspective, the optimal regression weights are (proportional to) the inverses of the variance, \(W_{i,t}=1/\mathrm{Var}(\varepsilon_{i,t})\). Since this is not known, we use the idiosyncratic volatility, or the estimated error variance of a 100-day rolling time-series regression of the returns of each stock against its market factor, \(W_{i,t}=100/\sum_{\tau=1}^{100}{e_{i,t-\tau}^2}\), for the fitted regression,

(2)#\[\begin{align} r_{i,t-\tau}&=a_i + r_{t-\tau}^Mb_i+e_{i,t-\tau},&\tau&=1,\ldots,100, \end{align}\]

where the market return is simply the market-cap weighted return of all assets in the estimation universe, \(r_{t}^M=\sum_{i\in\mathcal{I}_t^E}{\mathrm{MCap}_{i,t}r_{i,t}}\).

def statsmodels_regression(
    df: pd.DataFrame,
    market_caps: pd.DataFrame,
    industry_names: list[str],
    substyle_names: list[str],
) -> pd.DataFrame:
    df = df.dropna()
    df = df[df["estimation_universe"]]
    factor_names = ["Market", *industry_names, *substyle_names]
    weights = (1 / df["idio_vol"] ** 2).fillna(0.0)

    X = df[factor_names].copy().astype(np.float32)
    y = df["return"]

    # compute the adjustment for industry exposure summing to zero
    adj = market_caps[1:] / market_caps[0]
    X.loc[:, industry_names[1:]] -= adj.values * X.loc[:, industry_names[:1]].values
    X = X.drop(columns=industry_names[0])

    wls = WLS(
        endog=y,
        exog=X,
        weights=weights,
        missing="drop",
        hasconst=False,  # for r-squared calculation
    ).fit()
    sigma2_eps = DescrStatsW(wls.resid, weights=weights, ddof=0).var

    # if the mcap is zero, then t_stats are 0.0 and p_values are 1.0
    zero_mcap = [i.replace("trbc.", "") for i in market_caps.loc[market_caps == 0.0].index]
    wls.tvalues[zero_mcap] = 0.0
    wls.pvalues[zero_mcap] = 1.0

    wls_results = pd.concat(
        {
            "factor_returns": wls.params,
            "t_stats": wls.tvalues,
            "p_values": wls.pvalues,
            "r2": pd.Series([wls.rsquared], [None]),
            "sigma2": pd.Series(sigma2_eps, [None]),
        },
    )
    return wls_results

Extract the required regression input data to compute factor returns#

A. Exposures#

First, let’s get the exposures for this universe as a pandas DataFrame.

X = risk_model.exposures().to_pandas()
X.head()

	date	bayesid	market.Market	trbc.Energy	trbc.Financials	trbc.Healthcare	...	style.Size	style.Value	style.Growth	style.Volatility	style.Momentum	style.Dividend	style.Leverage
0	2024-08-25	IC0007D96F	1.0	0.0	1.0	0.0	...	0.000000	0.000000	0.000000	-0.004295	4.277344	0.000000	0.000000
1	2024-08-25	IC000B1557	1.0	1.0	0.0	0.0	...	0.194702	0.594238	-0.021225	-0.271973	0.373047	-0.303711	-0.020889
2	2024-08-25	IC0021AFB7	1.0	0.0	1.0	0.0	...	0.000000	0.000000	0.000000	0.130005	11.117188	0.000000	0.000000
3	2024-08-25	IC002CE8B9	1.0	0.0	0.0	1.0	...	8.109375	4.558594	-0.317627	0.375000	1.728516	0.000000	0.345947
4	2024-08-25	IC002DC646	1.0	0.0	1.0	0.0	...	0.000000	0.000000	0.000000	0.388916	15.921875	0.000000	0.000000

5 rows × 23 columns

Next, we’ll clean up our exposures dataframe slightly and massage it into the desired shape and structure.

# rename to be cleaner names corresponding to df
X.columns = [col.split('.')[-1] if '.' in col else col for col in X.columns]

X['date'] = pd.to_datetime(X["date"])
exposures = X.set_index(["date", "bayesid"])
exposures.head()

		Market	Energy	Basic Materials	Industrials	Consumer Cyclicals	Consumer Non-Cyclicals	Financials	Healthcare	Technology	Utilities	...	Institutions, Associations & Organizations	Government Activity	Academic & Educational Services	Size	Value	Growth	Volatility	Momentum	Dividend	Leverage
date	bayesid
2024-08-25	IC0007D96F	1.0	0.0	0.0	0.0	0.0	0.0	1.0	0.0	0.0	0.0	...	0.0	0.0	0.0	0.000000	0.000000	0.000000	-0.004295	4.277344	0.000000	0.000000
	IC000B1557	1.0	1.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	...	0.0	0.0	0.0	0.194702	0.594238	-0.021225	-0.271973	0.373047	-0.303711	-0.020889
	IC0021AFB7	1.0	0.0	0.0	0.0	0.0	0.0	1.0	0.0	0.0	0.0	...	0.0	0.0	0.0	0.000000	0.000000	0.000000	0.130005	11.117188	0.000000	0.000000
	IC002CE8B9	1.0	0.0	0.0	0.0	0.0	0.0	0.0	1.0	0.0	0.0	...	0.0	0.0	0.0	8.109375	4.558594	-0.317627	0.375000	1.728516	0.000000	0.345947
	IC002DC646	1.0	0.0	0.0	0.0	0.0	0.0	1.0	0.0	0.0	0.0	...	0.0	0.0	0.0	0.000000	0.000000	0.000000	0.388916	15.921875	0.000000	0.000000

5 rows × 21 columns

B. Asset Returns#

Next, we get the returns of our asset universe across dates. Note that some values are NaN since our estimation universe does not include all the securities in the larger universe across all dates.

returns = (
    risk_model.future_asset_returns()
    .to_pandas()
    .set_index("date")
    .stack(dropna=False)
    .rename("return")
    .to_frame()
)
returns.head()

		return
date
2024-08-26	IC0007D96F	0.000000
	IC000B1557	-0.005362
	IC0021AFB7	-0.002999
	IC002CE8B9	-0.027211
	IC002DC646	-0.005133

C. Idiosyncractic Volatility#

idio_vol = (
    risk_model.weights()
    .to_pandas()
    .set_index("date")
    .stack()
    .rename("idio_vol")
    .to_frame()
)
idio_vol.head()

		idio_vol
date
2024-08-26	IC0007D96F	0.002448
	IC000B1557	0.026666
	IC0021AFB7	0.007886
	IC002CE8B9	0.041966
	IC002DC646	0.007236

D. Estimation Universe#

estimation_universe = (
    risk_model.estimation_universe()
    .to_pandas()
    .set_index("date")
    .stack()
    .rename("estimation_universe")
    .astype('bool')
    .to_frame()
)
estimation_universe.head()

		estimation_universe
date
2024-08-26	IC0007D96F	True
	IC000B1557	True
	IC0021AFB7	True
	IC002CE8B9	False
	IC002DC646	True

Join all the regression components#

df_all = pd.concat(
    [exposures, returns, idio_vol, estimation_universe],
    axis=1
)

key = pd.MultiIndex.from_product(
    [sorted(X['date'].unique()), sorted(X['bayesid'].unique())],
)
df_all = df_all.reindex(key)
df_all[["estimation_universe"]] = df_all[["estimation_universe"]].fillna(False)
df_all.head()

		Market	Energy	Financials	Healthcare	...	Size	Value	Growth	Volatility	Momentum	Dividend	Leverage	return	idio_vol	estimation_universe
2024-08-25	IC0007D96F	1.0	0.0	1.0	0.0	...	0.000000	0.000000	0.000000	-0.004295	4.277344	0.000000	0.000000	NaN	NaN	False
	IC000B1557	1.0	1.0	0.0	0.0	...	0.194702	0.594238	-0.021225	-0.271973	0.373047	-0.303711	-0.020889	NaN	NaN	False
	IC0021AFB7	1.0	0.0	1.0	0.0	...	0.000000	0.000000	0.000000	0.130005	11.117188	0.000000	0.000000	NaN	NaN	False
	IC002CE8B9	1.0	0.0	0.0	1.0	...	8.109375	4.558594	-0.317627	0.375000	1.728516	0.000000	0.345947	NaN	NaN	False
	IC002DC646	1.0	0.0	1.0	0.0	...	0.000000	0.000000	0.000000	0.388916	15.921875	0.000000	0.000000	NaN	NaN	False

5 rows × 24 columns

Extract Aggregated Market Cap#

We also need market cap data, but this comes in the form of already industry-aggregated market caps, so cannot be merged with the rest of the input data as it does not span across assets. We only need the market caps corresponding to our industry exposures, so we can filter out any other columns.

# get the market caps
market_caps = (
    risk_model.market_caps()
    .to_pandas()
    .set_index("date")
    .reindex(df_all.index.get_level_values(0).unique())
)

# get all industry names and reformat strings from factor returns
industries = factor_returns.loc[:,'trbc'].columns.to_list()
industries_long = [f"trbc.{i}" for i in industries]

# filter market caps to industries
market_caps = market_caps[industries_long]

Drop weekly and forward fill#

Since factor returns tend to be stable over short periods of time, we only compute the regression every week on Wednesdays. Thus, we will drop all data on non-Wednesdays and forward fill.

# slice on reset day (every wednesday) and fill forward
forward_fill_cols = [*exposures.columns, "estimation_universe", "idio_vol"]

# fill returns
df_all.loc[df_all["estimation_universe"], "return"] = df_all.loc[
    df_all["estimation_universe"],
    "return",
].fillna(0.0)

df_all.head()

		Market	Energy	Financials	Healthcare	...	Size	Value	Growth	Volatility	Momentum	Dividend	Leverage	return	idio_vol	estimation_universe
2024-08-25	IC0007D96F	1.0	0.0	1.0	0.0	...	0.000000	0.000000	0.000000	-0.004295	4.277344	0.000000	0.000000	NaN	NaN	False
	IC000B1557	1.0	1.0	0.0	0.0	...	0.194702	0.594238	-0.021225	-0.271973	0.373047	-0.303711	-0.020889	NaN	NaN	False
	IC0021AFB7	1.0	0.0	1.0	0.0	...	0.000000	0.000000	0.000000	0.130005	11.117188	0.000000	0.000000	NaN	NaN	False
	IC002CE8B9	1.0	0.0	0.0	1.0	...	8.109375	4.558594	-0.317627	0.375000	1.728516	0.000000	0.345947	NaN	NaN	False
	IC002DC646	1.0	0.0	1.0	0.0	...	0.000000	0.000000	0.000000	0.388916	15.921875	0.000000	0.000000	NaN	NaN	False

5 rows × 24 columns

market_caps = market_caps.ffill()
market_caps.index.name = "date"

market_caps.head()

	trbc.Energy	trbc.Basic Materials	trbc.Industrials	trbc.Consumer Cyclicals	trbc.Consumer Non-Cyclicals	trbc.Financials	trbc.Healthcare	trbc.Technology	trbc.Utilities	trbc.Real Estate	trbc.Institutions, Associations & Organizations	trbc.Government Activity	trbc.Academic & Educational Services
date
2024-08-25	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
2024-08-26	2.426835e+12	1.416096e+12	5.001533e+12	7.000492e+12	4.370353e+12	1.222101e+13	6.946814e+12	2.392271e+13	1.378550e+12	1.582750e+12	0.0	184134800.0	2.553390e+10
2024-08-27	2.404425e+12	1.412721e+12	4.995433e+12	6.972647e+12	4.384789e+12	1.223004e+13	6.951653e+12	2.401241e+13	1.370135e+12	1.585482e+12	0.0	179809504.0	2.555385e+10
2024-08-28	2.387219e+12	1.401278e+12	4.989150e+12	6.889843e+12	4.389141e+12	1.221587e+13	6.950168e+12	2.372629e+13	1.368753e+12	1.580623e+12	0.0	169923104.0	2.553539e+10
2024-08-29	2.415747e+12	1.410031e+12	5.024398e+12	6.886805e+12	4.400771e+12	1.225904e+13	6.968942e+12	2.366477e+13	1.377284e+12	1.575535e+12	0.0	178573792.0	2.552819e+10

Run `statsmodels` regression to compute factor returns#

Now, we will take all the components we created above and run our manual regression function to compute factor returns.

# get trade days from the universe
trade_days = risk_model.universe().to_pandas().set_index("date").index

# get all style names factor returns
styles = factor_returns.loc[:,'style'].columns.to_list()

# perform regression on each date
manual_computation = {
    g: statsmodels_regression(df, mcap, industries, styles)
    for (g, df), (_, mcap) in zip(df_all.groupby(level=0), market_caps.iterrows())
    if g in trade_days
}

manual_computation = (
    pd.concat(manual_computation, axis=1, names=["date"])
    .T.reindex(trade_days)
    .shift(1)
    .tail(-1)
)

manual_computation.head()

	factor_returns										...	p_values								r2	sigma2
	Market	Basic Materials	Industrials	Consumer Cyclicals	Consumer Non-Cyclicals	Financials	Healthcare	Technology	Utilities	Real Estate	...	Academic & Educational Services	Size	Value	Growth	Volatility	Momentum	Dividend	Leverage	NaN	NaN
date
2024-08-27	-0.001146	0.000545	-0.001435	-0.001540	0.001061	0.001554	-0.000899	0.000959	-0.005658	0.001091	...	0.030251	1.830720e-02	0.000836	0.420312	3.909110e-03	2.630122e-01	8.798910e-02	1.170795e-01	0.045071	0.000020
2024-08-28	-0.004645	0.001472	0.001878	-0.000946	-0.000169	0.002164	0.000441	-0.001346	0.000131	0.003239	...	0.525912	1.008470e-10	0.000581	0.899202	2.825100e-99	1.425678e-51	2.079678e-18	7.012554e-01	0.221222	0.000023
2024-08-29	0.004473	-0.001839	0.002499	-0.004165	-0.003382	-0.003082	0.002344	0.001913	-0.001034	-0.004374	...	0.685452	2.597278e-08	0.973088	0.185428	1.651366e-02	3.251447e-12	5.867095e-08	6.912765e-03	0.107825	0.000026
2024-08-30	0.006078	-0.001179	0.003683	-0.002518	0.004085	-0.002175	0.001650	-0.000175	0.004160	0.002690	...	0.427787	1.032722e-25	0.000207	0.569683	2.253188e-74	1.442490e-84	1.251940e-10	8.500013e-03	0.312595	0.000024
2024-09-03	-0.019542	0.002842	-0.010038	0.006229	0.002907	0.010235	-0.003587	-0.003931	0.000904	0.005696	...	0.401246	1.228712e-77	0.613939	0.003953	0.000000e+00	2.435693e-159	2.804079e-25	3.318839e-08	0.555630	0.000076

5 rows × 62 columns

3. Compare results of extracted factor returns and `statsmodels` regression.#

In order to perform a complete comparison, we take our factor returns and check that all of the following match the statsmodels regression computation.

factor returns
t-stats
p-values

Compare factor returns#

# drop the first industry due to the constraint handling
first_ind = "trbc." + industries[0]
drop_cols = ["date", first_ind]

factor_returns_array = risk_model.fret()
factor_returns_array = factor_returns_array.drop(drop_cols).to_numpy()

tst.assert_array_almost_equal(factor_returns_array, manual_computation["factor_returns"].values)

Compare t-stats#

t_stats = risk_model.t_stats()
t_stats = t_stats.drop(drop_cols).to_numpy()

tst.assert_array_almost_equal(t_stats, manual_computation["t_stats"].values, decimal=3)

Compare p-values#

p_values = risk_model.p_values()
p_values = p_values.drop(drop_cols).to_numpy()

# compare the p-values in different ranges with different precisions (high for low p-values)
# this is due to numerical precision of the t-distribution CDF approximation
tst.assert_array_almost_equal(p_values.clip(None, 0.05), manual_computation["p_values"].clip(None, 0.05).values, decimal=4)
tst.assert_array_almost_equal(p_values.clip(0.05, 0.9), manual_computation["p_values"].clip(0.05, 0.9).values, decimal=3)
tst.assert_array_almost_equal(p_values.clip(0.9, None), manual_computation["p_values"].clip(0.9, None).values, decimal=2)

Factor Returns Computation Validation

Contents