ME_IA_ROE

My R Code

Load packages

options(scipen=999) # shut down scientific notation
library(tidyverse);library(Hmisc)
rm(list=ls());gc()

##           used  (Mb) gc trigger  (Mb) max used  (Mb)
## Ncells 2134453 114.0    4300762 229.7  2766532 147.8
## Vcells 3575815  27.3    8388608  64.0  7008183  53.5

Load data

CRSP

You will need to download the data from CRSP, to which WRDS has access. The data items used here are given below.

CRSP Permanent Number (PERMNO) (select PERMNO as company codes)
CRSP Permanent Company Number (permco)
Share Code (shrcd)
Exchange Code (exchcd)
Primary Exchange (primexch)
Price (prc)
Holding Period Return (ret)
Number of Shares Outstanding (shrout)
Delisting Return (dlret)

Note the difference between PERMNO and PERMCO, where PERMNO is security-specific and PERMCO is company-specific. For example, Google has two class of shares, therefore Google has one PERMCO but two PERMNO.

CCM

In addition to CRSP, I use the CRSP-COMPUSTAT-Merged link table (CCM) to connect CRSP to COMPUSTAT. CCM is accessible under CRSP in WRDS. I download the entire link table as it is not a big file. Even though constructing momemntum factor does not require combining CRSP and COMPUSTAT, I apply it to keep consistency.

Code

# Load data
crsp <- readRDS("Data/CRSP_2022.rds")
cmp_a <- readRDS("Data/COMPUSTAT_2022.rds")
cmp_q <- readRDS("Data/COMPUSTAT_Qtr_2022.rds")
ccm <- read.csv("Data/CCM.csv")

# convert all column names to lower case to keep consistent
names(crsp) <- tolower(names(crsp))
names(cmp_a) <- tolower(names(cmp_a))
names(cmp_q) <- tolower(names(cmp_q))
names(ccm) <- tolower(names(ccm))

Process CRSP

This section deals with the raw CRSP data. It includes converting date format, computing firm size, and filling missing returns with zeros.

One can argue about whether filling missing returns with zeros is reasonable, as there can be many different reasons for a return to be missing. If a stock is suspended from trading within month t while it is selected in your portfolio (long or short) at the end of month t-1, it will have a missing return but it takes a non-zero weight in your portfolio in month t. Without further information, taking a missing return as zero would make more sense than ignoring it in your portfolio in this case.

You can comment out the three-line code to ignore delisting returns, though it does not make any difference.

Code

######################
##   process crsp   ##
######################
df <- crsp

df$date.o <- df$date
df$date <- as.integer(substr(df$date, 1, 6))
df$year <- as.integer(substr(df$date, 1, 4))
df$month <- as.integer(substr(df$date, 5, 6))
# a minus sign before PRC means it is the average of bid and ask prices
df$prc <- abs(df$prc) # absolute value of price 
df$me <- df$prc * df$shrout # firm ME (ME--market equity)
df$ret <- as.numeric(df$ret) # NAs will be introduced by coercion
df$dlret <- as.numeric(df$dlret)

## Deal with Delist returns and NA returns
df$ret.o <- df$ret # keep the original return data
df[is.na(df$ret), 'ret'] <- 0
df[is.na(df$dlret), 'dlret'] <- 0
# Take the delist return into account
df$ret <- (1+df$ret)*(1+df$dlret) - 1

# lag me 
df <- df %>%
  group_by(permno) %>%
  arrange(date, .by_group=T) %>%
  mutate(
    me.lag = dplyr::lag(me, 1),
    date.dif = date - dplyr::lag(date,1)
  )
# df$me.lag[!df$date.dif %in% c(1,89) & !is.na(df$date.dif)] <- NA

# screening stocks following conventions
df <- df[(df$shrcd%in%c(10,11)) & (df$primexch%in%c('N','A','Q')),]

# portfolios are formed at the end of June in each year
df$year.ff <- df$year
df$year.ff[df$month<=6] <- df$year.ff[df$month<=6]-1

# save 
crsp.ready <- df

Process CCM and merge CRSP and CCM

This section merges CRSP with CCM. You can find more information about how to use CCM on WRDS: https://wrds-www.wharton.upenn.edu/pages/classroom/using-crspcompustat-merged-database/.

After combining the two datasets, I use gvkey as the firm id to construct portfolios. This procedure can help keep consistent with the other factor constructions.

Code

#####################
##   process ccm   ##
#####################
df <- ccm

