Beta Estimation

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 a useful computational concept: vectorized rolling-window estimation that computes the rolling betas of all stocks in closed form.

We use the following packages throughout this chapter:

library(tidyverse)
library(nanoparquet)
library(slider)

theme_set(theme_minimal())

Compared to previous chapters, we introduce slider (Vaughan 2021) for the rolling sums that we use in the closed-form beta estimation.

import polars as pl
import numpy as np
import pyfixest as pf

from plotnine import *
from mizani.formatters import percent_format

theme_set(theme_minimal())

Compared to previous chapters, we introduce pyfixest (Fischer 2024) for regression analysis.

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 <- read_parquet("data/crsp_monthly.parquet") |>
  select(permno, date, industry, ret_excess)

factors_ff3_monthly <- read_parquet("data/factors_ff3_monthly.parquet") |>
  select(date, mkt_excess)

crsp_monthly <- crsp_monthly |>
  left_join(factors_ff3_monthly, join_by(date))

crsp_monthly = (pl.read_parquet("data/crsp_monthly.parquet")
    .select("permno", "date", "industry", "ret_excess")
)

factors_ff3_monthly = (pl.read_parquet("data/factors_ff3_monthly.parquet")
    .select("date", "mkt_excess")
)

crsp_monthly = (crsp_monthly
    .join(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.

R provides a simple solution to estimate (linear) models with the function lm(). lm() 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, lm() 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 <- lm(
  "ret_excess ~ mkt_excess",
  data = crsp_monthly |>
    filter(permno == "14593")
)
coefficients <- summary(model_fit)$coefficients
coefficients

            Estimate Std. Error t value Pr(>|t|)
(Intercept)  0.00992    0.00488    2.03 4.27e-02
mkt_excess   1.37560    0.10761   12.78 8.89e-33

lm() returns an object of class lm which contains all information we usually care about with linear models. summary() 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\)).

Python provides a simple solution to estimate (linear) models with the function pf.feols(). 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, pf.feols() 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 = pf.feols(
    "ret_excess ~ mkt_excess",
    data=crsp_monthly.filter(pl.col("permno") == 14593)
)
coefficients = model_fit.tidy()
coefficients

             Estimate  Std. Error    t value  Pr(>|t|)      2.5%     97.5%
Coefficient                                                               
Intercept    0.009919    0.004884   2.031102  0.042747  0.000325  0.019513
mkt_excess   1.375598    0.107606  12.783617  0.000000  1.164207  1.586989

pf.feols() returns an object of class Feols, which contains all the information we usually care about with linear models. The tidy() method 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 the entire CRSP sample and perform rolling-window estimations for every stock. Fitting millions of individual regressions—one for each stock and month—would be prohibitively slow. Fortunately, for a single-regressor model such as the CAPM, the rolling OLS coefficients have a closed-form solution that we can compute from rolling sums. Writing \(x\) for the market excess return and \(y\) for the stock excess return, the slope and intercept over a window read

\[ \hat\beta_i = \frac{S_{xy}}{S_{xx}}, \qquad \hat\alpha_i = \bar y - \hat\beta_i \bar x, \tag{2}\]

where \(S_{xy} = \sum_t x_t y_t - \frac{1}{n}\sum_t x_t \sum_t y_t\) and \(S_{xx} = \sum_t x_t^2 - \frac{1}{n}\left(\sum_t x_t\right)^2\). Because these quantities only require the rolling sums \(\sum_t x_t\), \(\sum_t y_t\), \(\sum_t x_t^2\), \(\sum_t y_t^2\), and \(\sum_t x_t y_t\), we can compute the rolling betas—together with their \(t\)-statistics—for all stocks at once using vectorized rolling aggregations, without fitting a single regression. We refer to our blog post Fast, Vectorized Beta Estimation for a detailed derivation.

The following function implements this approach in three steps: (i) it computes the per-period cumulants for each stock, (ii) it rolls the sums over the requested look_back window, and (iii) it maps the aggregated sums into the closed-form estimates. We only keep windows with at least min_obs observations to avoid huge fluctuations if the time series is too short. As we demonstrate further below, we can apply the same function to daily returns data.

The slide_index_sum() function from the slider package computes the rolling sums over a monthly index in a straightforward manner.

