Archive

Archive for March, 2012

Gini Efficient Frontier

March 24, 2012 Leave a comment

David Varadi have recently wrote two posts about Gini Coefficient: I Dream of Gini, and Mean-Gini Optimization. I want to show how to use Gini risk measure to construct efficient frontier and compare it with alternative risk measures I discussed previously.

I will use Gini mean difference risk measure – the mean of the difference between every possible pair of returns to construct Mean-Gini Efficient Frontier. I will use methods presented in “The Generation of Mean Gini Efficient Sets” by J. Okunev (1991) paper to construct optimal portfolios.

Let x.i, i= 1,…,N be weights of instruments in the portfolio. Let us denote by r.it the return of i-th asset in the time period t for i= 1,…,N and t= 1,…,T. The portfolio’s Gini mean difference (page 5) can be written as:

\Gamma = \frac{1}{T^2}\sum_{j=1}^{T}\sum_{k>j}^{T}\left | Y_{j} - Y_{k} \right |

It can be formulated as a linear programming problem

\left | Y_{j} - Y_{k} \right | = a_{jk}^{+} - a_{jk}^{-} \newline \newline min \frac{1}{T^2}\sum_{j=1}^{T}\sum_{k>j}^{T}\left ( a_{jk}^{+} - a_{jk}^{-} \right ) \newline \sum_{i=1}^{N}x_{i}\left ( r_{ij} - r_{ik} \right ) - a_{jk}^{+} + a_{jk}^{-} = 0 \newline \left ( \begin{matrix} for & j=1 & to & T, & and & k>j\end{matrix} \right ) \newline a_{jk}^{+}, a_{jk}^{-}\geqslant 0

This linear programming problem can be easily implemented

min.gini.portfolio <- function
(
	ia,		# input assumptions
	constraints	# constraints
)
{
	n = ia$n
	nt = nrow(ia$hist.returns)
		
	# objective : Gini mean difference - the mean of the difference between every possible pair of returns
	#  1/(T^2) * [ SUM <over j = 1,...,T , k>j> a.long.jk + a.short.jk ]
	f.obj = c(rep(0, n), (1/(nt^2)) * rep(1, nt*(nt-1)))

	# adjust constraints, add a.long.jk , a.short.jk
	constraints = add.variables(nt*(nt-1), constraints, lb=0)

	# [ SUM <over i> x.i * (r.ij - r.ik) ] - a.long.jk + a.short.jk = 0
	# for each j = 1,...,T , k>j	
	a = matrix(0, n + nt*(nt-1), nt*(nt-1)/2)
		diag(a[(n+1) : (n + nt*(nt-1)/2), ]) = -1
		diag(a[(n+1+nt*(nt-1)/2) : (n + nt*(nt-1)), ]) = 1
	hist.returns = as.matrix(ia$hist.returns)
	
	i.start = 0
	for(t in 1:(nt-1)) {
		index = (i.start+1) : (i.start + nt -t)
		for(i in 1:n) {
			a[i, index] = ( hist.returns[t,i] - hist.returns[,i] ) [ (t+1) : nt ] 
		}
		i.start = i.start + nt -t		
	}
	constraints = add.constraints(a, 0, '=', constraints)
				
	# setup linear programming	
	f.con = constraints$A
	f.dir = c(rep('=', constraints$meq), rep('>=', len(constraints$b) - constraints$meq))
	f.rhs = constraints$b

	# find optimal solution
	x = NA
	sol = try(solve.LP.bounds('min', f.obj, t(f.con), f.dir, f.rhs, 
					lb = constraints$lb, ub = constraints$ub), TRUE)

	if(!inherits(sol, 'try-error')) {
		x = sol$solution[1:n]
	}

	return( x )
}

Let’s examine efficient frontiers computed under Gini and Standard deviation risk measures using sample historical input assumptions.

