Archive

Archive for the ‘Factors’ Category

Weekend Reading – Facebook’s P/E ratio

The Barron’s article Still Too Pricey by Andrew Bary looks at the share price of the Facebook and based on the P/E ration valuation metrics concludes that even at the current prices, stock is overvalued. I want to show how to do this type of fundamental analysis using the Systematic Investor Toolbox.

First let’s load historical prices and earnings per share (EPS) for Facebook and a few stocks in the technology sector: LinkedIn, Groupon, Apple, and Google.

###############################################################################
# Load Systematic Investor Toolbox (SIT)
# http://systematicinvestor.wordpress.com/systematic-investor-toolbox/
###############################################################################
setInternet2(TRUE)
con = gzcon(url('http://www.systematicportfolio.com/sit.gz', 'rb'))
    source(con)
close(con)

	#*****************************************************************
	# Load historical fundamental and pricing data
	#****************************************************************** 
	load.packages('quantmod') 
	tickers = spl('FB,LNKD,GRPN,AAPL,GOOG')
	tickers.temp = spl('NASDAQ:FB,NYSE:LNKD,NASDAQ:GRPN,NASDAQ:AAPL,NASDAQ:GOOG')
	
	# get fundamental data
	data.fund <- new.env()
	for(i in 1:len(tickers)) {
			cat(tickers[i],'\n')
			data.fund[[tickers[i]]] = fund.data(tickers.temp[i], 80)
	}
	
	# get pricing data
	data <- new.env()
	getSymbols(tickers, src = 'yahoo', from = '1970-01-01', env = data, auto.assign = T)
		for(i in ls(data)) data[[i]] = adjustOHLC(data[[i]], use.Adjusted=T)			

Next, let’s combine fundamental and pricing data and create P/E ratio for all stocks.

	#*****************************************************************
	# Combine fundamental and pricing data
	#****************************************************************** 				
	for(i in tickers) {
		fund = data.fund[[i]]
		fund.date = date.fund.data(fund)
			
		# Earnings per Share		
		EPS = 4 * get.fund.data('Diluted EPS from Total Operations', fund, fund.date)
		if(nrow(EPS) > 3)
			EPS = rbind(EPS[1:3], get.fund.data('Diluted EPS from Total Operations', fund, fund.date, is.12m.rolling=T)[-c(1:3)])
		
		# merge	
		data[[i]] = merge(data[[i]], EPS)
	}

	bt.prep(data, align='keep.all', dates='1995::')
	
	#*****************************************************************
	# Create PE
	#****************************************************************** 
	prices = data$prices
		prices = bt.apply.matrix(prices, function(x) ifna.prev(x))
		
	EPS =  bt.apply(data, function(x) ifna.prev(x[, 'EPS']))
	
	PE = ifna(prices / EPS, NA)
		PE[ abs(EPS) < 0.001 ] = NA	

Please note that for very small EPS, the P/E ratio will be very big; therefore, I set P/E to NA in such cases.

The hard part is done, not let’ plot P/E ratios for all companies.

    #*****************************************************************
    # Create Report
    #******************************************************************       
    plota.matplot(PE)

    plota.matplot(PE, type='b',pch=20, dates='2012::')
		
    plota.matplot(EPS)
	
    plota.matplot(prices)

P/E ratios for all companies:

P/E ratios for all companies in 2012:

Earnings per share (EPS) for all companies:

Prices for all companies:

From these charts I would say it is too early to decide if Facebook is overvalued based on historical P/E ratio basis only, because we only have 3 financial statements, not enough to make an informed conclusion. You might use project one year (FY1) and two year (FY2) earnings estimates to make a better decision.

What is interesting in these charts is how LinkedIn is managing to sustain its astronomical P/E ratio?

I have previously shown examples of how to get and use fundamental data. Here are links for your reference:

To view the complete source code for this example, please have a look at the fundamental.fb.test() function in fundamental.test.r at github.

Categories: Factors, R

More on Factor Attribution to improve performance of the 1-Month Reversal Strategy

In my last post, Factor Attribution to improve performance of the 1-Month Reversal Strategy, I discussed how Factor Attribution can be used to boost performance of the 1-Month Reversal Strategy. Today I want to dig a little dipper and examine this strategy for each sector and also run a sector-neutral back-test.