roll_capm_estimation <- function(data, look_back = 60, min_obs = 48) {
  cumulants <- data |>
    summarize(
      n = n(),
      sum_x = sum(mkt_excess),
      sum_y = sum(ret_excess),
      sum_xx = sum(mkt_excess^2),
      sum_yy = sum(ret_excess^2),
      sum_xy = sum(mkt_excess * ret_excess),
      .by = c(permno, date)
    ) |>
    arrange(permno, date) |>
    mutate(
      period = year(date) * 12 + month(date),
      .by = permno
    )

  rolling_sums <- cumulants |>
    mutate(
      across(
        c(n, sum_x, sum_y, sum_xx, sum_yy, sum_xy),
        \(col) slide_index_sum(col, period, before = look_back - 1)
      ),
      .by = permno
    ) |>
    filter(n >= min_obs)

  estimates <- rolling_sums |>
    mutate(
      s_xx = sum_xx - sum_x^2 / n,
      s_xy = sum_xy - sum_x * sum_y / n,
      s_yy = sum_yy - sum_y^2 / n,
      beta = s_xy / s_xx,
      alpha = (sum_y - beta * sum_x) / n,
      sigma2 = (s_yy - beta * s_xy) / (n - 2),
      t_alpha = alpha / sqrt(sigma2 * (1 / n + (sum_x / n)^2 / s_xx)),
      t_beta = beta / sqrt(sigma2 / s_xx)
    )

  bind_rows(
    estimates |>
      transmute(
        permno,
        date,
        coefficient = "alpha",
        estimate = alpha,
        t_statistic = t_alpha
      ),
    estimates |>
      transmute(
        permno,
        date,
        coefficient = "mkt_excess",
        estimate = beta,
        t_statistic = t_beta
      )
  ) |>
    arrange(permno, date, coefficient) |>
    as_tibble()
}

The rolling() method from polars computes the rolling sums over the monthly date index in a straightforward manner.

def roll_capm_estimation(data, look_back=60, min_obs=48):
    cumulants = (data
        .group_by("permno", "date")
        .agg(
            n=pl.len().cast(pl.Float64),
            sum_x=pl.col("mkt_excess").sum(),
            sum_y=pl.col("ret_excess").sum(),
            sum_xx=(pl.col("mkt_excess") ** 2).sum(),
            sum_yy=(pl.col("ret_excess") ** 2).sum(),
            sum_xy=(pl.col("mkt_excess") * pl.col("ret_excess")).sum(),
        )
        .sort("permno", "date")
    )

    rolling_sums = (cumulants
        .rolling(
            index_column="date",
            period=f"{look_back}mo",
            group_by="permno",
            closed="right",
        )
        .agg(
            pl.col("n", "sum_x", "sum_y", "sum_xx", "sum_yy", "sum_xy").sum(),
        )
        .filter(pl.col("n") >= min_obs)
    )

    estimates = (rolling_sums
        .with_columns(
            s_xx=pl.col("sum_xx") - pl.col("sum_x") ** 2 / pl.col("n"),
            s_xy=pl.col("sum_xy") - pl.col("sum_x") * pl.col("sum_y") / pl.col("n"),
            s_yy=pl.col("sum_yy") - pl.col("sum_y") ** 2 / pl.col("n"),
        )
        .with_columns(beta=pl.col("s_xy") / pl.col("s_xx"))
        .with_columns(
            alpha=(pl.col("sum_y") - pl.col("beta") * pl.col("sum_x")) / pl.col("n"),
            sigma2=(pl.col("s_yy") - pl.col("beta") * pl.col("s_xy"))
                / (pl.col("n") - 2),
        )
        .with_columns(
            t_alpha=pl.col("alpha") / (
                pl.col("sigma2")
                * (1 / pl.col("n")
                    + (pl.col("sum_x") / pl.col("n")) ** 2 / pl.col("s_xx"))
            ).sqrt(),
            t_beta=pl.col("beta") / (pl.col("sigma2") / pl.col("s_xx")).sqrt(),
        )
    )

    return (pl.concat([
            estimates.select(
                "permno", "date",
                coefficient=pl.lit("alpha"),
                estimate=pl.col("alpha"),
                t_statistic=pl.col("t_alpha"),
            ),
            estimates.select(
                "permno", "date",
                coefficient=pl.lit("mkt_excess"),
                estimate=pl.col("beta"),
                t_statistic=pl.col("t_beta"),
            ),
        ])
        .sort("permno", "date", "coefficient")
    )

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

