Archive

Archive for June, 2012

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

Volatility Position Sizing 2

I have discussed Volatility Position Sizing in the Volatility Position Sizing to improve Risk Adjusted Performance post using the Average True Range (ATR) as a measure of Volatility.

Today I want show how to use historical volatility to adjust portfolio leverage. Let’s start with Buy and Hold strategy using SPY and rescale it to the target volatility of 10%. I will use a 60 day rolling historical Volatility to adjust portfolio leverage in order to target 10% annual Volatility.

###############################################################################
# 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 = 'SPY'
	
	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)		
	bt.prep(data, align='keep.all', dates='1994::')
	
	#*****************************************************************
	# Buy and Hold
	#****************************************************************** 
	prices = data$prices
	models = list()
	
	data$weight[] = 1
	models$buy.hold = bt.run.share(data, clean.signal=T)

	#*****************************************************************
	# Buy and Hold with target 10% Volatility
	#****************************************************************** 
	ret.log = bt.apply.matrix(data$prices, ROC, type='continuous')
	hist.vol = sqrt(252) * bt.apply.matrix(ret.log, runSD, n = 60)

	data$weight[] = 0.1 / hist.vol
	models$buy.hold.volatility.weighted = bt.run.share(data, clean.signal=T)

	#*****************************************************************
	# Buy and Hold with target 10% Volatility and Max Total leverage 100%
	#****************************************************************** 		
	data$weight[] = 0.1 / hist.vol
		rs = rowSums(data$weight)
		data$weight[] = data$weight / iif(rs > 1, rs, 1) 			
	models$buy.hold.volatility.weighted.100 = bt.run.share(data, clean.signal=T)
		
	#*****************************************************************
	# Same, rebalanced Monthly
	#****************************************************************** 
	period.ends = endpoints(prices, 'months')
		period.ends = period.ends[period.ends > 0]
		
	data$weight[] = NA
	data$weight[period.ends,] = 0.1 / hist.vol[period.ends,]
		rs = rowSums(data$weight[period.ends,])
		data$weight[period.ends,] = data$weight[period.ends,] / iif(rs > 1, rs, 1) 			
	models$buy.hold.volatility.weighted.100.monthly = bt.run.share(data, clean.signal=T)
		
	#*****************************************************************
	# Create Report
	#****************************************************************** 
	# Plot perfromance
	plotbt(models, plotX = T, log = 'y', LeftMargin = 3)	    	
		mtext('Cumulative Performance', side = 2, line = 1)

	plotbt.custom.report.part2(rev(models))
	
	# Plot Portfolio Turnover for each strategy
	layout(1)
	barplot.with.labels(sapply(models, compute.turnover, data), 'Average Annual Portfolio Turnover', plotX = F, label='both')

All strategies do a good job of rescaling portfolio leverage to the target 10% volatility.

Next let’s investigate other volatility measures. The volatility function in the TTR package includes following volatility calculation types:

  • Historical volatility
  • Garman and Klass
  • Parkinson
  • Rogers and Satchell
  • Garman and Klass – Yang and Zhang
  • Yang and Zhang

We can easily modify the algorithm above to use these different volatility calculations. One of my favorite sides of R is that there is a function or package for anything you can imagine :)

	#*****************************************************************
	# Next let's examine other volatility measures
	#****************************************************************** 
	models = models[c('buy.hold' ,'buy.hold.volatility.weighted.100.monthly')]
		
	# TTR volatility calc types
	calc = c("close", "garman.klass", "parkinson", "rogers.satchell", "gk.yz", "yang.zhang")
	
	ohlc = OHLC(data$SPY)
	for(icalc in calc) {
		vol = volatility(ohlc, calc = icalc, n = 60, N = 252)
		
		data$weight[] = NA
		data$weight[period.ends,] = 0.1 / vol[period.ends,]
			rs = rowSums(data$weight[period.ends,])
			data$weight[period.ends,] = data$weight[period.ends,] / iif(rs > 1, rs, 1) 			
		models[[icalc]] = bt.run.share(data, clean.signal=T)
	}
	
	#*****************************************************************
	# Create Report
	#****************************************************************** 
	# Plot performance
	plotbt(models, plotX = T, log = 'y', LeftMargin = 3)	    	
		mtext('Cumulative Performance', side = 2, line = 1)
	
	plotbt.strategy.sidebyside(models)

The simple historical volatility does the best job at targeting volatility and controlling drawdowns.

Next let’s apply idea of Volatility Position Sizing to the strategy’s Equity Curve.

	#*****************************************************************
	# Volatility Position Sizing applied to MA cross-over strategy's Equity Curve
	#****************************************************************** 
	models = list()	
	
	sma.fast = SMA(prices, 50)
	sma.slow = SMA(prices, 200)
	weight = iif(sma.fast >= sma.slow, 1, -1)
	
	data$weight[] = weight
	models$ma.crossover = bt.run.share(data, clean.signal=T)
		
	#*****************************************************************
	# Target 10% Volatility
	#****************************************************************** 
	ret.log = bt.apply.matrix(models$ma.crossover$equity, ROC, type='continuous')
	hist.vol = sqrt(252) * bt.apply.matrix(ret.log, runSD, n = 60)
		
	data$weight[] = NA
		data$weight[period.ends,] = (0.1 / hist.vol[period.ends,]) * weight[period.ends,]
		# limit total leverage to 100%		
		rs = rowSums(data$weight[period.ends,])
		data$weight[period.ends,] = data$weight[period.ends,] / iif(abs(rs) > 1, abs(rs), 1) 			
	models$ma.crossover.volatility.weighted.100.monthly = bt.run.share(data, clean.signal=T)
	
	#*****************************************************************
	# Create Report
	#****************************************************************** 
	# Plot perfromance
	plotbt(models, plotX = T, log = 'y', LeftMargin = 3)	    	
		mtext('Cumulative Performance', side = 2, line = 1)
	
	plotbt.custom.report.part2(rev(models))