The initial steps to load historical prices, load Fama/French factors and run Factor Attribution are the same as in the original post, so I will not repeat them here. The complete source code for this example is located in the bt.fa.sector.one.month.test() function in bt.test.r at github for your reference. The new functionality to create sector and sector-neutral back-tests is located in the bt.make.quintiles.sector() function below:

bt.make.quintiles.sector <- function(
	sector,		# sector data
	position.score,	# position.score is a factor to form Quintiles sampled at the period.ends
	data,		# back-test object
	period.ends,	
	n.quantiles = 5,
	start.t = 2,	# first index at which to form Quintiles
	prefix = ''	
) 
{
	#*****************************************************************
	# Re-organize sectors into matrix, assume that sectors are constant in time
	#****************************************************************** 
	temp = factor(sector)
	sector.names = levels(temp)	
		n.sectors = len(sector.names)
	sectors = matrix(unclass(temp), nr=nrow(position.score), nc=ncol(position.score), byrow=T)
	
	#*****************************************************************
	# Create Quantiles
	#****************************************************************** 
	position.score = coredata(position.score)
	quantiles = weights = position.score * NA			
	
	for( s in 1:n.sectors) {
		for( t in start.t:nrow(weights) ) {
			index = sectors[t,] == s
			n = sum(index)
			
			# require at least 3 companies in each quantile
			if(n > 3*n.quantiles) {			
				factor = as.vector(position.score[t, index])
				ranking = ceiling(n.quantiles * rank(factor, na.last = 'keep','first') / count(factor))
			
				quantiles[t, index] = ranking
				weights[t, index] = 1/tapply(rep(1,n), ranking, sum)[ranking]			
			}
		}
	}
	quantiles = ifna(quantiles,0)
	
	#*****************************************************************
	# Create Q1-QN spread for each Sector
	#****************************************************************** 
	long = weights * NA
	short = weights * NA
	models = list()
	
	for( s in 1:n.sectors) {
		long[] = 0
		long[quantiles == 1 & sectors == s] = weights[quantiles == 1 & sectors == s]
		long = long / rowSums(long,na.rm=T)
		
		short[] = 0
		short[quantiles == n.quantiles & sectors == s] = weights[quantiles == n.quantiles & sectors == s]
		short = short / rowSums(short,na.rm=T)
		
		data$weight[] = NA
			data$weight[period.ends,] = long - short
		models[[ paste(prefix,'spread.',sector.names[s], sep='') ]]	= bt.run(data, silent = T)	
	}

	#*****************************************************************
	# Create Sector - Neutral Q1-QN spread
	#****************************************************************** 		
	long[] = 0
	long[quantiles == 1] = weights[quantiles == 1]
	long = long / rowSums(long,na.rm=T)
	
	short[] = 0
	short[quantiles == n.quantiles] = weights[quantiles == n.quantiles]
	short = short / rowSums(short,na.rm=T)
		
	data$weight[] = NA
		data$weight[period.ends,] = long - short
	models$spread.sn = bt.run(data, silent = T)	

	return(models)
}

I have compared below, the overall and sector-neutral versions of the strategy based on the 1 Month returns (one.month) and the Residuals from the Factor Attribution regression (last.e_s). In all cases the Residual strategy outperformed the base one and also the sector-neutral versions have better risk-adjusted coefficients compared to the overall strategy.

Next I looked at the each sector performance for both strategies and surprisingly found that Energy, Materials, and Utilities underperformed in both cases and Financials and Consumer Staples did very well in both cases. Looking at the sector back-test charts, I think the Momentum strategy that selects top few sectors each month would do very well. You can try implementing this idea as a homework.

To view the complete source code for this example, please have a look at the bt.fa.sector.one.month.test() function in bt.test.r at github.

Factor Attribution to improve performance of the 1-Month Reversal Strategy

Today I want to show how to use Factor Attribution to boost performance of the 1-Month Reversal Strategy. The Short-Term Residual Reversal by D. Blitz, J. Huij, S. Lansdorp, M. Verbeek (2011) paper presents the idea and discusses the results as applied to US stock market since 1929. To improve 1-Month Reversal Strategy performance authors investigate an alternative position ranking metric based on the residuals from the rolling Factor Attribution regression.