###############################################################################
# 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)

	#--------------------------------------------------------------------------
	# Create Efficient Frontier
	#--------------------------------------------------------------------------
	ia = aa.test.create.ia.rebal()
	n = ia$n		

	# 0 <= x.i <= 1 
	constraints = new.constraints(n, lb = 0, ub = 1)
		
	# SUM x.i = 1
	constraints = add.constraints(rep(1, n), 1, type = '=', constraints)		

	# create efficient frontier(s)
	ef.risk = portopt(ia, constraints, 50, 'Risk')
	ef.gini = portopt(ia, constraints, 50, 'GINI', min.gini.portfolio)
		
	#--------------------------------------------------------------------------
	# Create Plots
	#--------------------------------------------------------------------------
	layout( matrix(1:4, nrow = 2) )
	plot.ef(ia, list(ef.risk, ef.gini), portfolio.risk, F)	
	plot.ef(ia, list(ef.risk, ef.gini), portfolio.gini.coefficient, F)	

	plot.transition.map(ef.risk)
	plot.transition.map(ef.gini)

The Gini efficient frontier is almost identical to Standard deviation efficient frontier, labeled ‘Risk’. This is not a surprise because asset returns that are used in the sample input assumptions are well behaved. The Gini measure of risk would be most appropriate if asset returns contained large outliers.

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

Next I added Gini risk measure to the mix of Asset Allocation strategies that I examined in the Backtesting Asset Allocation portfolios post.

The Gini portfolios and Minimum Variance portfolios show very similar perfromance

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

Advertisement
Categories: Asset Allocation, Backtesting, Portfolio Construction, R, Risk Measures, Strategy, Trading Strategies

Backtesting Asset Allocation portfolios

March 19, 2012 9 comments

In the last post, Portfolio Optimization: Specify constraints with GNU MathProg language, Paolo and MC raised a question: “How would you construct an equal risk contribution portfolio?” Unfortunately, this problem cannot be expressed as a Linear or Quadratic Programming problem.

The outline for this post:

  • I will show how Equal Risk Contribution portfolio can be formulated and solved using a non-linear solver.
  • I will backtest Equal Risk Contribution portfolio and other Asset Allocation portfolios based on various risk measures I described in the Asset Allocation series of post.

Pat Burns wrote an excellent post: Unproxying weight constraints that explains Risk Contribution – partition the variance of a portfolio into pieces attributed to each asset. The Equal Risk Contribution portfolio is a portfolio that splits total portfolio risk equally among its assets. (The concept is similar to 1/N portfolio – a portfolio that splits total portfolio weight equally among its assets.)

