Beta Estimation

Note

You are reading Tidy Finance with Python. You can find the equivalent chapter for the sibling Tidy Finance with R here.

In this chapter, we introduce an important concept in financial economics: The exposure of an individual stock to changes in the market portfolio. According to the Capital Asset Pricing Model (CAPM) of Sharpe (1964), Lintner (1965), and Mossin (1966), cross-sectional variation in expected asset returns should be a function of the covariance between the excess return of the asset and the excess return on the market portfolio. The regression coefficient of excess market returns on excess stock returns is usually called the market beta. We show an estimation procedure for the market betas. We do not go into details about the foundations of market beta but simply refer to any treatment of the CAPM for further information. Instead, we provide details about all the functions that we use to compute the results. In particular, we leverage useful computational concepts: Rolling-window estimation and parallelization.

We use the following Python packages throughout this chapter:

import pandas as pd
import numpy as np
import statsmodels.formula.api as smf
import pyarrow.dataset as ds
import os

from plotnine import *
from mizani.formatters import percent_format
from joblib import Parallel, delayed, cpu_count
from itertools import product
from dateutil.relativedelta import relativedelta

Compared to previous chapters, we introduce statsmodels (Seabold and Perktold 2010) for regression analysis and for sliding-window regressions and joblib (Team 2023) for parallelization.

Estimating Beta Using Monthly Returns

The estimation procedure is based on a rolling-window estimation, where we may use either monthly or daily returns and different window lengths. First, let us start with loading the monthly CRSP data from our Parquet files introduced in Accessing and Managing Financial Data and WRDS, CRSP, and Compustat.

crsp_monthly = (pd.read_parquet("data-python/crsp_monthly.parquet")
    .get(["permno", "date", "industry", "ret_excess"])
    .dropna()
)

factors_ff3_monthly = (pd.read_parquet("data-python/factors_ff3_monthly.parquet")
    .get(["date", "mkt_excess"])
)

crsp_monthly = (crsp_monthly
    .merge(factors_ff3_monthly, how="left", on="date")
)

To estimate the CAPM regression coefficients

\[ r_{i, t} - r_{f, t} = \alpha_i + \beta_i(r_{m, t}-r_{f,t})+\varepsilon_{i, t}, \tag{1}\]

we regress stock excess returns ret_excess on excess returns of the market portfolio mkt_excess.

Python provides a simple solution to estimate (linear) models with the function smf.ols(). The function requires a formula as input that is specified in a compact symbolic form. An expression of the form y ~ model is interpreted as a specification that the response y is modeled by a linear predictor specified symbolically by model. Such a model consists of a series of terms separated by + operators. In addition to standard linear models, smf.ols() provides a lot of flexibility. You should check out the documentation for more information. To start, we restrict the data only to the time series of observations in CRSP that correspond to Apple’s stock (i.e., to permno 14593 for Apple) and compute \(\hat\alpha_i\) as well as \(\hat\beta_i\).

model_fit = smf.ols(
    formula="ret_excess ~ mkt_excess", data=crsp_monthly.query("permno == 14593")
).fit()
coefficients = model_fit.summary2().tables[1]
coefficients

	Coef.	Std.Err.	t	P>|t|	[0.025	0.975]
Intercept	0.009920	0.004884	2.031397	4.271697e-02	0.000327	0.019514
mkt_excess	1.375659	0.107611	12.783639	8.891466e-33	1.164259	1.587059

smf.ols() returns an object of class RegressionModel, which contains all the information we usually care about with linear models. summary2() returns information about the estimated parameters. The output above indicates that Apple moves excessively with the market as the estimated \(\hat\beta_i\) is above one (\(\hat\beta_i \approx 1.4\)).

Rolling-Window Estimation

After we estimated the regression coefficients on an example, we scale the estimation of \(\beta_i\) to a whole different level and perform rolling-window estimations for the entire CRSP sample. The following function implements the CAPM regression for a data frame (or a part thereof) containing at least min_obs observations to avoid huge fluctuations if the time series is too short. The function conveniently returns the regression results as a data frame, which ensures that our approach is scalable. If the min_obs-condition is violated, that is, the time series is too short, the function returns an empty data frame for consistency.