The base 1-Month Reversal Strategy forms portfolios each month by buying 20% of loosers and short selling 20% of winners from the S&P 500 index based on the prior 1-Month returns. The 1-Month Residual Reversal Strategy forms portfolios each month by buying 20% of loosers and short selling 20% of winners from the S&P 500 index based on the residuals from the rolling Factor Attribution regression. Following are the steps to form 1-Month Residual Reversal Strategy portfolio each month:

  • 1. for each company in the S&P 500 index, run a rolling Factor Attribution regression on the prior 36 months and compute residuals: e.1, e.2, …, e.t, …, e.T for t in 1:36
  • 2. alternative position ranking metric = e.T / standard deviation of (e)
  • 3. form portfolios by buying 20% of loosers and short selling 20% of winners from the S&P 500 index based on the alternative position ranking metric

Let’s start by loading historical prices for all companies in the S&P 500 and create SPY and Equal Weight benchmarks using the Systematic Investor Toolbox:

###############################################################################
# Load Systematic Investor Toolbox (SIT)
# http://systematicinvestor.wordpress.com/systematic-investor-toolbox/
###############################################################################
setInternet2(TRUE)
con = gzcon(url('http://www.systematicportfolio.com/sit.gz', 'rb'))
    source(con)
close(con)

	#*****************************************************************
	# Load historical data
	#****************************************************************** 
	load.packages('quantmod')	
	tickers = sp500.components()$tickers
	
	data <- new.env()
	getSymbols(tickers, src = 'yahoo', from = '1970-01-01', env = data, auto.assign = T)
		# remove companies with less than 5 years of data
		rm.index = which( sapply(ls(data), function(x) nrow(data[[x]])) < 1000 )	
		rm(list=names(rm.index), envir=data)
		
		for(i in ls(data)) data[[i]] = adjustOHLC(data[[i]], use.Adjusted=T)		
	bt.prep(data, align='keep.all', dates='1994::')
		tickers = data$symbolnames
	
	
	data.spy <- new.env()
	getSymbols('SPY', src = 'yahoo', from = '1970-01-01', env = data.spy, auto.assign = T)
	bt.prep(data.spy, align='keep.all', dates='1994::')
	
	#*****************************************************************
	# Code Strategies
	#****************************************************************** 
	prices = data$prices
		n = ncol(prices)
					
	#*****************************************************************
	# Setup monthly periods
	#****************************************************************** 
	periodicity = 'months'
	
	period.ends = endpoints(data$prices, periodicity)
		period.ends = period.ends[period.ends > 0]
	
	prices = prices[period.ends, ]		
		
	#*****************************************************************
	# Create Benchmarks, omit results for the first 36 months - to be consistent with Factor Attribution
	#****************************************************************** 	
	models = list()
	
	# SPY
	data.spy$weight[] = NA
		data.spy$weight[] = 1
		data.spy$weight[1:period.ends[36],] = NA
	models$spy = bt.run(data.spy)
	
	# Equal Weight
	data$weight[] = NA
		data$weight[period.ends,] = ntop(prices, n)
		data$weight[1:period.ends[36],] = NA		
	models$equal.weight = bt.run(data)

Next let’s run Factor Attribution on each the stocks in the S&P 500 index to determine it’s alternative position ranking metric. I will save both e.T and e.T / standard deviation of (e) metrics.

# function to compute additional custom stats for factor.rolling.regression
factor.rolling.regression.custom.stats <- function(x,y,fit) {
	n = len(y)
	e = y - x %*% fit$coefficients
	se = sd(e)
	return(c(e[n], e[n]/se))
}

	#*****************************************************************
	# Load factors and align them with prices
	#****************************************************************** 	
	# load Fama/French factors
	factors = get.fama.french.data('F-F_Research_Data_Factors', periodicity = periodicity,download = T, clean = F)
	
	# align monthly dates
	map = match(format(index(factors$data), '%Y%m'), format(index(prices), '%Y%m'))
		dates = index(factors$data)
		dates[!is.na(map)] = index(prices)[na.omit(map)]
	index(factors$data) = as.Date(dates)
		
	
	# add factors and align
	data.fa <- new.env()
		for(i in tickers) data.fa[[i]] = data[[i]][period.ends, ]
		data.fa$factors = factors$data / 100
	bt.prep(data.fa, align='remove.na')

	
	index = match( index(data.fa$prices), index(data$prices) )
		prices = data$prices[index, ]
		
	#*****************************************************************
	# Compute Factor Attribution for each ticker
	#****************************************************************** 	

	temp = NA * prices
	factors	= list()
		factors$last.e = temp
		factors$last.e_s = temp
	
	for(i in tickers) {
		cat(i, '\n')
		
		# Facto Loadings Regression
		obj = factor.rolling.regression(data.fa, i, 36, silent=T,
			factor.rolling.regression.custom.stats)

		for(j in 1:len(factors))		
			factors[[j]][,i] = obj$fl$custom[,j]
			
	}

	# add base strategy
	factors$one.month = coredata(prices / mlag(prices))