# "NU" "LU" "NR" "LC" "LX" "LD" "LN" "LS" "NP"
df <- df[df$linktype %in% c('LC','LU'), ]
df <- df[df$linkprim %in% c('P','C','J'),]
df <- df[c('lpermno','gvkey','cik','linkdt','linkenddt')]
colnames(df)[1] <- 'permno'
df$linkenddt <- as.integer(df$linkenddt)
df$linkenddt <- ifelse(is.na(df$linkenddt), Inf, df$linkenddt)

# merge crsp and ccm
df <- merge(crsp.ready, df, by='permno')
# the ccm link is only valid between linkdt and linkenddt
df <- df[df$date.o >= df$linkdt & df$date.o <= df$linkenddt, ]
df$linkdt <- NULL
df$linkenddt <- NULL

# sum up the size (me) of different permno under the same gvkey
df <- df %>%
  group_by(gvkey, date) %>%
  mutate(
    me.sum = sum(me),
    me.sum.lag = sum(me.lag)
  )

df[c("permco","cusip",'shrcd','shrcls','dlret','retx',"shrout")] <- NULL

# save
crsp_ccm <- df

June in year t (dataframe pt6) and December in year t-1

Code

##################
## FF Procedure ##
##################
# ME in June t-1
df <- crsp_ccm[crsp_ccm$month==6, c("gvkey","year","me.sum","primexch")]
df$me.jun <- df$me.sum
df <- df[!duplicated(df[c('gvkey','year')]),]
pt6 <- df

# need ME.Dec as to cut BEME
# ME in Dec t-1
df <- crsp_ccm[crsp_ccm$month==12, c("gvkey","year","me.sum")]
df <- df[!duplicated(df[c('gvkey','year')]),]
df$year <- df$year + 1
names(df)[3] <- 'me.dec'
pt6 <- merge(pt6, df, by=c('gvkey','year'), all.x = T)

Compustat Annual Data (dataframe cmp_a)

Code

#################################
##   process compustat annual  ##
#################################
df <- cmp_a

# index is from 1 to n
df <- df %>%
  group_by(gvkey) %>%
  arrange(fyear, .by_group=T) %>%
  mutate(
    index = 1:n(),
    year.dif = fyear - dplyr::lag(fyear,1)
  )
# remove the first two obs for each firm due the back-filling practice
# df <- df[df$index>2,] # remove the first two obs

# # Book value (annual) according to Fama-French (1993)
# # create shareholders' equity (seq)
# df$seq <- ifelse(is.na(df$seq), 0, df$seq)
# # create preferred stock (ps)
# # (pstkrv, pstkl, pstk, 0)
# df$ps <- ifelse(is.na(df$pstkrv), df$pstkl, df$pstkrv)
# df$ps <- ifelse(is.na(df$ps), df$pstk, df$ps)
# df$ps <- ifelse(is.na(df$ps), 0, df$ps)
# # fill NA of txditcq with 0
# df$txditc <- ifelse(is.na(df$txditc), 0, df$txditc)
# # create book equity (be)
# df$be <- df$seq + df$txditc - df$ps

# The easy way to get book equity
df$be <- df$ceq

# lagged at
df <- df %>%
  group_by(gvkey) %>%
  arrange(fyear, .by_group = T) %>%
  mutate(
    across(all_of(c('at')), ~dplyr::lag(.x, 1), .names = "{.col}.lag")
  )
df$at.lag[df$year.dif!=1&!is.na(df$year.dif)] <- NA

# RMW factor
# to be included, a firm needs a non-missing revt 
#   and non-missing value for one of ('cogs','xsga','xint')
# roe.ok == T means a firm is qualified to be included 
df$revt.ok <- ifelse(!is.na(df$revt), T, F) # T means revt is ok
df$roe.ok <- rowSums(!is.na(df[c('cogs','xsga','xint')]))>=1 # T means others are ok
df$roe.ok <- df$roe.ok + df$revt.ok == 2 # T mean both requirements are satisfied

df$revt[is.na(df$revt)] <- 0
df$cogs[is.na(df$cogs)] <- 0
df$xsga[is.na(df$xsga)] <- 0
df$xint[is.na(df$xint)] <- 0

## FF profitability: (revt-cogs-xsga-xint)/be
df$roe.ff <- (df$revt-df$cogs-df$xsga-df$xint)/df$be

########
df$dat <- (df$at - df$at.lag)
df$dat__at <- df$dat / df$at.lag

