Factor Returns Computation Validation

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:

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

  2. Manually compute factor returns with statsmodels.

  3. 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,
    CategoricalFilterSettings,
    ContinuousExposureGroupSettings,
    ExposureSettings,
    FactorRiskModelSettings,
    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-500-1y"),
    exposures=ExposureSettings(
        exposures=[
            ContinuousExposureGroupSettings(hierarchy="market"),
            CategoricalExposureGroupSettings(hierarchy="trbc"),
            ContinuousExposureGroupSettings(
                hierarchy="style", standardize_method="equal_weighted"
            ),
        ],
    ),
    modelconstruction=ModelConstructionSettings(
        weights="InvIdioVar",
        estimation_universe=UniverseSettings(
            dataset="Bayesline-US-500-1y", 
            categorical_filters=[
                CategoricalFilterSettings(hierarchy="estimation_universe")
            ],
        ),
        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-500-1y",
      "id_type": "bayesid",
      "calendar": {
        "dataset": "Bayesline-US-500-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",
      "estimation_universe": {
        "dataset": "Bayesline-US-500-1y",
        "id_type": "bayesid",
        "calendar": {
          "dataset": "Bayesline-US-500-1y",
          "filters": [
            [
              "XNYS"
            ]
          ]
        },
        "categorical_filters": [
          {
            "hierarchy": "estimation_universe",
            "include": "All",
            "exclude": []
          }
        ],
        "portfolio_filter": null,
        "mcap_filter": {
          "lower": 0.0,
          "upper": 1e20
        }
      },
      "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()

dataset = bln.equity.riskdatasets.load("Bayesline-US-500-1y")

dataset.describe().exposure_settings_menu.continuous_hierarchies

{'market': ['market'],
 'style': {'size': ['log_market_cap', 'log_total_assets'],
  'value': ['book_to_price'],
  'growth': ['price_to_earnings'],
  'volatility': ['sigma', 'sigma_eps', 'beta'],
  'momentum': ['mom6', 'mom12'],
  'dividend': ['dividend_yield'],
  'leverage': ['debt_to_assets', 'debt_to_equity']},
 'other': ['risk_free_rate', 'international_return']}

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.001238 -0.007247 -0.005044 -0.002368 -0.001060 -0.000460 0.000239 0.001188 0.001688 -0.009219 ... 0.0 0.0 0.0 0.000316 -0.000837 0.000967 -0.000333 0.003063 -0.000835 -0.000097
2024-08-28 -0.005025 -0.007031 -0.002325 0.003563 -0.000956 0.000907 0.003229 0.006363 -0.002291 0.000299 ... 0.0 0.0 0.0 -0.001170 0.001133 -0.001363 -0.005807 0.000949 0.000700 -0.000148
2024-08-29 0.004794 0.009394 0.001770 0.001462 -0.004767 -0.003243 0.001350 0.003728 -0.000522 0.001464 ... 0.0 0.0 0.0 -0.000339 0.000481 -0.000173 0.000332 -0.000209 -0.001748 -0.000502
2024-08-30 0.010128 -0.005914 0.002732 0.002031 -0.003292 0.001299 -0.001605 -0.002479 0.001657 -0.002147 ... 0.0 0.0 0.0 0.000358 0.000308 -0.001303 0.002947 0.000556 0.000532 0.000607
2024-09-03 -0.019321 -0.006828 -0.010553 -0.002509 -0.000859 0.008296 0.005331 0.007878 -0.003616 0.007325 ... 0.0 0.0 0.0 0.000555 -0.001163 -0.001286 -0.015321 -0.001910 0.001139 0.001130

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.Basic Materials trbc.Industrials trbc.Consumer Cyclicals trbc.Consumer Non-Cyclicals trbc.Financials trbc.Healthcare ... trbc.Institutions, Associations & Organizations trbc.Government Activity trbc.Academic & Educational Services style.Size style.Value style.Growth style.Volatility style.Momentum style.Dividend style.Leverage
0 2024-08-25 IC006CA2E0 1.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 ... 0.0 0.0 0.0 1.420898 0.277344 0.116577 -0.616699 0.073608 0.398193 -0.333008
1 2024-08-25 IC0121F541 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 -1.097656 -1.245117 0.749023 -0.600098 0.538574 -0.801270 1.960938
2 2024-08-25 IC012D373B 1.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 -0.143188 -0.257812 -0.162964 1.056641 0.922852 -0.485107 0.754395
3 2024-08-25 IC01430182 1.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 -0.520020 -0.968262 0.623535 0.841309 -2.435547 -1.234375 -1.688477
4 2024-08-25 IC015A481B 1.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 -0.522461 0.222046 0.645508 0.070251 1.075195 -1.151367 -0.738770

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 IC006CA2E0 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 1.420898 0.277344 0.116577 -0.616699 0.073608 0.398193 -0.333008
IC0121F541 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 ... 0.0 0.0 0.0 -1.097656 -1.245117 0.749023 -0.600098 0.538574 -0.801270 1.960938
IC012D373B 1.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 -0.143188 -0.257812 -0.162964 1.056641 0.922852 -0.485107 0.754395
IC01430182 1.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 -0.520020 -0.968262 0.623535 0.841309 -2.435547 -1.234375 -1.688477
IC015A481B 1.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 -0.522461 0.222046 0.645508 0.070251 1.075195 -1.151367 -0.738770

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 IC006CA2E0 0.003386
IC0121F541 0.003508
IC012D373B -0.009322
IC01430182 -0.001767
IC015A481B -0.001872

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 IC006CA2E0 0.009615
IC0121F541 0.012974
IC012D373B 0.017911
IC01430182 0.018917
IC015A481B 0.014382

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 IC006CA2E0 True
IC0121F541 True
IC012D373B True
IC01430182 True
IC015A481B 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 Basic Materials Industrials Consumer Cyclicals Consumer Non-Cyclicals Financials Healthcare Technology Utilities ... Size Value Growth Volatility Momentum Dividend Leverage return idio_vol estimation_universe
2024-08-25 IC006CA2E0 1.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 ... 1.420898 0.277344 0.116577 -0.616699 0.073608 0.398193 -0.333008 NaN NaN False
IC0121F541 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 ... -1.097656 -1.245117 0.749023 -0.600098 0.538574 -0.801270 1.960938 NaN NaN False
IC012D373B 1.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 ... -0.143188 -0.257812 -0.162964 1.056641 0.922852 -0.485107 0.754395 NaN NaN False
IC01430182 1.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 ... -0.520020 -0.968262 0.623535 0.841309 -2.435547 -1.234375 -1.688477 NaN NaN False
IC015A481B 1.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 ... -0.522461 0.222046 0.645508 0.070251 1.075195 -1.151367 -0.738770 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 Basic Materials Industrials Consumer Cyclicals Consumer Non-Cyclicals Financials Healthcare Technology Utilities ... Size Value Growth Volatility Momentum Dividend Leverage return idio_vol estimation_universe
2024-08-25 IC006CA2E0 1.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 ... 1.420898 0.277344 0.116577 -0.616699 0.073608 0.398193 -0.333008 NaN NaN False
IC0121F541 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 ... -1.097656 -1.245117 0.749023 -0.600098 0.538574 -0.801270 1.960938 NaN NaN False
IC012D373B 1.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 ... -0.143188 -0.257812 -0.162964 1.056641 0.922852 -0.485107 0.754395 NaN NaN False
IC01430182 1.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 ... -0.520020 -0.968262 0.623535 0.841309 -2.435547 -1.234375 -1.688477 NaN NaN False
IC015A481B 1.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 ... -0.522461 0.222046 0.645508 0.070251 1.075195 -1.151367 -0.738770 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 1.947083e+12 1.045178e+12 4.315087e+12 5.791688e+12 3.419057e+12 4.610744e+12 5.912425e+12 2.193800e+13 1.042832e+12 1.010499e+12 0.0 0.0 0.0
2024-08-27 1.929949e+12 1.043345e+12 4.320490e+12 5.770445e+12 3.433059e+12 4.620236e+12 5.920614e+12 2.202884e+13 1.035201e+12 1.013685e+12 0.0 0.0 0.0
2024-08-28 1.917057e+12 1.034918e+12 4.315773e+12 5.701664e+12 3.440224e+12 4.634233e+12 5.926554e+12 2.176459e+13 1.036469e+12 1.009837e+12 0.0 0.0 0.0
2024-08-29 1.938561e+12 1.040789e+12 4.352542e+12 5.700989e+12 3.449485e+12 4.662014e+12 5.940501e+12 2.167618e+13 1.041438e+12 1.005031e+12 0.0 0.0 0.0

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.001238 -0.005044 -0.002368 -0.001060 -0.000460 0.000239 0.001188 0.001688 -0.009219 0.001749 ... 1.0 0.460868 0.093702 0.080417 5.545968e-01 3.344099e-09 0.088992 0.828676 0.246861 0.000070
2024-08-28 -0.005025 -0.002325 0.003563 -0.000956 0.000907 0.003229 0.006363 -0.002291 0.000299 0.000556 ... 1.0 0.003226 0.012986 0.010407 1.307463e-24 3.088788e-02 0.120383 0.720238 0.390573 0.000060
2024-08-29 0.004794 0.001770 0.001462 -0.004767 -0.003243 0.001350 0.003728 -0.000522 0.001464 -0.006254 ... 1.0 0.632177 0.557501 0.852799 7.312268e-01 7.838325e-01 0.030701 0.502079 0.165913 0.000191
2024-08-30 0.010128 0.002732 0.002031 -0.003292 0.001299 -0.001605 -0.002479 0.001657 -0.002147 0.002580 ... 1.0 0.433436 0.561651 0.028482 1.084879e-06 2.769736e-01 0.308317 0.203011 0.516335 0.000079
2024-09-03 -0.019321 -0.010553 -0.002509 -0.000859 0.008296 0.005331 0.007878 -0.003616 0.007325 0.010604 ... 1.0 0.418282 0.145164 0.150278 2.347501e-51 1.674155e-02 0.146903 0.116333 0.668897 0.000180

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)