Fama–French 3-Factor Model

Fama and French (1993)

Paper

Fama, E. and K. French (1993), Common Risk Factors in the Returns on Stocks and Bonds, Journal of Financial Economics

Data

CRSP monthly stock file, msf_v2 CIZ format; Compustat annual fundamentals; CRSP/Compustat link table

Sample

NYSE, AMEX and NASDAQ common stocks, factors through June 2026

Author

Qingyi (Freda) Song Drechsler

Download notebook (.ipynb)

This set of Python code replicates the Fama French risk factors SMB and HML, in addition to the excess market risk factor. It utilizes CRSP data for pricing related items and Compustat data for fundamental data.

##########################################
# Fama French 3 Factors-CRSP New Format  #
# Qingyi (Freda) Song Drechsler          #
# Date: August 2022                      #
# Updated: October 2026                  #
##########################################

import pandas as pd
import numpy as np
import datetime as dt
import wrds
#import psycopg2 
import matplotlib.pyplot as plt
import plotly.graph_objects as go
import plotly.io as pio
from plotly.subplots import make_subplots
from dateutil.relativedelta import *
from pandas.tseries.offsets import *
from scipy import stats

pio.renderers.default = 'notebook_connected'

conn = wrds.Connection()

Loading library list...
Done

###################
# Compustat Block #
###################
comp = conn.raw_sql("""
                    select gvkey, datadate, at, pstkl, txditc,
                    pstkrv, seq, pstk
                    from comp.funda
                    where indfmt='INDL' 
                    and datafmt='STD'
                    and popsrc='D'
                    and consol='C'
                    and datadate >= '01/01/1959'
                    """, date_cols=['datadate'])

comp['year']=comp['datadate'].dt.year

# create preferrerd stock
comp['ps']=np.where(comp['pstkrv'].isnull(), comp['pstkl'], comp['pstkrv'])
comp['ps']=np.where(comp['ps'].isnull(),comp['pstk'], comp['ps'])
comp['ps']=np.where(comp['ps'].isnull(),0,comp['ps'])
comp['txditc']=comp['txditc'].fillna(0)

# create book equity
comp['be']=comp['seq']+comp['txditc']-comp['ps']
comp['be']=np.where(comp['be']>0, comp['be'], np.nan)

# number of years in Compustat
comp=comp.sort_values(by=['gvkey','datadate'])
comp['count']=comp.groupby(['gvkey']).cumcount()

comp=comp[['gvkey','datadate','year','be','count']]

###################
# CRSP Block      #
###################
# sql similar to crspmerge macro
# CRSP monthly update (crsp_m_stock) covers the sample through June 2026; the annual crsp library ends December 2025
crsp_m = conn.raw_sql("""
                      select a.permno, a.permco, a.mthcaldt, 
                      a.issuertype, a.securitytype, a.securitysubtype, a.sharetype, a.usincflg, 
                      a.primaryexch, a.conditionaltype, a.tradingstatusflg,
                      a.mthret, a.mthretx, a.shrout, a.mthprc
                      from crsp_m_stock.msf_v2 as a
                      where a.mthcaldt between '01/01/1959' and '06/30/2026'
                      """, date_cols=['mthcaldt']) 

crsp_m.shape

(4970181, 15)

Select common stock universe

• equivalent to legacy code shrcd = 10 or 11

crsp_m = crsp_m.loc[(crsp_m.sharetype=='NS') & \
                    (crsp_m.securitytype=='EQTY') & \
                    (crsp_m.securitysubtype=='COM') & \
                    (crsp_m.usincflg=='Y') & \
                    (crsp_m.issuertype.isin(['ACOR', 'CORP']))]
crsp_m.shape

(3627503, 15)

Select stocks traded on NYSE, AMEX and NASDAQ

• equivalent to legacy code exchcd = 1, 2 or 3

crsp_m = crsp_m.loc[(crsp_m.primaryexch.isin(['N', 'A', 'Q'])) & \
                   (crsp_m.conditionaltype =='RW') & \
                   (crsp_m.tradingstatusflg =='A')]
crsp_m.shape

(3565359, 15)

# change variable format to int
crsp_m[['permco','permno']]=crsp_m[['permco','permno']].astype(int)

# Line up date to be end of month
crsp_m['jdate']=crsp_m['mthcaldt']+MonthEnd(0)

No need to add delisting return in the new CIZ CRSP format

• last trading return and price prior to delisting are already part of the time series output