# we remove the obs with at.lag=0
df$dat__at[df$at.lag==0&!is.na(df$at.lag)] <- NA

## Why use the fyear + 1? Refer to FF (2015)
df$aval_year <- df$fyear + 1

# save
df.save <- df

# fundamental data (df)  is on fyear+1 and 
#        crsp data (pt6) is on June of each calendar year
df <- merge(pt6,
            df.save[c("gvkey","aval_year","index",
                      "be","roe.ff","roe.ok","dat__at",'sic')],
            by.x=c('gvkey','year'),
            by.y=c('gvkey','aval_year'), 
            all.x = T)
# Size cut uses me in June t, beme uses me in Dec t-1
df$me <- df$me.jun
df$beme <- df$be / df$me.dec
pt6_cmp <- df

Compustat Quarterly Data (dataframe cmp_q)

Code

##############################
##   process compustat qtr  ##
##############################
df <- cmp_q

# Book value (annual) according to Hou, Xue, and Zhang (2015)
# create shareholders' equity (seq)
# Depending on availability, we use 
# stockholders' equity (item SEQQ), 
df$se <- df$seqq
# or common equity (item CEQQ) plus the carrying value of preferred stock (item PSTKQ), 
df$se <- ifelse(is.na(df$se), df$ceqq + df$pstkq, df$se)
# or total assets (item ATQ) minus total liabilities (item LTQ) 
df$se <- ifelse(is.na(df$se), df$atq - df$ltq, df$se)
# in that order as shareholders' equity.
# We use redemption value (item PSTKRQ) if available, 
# or carrying value for the book value of preferred stock.
df$ps <- ifelse(is.na(df$pstkrq), df$pstkq, df$pstkrq)
# df$pstkq
# df$ps
# df$pstkrq

# In particular, book equity is 
# shareholders'equity, 
# plus balance-sheet deferred taxes and investment tax credit (item TXDITCQ) if available, 
# minus the book value of preferred stock.
df$txditcq[is.na(df$txditcq)] <- 0
#
df$be <- df$se + df$txditcq - df$ps

# lagged at
df <- df %>%
  group_by(gvkey) %>%
  arrange(datadate, .by_group = T) %>%
  mutate(
    across(all_of(c('be')), ~dplyr::lag(.x, 1), .names = "{.col}.lag")
  )

## ROE
df$ibq[is.na(df$ibq)] <- 0
df$roe <- df$ibq/df$be.lag
df$roe[df$be.lag<=0&!is.na(df$be.lag)] <- NA

#
df$date.rdq <- as.Date(paste0(df$rdq), '%Y%m%d')
df$date.fqr <- as.Date(paste0(df$datadate), '%Y%m%d')
df$date.dif <- df$date.rdq - df$date.fqr
df$roe[df$date.dif>190&!is.na(df$date.dif)] <- NA

# forward one month 
#   because if roe is available in June
#   I can only use it starting July, i.e., the end of June
df$date.rdq.plus1 <- df$date.rdq+months(1)
df$year <- substr(df$date.rdq.plus1,1,4)
df$month <- substr(df$date.rdq.plus1,6,7)
df$date <- as.integer(paste0(df$year, df$month))
#
df <- df[!is.na(df$date), c("gvkey","date","datadate","roe")]

## order by datadate
df <- df[order(df$datadate, decreasing=TRUE),]
## remove duplicated obs
df <- df[!duplicated(df[c('gvkey','date')]),]

# there should be zero duplications
sum(duplicated(df[c("gvkey","date")])) # 0

## [1] 0

df$datadate <- NULL
#
tmp <- crsp_ccm[c("gvkey",'date','primexch')]
tmp <- tmp[!duplicated(tmp),]

ptq <- merge(tmp, df, by=c('gvkey','date'), all.x = T)
ptq <- ptq %>% fill(roe)

Form portfolios and compute factor returns

This section cuts firms into 2 ME groups and 3 IA groups independently at the end of June in year t, and independently 3 ROE groups at the end of each month.
Note that Fama-French and Hou-Xue-Zhang always use NYSE cut-points to form portfolios, hence I use [primexch==“N”] when calculating the cut-points to take care of this convention.

Code

#########################################
## Form portfolios and compute returns ##
#########################################
year.st = 1972
# year.ed = 2021

# Financial firms and 
#  firms with negative book equity are excluded
fin.ex = T

# variables used
var.test.1 <- 'me'
var.test.2 <- 'dat__at'
var.test.3 <- 'roe'