Next let’s group stocks into Quantiles based on 1-Month Reversal factors and create reports.

	#*****************************************************************
	# Create Quantiles
	#****************************************************************** 
	quantiles = list()
	
	for(name in names(factors)) {
		cat(name, '\n')
		quantiles[[name]] = bt.make.quintiles(factors[[name]], data, index, start.t =  1+36, prefix=paste(name,'_',sep=''))
	}	

	#*****************************************************************
	# Create Report
	#****************************************************************** 					
	plotbt.custom.report.part1(quantiles$one.month$spread,quantiles$last.e$spread,quantiles$last.e_s$spread)
	
	plotbt.strategy.sidebyside(quantiles$one.month$spread,quantiles$last.e$spread,quantiles$last.e_s$spread)
		
	plotbt.custom.report.part1(quantiles$last.e_s)

The 1-Month Residual Reversal Strategy have done well over the last 10 years and handsomely outperformed the base 1-Month Reversal Strategy.

To view the complete source code for this example, please have a look at the bt.fa.one.month.test() function in bt.test.r at github.

1-Month Reversal Strategy

Today I want to show a simple example of the 1-Month Reversal Strategy. Each month we will buy 20% of loosers and short sell 20% of winners from the S&P 500 index. The loosers and winners are measured by prior 1-Month returns. I will use this post to set the stage for my next post that will show how Factor Attribution can boost performance of the 1-Month Reversal Strategy. Following is the references for my next post, in case you want to get a flavor, Short-Term Residual Reversal by D. Blitz, J. Huij, S. Lansdorp, M. Verbeek (2011) paper.

Let’s start by loading historical prices for all companies in the S&P 500 and create SPY and Equal Weight benchmarks using the Systematic Investor Toolbox:

###############################################################################
# Load Systematic Investor Toolbox (SIT)
# http://systematicinvestor.wordpress.com/systematic-investor-toolbox/
###############################################################################
setInternet2(TRUE)
con = gzcon(url('http://www.systematicportfolio.com/sit.gz', 'rb'))
    source(con)
close(con)

	#*****************************************************************
	# Load historical data
	#****************************************************************** 
	load.packages('quantmod')	
	tickers = sp500.components()$tickers
	
	data <- new.env()
	getSymbols(tickers, src = 'yahoo', from = '1970-01-01', env = data, auto.assign = T)
		# remove companies with less than 5 years of data
		rm.index = which( sapply(ls(data), function(x) nrow(data[[x]])) < 1000 )	
		rm(list=names(rm.index), envir=data)
		
		for(i in ls(data)) data[[i]] = adjustOHLC(data[[i]], use.Adjusted=T)		
	bt.prep(data, align='keep.all', dates='1994::')
		tickers = data$symbolnames
	
	
	data.spy <- new.env()
	getSymbols('SPY', src = 'yahoo', from = '1970-01-01', env = data.spy, auto.assign = T)
	bt.prep(data.spy, align='keep.all', dates='1994::')
	
	#*****************************************************************
	# Code Strategies
	#****************************************************************** 
	prices = data$prices
		n = ncol(prices)
					
	#*****************************************************************
	# Setup monthly periods
	#****************************************************************** 
	periodicity = 'months'
	
	period.ends = endpoints(data$prices, periodicity)
		period.ends = period.ends[period.ends > 0]
	
	prices = prices[period.ends, ]		
		
	#*****************************************************************
	# Create Benchmarks, omit results for the first 36 months - to be consistent with Factor Attribution
	#****************************************************************** 	
	models = list()
	
	# SPY
	data.spy$weight[] = NA
		data.spy$weight[] = 1
		data.spy$weight[1:period.ends[36],] = NA
	models$spy = bt.run(data.spy)
	
	# Equal Weight
	data$weight[] = NA
		data$weight[period.ends,] = ntop(prices, n)
		data$weight[1:period.ends[36],] = NA		
	models$equal.weight = bt.run(data)

