[R-sig-finance] function to calc Fama-French 3 factor model

Andrew West jgalt70 at yahoo.com
Thu Sep 8 15:43:10 CEST 2005


Thanks to Gabor and Dirk, I think I fixed some of the
worst problems, but then continued on to add new
output, perhaps intoroducing new problems.

I added an extra plot showing fits vs. actual returns
for the 2 models (CAPM & FamaFrench), and added a test
of AIC scores for the 2 models, resulting in a
statement of which is the better model for the stock.

For a good example of why using Fama French for COE
estimation might be a good idea, try:
getffBeta("DUK")

Notes of minor frustrations: 
1.Couldn't figure out how to extract a single AIC for
a WLE model, or pull an Adj. R2 from a WLE summary.
(could get AIC in best subsets form, but that's not
the same thing, so had to do AIC comparisons on LM
models)

2. Couldn't figure out how to convince a "pairs" plot
to share a window with another plot (tried par,
split.screen).

3. Was hoping that the strucchange package could
automatically show me time-varying coefficients, but
somehow the plot outputs of efp(type="ME" or "RE",
functional=NULL) didn't correspond to what I thought
it might, (rolling window coefficient estimates, which
I've done before but which is another lengthy
procedure).

4. It would be nice to offer a choice of more verbose
output, e.g. with model summaries. 

5. It would be nice to only download the Fama-French
data once per session.

I also haven't acquired/learned a substitute for
MSWord for code editing, thus the function paste below
may need reformatting to avoid inappropriate line
breaks.

Here's the code:

getffBeta=function(tool) {
require(MASS);
require(tseries);
require(wle)
require(zoo)
stock= (get.hist.quote(instrument = tool, quote
=c("Ad"),compression="m",origin=as.Date(0)));
spx=(get.hist.quote(instrument = "^gspc", quote
=c("Ad"),compression="m",origin=as.Date(0)));
today=Sys.Date();
seqmo=length(seq(as.Date("1926-07-01"),today,by="month"))-2;
url="http://mba.tuck.dartmouth.edu/pages/faculty/ken.french/ftp/F-F_Research_Data_Factors.zip";
destfile=tempfile();
dataff=download.file(url, destfile, mode='wb');
unzip=unz(destfile,"F-F_Research_Data_Factors.txt");
ff4 <- read.table(unzip, header=FALSE,
sep="",na.strings="NA", dec=".", strip.white=TRUE,
skip=4,nrows=seqmo)
ffdata <- ff4
attach(ffdata);
ffdata=ffdata/100
#find out starting year month for ffdata;
ffdatats=ts(ffdata, start=c(1926,7), frequency=12);
stock=ts(na.omit(coredata(stock)), start=
as.numeric(as.yearmon(as.Date(start(stock)[1]))),
frequency=12);
spx=ts(na.remove(spx), start=
as.numeric(as.yearmon(as.Date(start(spx)[1]))),
frequency=12);
combined= na.remove(ts.union(stock,spx));
combreturn= na.remove(diff(log(combined)));
combined=na.remove(ts.intersect(combreturn,ffdatats),
names=list("stockret", "spxret",
"dates","ffmktret","smb","hml","rf"));
combined[,1]=combined[,1]-combined[,7]
combined[,2]=combined[,2]-combined[,7]

simplereg=lm(combined[,1]~combined[,2]);
stockbeta=simplereg$coef[2];
textout=paste("CAPM beta=", stockbeta);
robustreg=wle.lm(combined[,1]~combined[,2]);
robstockbeta=robustreg$coef[2];
textout2=paste("robust stock beta=", robstockbeta);

plot(as.numeric(combined[,1])~as.numeric(combined[,2]),xlab="S&P",
ylab=tool)
title(main=textout, sub=list(textout2, col="red"))
abline(coef(simplereg));
abline(coef(robustreg),col="red")

print(textout);
print(textout2);
ffreg=wle.lm(combined[,1]~combined[,4]+combined[,5]+combined[,6]);
ffbeta=ffreg$coef[2];
ffsmb=ffreg$coef[3];
ffhml=ffreg$coef[4] ;
textout3=paste('robust ff beta=',ffbeta) ;
textout4=paste('robust ff smb=',ffsmb) ;
textout5=paste('robust ff hml=',ffhml) ;
print(textout3) ;
print(textout4) ;
print(textout5) ;
rf=5
smb=1
hml=2
rm=5
ffcoe=rf+ffbeta*rm+ffsmb*smb+ffhml*hml
robcoe=rf+robstockbeta*rm
regcoe=rf+stockbeta*rm
textout6=paste('regular coe=',regcoe,' robustcoe=',
robcoe)
textout7=paste('Fama-French coe=',ffcoe)
print(textout6)
print(textout7) ;
ffreg2=lm(combined[,1]~combined[,4]+combined[,5]+combined[,6]);

a1=as.vector(combined[,1])
b1=as.vector(ffreg2$fitted)
c1=as.vector(simplereg$fitted)
plot3= (data.frame(a1,b1,c1))
win.graph()
panel.cor <- function(x, y, digits=2, prefix="",
cex.cor)
      {
          usr <- par("usr"); on.exit(par(usr))
          par(usr = c(0, 1, 0, 1))
          r <- abs(cor(x, y))
          txt <- format(c(r, 0.123456789),
digits=digits)[1]
          txt <- paste(prefix, txt, sep="")
          if(missing(cex.cor)) cex <-
0.8/strwidth(txt)
          text(0.5, 0.5, txt, cex = cex * r)
      }
pairs(plot3, upper.panel=panel.cor,labels=c("stock
returns","FF fits","CAPM fits"))


ffreg2=lm(combined[,1]~combined[,4]+combined[,5]+combined[,6]);
if((AIC(simplereg))<(AIC(ffreg2))) "CAPM superior to
FF" else "FF superior to CAPM";
}



More information about the R-sig-finance mailing list