This code is written to replicate the characteristics-based benchmarks proposed by Daniel, Grinblatt, Titman and Wermers (1997), hence the short form DGTW benchmarks. In asset pricing literature, the DGTW benchmarks are widely used to measure the performance of institutional investors based on the portfolio holdings data.
As always, the first step is to import relevant Python packages and establish connection to WRDS server.
########################################### Characteristics-Based Benchmarks (V4) ## Qingyi (Freda) Song Drechsler ## Date: May 2018 ## Updated: April 2021 ## Updated: October 2026, CRSP CIZ format ## ## V3 Note: Updated book-to-market calc ## to handle division by zero problem ## caused by 0 marketcap ## V4 Note: Fixed one permno to 2 dgtw ## portfolio assignment problem ###########################################import pandas as pdimport numpy as npimport datetime as dtimport wrdsimport plotly.graph_objects as goimport plotly.io as piofrom dateutil.relativedelta import*from pandas.tseries.offsets import*pio.renderers.default ='notebook_connected'#################### Connect to WRDS ####################conn = wrds.Connection()
Loading library list...
Done
Pricing related information (return, price, shares outstanding) is obtained from CRSP using the following “raw_sql” step.
#################### CRSP Block ##################### sql similar to crspmerge macro# CRSP CIZ format: msf_v2 carries the security attributes on each monthly# record, so no separate names-file merge is needed.# Common stocks only: equivalent to legacy shrcd = 10 or 11.# No exchange restriction, as in the SAS version.crsp_m = conn.raw_sql(""" select distinct permno, permco, cusip, mthcaldt, siccd, primaryexch, mthret, mthvol, shrout, mthprc, mthcumfacpr, mthcumfacshr from crsp_m_stock.msf_v2 where mthcaldt between '01/01/1970' and '06/30/2026' and sharetype='NS' and securitytype='EQTY' and securitysubtype='COM' and usincflg='Y' and issuertype in ('ACOR', 'CORP') """, date_cols=['mthcaldt'],)# change variable format to intcrsp_m[['permco', 'permno']] = crsp_m[['permco', 'permno']].astype(int)# Line up date to be end of monthcrsp_m['jdate'] = crsp_m['mthcaldt'] + MonthEnd(0)# price adjustedcrsp_m['p'] = crsp_m['mthprc'].abs() / crsp_m['mthcumfacpr']# total shares out adjustedcrsp_m['tso'] = crsp_m['shrout'] * crsp_m['mthcumfacshr'] *1e3# market cap in $milcrsp_m['me'] = crsp_m['p'] * crsp_m['tso'] /1e6# sum of me across different permno belonging to same permco a given datecrsp_summe = ( crsp_m.groupby(['jdate', 'permco'])['me'] .sum() .reset_index() .rename(columns={'me': 'me_comp'}))crsp_m = pd.merge(crsp_m, crsp_summe, how='inner', on=['jdate', 'permco'])
Compustat Fundamental Annual contains the financial statement related information needed for calculating book value of equity.
#################### Compustat Block ####################comp = conn.raw_sql(""" select gvkey, datadate, cusip, sich, seq, pstkrv, pstkl, pstk, txdb, itcb from comp.funda where indfmt='INDL' and datafmt='STD' and popsrc='D' and consol='C' and datadate >= '01/01/1970' """, date_cols=['datadate'],)comp['year'] = comp['datadate'].dt.year# select sample where seq is greater than 0comp = comp.loc[comp['seq'] >0]# create preferrerd stock:# 1st choice: Preferred stock - Redemption Value# 2nd choice: Preferred stock - Liquidating Value# 3rd choice: Preferred stock - Carrying Value, Stock (Capital) - Totalcomp['pref'] = np.where( comp['pstkrv'].isnull(), comp['pstkl'], comp['pstkrv'])comp['pref'] = np.where(comp['pref'].isnull(), comp['pstk'], comp['pref'])comp['pref'] = np.where(comp['pref'].isnull(), 0, comp['pref'])# fill in missing values for deferred taxes and investment tax creditcomp['txdb'] = comp['txdb'].fillna(0)comp['itcb'] = comp['itcb'].fillna(0)# create book equity# Daniel and Titman (JF 1997):# BE = stockholders' equity + deferred taxes + investment tax credit# - Preferred Stockcomp['be'] = comp['seq'] + comp['txdb'] + comp['itcb'] - comp['pref']# keep only records with non-negative book equitycomp = comp.loc[comp['be'] >=0]comp = comp[['gvkey', 'datadate', 'year', 'be', 'sich']]
To link these two sources together, CRSP provides a handy linking table to join the key identifiers, permco and gvkey, from both databases.
########################## Add Historical PERMCO ##########################ccm = conn.raw_sql(""" select gvkey, lpermno as permno, lpermco as permco, linktype, linkprim, linkdt, linkenddt from crsp_m_ccm.ccmxpf_linktable where (linktype ='LU' or linktype='LC') and linkprim in ('P', 'C') and usedflag=1 """, date_cols=['linkdt', 'linkenddt'],)# if linkenddt is missing then set to today dateccm['linkenddt'] = ccm['linkenddt'].fillna(pd.to_datetime('today'))ccm1 = pd.merge(comp, ccm, how='left', on=['gvkey'])ccm1['jdate'] = ccm1['datadate'] + MonthEnd(0)ccm1['year'] = ccm1.datadate.dt.year# Impose date rangesccm2 = ccm1[ (ccm1['datadate'] >= ccm1['linkdt'])& (ccm1['datadate'] <= ccm1['linkenddt'])].drop(columns=['linktype', 'linkdt', 'linkenddt'])comp1 = pd.merge( ccm2, crsp_m[ ['permno', 'permco', 'mthcaldt', 'jdate', 'siccd', 'me', 'me_comp'] ], how='inner', on=['permco', 'permno', 'jdate'],)# Edit December 2020: correct the division by 0 problem caused by# missing market cap# Use a conditioning statement to address this problem# Old Code: comp3['bm']=comp3['be'].div(comp3['me_comp'])# New calculationcomp1['bm'] = np.where(comp1.me_comp >0, comp1.be / comp1.me_comp, np.nan)comp2 = comp1.sort_values( ['permno', 'year', 'datadate', 'linkprim', 'bm']).drop_duplicates()# pick max datadate for a given permno year combo# (firm changes fiscal period)maxdatadate = ( comp2.groupby(['permno', 'year'])['datadate'].max().reset_index())comp3 = pd.merge( comp2, maxdatadate, how='inner', on=['permno', 'year', 'datadate'])
Below is a somewhat clumsy function to assign stocks into Fama-French 48 industry classification based on the historical SIC codes. Please refer to Ken French’s site for a detailed text file describing the classification scheme.
The function “ffi48” is then applied to the data frame containing joint information of CRSP and Compustat created in the previous step to assign each stock to its proper Fama-French 48 Industry group.
# assign SIC codecomp4 = comp3# First use historical Compustat SIC Code# Then if missing use historical CRSP SIC Codecomp4['sic'] = np.where(comp4['sich'] >0, comp4['sich'], comp4['siccd'])# and adjust some SIC code to fit F&F 48 ind delineationcomp4['sic'] = np.where( (comp4['sic'].isin([3990, 9995, 9997]))& (comp4['siccd'] >0)& (comp4['sic'] != comp4['siccd']), comp4['siccd'], comp4['sic'],)comp4['sic'] = np.where( comp4['sic'].isin([3990, 3999]), 3991, comp4['sic'])comp4['sic'] = comp4.sic.astype(int)# assign the ffi48 function to comp4_sic = comp4['sic'].unique()_sicff = pd.DataFrame(_sic).rename(columns={0: 'sic'})_sicff = _sicff.apply(ffi48, axis=1)comp4 = pd.merge(comp4, _sicff, how='left', on=['sic'])# keep only records with non-missing bm and ffi48 classificationcomp4 = comp4[(comp4['bm'] != np.NaN) & (comp4['ffi48_desc'] !='')]comp4 = comp4.drop(['sich', 'siccd', 'datadate'], axis=1)comp4 = comp4.sort_values(['ffi48', 'year'])
Now that we have created an industry classification scheme, the rest are quite straightforward. AS DGTW benchmarks performance based on value, momentum and size of the stocks, we need to calculate the industry averages of these metrics.
Let’s start with the Book-to-Market ratio.
########################## Industry BM Average ########################### Calculate BM Industry Average Each Periodcomp4_tmp = comp4[(comp4['ffi48'] >0) & (comp4['bm'] >=0)]bm_ind = ( comp4_tmp.groupby(['ffi48', 'year'])['bm'] .mean() .reset_index() .rename(columns={'bm': 'bmind'}))# Calculate Long-Term Industry BtM Averagebm_ind['n'] = bm_ind.groupby(['ffi48'])['year'].cumcount()bm_ind['sumbm'] = bm_ind.groupby(['ffi48'])['bmind'].cumsum()bm_ind['bmavg'] = bm_ind['sumbm'] / (bm_ind['n'] +1)bm_ind = bm_ind.drop(['n', 'sumbm'], axis=1)# Adjust Firm-Specific BtM with Industry Averagescomp5 = pd.merge(comp4, bm_ind, how='left', on=['ffi48', 'year'])comp5['bm_adj'] = comp5['bm'] - comp5['bmavg']
Next is to create the momentum factors for the industry groups.
########################## Momentum Factor ########################### Create (12,1) Momentum Factor with at least 6 months of returns_tmp_crsp = ( crsp_m[['permno', 'mthcaldt', 'mthret', 'me', 'primaryexch']] .sort_values(['permno', 'mthcaldt']) .set_index('mthcaldt'))# replace missing return with 0_tmp_crsp['mthret'] = _tmp_crsp['mthret'].fillna(0)# CIZ returns include delisting returns: a -100% month gives log(0) = -inf,# which correctly yields a -100% cumulative returnwith np.errstate(divide='ignore'): _tmp_crsp['logret'] = np.log(1+ _tmp_crsp['mthret'])_tmp_cumret = ( _tmp_crsp.groupby(['permno'])['logret'] .rolling(12, min_periods=7) .sum())_tmp_cumret = _tmp_cumret.reset_index()_tmp_cumret['cumret'] = np.exp(_tmp_cumret['logret']) -1sizemom = pd.merge( _tmp_crsp.reset_index(), _tmp_cumret[['permno', 'mthcaldt', 'cumret']], how='left', on=['permno', 'mthcaldt'],)sizemom['mom'] = sizemom.groupby('permno')['cumret'].shift(1)sizemom = ( sizemom[sizemom['mthcaldt'].dt.month ==6] .drop(['logret', 'cumret'], axis=1) .rename(columns={'me': 'size'}))
Third is to create size portfolios based on the NYSE stocks size breakpoints.
########################## NYSE Size Breakpoint ########################### Get Size Breakpoints for NYSE firmssizemom = sizemom.sort_values(['mthcaldt', 'permno']).drop_duplicates()nyse = sizemom[sizemom['primaryexch'] =='N']nyse_break = ( nyse.groupby(['mthcaldt'])['size'] .describe(percentiles=[0.2, 0.4, 0.6, 0.8]) .reset_index())nyse_break = nyse_break[['mthcaldt', '20%', '40%', '60%', '80%']].rename( columns={'20%': 'dec20','40%': 'dec40','60%': 'dec60','80%': 'dec80', })sizemom = pd.merge(sizemom, nyse_break, how='left', on='mthcaldt')# Add NYSE Size Breakpoints to the Datadef size_group(row):if0<= row['size'] < row['dec20']: value =1elif row['size'] < row['dec40']: value =2elif row['size'] < row['dec60']: value =3elif row['size'] < row['dec80']: value =4elif row['size'] >= row['dec80']: value =5else: value = np.nanreturn valuesizemom['group'] = sizemom.apply(size_group, axis=1)sizemom['year'] = sizemom['mthcaldt'].dt.year -1sizemom = sizemom[ ['permno', 'mthcaldt', 'year', 'mom', 'group', 'size', 'mthret']]# Adjusted BtM from the calendar year preceding the formation datecomp6 = comp5[['gvkey', 'permno', 'year', 'bm_adj']]comp6 = pd.merge(comp6, sizemom, how='inner', on=['permno', 'year'])comp6 = comp6.dropna(subset=['size', 'mom', 'bm_adj', 'mthret'], how='any')
The final step involves calculating size, book-to-market, and momentum controlled value-weighted portfolio returns. I really wish Python can add a value-weighted option in the mean() function. How I miss SAS “proc means” statement’s weight=size capacity.
For those of you who have already created your DGTW excess return using the sample SAS code available on WRDS website Research Application section, you may want to compare the output of this Python code to that. When this comparison was run on the original (legacy CRSP) version of this code, the correlation between the two dgtw_vwret series was 99.4%, with p-value equals 0.0.
I also plot the corner portfolios below, as the growth of $1 invested in each from July 1971 through June 2026. The DGTW portfolio is ordered by Size, Book-to-Market and Momentum, so portfolio 111 holds the smallest stocks with the lowest book-to-market and momentum, and portfolio 555 the largest stocks with the highest.
# Growth of $1 in the corner DGTW benchmark portfolios# dgtw_port = size quintile + book-to-market quintile + momentum quintileCORNERS = [ ('111', 'Small, low B/M, low momentum', '#2a78d6'), ('115', 'Small, low B/M, high momentum', '#eb6834'), ('551', 'Big, high B/M, low momentum', '#1baf7a'), ('555', 'Big, high B/M, high momentum', '#eda100'),]_w = dgtw_vwret.pivot( index='mthcaldt', columns='dgtw_port', values='dgtw_vwret').sort_index()fig = go.Figure()for port, label, color in CORNERS: growth = (1+ _w[port].fillna(0)).cumprod() fig.add_trace( go.Scatter( x=growth.index, y=growth, name=f'{port}: {label}', mode='lines', line=dict(color=color, width=2), customdata=_w[port], hovertemplate='$%{y:,.2f} (month: %{customdata:+.1%})''<extra>'+ port +'</extra>', ) ) fig.add_annotation( x=growth.index[-1], y=np.log10(growth.iloc[-1]), text=f'<b>{port}</b> '+ (f'${growth.iloc[-1]:,.0f}'if growth.iloc[-1] >=10elsef'${growth.iloc[-1]:,.2f}' ), showarrow=False, xanchor='left', xshift=8, font=dict(size=12, color='#0b0b0b'), )_g = pd.concat([(1+ _w[p].fillna(0)).cumprod() for p, _, _ in CORNERS])ticks = [10**e for e inrange(-2, 8) if _g.min() /10<=10**e <= _g.max() *10]fig.update_layout( title=dict( text='DGTW Benchmark Portfolios: Growth of $1''<br><sup style="color:#52514e">Corner portfolios of the ''size / book-to-market / momentum sort, value-weighted, 'f'{_w.index[0]:%b %Y}\u2013{_w.index[-1]:%b %Y}</sup>', x=0, xanchor='left', ), template='plotly_white', height=540, margin=dict(l=70, r=110, t=90, b=130), 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=[_w.index[0], _w.index[-1]], ), yaxis=dict(type='log', title='Growth of $1 (log scale)', tickvals=ticks, ticktext=[f'${v:,.2g}'if v <1elsef'${v:,.0f}'for v in ticks], gridcolor='#ecebe7', zeroline=False, ),)fig.show()