Next let’s group stocks into Quantiles based on 1-Month returns and create back-test for each Quantile. I will rely on the code in the Volatility Quantiles post to create Quantiles.

	#*****************************************************************
	# Create Reversal Quantiles
	#****************************************************************** 
	n.quantiles = 5
	start.t = 1 + 36
	quantiles = weights = coredata(prices) * NA			
	
	one.month = coredata(prices / mlag(prices))
	
	for( t in start.t:nrow(weights) ) {
		factor = as.vector(one.month[t,])
		ranking = ceiling(n.quantiles * rank(factor, na.last = 'keep','first') / count(factor))
		
		quantiles[t,] = ranking
		weights[t,] = 1/tapply(rep(1,n), ranking, sum)[ranking]			
	}
	
	quantiles = ifna(quantiles,0)
	
	#*****************************************************************
	# Create backtest for each Quintile
	#****************************************************************** 
	temp = weights * NA

	for( i in 1:n.quantiles) {
		temp[] = 0
		temp[quantiles == i] = weights[quantiles == i]
	
		data$weight[] = NA
			data$weight[period.ends,] = temp
		models[[ paste('M1_Q',i,sep='') ]] = bt.run(data, silent = T)
	}

Finally, let’s construct Q1/Q5 spread and create summary performance report.

	#*****************************************************************
	# Create Q1-Q5 spread
	#****************************************************************** 
	temp[] = 0
	temp[quantiles == 1] = weights[quantiles == 1]
	temp[quantiles == n.quantiles] = -weights[quantiles == n.quantiles]
	
	data$weight[] = NA
		data$weight[period.ends,] = temp
	models$spread = bt.run(data, silent = T)	

	#*****************************************************************
	# Create Report
	#****************************************************************** 	
	plotbt.custom.report.part1(models)

	plotbt.custom.report.part1(models[spl('spy,equal.weight,spread')])

In the next post I will show how Factor Attribution can boost performance of the 1-Month Reversal Strategy using the methodology presented in the Short-Term Residual Reversal by D. Blitz, J. Huij, S. Lansdorp, M. Verbeek (2011) paper.

To view the complete source code for this example, please have a look at the bt.one.month.test() function in bt.test.r at github.

Example of Factor Attribution

In the prior post, Factor Attribution 2, I have shown how Factor Attribution can be applied to decompose fund’s returns in to Market, Capitalization, and Value factors, the “three-factor model” of Fama and French. Today, I want to show you a different application of Factor Attribution. First, let’s run Factor Attribution on each the stocks in the S&P 500 to determine it’s Value exposure. Next let’s group stocks into Quantiles based on Value exposure and create back-test for each Quantile. I will rely on the code in the Volatility Quantiles post to create Quantiles.

Let’s start by loading historical prices for all current components of the S&P 500 index.

###############################################################################
# Load Systematic Investor Toolbox (SIT)
# http://systematicinvestor.wordpress.com/systematic-investor-toolbox/
###############################################################################
setInternet2(TRUE)
con = gzcon(url('http://www.systematicportfolio.com/sit.gz', 'rb'))
    source(con)
close(con)

	#*****************************************************************
	# Load historical data
	#****************************************************************** 
	load.packages('quantmod')	
	tickers = sp500.components()$tickers
	
	data <- new.env()
	getSymbols(tickers, src = 'yahoo', from = '1970-01-01', env = data, auto.assign = T)
		# remove companies with less than 5 years of data
		rm.index = which( sapply(ls(data), function(x) nrow(data[[x]])) < 1000 )	
		rm(list=names(rm.index), envir=data)
		
		for(i in ls(data)) data[[i]] = adjustOHLC(data[[i]], use.Adjusted=T)		
	bt.prep(data, align='keep.all', dates='1994::')
		tickers = data$symbolnames
	
	
	data.spy <- new.env()
	getSymbols('SPY', src = 'yahoo', from = '1970-01-01', env = data.spy, auto.assign = T)
	bt.prep(data.spy, align='keep.all', dates='1994::')
	
	#*****************************************************************
	# Code Strategies
	#****************************************************************** 
	prices = data$prices
		nperiods = nrow(prices)
		n = ncol(prices)
			
	models = list()
	
	# SPY
	data.spy$weight[] = NA
		data.spy$weight[] = 1
	models$spy = bt.run(data.spy)
	
	# Equal Weight
	data$weight[] = NA
		data$weight[] = ntop(prices, n)
	models$equal.weight = bt.run(data)

