Portfolio sorts and Fama-MacBeth regressions

Python
Backtesting
Tutorial
Code
When only the top fifth of simulated stocks earns extra return, five group averages show it, and a straight line predicts more return at every step.
Published

September 2, 2026

Part 1 of 3 in a series on stock characteristics and factors: 1. Portfolio sorts and Fama-MacBeth regressions, 2. Testing a stock characteristic as a factor, 3. How to do factor investing.

Does a stock characteristic, such as the ratio of book value to market value, predict next month’s return? I compare two methods on simulated data in which I set the answer: a portfolio sort, which compares the average returns of groups of stocks ranked by the characteristic, and a regression, which fits a straight line through the stocks’ returns.

In the simulated data, only the top fifth of stocks by the characteristic earns extra return. The table of group averages shows it in the top group alone, 1.47% a month against 0.57% to 0.66% for the other four groups. The fitted line has a positive slope, so it predicts an increase at every step up the ranking, which is false.

Simulating the data

The simulated data is a monthly panel of 500 stocks from January 2006 to December 2025, with one row per stock and month-end:

  • Characteristic. One number per stock and month, which changes slowly from month to month. I standardise it each month, so 0 is that month’s average and 1 is one standard deviation above it, where the standard deviation measures how far values typically are from their average.
  • True average return. The return the simulation sets for a stock: 0.6% a month, plus 0.8 percentage points, the extra return, when the stock was in the top fifth of the characteristic the month before.
  • Random part. Each month, a stock’s return is its true average return plus a random part, which is between −12% and 12% in about two months out of three. So averages measured in the data differ slightly from the true average returns.

A stock’s return in month \(t+1\) is its true average return plus its random part,

\[r_{i,t+1} = 0.006 + 0.008 \times \mathbf{1}\{x_{i,t} \ge q_{80,t}\} + 0.12\,\varepsilon_{i,t+1},\]

where \(r_{i,t+1}\) is stock \(i\)’s return in month \(t+1\), \(x_{i,t}\) is its standardised characteristic in month \(t\), \(q_{80,t}\) is month \(t\)’s 80th percentile of the characteristic, the value that 80% of that month’s stocks are below, and \(\mathbf{1}\{\cdot\}\) is 1 when the condition holds and 0 otherwise. \(\varepsilon_{i,t+1}\) is a draw from the standard normal distribution, which is bell-shaped with an average of 0 and a standard deviation of 1, independent across stocks and months and of the characteristic. So the random part, \(0.12\,\varepsilon_{i,t+1}\), has a standard deviation of 12%, 15 times the extra return. The 0.6%, the 0.8 percentage points and the 12% are assumptions, and returns are in excess of the return on cash.

Each stock’s characteristic keeps most of its value from one month to the next,

\[c_{i,t} = 0.95\,c_{i,t-1} + \sqrt{1 - 0.95^2}\;u_{i,t},\]

where \(c_{i,t}\) is stock \(i\)’s characteristic in month \(t\) and \(u_{i,t}\) a standard normal draw, independent across stocks and months. The first month’s values are standard normal draws too. The characteristic keeps an assumed 95% of last month’s value, and the factor \(\sqrt{1 - 0.95^2}\) sizes the fresh draw so that the variance, the square of the standard deviation, stays at 1: the two independent parts contribute \(0.95^2 \times 1 + (1 - 0.95^2) \times 1 = 1\). In one month’s 500 stocks the standard deviation is only close to 1, so I standardise the characteristic each month,

\[x_{i,t} = \frac{c_{i,t} - \bar c_t}{s_t},\]

where \(\bar c_t\) and \(s_t\) are the average and the sample standard deviation across that month’s 500 stocks. Market cap is $5 billion times \(e^{s_{i,t}}\), where the size score \(s_{i,t}\) follows the same two formulas with draws of its own.

The simulation returns the panel: one row per stock and month-end, with the date, the ticker, the month’s return, and the characteristic and market cap measured at the month-end, so no information from the next month enters the characteristic used to form the groups.

import numpy as np
import pandas as pd

STOCKS = 500                        # how many stocks the simulated data holds
MONTHS = 240                        # twenty years of monthly observations
CHARACTERISTIC_PERSISTENCE = 0.95   # how much of last month's value persists into this month
RANDOM_SD = 0.12                    # standard deviation of a stock's random part in a month
BASE_RETURN = 0.006                 # true average return in a month, before any extra return
TOP_FIFTH_RETURN = 0.008            # extra return of the top fifth
GROUPS = 5                          # how many groups the sort forms