#######################
df <- pt6_cmp
# df <- df[df$index>=2,]

# Fama-French definition for 12 industries
if(fin.ex==T){
  # df <- df[df$sic<4900|df$sic>4949,]
  df <- df[df$sic<6000|df$sic>6999,]
}
df <- df[df$year>=year.st,]
# df <- df[df$year<=year.ed,]

# me cut
pt.1 <- df[!is.na(df[,var.test.1]),] 
# Rank var.1 annual
pt.1 <- pt.1 %>%
  group_by(year) %>%
  mutate(
    across(all_of(var.test.1), ~findInterval(.x, quantile(.x[primexch=="N"], probs=c(0.5)))+1, .names = "rank.var1")
  )

# ia cut
pt.2 <- df[!is.na(df[,var.test.2]),] 
# Rank var.2 annual
pt.2 <- pt.2 %>%
  group_by(year) %>%
  mutate(
    across(all_of(var.test.2), ~findInterval(.x, quantile(.x[primexch=="N"], probs=c(0.3,0.7)))+1, .names = "rank.var2")
  )

# roe cut
pt.3 <- ptq[!is.na(ptq[,var.test.3]),]
# Rank var.3 monthly!
pt.3 <- pt.3 %>%
  group_by(date) %>%
  mutate(
    across(all_of(var.test.3), ~findInterval(.x, quantile(.x[primexch=="N"], probs=c(0.3,0.7)))+1, .names = "rank.var3")
  )

# merge back to CRSP
df <- merge(crsp_ccm[c("gvkey","year.ff","date","ret","me.lag")], 
            pt.1[c('gvkey','year','rank.var1')],
            by.x=c('gvkey','year.ff'), by.y=c('gvkey','year'), all.x = T)

df <- merge(df[c("gvkey","year.ff","date","ret","me.lag",'rank.var1')], 
            pt.2[c('gvkey','year','rank.var2')],
            by.x=c('gvkey','year.ff'), by.y=c('gvkey','year'), all.x = T)

df <- merge(df[c("gvkey","date","ret","me.lag",'rank.var1','rank.var2')], 
            pt.3[c('gvkey','date','rank.var3')],
            by=c('gvkey','date'), all.x = T)

df <- df[!is.na(df$me.lag),]
df <- df[!is.na(df$rank.var1),]
df <- df[!is.na(df$rank.var2),]
df <- df[!is.na(df$rank.var3),]

# factor: 2-by-3 portfolio
pt.factor <- df %>%
  group_by(date, rank.var1, rank.var2, rank.var3) %>%
  summarise(
    # nstocks = n(),
    ret_vw = weighted.mean(ret, me.lag, na.rm = TRUE)*100,
    # ret_ew = mean(ret, na.rm = TRUE)*100,
    .groups='drop')
pt.factor$group <- paste(pt.factor$rank.var1, pt.factor$rank.var2, pt.factor$rank.var3)
pt.factor <- reshape2::dcast(pt.factor, date ~ group, value.var = 'ret_vw')

#
res <- pt.factor

res$size <- (res$`1 1 1`+res$`1 1 2`+res$`1 1 3` + 
               res$`1 2 1`+res$`1 2 2`+res$`1 2 3` + 
               res$`1 3 1`+res$`1 3 2`+res$`1 3 3`)/9 - 
  (res$`2 1 1`+res$`2 1 2`+res$`2 1 3` + 
     res$`2 2 1`+res$`2 2 2`+res$`2 2 3` + 
     res$`2 3 1`+res$`2 3 2`+res$`2 3 3`)/9

res$ia <- (res$`1 1 1`+res$`1 1 2`+res$`1 1 3` + 
             res$`2 1 1`+res$`2 1 2`+res$`2 1 3`)/6 - 
  (res$`1 3 1`+res$`1 3 2`+res$`1 3 3` + 
     res$`2 3 1`+res$`2 3 2`+res$`2 3 3`)/6

res$roe <- (res$`1 1 3`+res$`1 2 3`+res$`1 3 3` + 
              res$`2 1 3`+res$`2 2 3`+res$`2 3 3`)/6 - 
  (res$`1 1 1`+res$`1 2 1`+res$`1 3 1` + 
     res$`2 1 1`+res$`2 2 1`+res$`2 3 1`)/6

res <- res[order(res$date),]
# # save as csv
# write.csv(res,
#           paste0('hxz_factors.csv'),
#           row.names = F)

Verify the results