Next let’s run Factor Attribution on each the stocks in the S&P 500 to determine it’s Value exposure.

	#*****************************************************************
	# Compute Factor Attribution for each ticker
	#****************************************************************** 
	periodicity = 'weeks'
	
	# load Fama/French factors
	factors = get.fama.french.data('F-F_Research_Data_Factors', periodicity = periodicity,download = T, clean = F)
	
	period.ends = endpoints(data$prices, periodicity)
		period.ends = period.ends[period.ends > 0]
	
	# add factors and align
	data.fa <- new.env()
		for(i in tickers) data.fa[[i]] = data[[i]][period.ends,]
	data.fa$factors = factors$data / 100
	bt.prep(data.fa, align='remove.na')

	
	index = match( index(data.fa$prices), index(data$prices) )
	measure = data$prices[ index, ]	
	for(i in tickers) {
		cat(i, '\n')
		
		# Facto Loadings Regression
		obj = factor.rolling.regression(data.fa, i, 36, silent=T)
		
		measure[,i] = coredata(obj$fl$estimate$HML)
	}

Finally, let’s group stocks into Quantiles based on Value exposure and create back-test for each Quantile.

	#*****************************************************************
	# Create Value Quantiles
	#****************************************************************** 
	n.quantiles=5
	start.t = 1+36
	quantiles = weights = coredata(measure) * NA			
	
	for( t in start.t:nrow(weights) ) {
		factor = as.vector(coredata(measure[t,]))
		ranking = ceiling(n.quantiles * rank(factor, na.last = 'keep','first') / count(factor))
		
		quantiles[t,] = ranking
		weights[t,] = 1/tapply(rep(1,n), ranking, sum)[ranking]			
	}

	quantiles = ifna(quantiles,0)
	
	#*****************************************************************
	# Create backtest for each Quintile
	#****************************************************************** 
	for( i in 1:n.quantiles) {
		temp = weights * NA
			temp[] = 0
		temp[quantiles == i] = weights[quantiles == i]
	
		data$weight[] = NA
			data$weight[index,] = temp
		models[[ paste('Q',i,sep='_') ]] = bt.run(data, silent = T)
	}
	
	#*****************************************************************
	# Create Report
	#****************************************************************** 					
	plotbt.custom.report.part1(models)		
	
	plotbt.strategy.sidebyside(models)

There is no linear relationship between Value Quantiles and historical performance. I’m also suspecting that that implied Value exposure might be quite different than the real Price/Book ratio for each stock. Let me know what do you think about this approach.

In the next post I will show another example of Factor Attribution.

To view the complete source code for this example, please have a look at the bt.fa.value.quantiles.test() function in bt.test.r at github.

Categories: Backtesting, Factors, R

Factor Attribution 2

I want to continue with Factor Attribution theme that I presented in the Factor Attribution post. I have re-organized the code logic into the following 4 functions:

  • factor.rolling.regression – Factor Attribution over given rolling window
  • factor.rolling.regression.detail.plot – detail time-series plot and histogram for each factor
  • factor.rolling.regression.style.plot – historical style plot for selected 2 factors
  • factor.rolling.regression.bt.plot – compare fund’s performance with portfolios implied by Factor Attribution

Let’s first replicate style and performance charts from the Three Factor Rolling Regression Viewer at the mas financial tools web site.

###############################################################################
# Load Systematic Investor Toolbox (SIT)
# http://systematicinvestor.wordpress.com/systematic-investor-toolbox/
###############################################################################
setInternet2(TRUE)
con = gzcon(url('http://www.systematicportfolio.com/sit.gz', 'rb'))
    source(con)