def standardise(values):
    """Rescale each month to average 0 and standard deviation 1.

    One row is one month. Doing this month by month fixes what a value of 1 means
    without touching how the stocks are ranked against each other inside that month.
    """
    # keepdims=True keeps each month's average and standard deviation one column wide, so
    # every month is rescaled by its own two numbers
    return ((values - values.mean(axis=1, keepdims=True))
            / values.std(axis=1, ddof=1, keepdims=True))


def slow_moving_characteristic(generator, months, stocks):
    """One number per stock per month that mostly persists from the month before.

    A real characteristic changes slowly. A firm's ratio of book value to market value
    this month is close to last month's, so the loop keeps 95% of the previous value
    and adds a fresh random draw on top. The square root sizes that draw so the series
    keeps a standard deviation of one instead of growing over the sample.
    """
    fresh_share = np.sqrt(1 - CHARACTERISTIC_PERSISTENCE ** 2)
    values = np.empty((months, stocks))
    values[0] = generator.normal(0, 1, stocks)
    for month in range(1, months):
        carried_over = CHARACTERISTIC_PERSISTENCE * values[month - 1]
        values[month] = carried_over + fresh_share * generator.normal(0, 1, stocks)
    return standardise(values)

simulate_market builds the simulated data. generator.normal(0, RANDOM_SD, ...) draws the random part \(0.12\,\varepsilon\) directly, because its second argument is the standard deviation. month_cutoff is \(q_{80,t}\), the comparison characteristic >= month_cutoff is the indicator, and the extra return is added to the following month’s return. The size score is drawn after the returns and never enters them.

def simulate_market(seed):
    """Build the simulated data: one row per ticker and month-end, with the month's return and
    the characteristic and market cap measured at the month-end.

    A stock's true average return is BASE_RETURN, plus TOP_FIFTH_RETURN when the stock was in
    the top fifth of the characteristic in the month before.
    """
    generator = np.random.default_rng(seed)
    characteristic = slow_moving_characteristic(generator, MONTHS, STOCKS)
    # Each stock's random part, drawn separately for every stock and month
    stock_return = BASE_RETURN + generator.normal(0, RANDOM_SD, (MONTHS, STOCKS))
    # A size score, a second slow-moving number, drawn after the returns and with no effect on
    # them; market cap is $5 billion times e to the power of the score
    size_score = slow_moving_characteristic(generator, MONTHS, STOCKS)
    market_cap = 5e9 * np.exp(size_score)

    # keepdims=True keeps the cutoff one column wide, so each month is compared with its
    # own 80th percentile rather than with one number for the whole sample
    month_cutoff = np.percentile(characteristic, 80, axis=1, keepdims=True)
    extra_return = TOP_FIFTH_RETURN * (characteristic >= month_cutoff)

    # Row 1 onward gets the extra return from the row above, so this month's return uses
    # last month's characteristic and never its own.
    stock_return[1:] += extra_return[:-1]

    # One row per ticker and month-end: repeat gives each date once per ticker, tile gives the
    # tickers once per date, and ravel lays each months-by-stocks table out in that order
    dates = pd.date_range("2006-01-31", periods=MONTHS, freq="ME")    # the month-ends
    tickers = [f"S{number:03d}" for number in range(1, STOCKS + 1)]   # S001 to S500
    return pd.DataFrame({
        "date": dates.repeat(STOCKS),
        "ticker": np.tile(tickers, MONTHS),
        "return": stock_return.ravel(),
        "characteristic": characteristic.ravel(),
        "market_cap": market_cap.ravel(),
    })


# One fixed seed, so the page draws the same data every time it is rendered.
panel = simulate_market(7)

The block below shows the panel’s first three rows.

# The first three rows, with market cap in billions of dollars
first_rows = panel.head(3).copy()
first_rows["market_cap"] = first_rows["market_cap"] / 1e9
first_rows.round({"return": 4, "characteristic": 2, "market_cap": 1}).rename(
    columns={"market_cap": "market cap ($bn)"})
date ticker return characteristic market cap ($bn)
0 2006-01-31 S001 -0.2564 0.14 0.7
1 2006-01-31 S002 0.0126 0.46 20.1
2 2006-01-31 S003 -0.0217 -0.16 9.9

Each row’s next-month return is the same ticker’s return one calendar month later, matched by date, so a gap in the data never borrows a later month’s return. The final month, with no next month, is left out.

# Each ticker may appear once per month-end, and the month-ends, in date order, must follow
# each other
if panel.duplicated(["date", "ticker"]).any():
    raise ValueError("a ticker appears twice on one month-end")