def estimate_capm(data, min_obs=1):
    if data.shape[0] < min_obs:
        capm = pd.DataFrame()
    else:
        fit = smf.ols(formula="ret_excess ~ mkt_excess", data=data).fit()
        coefficients = fit.summary2().tables[1]

        capm = pd.DataFrame(
            {
                "coefficient": coefficients.index,
                "estimate": coefficients["Coef."],
                "t_statistic": coefficients["t"],
            }
        ).assign(
            coefficient=lambda x: np.where(
                x["coefficient"] == "Intercept", "alpha", x["coefficient"]
            )
        )

    return capm

Next, we define a function that does the rolling estimation. We use a simple for-loop to implement the sliding window estimation in a straightforward manner. The following function takes input data and slides across the date vector, considering only a total of look_back months. The function essentially performs three steps: (i) arrange all rows, (ii) compute betas by sliding across months, and (iii) return a tibble with dates and corresponding parameter estimates. As we demonstrate further below, we can also apply the same function to daily return data.

def roll_capm_estimation(data, look_back=60, min_obs=48):
    results = []
    dates = data["date"].sort_values().drop_duplicates()

    for i in range(look_back - 1, len(dates)):
        end_date = dates.iloc[i]
        start_date = end_date - relativedelta(months=look_back - 1)

        window_data = data.query("date >= @start_date & date <= @end_date")

        result = estimate_capm(window_data, min_obs=min_obs)
        result["date"] = np.max(window_data["date"])
        results.append(result)

    if results:
        rolling_capm_estimation = pd.concat(results, ignore_index=True)
    else:
        rolling_capm_estimation = pd.DataFrame()

    return rolling_capm_estimation

Before we approach the whole CRSP sample, let us focus on a couple of examples for well-known firms.

examples = pd.DataFrame({
    "permno": [14593, 10107, 93436, 17778],
    "company": ["Apple", "Microsoft", "Tesla", "Berkshire Hathaway"]
})

The main idea is to apply the function to each stock individually and then combine the results into a single data frame. First, we nest the data by permno. Nested data means we now have a list of permno with corresponding grouped time series data. We get one row of output for each unique combination of non-nested variables which is only permno in this case.

capm_examples_nested = (crsp_monthly
    .query("permno in @examples['permno']")
    .groupby("permno", group_keys=True)
)
capm_examples_nested

<pandas.core.groupby.generic.DataFrameGroupBy object at 0x107562b00>

Next, we want to apply the roll_capm_estimation() function to each stock. This situation is an ideal use case for apply(), which takes a list or vector as input and returns an object of the same length as the input. In our case, apply() returns a single data frame with a time series of beta estimates for each stock. Therefore, we use reset_index() to transform the list of outputs to a tidy data frame.

capm_examples = (capm_examples_nested
    .apply(lambda x: roll_capm_estimation(x), include_groups=False)
    .reset_index()
    .get(["permno", "date", "coefficient", "estimate", "t_statistic"])
)
capm_examples

	permno	date	coefficient	estimate	t_statistic
0	10107	1991-03-01	alpha	0.042753	2.805863
1	10107	1991-03-01	mkt_excess	1.414250	5.018872
2	10107	1991-04-01	alpha	0.038412	2.534832
3	10107	1991-04-01	mkt_excess	1.434267	5.115538
4	10107	1991-05-01	alpha	0.039112	2.583179
...	...	...	...	...	...
3013	93436	2024-10-01	mkt_excess	2.380725	5.489195
3014	93436	2024-11-01	alpha	0.038044	1.593911
3015	93436	2024-11-01	mkt_excess	2.451193	5.645665
3016	93436	2024-12-01	alpha	0.039501	1.654410
3017	93436	2024-12-01	mkt_excess	2.385599	5.496380

3018 rows × 5 columns

Figure 1 displays the resulting beta estimates, focusing exclusively on the coefficient fo "mkt_excess".