close(con)

	#*****************************************************************
	# Load historical data
	#****************************************************************** 
	load.packages('quantmod')	
	tickers = 'VISVX'

	periodicity = 'months'
		
	data <- new.env()
	getSymbols(tickers, src = 'yahoo', from = '1980-01-01', env = data, auto.assign = T)
	for(i in ls(data)) {
		temp = adjustOHLC(data[[i]], use.Adjusted=T)							
		
		period.ends = endpoints(temp, periodicity)
			period.ends = period.ends[period.ends > 0]

		if(periodicity == 'months') {
			# reformat date to match Fama French Data
			monthly.dates = as.Date(paste(format(index(temp)[period.ends], '%Y%m'),'01',sep=''), '%Y%m%d')
			data[[i]] = make.xts(coredata(temp[period.ends,]), monthly.dates)
		} else
			data[[i]] = temp[period.ends,]
	}
	data.fund = data[[tickers]]
	
	#*****************************************************************
	# Fama/French factors
	#****************************************************************** 
	factors = get.fama.french.data('F-F_Research_Data_Factors', periodicity = periodicity,download = T, clean = F)

	# add factors and align
	data <- new.env()
		data[[tickers]] = data.fund
	data$factors = factors$data / 100
	bt.prep(data, align='remove.na', dates='1994::')

	#*****************************************************************
	# Facto Loadings Regression
	#****************************************************************** 
	obj = factor.rolling.regression(data, tickers, 36)

	#*****************************************************************
	# Reports
	#****************************************************************** 
	factor.rolling.regression.detail.plot(obj)

	factor.rolling.regression.style.plot(obj)

	factor.rolling.regression.bt.plot(obj)

Next let’s add the Momentum factor from the Kenneth R French: Data Library and run Factor Attribution one more time.

	#*****************************************************************
	# Fama/French factors + Momentum
	#****************************************************************** 
	factors = get.fama.french.data('F-F_Research_Data_Factors', periodicity = periodicity,download = T, clean = F)

	factors.extra = get.fama.french.data('F-F_Momentum_Factor', periodicity = periodicity,download = T, clean = F)	
		factors$data = merge(factors$data, factors.extra$data) 
	
	# add factors and align
	data <- new.env()
		data[[tickers]] = data.fund
	data$factors = factors$data / 100
	bt.prep(data, align='remove.na', dates='1994::')

	#*****************************************************************
	# Facto Loadings Regression
	#****************************************************************** 
	obj = factor.rolling.regression(data, tickers, 36)

	#*****************************************************************
	# Reports
	#****************************************************************** 
	factor.rolling.regression.detail.plot(obj)

	factor.rolling.regression.style.plot(obj)

	factor.rolling.regression.bt.plot(obj)


To visualize style from the Momentum point of view, let’s create a style chart that shows fund’s attribution in the HML / Momentum space.

	factor.rolling.regression.style.plot(obj, xfactor='HML', yfactor='Mom')

I designed the Factor Attribution functions to take any user specified factors. This way you can easily run Factor Attribution on any combination of the historical factor returns from the Kenneth R French: Data Library Or use your own historical factor returns data.

To view the complete source code for this example, please have a look at the three.factor.rolling.regression() function in bt.test.r at github.

Factor Attribution

I came across a very descriptive visualization of the Factor Attribution that I will replicate today. There is the Three Factor Rolling Regression Viewer at the mas financial tools web site that performs rolling window Factor Analysis of the “three-factor model” of Fama and French. The factor returns are available from the Kenneth R French: Data Library. I recommend reading the Efficient Frontier: Rolling Your Own: Three-Factor Analysis by W. Bernstein for a step by step instructions.

Let’s start by loading historical returns for the Vanguard Small Cap Value Index (VISVX) and aligning them with Fama/French Monthly Factors. Please note I wrote a helper function, get.fama.french.data(), to simplify process of loading and analyzing factor data from the Kenneth R French: Data Library.

###############################################################################
# Load Systematic Investor Toolbox (SIT)
# http://systematicinvestor.wordpress.com/systematic-investor-toolbox/
###############################################################################
setInternet2(TRUE)
con = gzcon(url('http://www.systematicportfolio.com/sit.gz', 'rb'))
    source(con)
close(con)

	#*****************************************************************
	# Load historical data
	#****************************************************************** 
	load.packages('quantmod')	
	tickers = 'VISVX'
	
	data <- new.env()
	getSymbols(tickers, src = 'yahoo', from = '1980-01-01', env = data, auto.assign = T)
	for(i in ls(data)) {
		temp = adjustOHLC(data[[i]], use.Adjusted=T)							
		
		period.ends = endpoints(temp, 'months')
			period.ends = period.ends[period.ends > 0]

		# reformat date to match Fama French Data
		monthly.dates = as.Date(paste(format(index(temp)[period.ends], '%Y%m'),'01',sep=''), '%Y%m%d')
		data[[i]] = make.xts(coredata(temp[period.ends,]), monthly.dates)
	}
	
	# Fama/French factors
	factors = get.fama.french.data('F-F_Research_Data_Factors', periodicity = 'months',download = T, clean = F)
	
	# add factors and align
	data$factors = factors$data / 100
	bt.prep(data, align='remove.na', dates='1994::')