This section verifies the results. You need to download it from their website and save it in your “Data” folder.

Code

##########################
### Verify the results ###
##########################
# load HXZ factor
hxz <- read.csv("Data/HXZ_All_Factors.csv")
names(hxz) <- tolower(names(hxz))

# names(res)
res <- res[c("date","size","ia","roe")]
names(res)[2] <- 'me'
res$year <- as.integer(substr(res$date,1,4))
res$month <- as.integer(substr(res$date,5,6))
# res$date <- NULL

res <- merge(res, hxz, by=c('year','month'))
res <- res[order(res$date),]

Print summary and plot the cumulative returns

This section plots the cumulative returns of my result (with .rp) and factor returns from the authors’ websites. It helps visualize the consistency and it concludes my replication.

Code

#
res.sum <- data.frame(factor = c('r_me','r_ia','r_roe'),
                      corr = NA, mean.hxz = NA, mean.rp = NA,
                      t_stat.hxz = NA, t_stat.rp = NA)
#
for (i in 1:length(res.sum$factor)) {
  # i=1
  #
  i <- res.sum$factor[i]
  i.rp <- substr(i, 3, nchar(i))
  #
  res.sum$corr[res.sum$factor==i] <- cor(res[i], res[i.rp])
  
  tmp <- t.test(res[i])
  res.sum$mean.hxz[res.sum$factor==i] <- tmp$estimate
  res.sum$t_stat.hxz[res.sum$factor==i] <- tmp$statistic
  
  tmp <- t.test(res[i.rp])
  res.sum$mean.rp[res.sum$factor==i] <- tmp$estimate
  res.sum$t_stat.rp[res.sum$factor==i] <- tmp$statistic
  
  ############################
  ## Plot cumulative return ##
  ############################
  res.p <- res[c('date',i,i.rp)]
  # convert to R date format
  res.p$date <- as.Date(paste0(res.p$date,'01'), '%Y%m%d')
  res.p$date <- res.p$date + months(1) - 1
  
  # transform wide table into long table
  res.p <- reshape2::melt(res.p, id='date')
  names(res.p) <- c("Date","Factor","Return")
  # res.p <- res.p[order(res.p$Date),]
  
  # compute cumulative return
  res.p <- res.p %>%
    group_by(Factor) %>% 
    arrange(Date, .by_group = T) %>%
    mutate(
      Return.Cum = cumprod(1+Return/100)
    )
  
  # cut 1950-2021 into 4 subsamples
  res.p$Year <- as.integer(substr(res.p$Date,1,4))
  res.p$Subsample <- cut2(res.p$Year, g=4)
  
  # ggplot
  fig <- ggplot(data = res.p, aes(x=Date, y=Return.Cum, group=Factor)) +
    geom_line(aes(linetype=Factor, color=Factor)) +
    scale_linetype_manual(values=c('solid',"longdash","dotted","dashed","twodash"))+
    scale_color_manual(values=c('blue','brown','black','black','coral4'))+
    ylab('Portfolio worth if invest $1 at the beginning') +
    scale_x_date(date_labels = "%m-%Y", date_breaks = "4 year") +
    theme_bw() +
    theme(legend.title = element_blank()
    ) + 
    facet_wrap(vars(Subsample),
               scales = "free") + 
    scale_y_continuous() + 
    labs(title=paste('Cumulative return of factors from',min(res$date)))
  print(fig)
  # ggsave(paste0('cumulative_return_',i,'.pdf'), width = 12, height = 7)
  # ggsave(paste0('cumulative_return_',i,'.png'), width = 12, height = 7)
  
}

res.sum[,-1] <- round(res.sum[,-1], 2)
res.sum

##   factor corr mean.hxz mean.rp t_stat.hxz t_stat.rp
## 1   r_me 0.99     0.23    0.23       1.84      1.81
## 2   r_ia 0.95     0.32    0.35       4.18      4.19
## 3  r_roe 0.97     0.50    0.48       4.61      4.60

# write.csv(res.sum, 'factor_replicating_results_hxz.csv', row.names = F)

ME_IA_ROE_HXZ4

Kai Wang

2/22/2022

To Start

My R Code

Load packages

Load data

Process CRSP

Process CCM and merge CRSP and CCM

June in year t (dataframe pt6) and December in year t-1

Compustat Annual Data (dataframe cmp_a)

Compustat Quarterly Data (dataframe cmp_q)

Form portfolios and compute factor returns

Verify the results

Print summary and plot the cumulative returns