beta_examples_sub = capm_examples.merge(examples, how="left", on="permno").query(
    "coefficient == 'mkt_excess'"
)

beta_figure = (
    ggplot(
        beta_examples_sub,
        aes(x="date", y="estimate", color="company", linetype="company"),
    )
    + geom_line()
    + labs(
        x="",
        y="",
        color="",
        linetype="",
        title="Monthly beta estimates for example stocks using 5 years of data",
    )
    + scale_x_datetime(date_breaks="5 year", date_labels="%Y")
)
beta_figure.show()

Figure 1: The figure shows monthly beta estimates for example stocks using five years of data. The CAPM betas are estimated with monthly data and a rolling window of length five years based on adjusted excess returns from CRSP. We use market excess returns from Kenneth French data library.

Parallelized Rolling-Window Estimation

Even though we could now just apply the function using .groubby() on the whole CRSP sample, we advise against doing it as it is computationally quite expensive. Remember that we have to perform rolling-window estimations across all stocks and time periods. However, this estimation problem is an ideal scenario to employ the power of parallelization. Parallelization means that we split the tasks which perform rolling-window estimations across different workers (or cores on your local machine).

If you have a Windows or Mac machine, it makes most sense use the default parallelization backend of joblib, which means that separate Python processes are running in the background on the same machine to perform the individual jobs. If you check out the documentation of joblib.parallel_config(), you can also see other ways to resolve the parallelization in different environments. Note that we use availableCores() to determine the number of cores available for parallelization, but keep one core free for other tasks. Some machines might freeze if all cores are busy with Python jobs.

n_cores = cpu_count() - 1

Using eight cores, the estimation for our sample of around 25k stocks takes around 20 minutes. Of course, you can speed up things considerably by having more cores available to share the workload or by having more powerful cores. Instead of using .apply() on groups, we use Parallel() to execute multiple tasks concurrently and delayed() to wrap each function call, allowing the calls to be queued and distributed to worker processes rather than executed immediately.

crsp_monthly_nested = (crsp_monthly
    .groupby("permno", group_keys=False)
)

capm_monthly = pd.concat(
    Parallel(n_jobs=n_cores)(
        delayed(
            lambda name, group: roll_capm_estimation(group).assign(permno=name)
        )(
            name, group
        )
        for name, group in crsp_monthly_nested
    )
).get(["permno", "date", "coefficient", "estimate", "t_statistic"])
capm_monthly

	permno	date	coefficient	estimate	t_statistic
0	10001	1991-01-01	alpha	0.008016	1.222418
1	10001	1991-01-01	mkt_excess	0.098174	0.813336
2	10001	1991-02-01	alpha	0.007914	1.206700
3	10001	1991-02-01	mkt_excess	0.095534	0.791486
4	10001	1991-03-01	alpha	0.007375	1.122804
...	...	...	...	...	...
225	93436	2024-10-01	mkt_excess	2.380725	5.489195
226	93436	2024-11-01	alpha	0.038044	1.593911
227	93436	2024-11-01	mkt_excess	2.451193	5.645665
228	93436	2024-12-01	alpha	0.039501	1.654410
229	93436	2024-12-01	mkt_excess	2.385599	5.496380

4268612 rows × 5 columns

Estimating Beta Using Daily Returns

Before we provide some descriptive statistics of our beta estimates, we implement the estimation for the daily CRSP sample as well. Depending on the application, you might either use longer horizon beta estimates based on monthly data or shorter horizon estimates based on daily returns. As loading the full daily CRSP data requires relatively large amounts of memory, we split the beta estimation into smaller chunks. The logic follows the approach that we use to download the daily CRSP data (see WRDS, CRSP, and Compustat).

First, we load the daily Fama-French market excess returns and extract the vector of dates.

factors_ff3_daily = (pd.read_parquet("data-python/factors_ff3_daily.parquet")
    .get(["date", "mkt_excess"])
)