Next, let’s run Factor Attribution over all available data:

	#*****************************************************************
	# Facto Loadings Regression over whole period
	#****************************************************************** 
	prices = data$prices
		nperiods = nrow(prices)
		dates = index(data$prices)
	
	# compute simple returns	
	hist.returns = ROC(prices[,tickers], type = 'discrete')
		hist.returns = hist.returns - data$factors$RF
		colnames(hist.returns) = 'fund'
	hist.returns = cbind(hist.returns, data$factors$Mkt.RF,
					data$factors$SMB, data$factors$HML)

	fit.all = summary(lm(fund~Mkt.RF+SMB+HML, data=hist.returns))
		estimate.all = c(fit.all$coefficients[,'Estimate'], fit.all$r.squared)
		std.error.all = c(fit.all$coefficients[,'Std. Error'], NA)
Coefficients:
              Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.0006828  0.0012695  -0.538    0.591
Mkt.RF       0.9973980  0.0262881  37.941   <2e-16 ***
SMB          0.5478299  0.0364984  15.010   <2e-16 ***
HML          0.6316528  0.0367979  17.165   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Multiple R-squared: 0.9276,     Adjusted R-squared: 0.9262

All factor loadings are significant.

Next, let’s run a 36 month rolling window Factor Attribution:

	#*****************************************************************
	# Facto Loadings Regression over 36 Month window
	#****************************************************************** 							
	window.len = 36
	
	colnames = spl('alpha,MKT,SMB,HML,R2')
	estimate = make.xts(matrix(NA, nr = nperiods, len(colnames)), dates)
	colnames(estimate) = colnames
	std.error = estimate
	
	# main loop
	for( i in window.len:nperiods ) {
		window.index = (i - window.len + 1) : i
		
		fit = summary(lm(fund~Mkt.RF+SMB+HML, data=hist.returns[window.index,]))
		estimate[i,] = c(fit$coefficients[,'Estimate'], fit$r.squared)
		std.error[i,] = c(fit$coefficients[,'Std. Error'], NA)
		
		if( i %% 10 == 0) cat(i, '\n')
	}

Finally, let’s re-create the timeseries and histogram charts as presented by Three Factor Rolling Regression Viewer.

	#*****************************************************************
	# Reports
	#****************************************************************** 
	layout(matrix(1:10,nc=2,byrow=T))
	
	for(i in 1:5) {	
		#-------------------------------------------------------------------------
		# Time plot
		#-------------------------------------------------------------------------
		est = estimate[,i]
		est.std.error = ifna(std.error[,i], 0)
		
		plota(est, 
			ylim = range( c(
				range(est + est.std.error, na.rm=T),
				range(est - est.std.error, na.rm=T)		
				)))
	
		polygon(c(dates,rev(dates)), 
			c(coredata(est + est.std.error), 
			rev(coredata(est - est.std.error))), 
			border=NA, col=col.add.alpha('red',50))
	
		est = estimate.all[i]
		est.std.error = std.error.all[i]
	
		polygon(c(range(dates),rev(range(dates))), 
			c(rep(est + est.std.error,2),
			rep(est - est.std.error,2)),
			border=NA, col=col.add.alpha('blue',50))
		
		abline(h=0, col='blue', lty='dashed')
			
		abline(h=est, col='blue')
	
		plota.lines(estimate[,i], type='l', col='red')
		
		#-------------------------------------------------------------------------
		# Histogram
		#-------------------------------------------------------------------------
		par(mar = c(4,3,2,1))
		hist(estimate[,i], col='red', border='gray', las=1,
			xlab='', ylab='', main=colnames(estimate)[i])
			abline(v=estimate.all[i], col='blue', lwd=2)
	}

In the next post, I plan to replicate the Style charts and provide more examples.

To view the complete source code for this example, please have a look at the three.factor.rolling.regression() function in bt.test.r at github.

Categories: Factor Model, Factors, R
Follow

Get every new post delivered to your Inbox.

Join 246 other followers