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.
Multiple Factor Model Summary
In this post I want to summarize all the material I covered in the Multiple Factor Models series. The Multiple Factor Model can be used to decompose returns and calculate risk. Following are some examples of the Multiple Factor Models:
- The expected returns factor model: Commonality In The Determinants Of Expected Stock Returns by R. Haugen, N. Baker (1996)
- The expected returns factor model: CSFB Quantitative Research, Alpha Factor Framework on page 11, page 49 by P. N. Patel, S. Yao, R. Carlson, A. Banerji, J. Handelman
- The risk factor model: MSCI Barra United States Equity Multi-Factor Model, page 101
The factors in the model are usually created using pricing, fundamental, analyst estimates, and proprietary data. I will only show examples of factors using pricing and fundamental data because these infromation is readily available from Yahoo Fiance and ADVFN.
Following is a summary of all posts that I wrote about Multiple Factor Models:
- Multiple Factor Model – Fundamental Data – in this post I demonstrate how to get company’s Fundamental Data into R, create a simple factor, and run correlation analysis.
- Multiple Factor Model – Building Fundamental Factors – in this post I demonstrate how to build Fundamental factors described in the CSFB Alpha Factor Framework and compute quantiles spreads. For details of the CSFB Alpha Factor Framework please read CSFB Quantitative Research, Alpha Factor Framework on page 11, page 49 by P. N. Patel, S. Yao, R. Carlson, A. Banerji, J. Handelman.
- Multiple Factor Model – Building CSFB Factors – in this post I demonstrate how to build majority of factors described in the CSFB Alpha Factor Framework, run cross sectional regression to estimate factor loading, create and test Alpha model.
- Multiple Factor Model – Building Risk Model – in this post I demonstrate how to build a multiple factor risk model, compute factor covariance using shrinkage estimator, forecast stocks specific variances using GARCH(1,1).
- Portfolio Optimization – Why do we need a Risk Model – in this post I explain why do we need a risk model and demonstrate how it is used during portfolio construction process.
- Multiple Factor Model – Building 130/30 Index – in this post I demonstrate how to build 130/30 Index based on the CSFB Factors and the Risk Model we created previously. The 130/30: The New Long-Only (2008) by A. Lo, P. Patel paper presents a very detailed step by step guide to building 130/30 Index using average CSFB Factors as the alpha model and MSCI Barra Multi-Factor Risk model.
- Multiple Factor Model – Building 130/30 Index (Update) – in this post I demonstrate how to build Market-Neutral and Minimum Variance strategies and compare their performance to the 130/30 Index.
There is an excellent discussion of portfolio construction problems and possible solutions in the The top 7 portfolio optimization problems post by Pat Burns. I want to highlight two problems that are relevant to the Multiple Factor Models.
- Problem 3: The expected returns provided by Alpha model sometimes need to scaled or converted to be used in optimization. The Converting Scores into Alphas – A Barra Aegis Case Study (2010) by I. Gleiser, D McKenna paper provides a step by step guide.
- Problem 7: Factor alignment problem is present when different factors are used in the Alpha model and Risk model. An Empirical Case Study of Factor Alignment Problems using the United States Expected Returns (USER) Model – Axioma Research Paper (2011) by A. Saxena, R. Stubbs paper investigates this problem and proposes an alternative portfolio construction methodology.
Multiple Factor Model – Building 130/30 Index (Update)
This is just a quick update to my last post: Multiple Factor Model – Building 130/30 Index. I want to introduce Market-Neutral and Minimum Variance strategies and compare their performance to the 130/30 Index.
I have updated the source code of the fm.long.short.test() function in factor.model.test.r at github to include Market-Neutral and Minimum Variance strategies.
Let’s start by creating a Minimum Variance strategy. I will use the same framework as in the Multiple Factor Model – Building 130/30 Index post. At each month end, I will solve for a portfolio that has minimum variance, ignoring all alpha information.
#*****************************************************************
# Construct LONG ONLY minimum variance portfolio using the multiple factor risk model
#******************************************************************
weights$long.min.var.alpha = weight
for(t in 36:nperiods) {
#--------------------------------------------------------------------------
# Create constraints
#--------------------------------------------------------------------------
# set min/max wgts for individual stocks: 0 =< x <= 10/100
constraints = new.constraints(n, lb = 0, ub = 10/100)
# wgts must sum to 1 (fully invested)
constraints = add.constraints(rep(1,n), type = '=', b = 1, constraints)
#--------------------------------------------------------------------------
# Create factor exposures constraints
#--------------------------------------------------------------------------
# adjust prior constraints, add factor exposures
constraints = add.variables(nfactors, constraints)
# BX - X1 = 0
constraints = add.constraints(rbind(ifna(factor.exposures[t,,], 0), -diag(nfactors)), rep(0, nfactors), type = '=', constraints)
#--------------------------------------------------------------------------
# Create Covariance matrix
# [Qu 0]
# [ 0 Qf]
#--------------------------------------------------------------------------
temp = diag(n)
diag(temp) = ifna(specific.variance[t,], mean(coredata(specific.variance[t,]), na.rm=T))^2
cov.temp = diag(n + nfactors)
cov.temp[1:n,1:n] = temp
cov.temp[(n+1):(n+nfactors),(n+1):(n+nfactors)] = factor.covariance[t,,]
#--------------------------------------------------------------------------
# Setup optimizations
#--------------------------------------------------------------------------
# set expected return
alpha = factors.avg$AVG[t,] / 5
expected.return = c(ifna(coredata(alpha),0), rep(0, nfactors))
# remove companies that have no beta from optimization
index = which(is.na(beta[t,]))
if( len(index) > 0) {
constraints$ub[index] = 0
constraints$lb[index] = 0
}
# find solution
sol = solve.QP.bounds(Dmat = cov.temp, dvec = 0 * expected.return,
Amat = constraints$A, bvec = constraints$b,
meq = constraints$meq, lb = constraints$lb, ub = constraints$ub)
weights$long.min.var.alpha[t,] = sol$solution[1:n]
cat(t, '\n')
}
Next, let’s construct Market-Neutral portfolio. I will restrict portfolio weights to be +/- 10% and will use portfolio construction technique that I documented in the 130/30 Portfolio Construction post.
#*****************************************************************
# Construct Market-Neutral portfolio 100:100 with beta=0 using the multiple factor risk model
# based on the examples in the aa.long.short.test functions
#******************************************************************
weights$market.neutral.alpha = weight
for(t in 36:nperiods) {
#--------------------------------------------------------------------------
# Split x into x.long and x.short, x_long and x_short >= 0
# SUM(x.long) - SUM(x.short) = 0
#--------------------------------------------------------------------------
# x.long and x.short >= 0
# x.long <= 0.1
# x.short <= 0.1
constraints = new.constraints(2*n, lb = 0, ub = c(rep(0.1,n), rep(0.1,n)))
# SUM (x.long - x.short) = 0
constraints = add.constraints(c(rep(1,n), -rep(1,n)), 0, type = '=', constraints)
# SUM (x.long + x.short) = 2
constraints = add.constraints(c(rep(1,n), rep(1,n)), 2, type = '=', constraints)
#--------------------------------------------------------------------------
# beta of portfolio is the weighted average of the individual asset betas
# http://www.duke.edu/~charvey/Classes/ba350/riskman/riskman.htm
#--------------------------------------------------------------------------
temp = ifna(as.vector(beta[t,]),0)
constraints = add.constraints(c(temp, -temp), type = '=', b = 0, constraints)
#--------------------------------------------------------------------------
# Create factor exposures constraints
#--------------------------------------------------------------------------
# adjust prior constraints, add factor exposures
constraints = add.variables(nfactors, constraints)
# BX - X1 = 0
temp = ifna(factor.exposures[t,,], 0)
constraints = add.constraints(rbind(temp, -temp, -diag(nfactors)), rep(0, nfactors), type = '=', constraints)
#--------------------------------------------------------------------------
# Add binary constraints
#--------------------------------------------------------------------------
# adjust prior constraints: add b.i
constraints = add.variables(n, constraints)
# index of binary variables b.i
constraints$binary.index = (2*n + nfactors + 1):(3*n + nfactors)
# binary variable b.i : x.long < b, x.short < (1 - b)
# x.long < b
constraints = add.constraints(rbind(diag(n), 0*diag(n), matrix(0,nfactors,n), -diag(n)), rep(0, n), type = '<=', constraints)
# x.short < (1 - b)
constraints = add.constraints(rbind(0*diag(n), diag(n), matrix(0,nfactors,n), diag(n)), rep(1, n), type = '<=', constraints)
#--------------------------------------------------------------------------
# set expected return
#--------------------------------------------------------------------------
# set expected return
alpha = factors.avg$AVG[t,] / 5
temp = ifna(coredata(alpha),0)
expected.return = c(temp, -temp, rep(0, nfactors), rep(0, n))
#--------------------------------------------------------------------------
# Create Covariance matrix
# [Qu 0]
# [ 0 Qf]
#--------------------------------------------------------------------------
temp = diag(n)
diag(temp) = ifna(specific.variance[t,], mean(coredata(specific.variance[t,]), na.rm=T))^2
# | cov -cov |
# |-cov cov |
temp = cbind( rbind(temp, -temp), rbind(-temp, temp) )
cov.temp = 0*diag(2*n + nfactors + n)
cov.temp[1:(2*n),1:(2*n)] = temp
cov.temp[(2*n+1):(2*n+nfactors),(2*n+1):(2*n+nfactors)] = factor.covariance[t,,]
#--------------------------------------------------------------------------
# Adjust Covariance matrix
#--------------------------------------------------------------------------
if(!is.positive.definite(cov.temp)) {
cov.temp <- make.positive.definite(cov.temp, 0.000000001)
}
#--------------------------------------------------------------------------
# page 9, Risk: We use the Barra default setting, risk aversion value of 0.0075, and
# AS-CF risk aversion ratio of 1.
#
# The Effects of Risk Aversion on Optimization (2010) by S. Liu, R. Xu
# page 4/5
#--------------------------------------------------------------------------
risk.aversion = 0.0075
# remove companies that have no beta from optimization
index = which(is.na(beta[t,]))
if( len(index) > 0) {
constraints$ub[index] = 0
constraints$lb[index] = 0
constraints$ub[2*index] = 0
constraints$lb[2*index] = 0
}
# find solution
sol = solve.QP.bounds(Dmat = 2* risk.aversion * cov.temp, dvec = expected.return,
Amat = constraints$A, bvec = constraints$b,
meq = constraints$meq, lb = constraints$lb, ub = constraints$ub,
binary.vec = constraints$binary.index)
x = sol$solution[1:n] - sol$solution[(n+1):(2*n)]
weights$market.neutral.alpha[t,] = x
cat(t, '\n')
}
Next, let’s re-create all summary charts from the Multiple Factor Model – Building 130/30 Index post.
Please note that we are not comparing apples to apples here. For example, Market-Neutral strategy has beta set equal to zero while 130/30 Index has beta set equal to one. But at the same time there are a few interesting observations:
- Minimum Variance strategy has the smallest drawdown among all long strategies, but its performance is lacking
- Market-Neutral strategy does surprisingly well in this environment and has the best Sharpe ratio
- Market-Neutral strategy also has the highest portfolio turnover; hence, high trading costs will make this strategy less attractive
To view the complete source code for this example, please have a look at the fm.long.short.test() function in factor.model.test.r at github.
Multiple Factor Model – Building 130/30 Index
Nico brought to my attention the 130/30: The New Long-Only (2008) by A. Lo, P. Patel paper in his comment to the Multiple Factor Model – Building CSFB Factors post. This paper presents a very detailed step by step guide to building 130/30 Index using average CSFB Factors as the alpha model and MSCI Barra Multi-Factor Risk model. Today, I want to adapt this methodology and to show how to build 130/30 Index based on the CSFB Factors we created in the Multiple Factor Model – Building CSFB Factors post and the Risk Model we created in the Multiple Factor Model – Building Risk Model post.
Let’s start by loading the CSFB factors that we saved at the end of the Multiple Factor Model – Building CSFB Factors post. [If you are missing data.factors.Rdata file, please execute fm.all.factor.test() function first to create and save CSFB factors.] Next, let’s load the multiple factor risk model we saved at the end of the Multiple Factor Model – Building Risk Model post. [If you are missing risk.model.Rdata file, please execute fm.risk.model.test() function first to create and save multiple factor risk model.] Next, I will compute betas over a two-year rolling window.
###############################################################################
# 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 data
#******************************************************************
load.packages('quantmod')
# Load CSFB factor data that we saved at the end of the fm.all.factor.test function
load(file='data.factors.Rdata')
nperiods = nrow(next.month.ret)
tickers = colnames(next.month.ret)
n = len(tickers)
# Load multiple factor risk model data that we saved at the end of the fm.risk.model.test function
load(file='risk.model.Rdata')
factor.exposures = all.data[,,-1]
factor.names = dimnames(factor.exposures)[[3]]
nfactors = len(factor.names)
#*****************************************************************
# Compute Betas: b = cov(r,m) / var(m)
# The betas are measured on a two-year rolling window
# http://en.wikipedia.org/wiki/Beta_(finance)
#******************************************************************
ret = mlag(next.month.ret)
beta = ret * NA
# 1/n benchmark portfolio
benchmark = ntop(ret, n)
benchmark.ret = rowSums(benchmark * ret, na.rm=T)
# estimate betas
for(t in 24:nperiods) {
t.index = (t-23):t
benchmark.var = var( benchmark.ret[t.index], na.rm=T )
t.count = count(ret[t.index, ])
t.cov = cov( ifna(ret[t.index,], 0), benchmark.ret[t.index], use='complete.obs' )
# require at least 20 months of history
beta[t,] = iif(t.count > 20, t.cov/benchmark.var, NA)
}
To construct efficient portfolio each month, I will solve mean-variance portfolio optimization problem of the following form:
Max weight * return - risk.aversion * weight * Covariance * weight Sum weight = 1 Sum weight * beta = 1 0 <= weight <= 0.1
Please read The Effects of Risk Aversion on Optimization (2010) by S. Liu, R. Xu paper for the detailed discussion of this optimization problem.
Please note that I will use risk.aversion = 0.0075 and portfolio beta = 1 as discussed on the page 22 of the 130/30: The New Long-Only (2008) by A. Lo, P. Patel paper. To model portfolio beta constraint, I will use the fact that portfolio beta is equal to the weighted average of the individual asset betas.
#*****************************************************************
# Construct LONG ONLY portfolio using the multiple factor risk model
#******************************************************************
load.packages('quadprog,corpcor,kernlab')
weight = NA * next.month.ret
weights = list()
weights$benchmark = ntop(beta, n)
weights$long.alpha = weight
for(t in 36:nperiods) {
#--------------------------------------------------------------------------
# Create constraints
#--------------------------------------------------------------------------
# set min/max wgts for individual stocks: 0 =< x <= 10/100
constraints = new.constraints(n, lb = 0, ub = 10/100)
# wgts must sum to 1 (fully invested)
constraints = add.constraints(rep(1,n), type = '=', b = 1, constraints)
#--------------------------------------------------------------------------
# beta of portfolio is the weighted average of the individual asset betas
# http://www.duke.edu/~charvey/Classes/ba350/riskman/riskman.htm
#--------------------------------------------------------------------------
constraints = add.constraints(ifna(as.vector(beta[t,]),0), type = '=', b = 1, constraints)
#--------------------------------------------------------------------------
# Create factor exposures constraints
#--------------------------------------------------------------------------
# adjust prior constraints, add factor exposures
constraints = add.variables(nfactors, constraints)
# BX - X1 = 0
constraints = add.constraints(rbind(ifna(factor.exposures[t,,], 0), -diag(nfactors)), rep(0, nfactors), type = '=', constraints)
#--------------------------------------------------------------------------
# Create Covariance matrix
# [Qu 0]
# [ 0 Qf]
#--------------------------------------------------------------------------
temp = diag(n)
diag(temp) = ifna(specific.variance[t,], mean(coredata(specific.variance[t,]), na.rm=T))^2
cov.temp = diag(n + nfactors)
cov.temp[1:n,1:n] = temp
cov.temp[(n+1):(n+nfactors),(n+1):(n+nfactors)] = factor.covariance[t,,]
#--------------------------------------------------------------------------
# create input assumptions
#--------------------------------------------------------------------------
ia = list()
ia$n = nrow(cov.temp)
ia$annual.factor = 12
ia$symbols = c(tickers, factor.names)
ia$cov = cov.temp
#--------------------------------------------------------------------------
# page 9, Risk: We use the Barra default setting, risk aversion value of 0.0075, and
# AS-CF risk aversion ratio of 1.
#
# The Effects of Risk Aversion on Optimization (2010) by S. Liu, R. Xu
# page 4/5
#--------------------------------------------------------------------------
risk.aversion = 0.0075
ia$cov.temp = ia$cov
# set expected return
alpha = factors.avg$AVG[t,] / 5
ia$expected.return = c(ifna(coredata(alpha),0), rep(0, nfactors))
# remove companies that have no beta from optimization
index = which(is.na(beta[t,]))
if( len(index) > 0) {
constraints$ub[index] = 0
constraints$lb[index] = 0
}
# find solution
sol = solve.QP.bounds(Dmat = 2* risk.aversion * ia$cov.temp, dvec = ia$expected.return,
Amat = constraints$A, bvec = constraints$b,
meq = constraints$meq, lb = constraints$lb, ub = constraints$ub)
weights$long.alpha[t,] = sol$solution[1:n]
}
Next, let’s construct 130/30 portfolio. I will restrict portfolio weights to be +/- 10% and will use portfolio construction technique that I documented in the 130/30 Portfolio Construction post.
#*****************************************************************
# Construct Long/Short 130:30 portfolio using the multiple factor risk model
# based on the examples in the aa.long.short.test functions
#******************************************************************
weights$long.short.alpha = weight
for(t in 36:nperiods) {
#--------------------------------------------------------------------------
# Create constraints
#--------------------------------------------------------------------------
# set min/max wgts for individual stocks: -10/100 =< x <= 10/100
constraints = new.constraints(n, lb = -10/100, ub = 10/100)
# wgts must sum to 1 (fully invested)
constraints = add.constraints(rep(1,n), type = '=', b = 1, constraints)
#--------------------------------------------------------------------------
# beta of portfolio is the weighted average of the individual asset betas
# http://www.duke.edu/~charvey/Classes/ba350/riskman/riskman.htm
#--------------------------------------------------------------------------
constraints = add.constraints(ifna(as.vector(beta[t,]),0), type = '=', b = 1, constraints)
#--------------------------------------------------------------------------
# Create factor exposures constraints
#--------------------------------------------------------------------------
# adjust prior constraints, add factor exposures
constraints = add.variables(nfactors, constraints)
# BX - X1 = 0
constraints = add.constraints(rbind(ifna(factor.exposures[t,,], 0), -diag(nfactors)), rep(0, nfactors), type = '=', constraints)
#--------------------------------------------------------------------------
# Create 130:30
# -v.i <= x.i <= v.i, v.i>0, SUM(v.i) = 1.6
#--------------------------------------------------------------------------
# adjust prior constraints, add v.i
constraints = add.variables(n, constraints)
# -v.i <= x.i <= v.i
# x.i + v.i >= 0
constraints = add.constraints(rbind(diag(n), matrix(0,nfactors,n) ,diag(n)), rep(0, n), type = '>=', constraints)
# x.i - v.i <= 0
constraints = add.constraints(rbind(diag(n), matrix(0,nfactors,n), -diag(n)), rep(0, n), type = '<=', constraints)
# SUM(v.i) = 1.6
constraints = add.constraints(c(rep(0, n), rep(0, nfactors), rep(1, n)), 1.6, type = '=', constraints)
#--------------------------------------------------------------------------
# Create Covariance matrix
# [Qu 0]
# [ 0 Qf]
#--------------------------------------------------------------------------
temp = diag(n)
diag(temp) = ifna(specific.variance[t,], mean(coredata(specific.variance[t,]), na.rm=T))^2
cov.temp = 0*diag(n + nfactors + n)
cov.temp[1:n,1:n] = temp
cov.temp[(n+1):(n+nfactors),(n+1):(n+nfactors)] = factor.covariance[t,,]
#--------------------------------------------------------------------------
# create input assumptions
#--------------------------------------------------------------------------
ia = list()
ia$n = nrow(cov.temp)
ia$annual.factor = 12
ia$symbols = c(tickers, factor.names, tickers)
ia$cov = cov.temp
#--------------------------------------------------------------------------
# page 9, Risk: We use the Barra default setting, risk aversion value of 0.0075, and
# AS-CF risk aversion ratio of 1.
#
# The Effects of Risk Aversion on Optimization (2010) by S. Liu, R. Xu
# page 4/5
#--------------------------------------------------------------------------
risk.aversion = 0.0075
ia$cov.temp = ia$cov
# set expected return
alpha = factors.avg$AVG[t,] / 5
ia$expected.return = c(ifna(coredata(alpha),0), rep(0, nfactors), rep(0, n))
# remove companies that have no beta from optimization
index = which(is.na(beta[t,]))
if( len(index) > 0) {
constraints$ub[index] = 0
constraints$lb[index] = 0
}
# find solution
sol = solve.QP.bounds(Dmat = 2* risk.aversion * ia$cov.temp, dvec = ia$expected.return,
Amat = constraints$A, bvec = constraints$b,
meq = constraints$meq, lb = constraints$lb, ub = constraints$ub)
weights$long.short.alpha[t,] = sol$solution[1:n]
}
At this point, we created monthly long-only and 130:30 portfolios. Let’s examine their transition maps.
#***************************************************************** # Plot Transition Maps #****************************************************************** layout(1:3) for(i in names(weights)) plotbt.transition.map(weights[[i]], i)
The portfolio weights for the long-only portfolio (long.alpha) sum up to 100% and for the 130:30 portfolio (long.short.alpha) is 130% long and 30% short as expected.
Next let’s create trading strategies based on these portfolios and check their performance.
#*****************************************************************
# Create strategies
#******************************************************************
prices = data$prices
prices = bt.apply.matrix(prices, function(x) ifna.prev(x))
# find month ends
month.ends = endpoints(prices, 'months')
# create strategies that invest in each qutile
models = list()
for(i in names(weights)) {
data$weight[] = NA
data$weight[month.ends,] = weights[[i]]
capital = 100000
data$weight[] = (capital / prices) * (data$weight)
models[[i]] = bt.run(data, type='share', capital=capital)
}
#*****************************************************************
# Create Report
#******************************************************************
models = rev(models)
plotbt.custom.report.part1(models, dates='1998::')
plotbt.custom.report.part2(models)
The 130:30 portfolio works best, producing high returns with smaller drawdowns than long-only and benchmark (1/N) portfolios.
The note of caution: the above results are based on $0 commission rate (i.e. trading is free). Also, I’m using the current Dow Jones index components through out the whole history; hence, introducing survivorship bias (i.e. Dow Jones index changed its composition a few times in the last 20 years).
To see how much of the problem is my assumption of $0 commission rate, let’s have a look at the annual portfolio turnovers.
# Plot Portfolio Turnover for each strategy layout(1) barplot.with.labels(sapply(models, compute.turnover, data), 'Average Annual Portfolio Turnover')
The 130:30 portfolio has a pretty high portfolio turnover; therefore, it will not perform as well in the real life as in our backtest.
To view the complete source code for this example, please have a look at the fm.long.short.test() function in factor.model.test.r at github.
Systematic Investor 