# calculate market equity
crsp = crsp_m.copy()

crsp['mthret']=crsp['mthret'].fillna(0)
crsp['mthretx']=crsp['mthretx'].fillna(0)
crsp['me']=crsp['mthprc']*crsp['shrout'] 
crsp=crsp.drop(['mthprc','shrout'], axis=1)
crsp=crsp.sort_values(by=['jdate','permco','me'])

### Aggregate Market Cap ###
# sum of me across different permno belonging to same permco a given date
crsp_summe = crsp.groupby(['jdate','permco'])['me'].sum().reset_index()

# largest mktcap within a permco/date
crsp_maxme = crsp.groupby(['jdate','permco'])['me'].max().reset_index()

# join by jdate/maxme to find the permno
crsp1=pd.merge(crsp, crsp_maxme, how='inner', on=['jdate','permco','me'])

# drop me column and replace with the sum me
crsp1=crsp1.drop(['me'], axis=1)

# join with sum of me to get the correct market cap info
crsp2=pd.merge(crsp1, crsp_summe, how='inner', on=['jdate','permco'])

# sort by permno and date and also drop duplicates
crsp2=crsp2.sort_values(by=['permno','jdate']).drop_duplicates()

# keep December market cap
crsp2['year']=crsp2['jdate'].dt.year
crsp2['month']=crsp2['jdate'].dt.month
decme=crsp2[crsp2['month']==12]
decme=decme[['permno','mthcaldt','jdate','me','year']].rename(columns={'me':'dec_me'})

### July to June dates
crsp2['ffdate']=crsp2['jdate']+MonthEnd(-6)
crsp2['ffyear']=crsp2['ffdate'].dt.year
crsp2['ffmonth']=crsp2['ffdate'].dt.month
crsp2['1+retx']=1+crsp2['mthretx']
crsp2=crsp2.sort_values(by=['permno','mthcaldt'])

# cumret by stock
crsp2['cumretx']=crsp2.groupby(['permno','ffyear'])['1+retx'].cumprod()

# lag cumret
crsp2['lcumretx']=crsp2.groupby(['permno'])['cumretx'].shift(1)

# lag market cap
crsp2['lme']=crsp2.groupby(['permno'])['me'].shift(1)

# if first permno then use me/(1+retx) to replace the missing value
crsp2['count']=crsp2.groupby(['permno']).cumcount()
crsp2['lme']=np.where(crsp2['count']==0, crsp2['me']/crsp2['1+retx'], crsp2['lme'])

# baseline me
mebase=crsp2[crsp2['ffmonth']==1][['permno','ffyear', 'lme']].rename(columns={'lme':'mebase'})

# merge result back together
crsp3=pd.merge(crsp2, mebase, how='left', on=['permno','ffyear'])
crsp3['wt']=np.where(crsp3['ffmonth']==1, crsp3['lme'], crsp3['mebase']*crsp3['lcumretx'])

decme['year']=decme['year']+1
decme=decme[['permno','year','dec_me']]

# Info as of June
crsp3_jun = crsp3[crsp3['month']==6]

crsp_jun = pd.merge(crsp3_jun, decme, how='inner', on=['permno','year'])
crsp_jun=crsp_jun[['permno','mthcaldt', 'jdate', 'sharetype', 'securitytype', 'securitysubtype', 'usincflg', 'issuertype', \
                   'primaryexch', 'conditionaltype', 'tradingstatusflg', \
                   'mthret','me','wt','cumretx','mebase','lme','dec_me']]
crsp_jun=crsp_jun.sort_values(by=['permno','jdate']).drop_duplicates()

#######################
# CCM Block           #
#######################
ccm=conn.raw_sql("""
                  select gvkey, lpermno as permno, linktype, linkprim, 
                  linkdt, linkenddt
                  from crsp_m_ccm.ccmxpf_linktable
                  where substr(linktype,1,1)='L'
                  and (linkprim ='C' or linkprim='P')
                  """, date_cols=['linkdt', 'linkenddt'])

# if linkenddt is missing then set to today date
ccm['linkenddt']=ccm['linkenddt'].fillna(pd.to_datetime('today'))

ccm1=pd.merge(comp[['gvkey','datadate','be', 'count']],ccm,how='left',on=['gvkey'])
ccm1['yearend']=ccm1['datadate']+YearEnd(0)
ccm1['jdate']=ccm1['yearend']+MonthEnd(6)