We then create a connection to the daily CRSP data, but we don’t load the whole table into our memory. We only extract all distinct permno because we loop the beta estimation over batches of stocks with size 500. To estimate the CAPM over a consistent lookback window while accommodating different return frequencies, we adjust the minimum required number of observations accordingly. Specifically, we require at least 1,000 daily returns over a five‑year period for a valid estimation. This threshold is consistent with the monthly requirement of 48 observations out of 60 months, given that there are roughly 252 trading days in a year.

crsp_daily_db = ds.dataset(
    "data-python/crsp_daily", format="parquet", partitioning="hive"
)

permnos = [
    int(folder.split("=")[1])
    for folder in os.listdir("data-python/crsp_daily")
    if folder.startswith("permno=")
]

batch_size = 500
batches = np.ceil(len(permnos) / batch_size).astype(int)
min_obs = 1_000

We then proceed to perform the same steps as with the monthly CRSP data, just in batches: Load in daily returns, nest the data by stock, and parallelize the beta estimation across stocks. Note that we also convert the daily date to the beginning of the month so that we can still look back over 60 months and get one beta estimate per month, even though we are using daily data.

capm_daily = []

for j in range(1, batches + 1):
    permno_batch = permnos[((j - 1) * batch_size) : (min(j * batch_size, len(permnos)))]

    crsp_daily_sub = (crsp_daily_db
        .scanner(
            columns=["permno", "date", "ret_excess"],
            filter=ds.field("permno").isin(permno_batch),
        )
        .to_table()
        .to_pandas()
    )

    crsp_daily_sub_nested = (crsp_daily_sub
        .merge(factors_ff3_daily, how="inner", on="date")
        .assign(date=lambda x: x["date"].dt.to_period("M").dt.to_timestamp())
        .groupby("permno", group_keys=False)
    )

    results = Parallel(n_jobs=n_cores)(
        delayed(
            lambda name, group: roll_capm_estimation(group, min_obs=min_obs).assign(
                permno=name
            )
        )(name, group)
        for name, group in crsp_daily_sub_nested
    )

    if results:
        capm_daily_sub = pd.concat(results).get(
            ["permno", "date", "coefficient", "estimate", "t_statistic"]
        )
        capm_daily.append(capm_daily_sub)
    else:
        print(f"Warning: Batch {j} produced no results (insufficient data)")

    print(f"Batch {j} out of {batches} done ({(j / batches) * 100:.2f}%)\n")

capm_daily = pd.concat(capm_daily)

Comparing Beta Estimates

What is a typical value for stock betas? First, let us extract the relevant estimates from our CAPM results based on monthly returns.

beta_monthly = (capm_monthly
    .query("coefficient == 'mkt_excess'")
    .get(["permno", "date", "estimate"])
    .rename(columns={"estimate": "beta"})
    .assign(return_type="monthly")
)

To get some feeling, we illustrate the dispersion of the estimated \(\hat\beta_i\) across different industries and across time below. Figure 2 shows that typical business models across industries imply different exposure to the general market economy. However, there are barely any firms that exhibit a negative exposure to the market factor.

beta_industries = (beta_monthly
    .merge(crsp_monthly, how="inner", on=["permno", "date"])
    .dropna(subset="beta")
    .groupby(["industry","permno"])["beta"]
    .aggregate("mean")
    .reset_index()
)

industry_order = (beta_industries
    .groupby("industry")["beta"]
    .aggregate("median")
    .sort_values()
    .index.tolist()
)

beta_industries_figure = (
    ggplot(
        beta_industries, 
        aes(x="industry", y="beta")
    )
    + geom_boxplot()
    + coord_flip()
    + labs(
        x="",
        y="", 
        title="Firm-specific beta distributions by industry"
        )
    + scale_x_discrete(limits=industry_order)
)
beta_industries_figure.show()

Figure 2: The box plots show the average firm-specific beta estimates by industry.