month_ends = pd.DatetimeIndex(panel["date"].unique()).sort_values()
if not month_ends.equals(pd.date_range(month_ends.min(), month_ends.max(), freq="ME")):
    raise ValueError("the month-ends have a gap")
# Moving each row's date back one month-end and merging on date and ticker attaches every
# return to the row of the month before, as that row's next-month return
following = panel[["date", "ticker", "return"]].rename(columns={"return": "next_return"})
following["date"] = following["date"] - pd.offsets.MonthEnd(1)
panel = panel.merge(following, on=["date", "ticker"], how="left")
# The last month has no next month, so it is left out on purpose; any other missing
# next-month return stops the calculation
panel = panel.loc[panel["date"] < panel["date"].max()].reset_index(drop=True)
missing = panel["next_return"].isna()
if missing.any():
    first = panel.loc[missing].iloc[0]
    raise ValueError(f"no next-month return for {first['ticker']} after {first['date']:%B %Y}")
print(f"{len(panel):,} rows: {panel['ticker'].nunique()} tickers over "
      f"{panel['date'].nunique()} month-ends that have a next month")
119,500 rows: 500 tickers over 239 month-ends that have a next month

Sorting on one characteristic

To measure the average return at each level of the characteristic, I rank the stocks each month by their characteristic and cut the ranking into five groups of 100 at the 20th, 40th, 60th and 80th percentiles, the breakpoints. The 80th percentile is the value that 80% of the month’s stocks are below, so group 5 holds exactly the stocks with the extra return. Following Bali, Engle and Murray (2016), group \(g\)’s measured average return is

\[\bar R_g = \frac{1}{239}\sum_{t} \frac{1}{100}\sum_{i \in G_{g,t}} r_{i,t+1},\]

where \(r_{i,t+1}\) is stock \(i\)’s return in month \(t+1\) and \(G_{g,t}\) the set of stocks in group \(g\) in month \(t\): the inner sum averages next month’s return across the group’s 100 stocks, and the outer sum averages those monthly averages over the 239 months. The block below builds the groups and the table.

def five_groups(month_values):
    """Label one month's stocks 0 to 4, from the lowest fifth of the characteristic up."""
    # qcut cuts at the 20th, 40th, 60th and 80th percentiles of the values it is given
    return pd.qcut(month_values, GROUPS, labels=False)


# Each month's groups: transform runs five_groups within each month-end and returns the labels
# in the original row order; qcut labels from 0, so + 1 numbers the groups 1 to 5
panel["group"] = panel.groupby("date")["characteristic"].transform(five_groups) + 1
# The average next-month return of each group in each month, across the group's stocks
monthly_group_returns = panel.groupby(["date", "group"])["next_return"].mean()
# Each group's measured average return, across the months
measured_average_return = monthly_group_returns.groupby("group").mean()
average_characteristic = panel.groupby("group")["characteristic"].mean()
# The true average return set in the simulation: BASE_RETURN in every group, plus TOP_FIFTH_RETURN
# in group 5, which holds exactly the stocks in the top fifth
true_average_return = pd.Series(BASE_RETURN, index=measured_average_return.index)
true_average_return.loc[GROUPS] = BASE_RETURN + TOP_FIFTH_RETURN

# The textbook layout: .T turns the table so that each group is a column. In a notebook, a table
# on the last line of a block is shown as a formatted table.
GROUP_LABELS = ["1 (low)", "2", "3", "4", "5 (high)"]
table = pd.DataFrame({"true average return, % a month": true_average_return * 100,
                      "measured average return, % a month": measured_average_return * 100}).T
table.columns = pd.Index(GROUP_LABELS, name="characteristic group")
print("stocks in each group each month:", panel.groupby(["date", "group"]).size().unique())
table.round(2)
stocks in each group each month: [100]
characteristic group 1 (low) 2 3 4 5 (high)
true average return, % a month 0.60 0.60 0.60 0.60 1.40
measured average return, % a month 0.64 0.61 0.66 0.57 1.47

The table shows the extra return in group 5 only. Averaging over 100 stocks and 239 months cancels most of the random part, so every measured average return is within 0.07 percentage points of its true average return.

Sorting on two characteristics

A two-way sort compares the groups of one characteristic within the groups of a second, as in a table of size crossed with book-to-market. Here the second is market cap, which changes slowly and has no effect on returns. Each month, I sort all stocks into five market-cap groups, separately from the characteristic sort, and cross the two labels, without sorting the characteristic again inside each market-cap group. So the number of stocks in each of the 25 groups varies around 20. Group \((s, g)\)’s measured average return is