# set link date bounds
ccm2=ccm1[(ccm1['jdate']>=ccm1['linkdt'])&(ccm1['jdate']<=ccm1['linkenddt'])]
ccm2=ccm2[['gvkey','permno','datadate','yearend', 'jdate','be', 'count']]

# link comp and crsp
ccm_jun=pd.merge(crsp_jun, ccm2, how='inner', on=['permno', 'jdate'])
ccm_jun['beme']=ccm_jun['be']*1000/ccm_jun['dec_me']

# select NYSE stocks for bucket breakdown
# legacy data format: exchcd = 1 and positive beme and positive me and shrcd in (10,11) and at least 2 years in comp
# new CIZ format: primaryexch == 'N', positive beme, positive me, at least 2 years in comp
# shrcd in 10 and 11 is already handled in the code earlier

nyse=ccm_jun[(ccm_jun['primaryexch']=='N') & (ccm_jun['beme']>0) & (ccm_jun['me']>0) & \
             (ccm_jun['count']>=1)]

# size breakdown
nyse_sz=nyse.groupby(['jdate'])['me'].median().to_frame().reset_index().rename(columns={'me':'sizemedn'})

# beme breakdown
nyse_bm=nyse.groupby(['jdate'])['beme'].describe(percentiles=[0.3, 0.7]).reset_index()
nyse_bm=nyse_bm[['jdate','30%','70%']].rename(columns={'30%':'bm30', '70%':'bm70'})

nyse_breaks = pd.merge(nyse_sz, nyse_bm, how='inner', on=['jdate'])

# join back size and beme breakdown
ccm1_jun = pd.merge(ccm_jun, nyse_breaks, how='left', on=['jdate'])

# function to assign sz and bm bucket
def sz_bucket(row):
    if row['me']==np.nan:
        value=''
    elif row['me']<=row['sizemedn']:
        value='S'
    else:
        value='B'
    return value

def bm_bucket(row):
    if 0<=row['beme']<=row['bm30']:
        value = 'L'
    elif row['beme']<=row['bm70']:
        value='M'
    elif row['beme']>row['bm70']:
        value='H'
    else:
        value=''
    return value

# assign size portfolio
ccm1_jun['szport']=np.where((ccm1_jun['beme']>0)&(ccm1_jun['me']>0)&(ccm1_jun['count']>=1), ccm1_jun.apply(sz_bucket, axis=1), '')

# assign book-to-market portfolio
ccm1_jun['bmport']=np.where((ccm1_jun['beme']>0)&(ccm1_jun['me']>0)&(ccm1_jun['count']>=1), ccm1_jun.apply(bm_bucket, axis=1), '')

# create positivebmeme and nonmissport variable
ccm1_jun['posbm']=np.where((ccm1_jun['beme']>0)&(ccm1_jun['me']>0)&(ccm1_jun['count']>=1), 1, 0)
ccm1_jun['nonmissport']=np.where((ccm1_jun['bmport']!=''), 1, 0)

# store portfolio assignment as of June

june=ccm1_jun[['permno','mthcaldt', 'jdate', 'bmport','szport','posbm','nonmissport']].copy()
june['ffyear']=june['jdate'].dt.year

# merge back with monthly records
crsp3 = crsp3[['mthcaldt','permno', 'sharetype', 'securitytype', 'securitysubtype', 'usincflg', 'issuertype', \
               'primaryexch', 'conditionaltype', 'tradingstatusflg', \
               'mthret', 'me','wt','cumretx','ffyear','jdate']]
ccm3=pd.merge(crsp3, 
        june[['permno','ffyear','szport','bmport','posbm','nonmissport']], how='left', on=['permno','ffyear'])

# keeping only records that meet the criteria
ccm4=ccm3[(ccm3['wt']>0)& (ccm3['posbm']==1) & (ccm3['nonmissport']==1)]

############################
# Form Fama French Factors #
############################

# function to calculate value weighted return
def wavg(group, avg_name, weight_name):
    d = group[avg_name]
    w = group[weight_name]
    try:
        return (d * w).sum() / w.sum()
    except ZeroDivisionError:
        return np.nan

# value-weigthed return
vwret=ccm4.groupby(['jdate','szport','bmport']).apply(wavg, 'mthret','wt').to_frame().reset_index().rename(columns={0: 'vwret'})
vwret['sbport']=vwret['szport']+vwret['bmport']