examples <- tibble(
  permno = c(14593, 10107, 93436, 17778),
  company = c("Apple", "Microsoft", "Tesla", "Berkshire Hathaway")
)

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

Because roll_capm_estimation() groups by permno internally, we do not have to loop over stocks or nest the data. We simply restrict the sample to the example stocks and hand the entire frame to the function, which returns a tidy data frame with a time series of beta estimates for each stock.

capm_examples <- crsp_monthly |>
  filter(permno %in% examples$permno) |>
  roll_capm_estimation()
capm_examples

# A tibble: 3,114 × 5
  permno date       coefficient estimate t_statistic
   <dbl> <date>     <chr>          <dbl>       <dbl>
1  10107 1990-03-01 alpha         0.0417        2.31
2  10107 1990-03-01 mkt_excess    1.40          4.17
3  10107 1990-04-01 alpha         0.0427        2.41
4  10107 1990-04-01 mkt_excess    1.39          4.20
5  10107 1990-05-01 alpha         0.0443        2.53
# ℹ 3,109 more rows

capm_examples = roll_capm_estimation(
    crsp_monthly.filter(pl.col("permno").is_in(examples["permno"]))
)
capm_examples

shape: (3_114, 5)

permno	date	coefficient	estimate	t_statistic
f64	date	str	f64	f64
10107.0	1990-03-01	"alpha"	0.041709	2.309559
10107.0	1990-03-01	"mkt_excess"	1.397341	4.173565
10107.0	1990-04-01	"alpha"	0.042693	2.413745
10107.0	1990-04-01	"mkt_excess"	1.385035	4.197367
10107.0	1990-05-01	"alpha"	0.044259	2.533768
…	…	…	…	…
93436.0	2024-10-01	"mkt_excess"	2.380626	5.488947
93436.0	2024-11-01	"alpha"	0.038041	1.593754
93436.0	2024-11-01	"mkt_excess"	2.451092	5.645404
93436.0	2024-12-01	"alpha"	0.039498	1.654254
93436.0	2024-12-01	"mkt_excess"	2.385498	5.496126

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

beta_examples <- capm_examples |>
  left_join(examples, join_by(permno)) |>
  filter(coefficient == "mkt_excess")

beta_examples |>
  ggplot(aes(x = date, y = estimate, color = company, linetype = company)) +
  geom_line() +
  labs(
    x = NULL,
    y = NULL,
    color = NULL,
    linetype = NULL,
    title = "Monthly beta estimates for example stocks using 5 years of data"
  )

Title: Monthly beta estimates for example stocks using five years of data. The figure shows a time series of beta estimates based on five years of monthly data for Apple, Berkshire Hathaway, Microsoft, and Tesla. The estimated betas vary over time and across stocks but are always positive for each stock. — 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.

beta_examples = (capm_examples
    .join(examples, how="left", on="permno")
    .filter(pl.col("coefficient") == "mkt_excess")
)