Risk Contributions (risk fractions) can be expressed in terms of portfolio weights (w) and covariance matrix (V):
f=\frac{w*Vw}{w'Vw}

Our objective is to find portfolio weights (w) such that Risk Contributions are equal for all assets. This objective function can be easily coded in R:

	risk.contribution = w * (cov %*% w)
	sum( abs(risk.contribution - mean(risk.contribution)) )

I recommend following references for a detailed discussion of Risk Contributions:

I will use a Nonlinear programming solver, Rdonlp2, which is based on donlp2 routine developed and copyright by Prof. Dr. Peter Spellucci to solve for Equal Risk Contribution portfolio. [Please note that following code might not properly execute on your computer because Rdonlp2 package is required and not available on CRAN]

#--------------------------------------------------------------------------
# Equal Risk Contribution portfolio
#--------------------------------------------------------------------------
ia = aa.test.create.ia()
n = ia$n		

# 0 <= x.i <= 1
constraints = new.constraints(n, lb = 0, ub = 1)

# SUM x.i = 1
constraints = add.constraints(rep(1, n), 1, type = '=', constraints)		

# find Equal Risk Contribution portfolio 
w = find.erc.portfolio(ia, constraints)	

# compute Risk Contributions 	
risk.contributions = portfolio.risk.contribution(w, ia)

Next, I want to expand on the Backtesting Minimum Variance portfolios post to include Equal Risk Contribution portfolio and and other Asset Allocation portfolios based on various risk measures I described in the Asset Allocation series of post.

###############################################################################
# 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,quadprog,corpcor,lpSolve')
	tickers = spl('SPY,QQQ,EEM,IWM,EFA,TLT,IYR,GLD')

	data <- new.env()
	getSymbols(tickers, src = 'yahoo', from = '1980-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='1990::2011')
 
	#*****************************************************************
	# Code Strategies
	#****************************************************************** 
	prices = data$prices   
	n = ncol(prices)
	
	# find week ends
	period.ends = endpoints(prices, 'weeks')
		period.ends = period.ends[period.ends > 0]

	#*****************************************************************
	# Create Constraints
	#*****************************************************************
	constraints = new.constraints(n, lb = 0, ub = 1)
	
	# SUM x.i = 1
	constraints = add.constraints(rep(1, n), 1, type = '=', constraints)		

	#*****************************************************************
	# Create Portfolios
	#*****************************************************************			
	ret = prices / mlag(prices) - 1
	start.i = which(period.ends >= (63 + 1))[1]
	
	weight = NA * prices[period.ends,]
	weights = list()
		# Equal Weight 1/N Benchmark
		weights$equal.weight = weight
			weights$equal.weight[] = ntop(prices[period.ends,], n)	
			weights$equal.weight[1:start.i,] = NA
			
		weights$min.var = weight
		weights$min.maxloss = weight
		weights$min.mad = weight
		weights$min.cvar = weight
		weights$min.cdar = weight
		weights$min.cor.insteadof.cov = weight
		weights$min.mad.downside = weight
		weights$min.risk.downside = weight
		
		# following optimizations use a non-linear solver
		weights$erc = weight		
		weights$min.avgcor = weight		
		
	risk.contributions = list()	
		risk.contributions$erc = weight		
		
	# construct portfolios
	for( j in start.i:len(period.ends) ) {
		i = period.ends[j]
		
		# one quarter = 63 days
		hist = ret[ (i- 63 +1):i, ]
		
		# create historical input assumptions
		ia = create.historical.ia(hist, 252)
			s0 = apply(coredata(hist),2,sd)		
			ia$correlation = cor(coredata(hist), use='complete.obs',method='pearson')
			ia$cov = ia$correlation * (s0 %*% t(s0))
			
		# construct portfolios based on various risk measures
		weights$min.var[j,] = min.risk.portfolio(ia, constraints)
		weights$min.maxloss[j,] = min.maxloss.portfolio(ia, constraints)
		weights$min.mad[j,] = min.mad.portfolio(ia, constraints)
		weights$min.cvar[j,] = min.cvar.portfolio(ia, constraints)
		weights$min.cdar[j,] = min.cdar.portfolio(ia, constraints)
		weights$min.cor.insteadof.cov[j,] = min.cor.insteadof.cov.portfolio(ia, constraints)
		weights$min.mad.downside[j,] = min.mad.downside.portfolio(ia, constraints)
		weights$min.risk.downside[j,] = min.risk.downside.portfolio(ia, constraints)

		# following optimizations use a non-linear solver		
		constraints$x0 = weights$erc[(j-1),]
		weights$erc[j,] = find.erc.portfolio(ia, constraints)		
		
		constraints$x0 = weights$min.avgcor[(j-1),]
		weights$min.avgcor[j,] = min.avgcor.portfolio(ia, constraints)						
		
		risk.contributions$erc[j,] = portfolio.risk.contribution(weights$erc[j,], ia)
	}

Next let’s backtest these portfolios and create summary statistics:

	#*****************************************************************
	# Create strategies
	#****************************************************************** 		
	models = list()
	for(i in names(weights)) {
		data$weight[] = NA
			data$weight[period.ends,] = weights[[i]]	
		models[[i]] = bt.run.share(data, clean.signal = F)
	}
		
	#*****************************************************************
	# Create Report
	#****************************************************************** 
	models = rev(models)

	# Plot perfromance
	plotbt(models, plotX = T, log = 'y', LeftMargin = 3)	    	
		mtext('Cumulative Performance', side = 2, line = 1)

	# Plot Strategy Statistics  Side by Side
	plotbt.strategy.sidebyside(models)
	
	# Plot transition maps
	layout(1:len(models))
	for(m in names(models)) {
		plotbt.transition.map(models[[m]]$weight, name=m)
			legend('topright', legend = m, bty = 'n')
	}

	# Plot risk contributions
	layout(1:len(risk.contributions))
	for(m in names(risk.contributions)) {
		plotbt.transition.map(risk.contributions[[m]], name=paste('Risk Contributions',m))
			legend('topright', legend = m, bty = 'n')
	}

	# Compute portfolio concentration and turnover stats based on the
	# On the property of equally-weighted risk contributions portfolios by S. Maillard, 
	# T. Roncalli and J. Teiletche (2008), page 22
	# http://www.thierry-roncalli.com/download/erc.pdf
	out = compute.stats( rev(weights),
		list(Gini=function(w) mean(portfolio.concentration.gini.coefficient(w), na.rm=T),
			Herfindahl=function(w) mean(portfolio.concentration.herfindahl.index(w), na.rm=T),
			Turnover=function(w) 52 * mean(portfolio.turnover(w), na.rm=T)
			)
		)
	
	out[] = plota.format(100 * out, 1, '', '%')
	plot.table(t(out))

The minimum variance (min.risk) portfolio performed very well during that period with 10.5% CAGR and 14% maximum drawdown. The Equal Risk Contribution portfolio (find.erc) also fares well with 10.5% CAGR and 19% maximum drawdown. The 1/N portfolio (equal.weight) is the worst strategy with 7.8% CAGR and 45% maximum drawdown.

One interesting way to modify this strategy is to consider different measures of volatility used to construct a covariance matrix. For example TTR package provides functions for the Garman Klass – Yang Zhang and the Yang Zhang volatility estimation methods. For more details, please have a look at the Different Volatility Measures Effect on Daily MR by Quantum Financier post.

Inspired by the I Dream of Gini by David Varadi, I will show how to create Gini efficient frontier in the next post.

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

Categories: Asset Allocation, Backtesting, Portfolio Construction, R, Risk Measures

Portfolio Optimization: Specify constraints with GNU MathProg language

March 14, 2012 4 comments

I have previously described a few examples of portfolio construction:

I created a number of helper functions to simplify process of making the constraints( i.e. minimum / maximum investment constraints, fully invested constraint – weights must sum to 1, and etc.)

  • new.constraints
  • add.constraints
  • add.variables

However, even with help of these functions, the process of describing the constraints is not simple and user-friendly. Fortunately, there is an alternative way to specify linear constraints using GNU MathProg language. MathProg resembles a subset of AMPL. To find more about GNU MathProg language, I recommend reading following resources:

Let’s start by solving a simple portfolio construction problem using helper functions to specify the constraints.

###############################################################################
# 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 packages
	#****************************************************************** 
	load.packages('quantmod,quadprog,corpcor')
				
	#--------------------------------------------------------------------------
	# Create historical input assumptions
	#--------------------------------------------------------------------------
	tickers = dow.jones.components()
	ia = aa.test.create.ia.custom(tickers, dates = '2000::2010')	
	
	#--------------------------------------------------------------------------
	# Create Constraints & Solve QP problem
	#--------------------------------------------------------------------------
	n = ia$n		

	# 0 <= x.i <= 1 
	constraints = new.constraints(n, lb = 0, ub = 1)

	# SUM x.i = 1
	constraints = add.constraints(rep(1, n), 1, type = '=', constraints)		
	
	# Solve QP problem
	x = min.var.portfolio(ia, constraints)	
	
	# plot weights
	barplot(100*x, las=2, main='Minimum Variance Portfolio')

Now let’s create a GNU MathProg language model that will impose the same constraints. Please copy and save the following model description in the “model1.mod” file:

###############################################################################
set SYMBOLS ;

# set min/max weights for individual stocks
var weight{i in SYMBOLS} >= 0, <= 1 ;
 
# objective function, NOT USED
minimize alpha : sum{i in SYMBOLS} weight[i] ;
 
# weights must sum to 1 (fully invested)
subject to fully_invested : sum{i in SYMBOLS} weight[i] = 1 ;

data;

set SYMBOLS :=  AA AXP BA BAC CAT CSCO CVX DD DIS GE HD HPQ IBM INTC JNJ JPM KFT KO MCD MMM MRK MSFT PFE PG T TRV UTX VZ WMT XOM ;
###############################################################################

Next, let’s use this model to find minimum variance portfolio.

	#*****************************************************************
	# Load packages
	#****************************************************************** 
	# load Rglpk to read GNU MathProg files
	load.packages('Rglpk')

	#--------------------------------------------------------------------------
	# Read GNU MathProg model/Setup constraints/Solve QP problem
	#--------------------------------------------------------------------------	
	model.file = 'model1.mod'

	# read model
	model = Rglpk.read.model(model.file,type = 'MathProg') 	

	# convert GNU MathProg model to constraint used in solve.QP
	constraints = Rglpk.create.constraints(model)$constraints	
		
	# Solve QP problem
	x = min.var.portfolio(ia, constraints)	
	
	# plot weights
	barplot(100*x, las=2, main='Minimum Variance Portfolio using GNU MathProg model')

Next, let’s describe the problem from the Minimum Investment and Number of Assets Portfolio Cardinality Constraints post. Please copy and save the following model description in the “model2.mod” file:

###############################################################################
set SYMBOLS ;

# set min/max weights for individual stocks
var weight{i in SYMBOLS} >= 0, <= 1 ;
 
# add binary, 1 if held, 0 if not held
var held{SYMBOLS} binary;

# objective function, NOT USED
minimize alpha : sum{i in SYMBOLS} weight[i] ;
 
# weights must sum to 1 (fully invested)
subject to fully_invested : sum{i in SYMBOLS} weight[i] = 1 ;

# min weight constraint for individual asset
subject to MinWgt {i in SYMBOLS} : weight[i] >= 0.025 * held[i];
 
# max weight constraint for individual asset
subject to MaxWgt {i in SYMBOLS} : weight[i] <= .20 * held[i] ;

# number of stocks in portfolio
subject to MaxAssetsLB : 0 <= sum {i in SYMBOLS} held[i] ;
subject to MaxAssetsUB : sum {i in SYMBOLS} held[i] <= 6 ;
 
data;

set SYMBOLS :=  AA AXP BA BAC CAT CSCO CVX DD DIS GE HD HPQ IBM INTC JNJ JPM KFT KO MCD MMM MRK MSFT PFE PG T TRV UTX VZ WMT XOM ;
###############################################################################

Next, let’s use this model to find minimum variance portfolio.

	#--------------------------------------------------------------------------
	# Read GNU MathProg model/Setup constraints/Solve QP problem
	#--------------------------------------------------------------------------	
	model.file = 'model2.mod'

	# read model
	model = Rglpk.read.model(model.file,type = 'MathProg') 	

	# convert GNU MathProg model to constraint used in solve.QP
	constraints = Rglpk.create.constraints(model)$constraints	
		
	# Solve QP problem
	x = min.var.portfolio(ia, constraints)	
	
	# plot weights
	barplot(100*x, las=2, main='Minimum Variance Portfolio using GNU MathProg model \n with Minimum Investment and Number of Assets Constraints')

I described another interesting portfolio construction problem in the 130/30 Portfolio Construction post. Please copy and save the following model description in the “model3.mod” file:

###############################################################################
set SYMBOLS ;

# set min/max weights for individual stocks
var long {i in SYMBOLS} >= 0, <= 0.8 ;
var short{i in SYMBOLS} >= 0, <= 0.5 ;
 
# add binary, 1 if long, 0 if short
var islong{SYMBOLS} binary;

# objective function, NOT USED
minimize alpha : sum{i in SYMBOLS} long[i] ;
 
# weights must sum to 1 (fully invested)
subject to fully_invested : sum{i in SYMBOLS} (long[i] - short[i]) = 1 ;

# leverage is 1.6 = longs + shorts
subject to leverage : sum{i in SYMBOLS} (long[i] + short[i]) = 1.6 ;

# force long and short to be mutually exclusive (only one of them is greater then 0 for each i)
subject to long_flag  {i in SYMBOLS} : long[i]  <= islong[i] ;
subject to short_flag {i in SYMBOLS} : short[i] <= (1 - islong[i]) ;
 
data;

set SYMBOLS :=  AA AXP BA BAC CAT CSCO CVX DD DIS GE HD HPQ IBM INTC JNJ JPM KFT KO MCD MMM MRK MSFT PFE PG T TRV UTX VZ WMT XOM ;
###############################################################################

Next, let’s use this model to find minimum variance portfolio.

	#--------------------------------------------------------------------------
	# Read GNU MathProg model/Setup constraints/Solve QP problem
	#--------------------------------------------------------------------------	
	model.file = 'model3.mod'

	# read model
	model = Rglpk.read.model(model.file,type = 'MathProg') 	

	# convert GNU MathProg model to constraint used in solve.QP
	constraints = Rglpk.create.constraints(model)$constraints	
		
	# Solve QP problem, modify Input Assumptions to include short positions
	x = min.var.portfolio(aa.test.ia.add.short(ia), constraints)	
	
	# Compute total weight = longs - short
	x = x[1:ia$n] - x[-c(1:ia$n)]

	# plot weights
	barplot(100*x, las=2, main='Minimum Variance Portfolio using GNU MathProg model \n with 130:30 Constraints')

Another interesting portfolio construction problem is limiting portfolio turnover, or limiting minimum trade size and number of trades. Following model is restricting the trade size to be between 5% and 20% and no more than 8 trades. Please copy and save the following model description in the “model4.mod” file:

###############################################################################
set SYMBOLS ;

param CurWgt{SYMBOLS} ;

# set min/max weights for individual stocks
var weight{i in SYMBOLS} >= 0, <= 1 ;

# TradePos[i] - TradeNeg[i] = CurWgt[i] - weight[i]
var TradePos{i in SYMBOLS} >= 0 ;
var TradeNeg{i in SYMBOLS} >= 0 ;

# Only one of TradePos or TradeNeg is > 0
var TradeFlag{SYMBOLS} binary;

# add binary, 1 if traded, 0 if not traded
var trade{SYMBOLS} binary;

# objective function, NOT USED
minimize alpha : sum{i in SYMBOLS} weight[i] ;
 
# weights must sum to 1 (fully invested)
subject to fully_invested : sum{i in SYMBOLS} weight[i] = 1 ;

# setup Trades for individual asset
subject to TradeRange {i in SYMBOLS} : (CurWgt[i] - weight[i]) = (TradePos[i] - TradeNeg[i]) ;

# Only one of TradePos or TradeNeg is > 0
subject to TradeFlagPos {i in SYMBOLS} : TradePos[i] <= 100 * TradeFlag[i];
subject to TradeFlagNeg {i in SYMBOLS} : TradeNeg[i] <= 100 * (1 - TradeFlag[i]);

# min trade size constraint for individual asset
subject to MinTradeSize {i in SYMBOLS} : (TradePos[i] + TradeNeg[i]) >= 0.05 * trade[i];
subject to MaxTradeSize {i in SYMBOLS} : (TradePos[i] + TradeNeg[i]) <= .20 * trade[i] ; 

# number of trades in portfolio
subject to MaxTrade : sum {i in SYMBOLS} trade[i] <= 8 ;
 
data;

set SYMBOLS :=  AA AXP BA BAC CAT CSCO CVX DD DIS GE ;

param : CurWgt:=
AA	0.1
AXP	0.1
BA	0.1
BAC	0.1
CAT	0.1
CSCO	0.1
CVX	0.1
DD	0.1
DIS	0.1
GE	0.1
; 
###############################################################################

Next, let’s use this model to find minimum variance portfolio.

	#--------------------------------------------------------------------------
	# Read GNU MathProg model/Setup constraints/Solve QP problem
	#--------------------------------------------------------------------------	
	model.file = 'model4.mod'

	# reduce problem size
	ia = aa.test.create.ia.custom(tickers[1:10], dates = '2000::2010')


	# read model
	model = Rglpk.read.model(model.file,type = 'MathProg') 	

	# convert GNU MathProg model to constraint used in solve.QP
	constraints = Rglpk.create.constraints(model)$constraints	
		
	# Solve QP problem
	x = min.var.portfolio(ia, constraints)	
	
	# plot weights
	barplot(100*x, las=2, main='Minimum Variance Portfolio using GNU MathProg model \n with Turnover Constraints')

BA and GE are held constant at 10% and other 8 stocks are traded with trade size at least 5% and no more than 20%.

Please let me know what other type of constraints you like to impose during portfolio construction process.

To view the complete source code for this example, please have a look at the portopt.mathprog.test() function in aa.gmpl.r at github.

Categories: Portfolio Construction, R

Multiple Factor Model Summary

March 11, 2012 6 comments

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 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:

  1. 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.
  2. 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.
  3. 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.
  4. 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).
  5. 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.
  6. 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.
  7. 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.

Categories: Factor Model, Factors

Multiple Factor Model – Building 130/30 Index (Update)

March 10, 2012 2 comments

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.

Categories: Backtesting, Factor Model, Factors, Portfolio Construction, Strategy, Trading Strategies