# firm count
vwret_n=ccm4.groupby(['jdate','szport','bmport'])['mthret'].count().reset_index().rename(columns={'mthret':'n_firms'})
vwret_n['sbport']=vwret_n['szport']+vwret_n['bmport']

# tranpose
ff_factors=vwret.pivot(index='jdate', columns='sbport', values='vwret').reset_index()
ff_nfirms=vwret_n.pivot(index='jdate', columns='sbport', values='n_firms').reset_index()

# create SMB and HML factors
ff_factors['WH']=(ff_factors['BH']+ff_factors['SH'])/2
ff_factors['WL']=(ff_factors['BL']+ff_factors['SL'])/2
ff_factors['WHML'] = ff_factors['WH']-ff_factors['WL']

ff_factors['WB']=(ff_factors['BL']+ff_factors['BM']+ff_factors['BH'])/3
ff_factors['WS']=(ff_factors['SL']+ff_factors['SM']+ff_factors['SH'])/3
ff_factors['WSMB'] = ff_factors['WS']-ff_factors['WB']
ff_factors=ff_factors.rename(columns={'jdate':'date'})

# n firm count
ff_nfirms['H']=ff_nfirms['SH']+ff_nfirms['BH']
ff_nfirms['L']=ff_nfirms['SL']+ff_nfirms['BL']
ff_nfirms['HML']=ff_nfirms['H']+ff_nfirms['L']

ff_nfirms['B']=ff_nfirms['BL']+ff_nfirms['BM']+ff_nfirms['BH']
ff_nfirms['S']=ff_nfirms['SL']+ff_nfirms['SM']+ff_nfirms['SH']
ff_nfirms['SMB']=ff_nfirms['B']+ff_nfirms['S']
ff_nfirms['TOTAL']=ff_nfirms['SMB']
ff_nfirms=ff_nfirms.rename(columns={'jdate':'date'})

###################
# Compare With FF #
###################
_ff = conn.get_table(library='ff', table='factors_monthly')
_ff=_ff[['date','smb','hml']]
_ff[['smb','hml']]=_ff[['smb','hml']].astype(float)
_ff['date']=_ff['date']+MonthEnd(0)

_ffcomp = pd.merge(_ff, ff_factors[['date','WSMB','WHML']], how='inner', on='date')
_ffcomp70=_ffcomp[_ffcomp['date']>='01/01/1970']


print(stats.pearsonr(_ffcomp70['smb'], _ffcomp70['WSMB']))
print(stats.pearsonr(_ffcomp70['hml'], _ffcomp70['WHML']))

PearsonRResult(statistic=0.9955581642809235, pvalue=0.0)
PearsonRResult(statistic=0.9835819838496571, pvalue=0.0)

_ffcomp.set_index('date', inplace=True)
_ffcomp.tail(2)

	smb	hml	WSMB	WHML
date
2026-05-31	-0.0159	-0.0231	-0.007666	-0.026781
2026-06-30	0.0339	0.0358	0.041610	0.024341

Comparison charts

Monthly factor returns, replicated from CRSP CIZ (blue) against the Ken French data library factors on WRDS (orange). Hover over the chart for both returns in a given month; drag across a period to zoom in.

#################################
# Monthly Factor Returns        #
#################################
fig = make_subplots(rows=2, cols=1, shared_xaxes=True, vertical_spacing=0.14,
                    subplot_titles=('SMB', 'HML'))
for row, (ff_col, w_col, name) in enumerate([('smb','WSMB','SMB'), ('hml','WHML','HML')], start=1):
    # Ken French drawn first, replicated on top, so differences show in both colors
    for col, label, color, width, opacity in [(ff_col, 'Ken French', '#eb6834', 1.6, 1.0),
                                              (w_col, 'Replicated (CRSP CIZ)', '#2a78d6', 1.2, 0.85)]:
        fig.add_trace(go.Scatter(x=_ffcomp.index, y=_ffcomp[col], name=label, mode='lines', opacity=opacity,
                                 line=dict(color=color, width=width), legendgroup=label, showlegend=(row == 1),
                                 hovertemplate='%{y:+.2%}<extra>' + f'{label}: {col}' + '</extra>'), row=row, col=1)
    fig.add_hline(y=0, line=dict(color='#8a8984', width=1), row=row, col=1)
    fig.update_yaxes(title_text=f'{name} monthly return', tickformat='.0%', gridcolor='#ecebe7',
                     zeroline=False, row=row, col=1)