Next, we illustrate the time-variation in the cross-section of estimated betas. Figure 3 shows the monthly deciles of estimated betas (based on monthly data) and indicates an interesting pattern: First, betas seem to vary over time in the sense that during some periods, there is a clear trend across all deciles. Second, the sample exhibits periods where the dispersion across stocks increases in the sense that the lower decile decreases and the upper decile increases, which indicates that for some stocks the correlation with the market increases while for others it decreases. Note also here: stocks with negative betas are a rare exception.

beta_quantiles = (
    beta_monthly.groupby("date")["beta"]
    .quantile(q=np.arange(0.1, 1.0, 0.1))
    .reset_index()
    .rename(columns={"level_1": "quantile"})
    .assign(quantile=lambda x: (x["quantile"] * 100).astype(int))
    .dropna()
)

linetypes = ["-", "--", "-.", ":"]
n_quantiles = beta_quantiles["quantile"].nunique()

beta_quantiles_figure = (
    ggplot(
        beta_quantiles,
        aes(x="date", y="beta", color="factor(quantile)", linetype="factor(quantile)"),
    )
    + geom_line()
    + labs(
        x="", y="", color="", linetype="", title="Monthly deciles of estimated betas"
    )
    + scale_x_datetime(date_breaks="5 year", date_labels="%Y")
    + scale_linetype_manual(
        values=[linetypes[l % len(linetypes)] for l in range(n_quantiles)]
    )
)
beta_quantiles_figure.show()

Figure 3: The figure shows monthly deciles of estimated betas. Each line corresponds to the monthly cross-sectional quantile of the estimated CAPM beta.

To compare the difference between daily and monthly data, we combine beta estimates to a single table.

beta_daily = (capm_daily
    .query("coefficient == 'mkt_excess'")
    .get(["permno", "date", "estimate"])
    .rename(columns={"estimate": "beta"})
    .assign(return_type="daily")
)

beta = pd.concat([beta_monthly, beta_daily], ignore_index=True)

To compare the difference between daily and monthly data, we combine beta estimates to a single table. Then, we use the table to plot a comparison of beta estimates for our example stocks in Figure 4.

beta_comparison = beta.merge(examples, how="inner", on="permno")

beta_comparison_figure = (
    ggplot(
        beta_comparison,
        aes(x="date", y="beta", color="return_type", linetype="return_type"),
    )
    + geom_line()
    + facet_wrap("~company", ncol=1)
    + labs(
        x="",
        y="",
        color="",
        linetype="",
        title="Comparison of beta estimates using monthly and daily data",
    )
    + scale_x_datetime(date_breaks="10 years", date_labels="%Y")
    + theme(figure_size=(6.4, 6.4))
)
beta_comparison_figure.show()

Figure 4: The figure shows the comparison of beta estimates using monthly and daily data. CAPM betas are computed using five years of monthly or three months of daily data. The two lines show the monthly estimates based on a rolling window for few exemplary stocks.

The estimates in Figure 4 look as expected. As you can see, it really depends on the data frequency how your beta estimates turn out because the estimates based on daily data are much smoother due to the higher number of observations in each regression.

Finally, we write the estimates to our local folder so that we can use them in later chapters.

beta.to_parquet("data-python/beta.parquet")

Whenever you perform some kind of estimation, it also makes sense to do rough plausibility tests. A possible check is to plot the share of stocks with beta estimates over time. This descriptive helps us discover potential errors in our data preparation or estimation procedure. For instance, suppose there was a gap in our output where we do not have any betas. In this case, we would have to go back and check all previous steps to find out what went wrong.

return_types = pd.DataFrame({"return_type": ["monthly", "daily"]})

beta_coverage = (crsp_monthly
    .merge(return_types, how="cross")
    .merge(beta, on=["permno", "date", "return_type"], how="left")
    .groupby(["date", "return_type"], as_index=False)
    .apply(lambda x: pd.Series({"share": x["beta"].notna().sum() / len(x)}))
)

beta_coverage_figure = (
    ggplot(
        beta_coverage,
        aes(x="date", y="share", color="return_type", linetype="return_type"),
    )
    + geom_line()
    + labs(
        x="",
        y="",
        color="",
        linetype="",
        title="End-of-month share of securities with beta estimates",
    )
    + scale_y_continuous(labels=percent_format())
    + scale_x_datetime(date_breaks="10 year", date_labels="%Y")
)
beta_coverage_figure.show()