beta_figure = (
    ggplot(
        beta_examples,
        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_date(date_breaks="5 year", date_labels="%Y")
)
beta_figure.show()

Estimating Beta for the Whole Sample

Because roll_capm_estimation() estimates the betas in closed form and in a single vectorized pass, we can now apply it directly to the whole CRSP sample. The estimation for our sample of around 25k stocks finishes in a few seconds on a typical laptop—no explicit loops over stocks and no parallelization required.

Note

We could distribute the estimation across several cores, but for the closed-form CAPM estimation the parallelization overhead is not worth it. We revisit parallelization in Size Sorts and p-Hacking, where repeatedly sorting portfolios across many specifications makes it pay off.

capm_monthly <- roll_capm_estimation(crsp_monthly)
capm_monthly

# A tibble: 4,665,538 × 5
  permno date       coefficient estimate t_statistic
   <dbl> <date>     <chr>          <dbl>       <dbl>
1  10001 1990-01-01 alpha         0.0111       1.38 
2  10001 1990-01-01 mkt_excess    0.0983       0.677
3  10001 1990-02-01 alpha         0.0106       1.34 
4  10001 1990-02-01 mkt_excess    0.0976       0.678
5  10001 1990-03-01 alpha         0.0105       1.35 
# ℹ 4,665,533 more rows

Instead of implementing the rolling-window estimation by hand, you can also use the estimate_betas() function from the tidyfinance package, which is built precisely for this task:

library(tidyfinance)

estimate_betas(
  data = crsp_monthly,
  model = "ret_excess ~ mkt_excess",
  lookback = months(60),
  min_obs = 48
)

capm_monthly = roll_capm_estimation(crsp_monthly)
capm_monthly

shape: (4_665_538, 5)

permno	date	coefficient	estimate	t_statistic
f64	date	str	f64	f64
10001.0	1990-01-01	"alpha"	0.011065	1.378016
10001.0	1990-01-01	"mkt_excess"	0.098333	0.677007
10001.0	1990-02-01	"alpha"	0.010577	1.342069
10001.0	1990-02-01	"mkt_excess"	0.097602	0.677908
10001.0	1990-03-01	"alpha"	0.010462	1.353963
…	…	…	…	…
93436.0	2024-10-01	"mkt_excess"	2.380626	5.488947
93436.0	2024-11-01	"alpha"	0.038041	1.593754
93436.0	2024-11-01	"mkt_excess"	2.451092	5.645404
93436.0	2024-12-01	"alpha"	0.039498	1.654254
93436.0	2024-12-01	"mkt_excess"	2.385498	5.496126

Instead of implementing the rolling-window estimation by hand, you can also use the estimate_betas() function from the tidyfinance package, which is built precisely for this task:

import tidyfinance as tf

tf.estimate_betas(
    data=crsp_monthly,
    model="ret_excess ~ mkt_excess",
    lookback="60mo",
    min_obs=48
)

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. Because roll_capm_estimation() estimates all betas in a single vectorized pass, we can apply it directly to the full daily sample—just as we did for the monthly data, without splitting the estimation into batches.

First, we load the daily Fama-French market excess returns.

factors_ff3_daily <- read_parquet("data/factors_ff3_daily.parquet") |>
  select(date, mkt_excess)

factors_ff3_daily = (pl.read_parquet("data/factors_ff3_daily.parquet")
    .select("date", "mkt_excess")
)

We then load the daily CRSP returns. 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.

permnos <- list.dirs(
  "data/crsp_daily",
  full.names = FALSE,
  recursive = FALSE
) |>
  str_remove("permno=") |>
  as.integer()

crsp_daily <- permnos |>
  map(\(p) {
    paste0("data/crsp_daily/permno=", p) |>
      list.files(full.names = TRUE) |>
      map(read_parquet) |>
      list_rbind() |>
      mutate(permno = p)
  }) |>
  list_rbind() |>
  select(permno, date, ret_excess)

min_obs <- 1000

crsp_daily = (
    pl.scan_parquet("data/crsp_daily", hive_partitioning=True)
    .select("permno", "date", "ret_excess")
    .collect()
    # hive partition keys are read as integers; align with crsp_monthly
    .with_columns(permno=pl.col("permno").cast(pl.Float64))
)

min_obs = 1_000

We then proceed exactly as with the monthly CRSP data: we join the daily market excess returns, truncate the daily dates to the beginning of the month so that we can still look back over 60 months and get one beta estimate per month, and hand the data to roll_capm_estimation(), which estimates the betas for all stocks at once.

capm_daily <- crsp_daily |>
  inner_join(factors_ff3_daily, join_by(date)) |>
  mutate(date = floor_date(date, "month")) |>
  roll_capm_estimation(min_obs = min_obs)

capm_daily = (crsp_daily
    .join(factors_ff3_daily, how="inner", on="date")
    .with_columns(date=pl.col("date").dt.truncate("1mo"))
    .pipe(roll_capm_estimation, min_obs=min_obs)
)

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 |>
  filter(coefficient == "mkt_excess") |>
  select(permno, date, beta = estimate) |>
  mutate(return_type = "monthly")

beta_monthly = (capm_monthly
    .filter(pl.col("coefficient") == "mkt_excess")
    .select("permno", "date", "estimate")
    .rename({"estimate": "beta"})
    .with_columns(return_type=pl.lit("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.

crsp_monthly |>
  left_join(beta_monthly, join_by(permno, date)) |>
  drop_na(beta) |>
  group_by(industry, permno) |>
  summarize(beta = mean(beta), .groups = "drop") |>
  ggplot(aes(x = reorder(industry, beta, FUN = median), y = beta)) +
  geom_boxplot() +
  coord_flip() +
  labs(
    x = NULL,
    y = NULL,
    title = "Firm-specific beta distributions by industry"
  )

Title: Firm-specific beta distributions by industry. The figure shows box plots for each industry. Firms with the highest average CAPM beta belong to the public administration industry. Firms from the utility sector have the lowest average CAPM beta. The figure indicates very few outliers with negative CAPM betas. The large majority of all stocks has CAPM betas between 0.5 and 1.5. — Figure 2: The box plots show the average firm-specific beta estimates by industry.

beta_industries = (beta_monthly
    .join(crsp_monthly, how="inner", on=["permno", "date"])
    .drop_nulls("beta")
    .group_by("industry", "permno")
    .agg(beta=pl.col("beta").mean())
)

industry_order = (beta_industries
    .group_by("industry")
    .agg(beta=pl.col("beta").median())
    .sort("beta")
    ["industry"].to_list()
)

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()

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_monthly |>
  group_by(date) |>
  reframe(
    x = quantile(beta, seq(0.1, 0.9, 0.1)),
    quantile = 100 * seq(0.1, 0.9, 0.1)
  ) |>
  ggplot(aes(
    x = date,
    y = x,
    color = as_factor(quantile),
    linetype = as_factor(quantile)
  )) +
  geom_line() +
  labs(
    x = NULL,
    y = NULL,
    color = NULL,
    linetype = NULL,
    title = "Monthly deciles of estimated betas",
  )

Title: Monthly deciles of estimated betas. The figure shows time series of deciles of estimated betas to illustrate the distribution of betas over time. The top ten percent quantile on average is around two but varies substantially over time. The lowest ten percent quantile is around 0.4 on average but is highly correlated with the top quantile such that in general CAPM market betas seem to go up and down jointly. — Figure 3: The figure shows monthly deciles of estimated betas. Each line corresponds to the monthly cross-sectional quantile of the estimated CAPM beta.

quantiles = np.arange(0.1, 1.0, 0.1)

beta_quantiles = (
    beta_monthly
    .group_by("date")
    .agg(
        pl.col("beta").quantile(q, interpolation="linear").alias(str(int(round(q * 100))))
        for q in quantiles
    )
    .unpivot(index="date", variable_name="quantile", value_name="beta")
    .with_columns(pl.col("quantile").cast(pl.Int64))
    .drop_nulls()
)

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

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_date(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()

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

beta_daily <- capm_daily |>
  filter(coefficient == "mkt_excess") |>
  select(permno, date, beta = estimate) |>
  mutate(return_type = "daily")

beta <- bind_rows(beta_monthly, beta_daily)

beta_daily = (capm_daily
    .filter(pl.col("coefficient") == "mkt_excess")
    .select("permno", "date", "estimate")
    .rename({"estimate": "beta"})
    .with_columns(return_type=pl.lit("daily"))
)

beta = pl.concat([beta_monthly, beta_daily])

Then, we use the table to plot a comparison of beta estimates for our example stocks in Figure 4.

beta |>
  inner_join(examples, join_by(permno)) |>
  ggplot(aes(
    x = date,
    y = beta,
    color = return_type,
    linetype = return_type
  )) +
  geom_line() +
  facet_wrap(~company, ncol = 1) +
  labs(
    x = NULL,
    y = NULL,
    color = NULL,
    linetype = NULL,
    title = "Comparison of beta estimates using monthly and daily data"
  )

Title: Comparison of beta estimates using monthly and daily data. The figure shows a time series of beta estimates using five years of monthly versus daily data for Apple, Berkshire Hathaway, Microsoft, and Tesla. The estimates based on monthly data are smooth relative to the estimates based on daily data. However, the general trend and level is similar, irrespective of the choice of frequency. — 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 daily data. The two lines show the monthly estimates based on a rolling window for few exemplary stocks.

beta_comparison = beta.join(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_date(date_breaks="10 years", date_labels="%Y")
    + theme(figure_size=(6.4, 6.4))
)
beta_comparison_figure.show()

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 such that we can use them in later chapters.

write_parquet(beta, "data/beta.parquet")

beta.write_parquet("data/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.

beta_coverage <- crossing(
  crsp_monthly,
  tibble(return_type = c("monthly", "daily"))
) |>
  left_join(beta, join_by(permno, date, return_type)) |>
  group_by(date, return_type) |>
  summarize(share = sum(!is.na(beta)) / n(), .groups = "drop")

beta_coverage |>
  ggplot(
    aes(x = date, y = share, color = return_type, linetype = return_type)
  ) +
  geom_line() +
  scale_y_continuous(labels = scales::percent) +
  labs(
    x = NULL,
    y = NULL,
    color = NULL,
    linetype = NULL,
    title = "End-of-month share of securities with beta estimates"
  ) +
  coord_cartesian(ylim = c(0, 1))

Title: End-of-month share of securities with beta estimates. The figure shows two time series with end-of-year shares of securities with beta estimates using five years of monthly or daily data. There is almost no missing data for the estimates based on daily data. For the beta estimates based on monthly data, around 75 percent of all stock-month combinations provide sufficient long historical periods to estimate the beta. — 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 daily data.

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

beta_coverage = (
    crsp_monthly.join(return_types, how="cross")
    .join(beta, on=["permno", "date", "return_type"], how="left")
    .group_by("date", "return_type")
    .agg(share=pl.col("beta").is_not_null().mean())
)

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_date(date_breaks="10 year", date_labels="%Y")
)
beta_coverage_figure.show()

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 |>
  group_by(return_type) |>
  summarize(
    mean = mean(beta),
    sd = sd(beta),
    min = min(beta),
    q05 = quantile(beta, 0.05),
    q50 = quantile(beta, 0.50),
    q95 = quantile(beta, 0.95),
    max = max(beta),
    n = n()
  )

# A tibble: 2 × 9
  return_type  mean    sd    min    q05   q50   q95   max       n
  <chr>       <dbl> <dbl>  <dbl>  <dbl> <dbl> <dbl> <dbl>   <int>
1 daily       0.794 0.495  -3.67 0.0857 0.755  1.66  4.97 2354575
2 monthly     1.11  0.714 -13.1  0.132  1.04   2.33 11.7  2332769

(beta
    .group_by("return_type")
    .agg(
        count=pl.len(),
        mean=pl.col("beta").mean(),
        std=pl.col("beta").std(),
        min=pl.col("beta").min(),
        q05=pl.col("beta").quantile(0.05, interpolation="linear"),
        q50=pl.col("beta").quantile(0.50, interpolation="linear"),
        q95=pl.col("beta").quantile(0.95, interpolation="linear"),
        max=pl.col("beta").max(),
    )
    .sort("return_type")
    .with_columns(pl.col(pl.Float64).round(2))
)

shape: (2, 9)

return_type	count	mean	std	min	q05	q50	q95	max
str	u32	f64	f64	f64	f64	f64	f64	f64
"daily"	2354575	0.79	0.49	-3.67	0.09	0.75	1.66	4.97
"monthly"	2332769	1.11	0.71	-13.05	0.13	1.04	2.33	11.72

The summary statistics indicate that estimates based on daily returns are, on average, lower and less variable than those derived from monthly returns.

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

beta |>
  pivot_wider(names_from = return_type, values_from = beta) |>
  select(monthly, daily) |>
  cor(use = "complete.obs")

        monthly daily
monthly   1.000 0.618
daily     0.618 1.000

(beta
    .pivot(index=["permno", "date"], on="return_type", values="beta")
    .drop_nulls(["monthly", "daily"])
    .select(correlation=pl.corr("monthly", "daily").round(2))
)

shape: (1, 1)

correlation
f64
0.62

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 with a vectorized rolling-window approach that computes the closed-form OLS coefficients for all stocks at once, avoiding millions of individual regressions.
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

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?
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?
Instead of using the closed-form estimation, estimate the rolling betas for a subset of 100 permnos of your choice by fitting individual regressions with lm() (R) or pf.feols() (Python) in a rolling window. Verify that you get the same results as with the vectorized code from above.
Filter out the stocks with negative betas. Do these stocks frequently exhibit negative betas, or do they resemble estimation errors?
Compute beta estimates for multi-factor models such as the Fama-French three-factor model by writing a rolling estimation that fits the model per stock with lm() (R) or pf.feols() (Python). In particular, your regression should support the model \[ r_{i, t} - r_{f, t} = \alpha_i + \sum\limits_{j=1}^k\beta_{i,j}(r_{j, t}-r_{f,t})+\varepsilon_{i, t} \tag{3}\] where \(r_{j, t}\) are the \(k\) factor returns. Thus, for the three-factor model, 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

Fischer, Alexander. 2024. PyFixest: Fast High-Dimensional Fixed Effects Regression in Python. Https://pypi.org/project/pyfixest/.

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.

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.

Vaughan, Davis. 2021. slider: Sliding window functions. https://CRAN.R-project.org/package=slider.