\[\bar R_{s,g} = \frac{1}{239}\sum_{t} \frac{1}{N_{s,g,t}}\sum_{i \in G_{s,g,t}} r_{i,t+1},\]

where \(G_{s,g,t}\) is the set of stocks in size group \(s\) and characteristic group \(g\) in month \(t\), and \(N_{s,g,t}\) their number. The table below gives each group’s measured average return in % a month. Each row is a size group and each column a characteristic group, so the top-right cell holds the small stocks with the highest characteristic.

# A second sort, on market cap, of all stocks each month, separate from the characteristic sort;
# grouping by both labels crosses the two sorts
panel["size_group"] = panel.groupby("date")["market_cap"].transform(five_groups) + 1
monthly_cell_returns = panel.groupby(["date", "size_group", "group"])["next_return"].mean()
# Each cell's measured average return across the months; unstack turns the characteristic
# groups into columns
two_way = monthly_cell_returns.groupby(["size_group", "group"]).mean().unstack() * 100
two_way.index = pd.Index(["1 (small)", "2", "3", "4", "5 (big)"], name="size group")
two_way.columns = pd.Index(GROUP_LABELS, name="characteristic group")
two_way.round(2)
characteristic group 1 (low) 2 3 4 5 (high)
size group
1 (small) 0.56 0.71 0.45 0.84 1.49
2 0.63 0.72 0.52 0.76 1.39
3 0.65 0.47 0.87 0.37 1.44
4 0.84 0.60 0.59 0.21 1.49
5 (big) 0.46 0.49 0.77 0.65 1.50

The top-right cell’s measured average return is 1.49% a month, and in every row, column 5 has the highest measured average return. So the extra return appears among small and big stocks alike. With fewer stocks, the cells differ more by chance than the five groups.

Estimating the relationship

The second method, a Fama-MacBeth regression (Fama and MacBeth, 1973), summarises the relation with a straight line. Each month, I fit a line that predicts next month’s return from this month’s characteristic, using that month’s 500 stocks,

\[\hat r_{i,t+1} = a_t + b_t\,x_{i,t},\]

where \(x_{i,t}\) is stock \(i\)’s standardised characteristic in month \(t\), \(\hat r_{i,t+1}\) its predicted return, the slope \(b_t\) the increase in predicted return for one more standard deviation of the characteristic, and the intercept \(a_t\) the predicted return at a characteristic of 0. Ordinary least squares (OLS) chooses \(a_t\) and \(b_t\) to make the sum of squared vertical distances between the stocks and the line as small as possible,

\[\sum_i \left(r_{i,t+1} - a_t - b_t\,x_{i,t}\right)^2.\]

The Fama-MacBeth estimates average the 239 monthly estimates,

\[a = \frac{1}{239}\sum_t a_t, \qquad b = \frac{1}{239}\sum_t b_t.\]

import statsmodels.api as sm


def fama_macbeth(frame):
    """Average of the monthly OLS coefficients of next month's return on the characteristic.

    Each month's regression fits a straight line through that month's stocks. Its slope is
    the increase in predicted return for one more standard deviation of the characteristic,
    and its intercept is the predicted return at a characteristic of 0. Fama and MacBeth (1973) average these
    monthly coefficients.
    """
    monthly_coefficients = []
    for date, month_rows in frame.groupby("date"):
        # add_constant adds a column of ones, whose coefficient is the intercept
        explanatory = sm.add_constant(month_rows[["characteristic"]])
        fitted = sm.OLS(month_rows["next_return"], explanatory).fit()
        monthly_coefficients.append(fitted.params)
    # One row per month and one column per coefficient; mean averages each column.
    return pd.DataFrame(monthly_coefficients).mean()


coefficients = fama_macbeth(panel)
slope = coefficients.loc["characteristic"]
print(f"average monthly slope: {slope * 100:.3f}% a month per standard deviation")
average monthly slope: 0.229% a month per standard deviation

The average slope is 0.229% a month per standard deviation, positive although only group 5 has extra return. At group \(g\)’s average characteristic \(\bar x_g\), the fitted line predicts \(\hat R_g = a + b\,\bar x_g\), so between two adjacent groups it predicts an increase of

\[\hat R_{g+1} - \hat R_g = b\,(\bar x_{g+1} - \bar x_g),\]

which is positive at every step, because \(b\) is positive and each group’s average characteristic is above the one before. The chart below draws each group’s measured average return as a bar at \(\bar x_g\), with the fitted line.

import matplotlib.pyplot as plt