Figure 5: The figure shows end-of-month share of securities with beta estimates. The two lines show the share of securities with beta estimates using five years of monthly or three months of daily data.

Figure 5 shows no issues, as the two coverage lines track each other closely, so we can proceed to the next check.

We also encourage everyone to always look at the distributional summary statistics of variables. You can easily spot outliers or weird distributions when looking at such tables.

(beta
    .groupby("return_type")["beta"]
    .describe()
    .round(2)
)

	count	mean	std	min	25%	50%	75%	max
return_type
daily	2170727.0	0.8	0.49	-3.67	0.42	0.76	1.12	4.74
monthly	2134306.0	1.1	0.70	-9.02	0.64	1.04	1.47	10.35

The summary statistics also look plausible for the two estimation procedures.

Finally, since we have two different estimators for the same theoretical object, we expect the estimators to be at least positively correlated (although not perfectly as the estimators are based on different sample periods and frequencies).

(beta
    .pivot_table(index=["permno", "date"], columns="return_type", values="beta")
    .reset_index()
    .get(["monthly", "daily"])
    .corr()
    .round(2)
)

return_type	monthly	daily
return_type
monthly	1.00	0.62
daily	0.62	1.00

Indeed, we find a positive correlation between our beta estimates. In the subsequent chapters, we mainly use the estimates based on monthly data, as most readers should be able to replicate them due to potential memory limitations that might arise with the daily data.

Key Takeaways

• CAPM betas can be estimated using rolling-window estimation and processed in parallel via joblib.
• Both monthly and daily return data can be used to estimate betas with different frequencies and window lengths, depending on the application.
• Summary statistics, visualization, and plausibility checks help to validate beta estimates across time and industries.

Exercises

1. Compute beta estimates based on monthly data using one, three, and five years of data and impose a minimum number of observations of 10, 28, and 48 months with return data, respectively. How strongly correlated are the estimated betas?
2. Compute beta estimates based on monthly data using five years of data and impose different numbers of minimum observations. How does the share of permno-date observations with successful beta estimates vary across the different requirements? Do you find a high correlation across the estimated betas?
3. Instead of using joblib, perform the beta estimation in a loop (using either monthly or daily data) for a subset of 100 permnos of your choice. Verify that you get the same results as with the parallelized code from above.
4. Filter out the stocks with negative betas. Do these stocks frequently exhibit negative betas, or do they resemble estimation errors?
5. Compute beta estimates for multi-factor models such as the Fama-French three-factor model by extending the estimate_capm() function with a model parameter. In particular, your regression should support the model \[ r_{i, t} - r_{f, t} = \alpha_i + \sum\limits_{j=1}^k\beta_{i,k}(r_{j, t}-r_{f,t})+\varepsilon_{i, t} \tag{2}\] where \(r_{i, t}\) are the \(k\) factor returns. Thus, you estimate four parameters (\(\alpha_i\) and the slope coefficients). Provide some summary statistics of the cross-section of firms and their exposure to the different factors.

References

Lintner, John. 1965. “Security prices, risk, and maximal gains from diversification.” The Journal of Finance 20 (4): 587–615. https://doi.org/10.1111/j.1540-6261.1965.tb02930.x.

Mossin, Jan. 1966. “Equilibrium in a capital asset market.” Econometrica 34 (4): 768–83. https://doi.org/10.2307/1910098.

Seabold, Skipper, and Josef Perktold. 2010. “Statsmodels: Econometric and Statistical Modeling with Python.” 9th Python in Science Conference.

Sharpe, William F. 1964. “Capital asset prices: A theory of market equilibrium under conditions of risk .” The Journal of Finance 19 (3): 425–42. https://doi.org/10.1111/j.1540-6261.1964.tb02865.x.

Team, Joblib Development. 2023. Joblib: Running Python Functions as Pipeline Jobs. Https://joblib.readthedocs.io/.