Archive
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)
# https://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)
# https://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.
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)
# https://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)
# https://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&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.
Systematic Investor 