The Volatility Position Sizing does keep strategy’s volatility close to the target 10% volatility.

In the final example, I will apply Volatility Position Sizing to the Timing strategy developed by M. Faber

	#*****************************************************************
	# Apply Volatility Position Sizing Timing stretegy by M. Faber
	#****************************************************************** 
	tickers = spl('SPY,QQQ,EEM,IWM,EFA,TLT,IYR,GLD')
	
	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)		
	bt.prep(data, align='remove.na', dates='1994::')
	
	#*****************************************************************
	# Code Strategies
	#****************************************************************** 
	prices = data$prices   
		n = ncol(prices)
	models = list()

	period.ends = endpoints(prices, 'months')
		period.ends = period.ends[period.ends > 0]
		
	#*****************************************************************
	# Equal Weight
	#****************************************************************** 
	data$weight[] = NA
		data$weight[period.ends,] = ntop(prices[period.ends,], n)
		data$weight[1:200,] = NA
	models$equal.weight = bt.run.share(data, clean.signal=F)
				
	#*****************************************************************
	# Timing by M. Faber
	#****************************************************************** 
	sma = bt.apply.matrix(prices, SMA, 200)
	
	weight = ntop(prices, n) * (prices > sma)
	data$weight[] = NA
		data$weight[period.ends,] = weight[period.ends,]
	models$timing = bt.run.share(data, clean.signal=F)
		
	#*****************************************************************
	# Timing with target 10% Volatility
	#****************************************************************** 
	ret.log = bt.apply.matrix(models$timing$equity, ROC, type='continuous')
	hist.vol = bt.apply.matrix(ret.log, runSD, n = 60)
		hist.vol = sqrt(252) * as.vector(hist.vol)
	
	data$weight[] = NA
		data$weight[period.ends,] = (0.1 / hist.vol[period.ends]) * weight[period.ends,]
		rs = rowSums(data$weight)
		data$weight[] = data$weight / iif(rs > 1, rs, 1) 				
		data$weight[1:200,] = NA
	models$timing.volatility.weighted.100.monthly = bt.run.share(data, clean.signal=T)
	
	#*****************************************************************
	# Create Report
	#****************************************************************** 
	# Plot perfromance
	plotbt(models, plotX = T, log = 'y', LeftMargin = 3)	    	
		mtext('Cumulative Performance', side = 2, line = 1)
	
	plotbt.custom.report.part2(rev(models))

Volatility Position Sizing is one of many Position Sizing algorithms that can be part of your money management rules. Let me what Position Sizing scheme works best for you.

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

Volatility Quantiles

Today I want to examine the performance of stocks in the S&P 500 grouped into Quantiles based on one year historical Volatility. The idea is very simple: each week we will form Volatility Quantiles portfolios by grouping stocks in the S&P 500 into Quantiles using one year historical Volatility. Next we will backtest each portfolio and check if low historical volatility corresponds to the low realized volatility.

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/
###############################################################################
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::')
	
	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, 500)
	models$equal.weight = bt.run(data)	

Next let’s divide stocks in the S&P 500 into Quantiles using one year historical Volatility and create backtest for each quantile.

	#*****************************************************************
	# Create Quantiles based on the historical one year volatility 
	#****************************************************************** 
	# setup re-balancing periods
	period.ends = endpoints(prices, 'weeks')
		period.ends = period.ends[period.ends > 0]

	# compute historical one year volatility
	p = bt.apply.matrix(coredata(prices), ifna.prev)	
	ret = p / mlag(p) - 1		
	sd252 = bt.apply.matrix(ret, runSD, 252)		
		
	# split stocks in the S&amp;P 500 into Quantiles using one year historical Volatility
	n.quantiles=5
	start.t = which(period.ends >= (252+2))[1]	
	quantiles = weights = p * NA
				
	for( t in start.t:len(period.ends) ) {
		i = period.ends[t]

		factor = sd252[i,]
		ranking = ceiling(n.quantiles * rank(factor, na.last = 'keep','first') / count(factor))
	
		quantiles[i,] = ranking
		weights[i,] = 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[period.ends,] = 0
		temp[quantiles == i] = weights[quantiles == i]
	
		data$weight[] = NA
			data$weight[] = temp
		models[[ paste('Q',i,sep='_') ]] = bt.run(data, silent = T)
	}

Finally let’s plot historical equity curves for each Volatility Quantile and create a summary table.

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

Our thesis is true for stocks in the S&P 500 index: low historical volatility leads to low realized volatility. There is a one-to-one correspondence between historical and realized volatility quantiles.

Please note that performance numbers have to be taken with a grain of salt because I used current components in the S&P 500 index; hence, introducing survivorship bias.

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

Follow

Get every new post delivered to your Inbox.

Join 242 other followers