fig.update_xaxes(showgrid=False, ticks='outside', tickcolor='#c3c2b7', linecolor='#c3c2b7',
                 range=[_ffcomp.index[0], _ffcomp.index[-1]])
fig.update_annotations(font=dict(size=14, color='#0b0b0b'), x=0, xanchor='left', yshift=6)
fig.update_layout(
    title=dict(text=f'Comparison of Results<br><sup style="color:#52514e">Monthly SMB and HML, replicated vs. Ken French, '
                    f'{_ffcomp.index[0]:%b %Y}–{_ffcomp.index[-1]:%b %Y}</sup>', x=0, xanchor='left'),
    template='plotly_white', height=760, margin=dict(l=80, r=40, t=100, b=90),
    font=dict(family='Inter, system-ui, sans-serif', size=13, color='#0b0b0b'),
    paper_bgcolor='#fcfcfb', plot_bgcolor='#fcfcfb', hovermode='x unified',
    legend=dict(orientation='h', x=0, y=-0.07, yanchor='top'),
)
fig.show()

Growth of $1 invested in each factor, replicated from CRSP CIZ (blue) against the Ken French data library factors on WRDS (orange), on a log scale. Hover over a line for the month, the value of $1 and the two monthly factor returns.

#################################
# Chart Style                   #
#################################
COLORS = {'ff': '#eb6834', 'w': '#2a78d6'}

def log_ticks(series):
    # $ ticks at 1-2-5 steps covering the plotted range
    lo, hi = min(s.min() for s in series), max(s.max() for s in series)
    steps = (1, 1.5, 2, 3, 4, 5, 7) if hi / lo < 20 else (1, 2, 5)   # finer steps for a narrow range
    vals = [m * 10**e for e in range(-2, 8) for m in steps]
    vals = [v for v in vals if lo * 0.9 <= v <= hi * 1.1]
    return dict(tickvals=vals, ticktext=[f'${v:,.0f}' if v >= 1 and v == int(v) else f'${v:g}' for v in vals])

def factor_chart(ff_col, w_col, name):
    df = _ffcomp[[ff_col, w_col]].dropna()
    growth = {c: (1 + df[c]).cumprod() for c in (w_col, ff_col)}
    corr = stats.pearsonr(_ffcomp70[ff_col], _ffcomp70[w_col])[0]
    labels = {w_col: f'Replicated {name} (CRSP CIZ)', ff_col: f'Ken French {name}'}
    colors = {w_col: COLORS['w'], ff_col: COLORS['ff']}

    fig = go.Figure()
    for c in (w_col, ff_col):
        fig.add_trace(go.Scatter(x=df.index, y=growth[c], name=labels[c], mode='lines',
                                 line=dict(color=colors[c], width=2), customdata=df[c],
                                 hovertemplate='$%{y:,.2f} (month: %{customdata:+.2%})<extra>' + labels[c] + '</extra>'))
        fig.add_annotation(x=df.index[-1], y=np.log10(growth[c].iloc[-1]),
                           text=f'<b>{labels[c].split(" (")[0]}</b><br>${growth[c].iloc[-1]:,.2f}',
                           showarrow=False, xanchor='left', xshift=8, align='left', font=dict(color='#0b0b0b'))
    fig.add_hline(y=1, line=dict(color='#8a8984', width=1, dash='dot'))
    fig.update_layout(
        title=dict(text=f'{name}: Replicated vs. Ken French<br><sup style="color:#52514e">'
                        f'{df.index[0]:%b %Y}–{df.index[-1]:%b %Y}; correlation of monthly returns since 1970: {corr:.3f}</sup>',
                   x=0, xanchor='left'),
        template='plotly_white', height=520, margin=dict(l=70, r=190, t=80, b=90),
        font=dict(family='Inter, system-ui, sans-serif', size=13, color='#0b0b0b'),
        paper_bgcolor='#fcfcfb', plot_bgcolor='#fcfcfb', hovermode='x unified',
        legend=dict(orientation='h', x=0, y=-0.12, yanchor='top'),
        xaxis=dict(showgrid=False, ticks='outside', tickcolor='#c3c2b7', linecolor='#c3c2b7',
                   range=[df.index[0], df.index[-1]]),
        yaxis=dict(type='log', title='Growth of $1 (log scale)', gridcolor='#ecebe7', zeroline=False,
                   **log_ticks(growth.values())),
    )
    return fig

factor_chart('smb', 'WSMB', 'SMB').show()

factor_chart('hml', 'WHML', 'HML').show()