RED, TEAL, GREY = "#C0392B", "#17868A", "#888888"
# The site's chart colours: a cream background, a light grid and a pale frame, so every chart
# on the site has the same look.
BG, INK, GRID, SPINE, AXIS_TEXT = "#FCEFE3", "#1f1f1f", "#EADCCC", "#D5C6B4", "#4a4a4a"


def style_chart(figure, axis):
    """Give a finished chart the site's cream background, light grid and pale frame."""
    figure.patch.set_facecolor(BG)
    axis.set_facecolor(BG)
    axis.grid(True, axis="y", color=GRID, lw=1.0)
    axis.set_axisbelow(True)
    for side in ("top", "right"):
        axis.spines[side].set_visible(False)
    for side in ("left", "bottom"):
        axis.spines[side].set_color(SPINE)
    axis.tick_params(axis="both", length=0, colors=INK, pad=6)
    axis.xaxis.label.set_color(AXIS_TEXT)
    axis.yaxis.label.set_color(AXIS_TEXT)
    axis.title.set_color(INK)
    # Each legend entry in the colour of what it labels: a line's colour or a bar's fill.
    # get_legend() returns None when the chart has no legend.
    legend = axis.get_legend()
    if legend is not None:
        for text, handle in zip(legend.get_texts(), legend.legend_handles):
            if hasattr(handle, "get_color"):
                text.set_color(handle.get_color())
            else:
                text.set_color(handle.get_facecolor())


# figsize is the figure's width and height in inches
fig, axis = plt.subplots(figsize=(7.6, 4.4))
positions = average_characteristic.to_numpy()
heights = measured_average_return.to_numpy() * 100
# width is each bar's width in units of the characteristic; label is its legend text
axis.bar(positions, heights, width=0.3, color=TEAL, label="measured average return")
for position, height in zip(positions, heights):
    # Each bar's value, in white just inside the bar's top: ha centres it, va puts its top there
    axis.text(position, height - 0.04, f"{height:.2f}%", ha="center", va="top", color="white")
# The fitted line at each group's average characteristic, a + b times that characteristic
fitted_line = (coefficients.loc["const"] + slope * average_characteristic) * 100
# lw is the line width, and marker="o" draws a circle at each group's average characteristic
axis.plot(positions, fitted_line, color=RED, lw=2.5, marker="o",
          label="fitted line from the regression")
# Group 3 averages a fraction below zero, which would print as -0.00
tick_labels = [f"group {group}\n{value:.2f}".replace("-0.00", "0.00")
               for group, value in average_characteristic.items()]
axis.set_xticks(positions, tick_labels)
axis.set_xlabel("group's average characteristic, standardised")
axis.set_ylabel("average next-month return, %")
axis.set_ylim(0, 1.7)
axis.set_title("All 239 months: group averages and the fitted straight line")
axis.legend(frameon=False, loc="upper left")   # frameon=False: no box around the legend
style_chart(fig, axis)

plt.tight_layout()   # fits the labels inside the figure
# dpi sets the saved image's pixels per inch, and bbox_inches="tight" trims the empty margin
plt.savefig("ps_nonlinear.png", dpi=140, bbox_inches="tight", facecolor=BG)
plt.show()

Five teal bars, the measured average returns of five groups of simulated stocks, at 0.64%, 0.61%, 0.66%, 0.57% and 1.47% a month, and a red fitted straight line that increases across all five groups, below the bars of groups 1 and 5 and above the bars of groups 2, 3 and 4.

The red line is below the bars of groups 1 and 5 and above the bars of groups 2, 3 and 4. Group 5’s extra return makes the slope positive, and because the line has one slope for all five groups, it predicts a higher return for group 4 than for group 1, while the table of group averages shows the extra return in group 5 alone.

Limitations

The straight line comes from the specification used here. A regression with a variable for each of groups 2 to 5, which is 1 for the group’s stocks and 0 otherwise, would show the step at group 5.

The sort has its own specification. Here the extra return starts exactly at a breakpoint, the 80th percentile. If only the top tenth of stocks had it, group 5 would mix 50 stocks with the extra return and 50 without, and the table would not show where inside group 5 the extra return is.

This post computes no standard errors. On real data, judging whether two measured average returns differ by more than chance needs the standard error of their difference.

Conclusion

The takeaway is that a positive slope does not show that the true average return increases at every step up the ranking. In the simulated data, all the extra return is in the top fifth, and the tables of group averages show where it is.

A two-way table is also a first look at combining characteristics, as a multifactor portfolio does. Here market cap has no effect on returns, so the table shows only that the characteristic’s extra return holds at every size.