Skip to content

Sharpe's algorithm for portfolio improvement

22 messages · Marc Delvaux, Enrico Schumann, Alexios Ghalanos +5 more

#
An attractive sounding algorithm for maximizing the expected utility of 
of a portfolio was proposed by William F. Sharpe in "An algorithm for 
portfolio improvement," Advances in Mathematical Programming and 
Financial Planning, 1987, pp. 155-170 and summarized by the same author 
in "Expected utility asset allocation," Financial Analysts Journal, vol. 
63, no. 5 (Sep.-Oct., 2007), pp. 18-30.

Has this algorithm been implemented in R?

If not, is there a substitute that is likely to work well for a 
user-specified concave utility function?  I've tried optim, DEoptim, and 
Rgenoud but encountered difficulty in getting them to converge for a 
long-only portfolio of around 30 assets.

Best regards,
John
#
Hi John,

what kind of "difficulty" did you encounter? If you would give more details
on what you tried, and how, then people should be better able to help you. 

I don't know the paper you mentioned, but I know a paper of W. Sharpe in
which he suggests to do repeated zero-sum changes to the portfolio, like
"increase one asset by x%, and decrease another one by x%". If that is what
you mean, this can be done with a local search. (But actually, other
functions like DEoptim should work just as well. The DEoptim package even
comes with a vignette that adresses portfolio optimisation.) 

An example for a local search procedure is given in one of the vignettes of
the NMOF package (which is available from
https://r-forge.r-project.org/R/?group_id=1128 ), even though I am not sure
how self-explanatory the vignette is.


Regards,
Enrico
#
Enrico Schumann wrote:
Thank you, Enrico, for your prompt and helpful response.  Focusing on 
Sharpe's algorithm, I originally omitted specifics about difficulties I 
have encountered with packages based on other algorithms.  However, in 
case it is of interest, I will now outline my project and difficulties. 
     I have about 120 monthly observations on gross real returns for 
about 30 assets. [Trying to mitigate estimation error, I've shrunk 
observed returns toward a common value as proposed by Philippe Jorion, 
"Bayes-Stein Estimation for Portfolio Analysis",
Journal of Financial and Quantitative Analysis, v. 21, n. 3 (Sept., 
1986), pp. 279-92.]   I initially specified the utility function as U = 
R/(0.5 + R) where R is R is the gross return on a portfolio and 
specified long-only constraints. The restriction that portfolio shares 
sum to 1 is handled by specifying n-1 variables (where n is the nubmer 
of assets) and making the share asset n be 1 minus the sum of other 
shares. If I apply DEoptim to a sufficiently small subset of the assets, 
it converges and  selects a plausible portfolio.  Yet if I ask DEoptim 
to analyze as many as 30 assets, it fails to converge in any number of 
iterations that I've tried.  Given the same problem, Rgenoud converged 
after 44 hours (more than 2 million generation, if memory serves).  I 
subsequently tried changing the utility function to ln(R) and asking 
Rgenoud to maximize that.  Thus far it has run for over 1.5 million 
generations without converging.  The portfolio shares have barely 
changed over many recent generations.  Perhaps I could just relax the 
default convergence criteria and declare the problem solved for 
practical purposes.  However, that might result in mistaking a local for 
a global maximum.  These experiences may just indicate that a 30 asset 
portfolio is hard to optimize using general purpose algorithms. Kris 
Boudt et al. note that "portfolio problems encountered in practice may 
require days to optimize on a personal computer" ("Large-scale portfolio 
optimization with DEoptim," p. 1). These experiences motivated my 
interest in trying an algorithm, such as that of Sharpe, designed 
specifically for portfolio optimization.
The algorithm you describe sounds very much like that covered in the 
articles I cited in my previous note.  It's probably the same algorithm.

(But actually, other
Perhaps I should just be more patient with DEoptim or buy a faster computer!
Thank you for the NMOF reference.  I've printed "Examples and Extensions 
for the NMOF package" and tried the command
install.packages("NMOF", repos = "http://R-Forge.R-project.org")
That command elicited the following messages:
Warning in install.packages("NMOF", repos = 
"http://R-Forge.R-project.org") :
   argument 'lib' is missing: using 
'/home/john/R/i486-pc-linux-gnu-library/2.8'
trying URL 'http://R-Forge.R-project.org/src/contrib/NMOF_0.14-2.tar.gz'
Content type 'application/x-gzip' length 1352123 bytes (1.3 Mb)
opened URL
==================================================
downloaded 1.3 Mb

* Installing *source* package 'NMOF' ...
** R
** data
**  moving datasets to lazyload DB
** inst
** preparing package for lazy loading
** help
  >>> Building/Updating help pages for package 'NMOF'
      Formats: text html latex example
   DEopt                             text    html    latex   example
   EuropeanCall                      text    html    latex   example
   GAopt                             text    html    latex   example
   LSopt                             text    html    latex   example
   MA                                text    html    latex   example
   NMOF-package                      text    html    latex   example
   NS                                text    html    latex   example
   NSf                               text    html    latex   example
   PSopt                             text    html    latex   example
   TAopt                             text    html    latex   example
   bundData                          text    html    latex   example
   callHestoncf                      text    html    latex   example
   fundData                          text    html    latex   example
   gridSearch                        text    html    latex   example
   qTable                            text    html    latex   example

ERROR in file 'repairMatrix.Rd': Environment (text enclosed in {}) found 
in \title{...}.
The title must be plain text!

ERROR: building help failed for package 'NMOF'
** Removing '/home/john/R/i486-pc-linux-gnu-library/2.8/NMOF'

The downloaded packages are in
	/tmp/RtmpNAuDvf/downloaded_packages
Warning message:
In install.packages("NMOF", repos = "http://R-Forge.R-project.org") :
   installation of package 'NMOF' had non-zero exit status

***************************************

Any suggestions for how to successfully install NMOF would be greatly 
appreciated.

Best regards,
John

  
    
#
I've found the cmaes (Covariance Matrix Adapting Evolutionary Strategy) 
solver on CRAN to be quite robust to some nonconvex portfolio
problems which required a global optimizer. Like other such global 
solvers which do not accept explicit constraint functions, you would 
need to create a sort of penalty barrier function i.e

-obj(x) + 100*(1-sum(x))^2

For state of the art global solvers you'll probably need to look outside 
R to a commercial implementation or try submitting your problem to the 
NEOS server (http://www.neos-server.org/neos/solvers/index.html...see in 
particular the Branch and Bound type solver BARON).

Regards,

Alexios
On 01/08/2011 22:07, John P. Burkett wrote:
#
On Mon, 2011-08-01 at 17:07 -0400, John P. Burkett wrote:
As one of the authors of the paper in question, I'll note that I was
talking about portfolios with hundreds or thousands of assets, and/or
with many rebalancing periods or observations.  Monthly series of 120
observations each shouldn't take millions of iterations to converge.  It
is possible that the starting set is not feasible, given the way your
full investment constraint is being applied.

If your objective function is truly a smooth convex function, then
optim() or some other linear or conical solver should be sufficient.  If
the function is not smooth (and real portfolio objectives and
constraints rarely are), then a random portfolios or a global optimizer
will be required.  

Perhaps you should start with a couple hundred random portfolio survey
of your feasible space?

This may be accomplished in R with Pat Burns' excellent Portfolio Probe
(commercial non-free), or with PortfolioAnalytics (open source/free, on
R-Forge).

PortfolioAnalytics will also allow you to use DEoptim after using random
portfolios to get close to the global minima.  

See out R/Finance seminar presentation here:
http://www.rinfinance.com/agenda/2010/Carl+Peterson+Boudt_Tutorial.pdf
and the code to replicate here:
https://r-forge.r-project.org/scm/viewvc.php/*checkout*/pkg/PortfolioAnalytics/sandbox/script.workshop2010.R?root=returnanalytics
PortfolioAnalytics does some of the heavy lifting to set good defaults
and give you a good set of starting population values to speed
convergence.

Also, Pat Burns' Portfolio Probe blog is here:
http://www.portfolioprobe.com/blog/

- Brian
#
alexios wrote:
Alexios,
Thank you for your suggestions.  CMAES looks promising for some 
applications. However, for maximizing expected utility or minimizing -1* 
(expected utility), detection of CMAES convergence could be tricky.  The 
reference manual, in the section on the cma_es function, states that "no 
sophisticated convergence detection is included" (p. 3)  If the user 
specifies a value for stopfitness, the computation will halt once the 
objective function value "is smaller than or equal to stopfitness" (p. 
3).  "This is the only way for CMA-ES to 'converge'" (p. 3). If we knew 
the minimum attainable value for -1*(expected utility), we could specify 
stopfitness as a number slightly greater than that value. 
Unfortunately, in portfolio optimization problems we seldom know the 
minimum attainable value of -1*(expected utility) until we have found 
the optimal portfolio.
Thanks for drawing my attention to BARON.  It looks good; when I've 
learned enough about GAMS format, I'll give it a try.

Best regards,
John

  
    
#
Brian G. Peterson wrote:
Thanks, Brian.  I'll double check my code to see if it's given DEoptim 
an impossible task.
These are very helpful leads.  I'll pursue them tomorrow morning. 
Thanks again!

Best regards,
John

  
    
#
With 120 observations and 30 assets and a smooth objective function and
relatively simple constraints, the running time should be *seconds*. One
potential problem with DE is how the constraints are implemented, more
specifically the lower/upper limits. When DE creates new solutions, it does
not observe any constraints except for such that one "implicitly" adds (eg,
like you did with the budget constraint). Thus, one typically repairs the
constraints after new solutions are created, or penalises solutions that
violate constraints. If you penalise too strongly, the algorithm will in
effect throw away many candidate solutions, and so the algorithm may indeed
have a hard time to find a solution, since in portfolio optimisation
problems the solution is often on the boundary. 

Regarding the installation error: the \title in one Rd-file uses quotation
marks (via \sQuote), which should be supported (according to R-exts since R
2.12.0). R CMD check never complained to me on Win XP or Ubuntu. Could you
please send me (off-list) the R-version/OS you work on? But I will remove
the quotation marks later today, so the new version should be available from
R-Forge tomorrow. Please let me know if this doesn't work either.

Regards,
Enrico
#
Another alternative to get the budget
constraint is to let the optimizer see
unconstrained solutions but then scale
the weights before evaluating the utility.

This has the possibility of confusing some
optimizers, but can work pretty well.
On 02/08/2011 03:24, John P. Burkett wrote:

  
    
#
n Tue, 2011-08-02 at 07:31 +0200, Enrico Schumann wrote:
This isn't *quite* correct.  DE will observe your lower and upper bounds
on each objective parameter (asset weight), but will not observe some
additional overall constraint such as a full investment constraint.
This mirrors Pat Burns' point about letting the optimizer do its thing,
and then 'adjusting' the parameter set inside your objective function to
apply the constraint.  For many portfolios, this works quite well, but
for those where it doesn't work, a penalty function needs to be applied.
I would put this difficulty a little differently.  If the constraint we
are concerned with is the full investment constraint, then the issue is
one of how large the feasible adjusted constrained space (portfolios
with weight bounds min_w and max_w where sum(w)==1) in comparison to the
total space the optimizer will search (all portfolios comprised of
assets in any combination of weights min_w to max_w).

If the feasible space is large in comparison to the total space, a
global optimizer will locate the solution.  If the feasible space is not
large in comparison to the total space, and a large penalty function is
used to tell the optimizer that it is not close, then it could take a
long time to find a solution.

There are several approaches to addressing this.  One is to graduate
your penalty response, for example penalizing with some multiplier
(which should be scaled somewhat larger than a 'good' value of the
objective function, so for global minima at -1, perhaps a multiplier of
10 or 100 will serve) the absolute value of the amount the sum of the
weights exceeds 1.  You can make this work even better by applying a
fraction of the penalty for any sum of weights that is less than one,
because this solution is likely closer to the feasible space.

I find that it is best to pre-seed the optimizer with a population of
candidate portfolios that are all inside the feasible space, so that
even random perturbations are still closer to the feasible space.  In
PortfolioAnalytics, we pre-seed be default, and use a full-investment
constraint adjustment by default (with a graduated penalty function as
another supported alternative), and of course all these options are
changeable by the user.
 
I haven't been able to locate the original Sharpe paper, but I suspect
he utilized a smooth objective function given the date of the paper
alone.  In that case, as I said previously, you don't need and *won't
benefit from* a global optimizer.  Some other author has probably
already cited Sharpe's paper and told you how to use an exact or
smooth-surface optimization technique. Of course, as soon as you want to
add another layer to your objective, or more complex constraints, you're
back in the land of global optimization.

Regards,

   - Brian
#
[comments inline]
Sorry if there was a misunderstanding: when I referred to 'DE', I meant
Differential Evolution in general, not 'DEoptim'. I don't know how the
constraints are incorporated in 'DEoptim'.
Just for the record: actually we have even two possibilities: really repair,
or just map the solutions into feasible ones. In the latter case, we go on
evolving the solutions that violate the constraints.
if it works properly, yes :)
Yes, the important point, I think, is that in such cases the feedback from
the penalty need be constructive, ie, "how much do I violate a constraint?"
(if I want to push some object from the table to the floor, then feedback
like 'the object is still on the table' is less useful than 'the object is
now 3cm off the edge of the table')

But you can also directly use the constraint in the creation of the new
solutions, which is what Sharpe suggested (and what works quite well, and
not just for smooth functions/constraints): if you only add zero-sum changes
to a feasible portfolio, it will remain feasible with respect to the budget
constraint. If you want min/max-holding sizes, choose asset weights such
that the min/max-constraints remain unviolated. This is straightforward in
methods like Simulated Annealing/Threshold Accepting, but somewhat more
difficult in 'DE'. (Which does not mean that 'DE' is not as good as the
other methods, just that the efficient approaches possibly differ, depending
on the method.)
best regards,
Enrico
#
Patrick Burns wrote:
Thank you, Patrick, for this suggestion.  I tried implementing it in 
DEoptim by writing the function to be minimized as shown below.  This 
code is intended to handle a long-only portfolio of 31 assets.  The 
variable x is meant to be a vector representing the shares of the assets 
in the portfolio. The Data matrix has dimensions 20x31. Each row 
corresponds a time period (six months for all but the first period, 
which is only four months long).  Each column corresponds to an asset. 
The elements of the matrix are gross returns. (The gross returns for the 
short first period are raised to the 3/2 power to make them comparable 
to those for the other periods.)  My attempt to rescale the variables to 
assure that they add up to one can be found 5 lines from the bottom of 
the code below.

Neu31<- function(x){
x1 <- x[1]
x2 <- x[2]
x3 <- x[3]
x4 <- x[4]
x5 <- x[5]
x6 <- x[6]
x7 <- x[7]
x8 <- x[8]
x9 <- x[9]
x10 <- x[10]
x11 <- x[11]
x12 <- x[12]
x13 <- x[13]
x14 <- x[14]
x15 <- x[15]
x16 <- x[16]
x17 <- x[17]
x18 <- x[18]
x19 <- x[19]
x20 <- x[20]
x21 <- x[21]
x22 <- x[22]
x23 <- x[23]
x24 <- x[24]
x25 <- x[25]
x26 <- x[26]
x27 <- x[27]
x28 <- x[28]
x29 <- x[29]
x30 <- x[30]
x31 <- x[31]
if (sum(x) == 0) {x <- x + 1e-2}
x = x/sum(x)
grp <- Data%*%x
  u <- grp/(.5 + grp)
objfun <- -sum(u)/nr
return(objfun)
}

To create initial shares summing to one, I used the following code:
npar = 31
NP = 10*npar
rum <- matrix(runif(npar*NP), ncol=npar, nrow=NP)
rowsums <- rowSums(rum)
ip <- rum/rowsums

Finally I invoked DEoptim as follows:
out <- DEoptim(Neu31, lower, upper, control=list(NP=NP, initial=ip, 
itermax=5000, F=.8))

That started the iterations, which progressed plausibly until getting 
stuck at iteration 4055.  At that point my computer became unresponsive 
and remained so until I rebooted.

Any suggestions for improving the code to obtain convergence would be 
greatly appreciated.

Best regards,
John

  
    
#
Enrico Schumann wrote:

            
Thanks, Enrico.  I'm interested in using the adding-up constraint in the 
creation of new solutions, via methods such as simulated annealing or 
threshold accepting.  I imagine that the DMOF package is useful for 
implementing that approach.  Are there also other packages that I should 
consider?

Best regards,
John
#
[see below]
There are many packages in R that handle optimisation, but I am not really
familiar with most if them. The task view for optimisation was already
suggested to you; you may also want to browse R-Forge/RForge.

I cannot comment on the model that you gave your other e-mail -- even
though, this objective function will essentially pick portfolios of just
very few assets with very high returns; I hope you have faith that your
return scenarios are representative of what will happen in the future.

Here is a complete example with TAopt from the NMOF package; it is mostly
copy--pasted from a package vignette. Since I don't have your data, I create
random data and put all data into one list which I later pass to the
function TAopt.

# --- #
require("NMOF")

# the neighbourhood: see the package vignette on portfolio optimisation
resample <- function(x, ...) x[sample.int(length(x), ...)] # see ?sample
neighbourU <- function(sol, data){
    wn <- sol$w
    toSell <- wn > data$winf
    toBuy  <- wn < data$wsup
    i <- resample(which(toSell), size = 1L)
    j <- resample(which(toBuy), size = 1L)
    eps <- runif(1) * data$eps
    eps <- min(wn[i] - data$winf, data$wsup - wn[j], eps)
    wn[i] <- wn[i] - eps
    wn[j] <- wn[j] + eps
    Rw <- sol$Rw + data$R[,c(i,j)] %*% c(-eps,eps)
    list(w = wn, Rw = Rw)
}

# random data
na <- 31L   # number of assets
nr <- 20L   # number of return obs
randomData <- rnorm(nr*na)*0.01  # random data: a plain matrix ;)
dim(randomData) <- c(nr, na)

# collect all in 'data': eps is stepsize for TA, winf/wsup are min/max
holding sizes

data <- list(R = randomData, na = na, nr = nr,  eps = 0.5/100, winf = 0,
wsup = 1)

# --- #

The neighbourhood may look a bit complicated, but it uses the following
fact: in a given iteration, your objective function can be split into two
pieces: first, you compute portfolio returns Rw <- Data %*% w; then you
compute your actual objective for these returns, say, fun(Rw). The
multiplication will typically take much of the computing time (in particular
if the actual objective is quite cheap); see the following test. But the
neighbourhood function changes only two weights: increase one, decrease one.
So you can update the multiplication. In your case, this may not make much
difference >> but try with 100 or more return observations

Thus, the representation of a solution is a list with two components: the
weight vector w, and the return associated with this portfolio Rw (ie, Rw <-
Data %*% w). In the objective function, you can compute anything from sol$w
and sol$Rw.


# ---
# create a feasible random solution 
w0 <- runif(data$na); w0 <- w0/sum(w0)
x0 <- list(w = w0, Rw = randomData %*% w0)
system.time(for (i in 1:100000) Rw <- randomData %*% w0)
system.time(for (i in 1:100000) ignore <- sum(Rw/(0.5 + Rw)))

# test the neighbourhood
x1 <- neighbourU(x0, data)
all.equal(data$R %*% x1$w, x1$Rw)
barplot(x0$w - x1$w)  # difference between x0 and x1

# the objective function
Neu31 <- function(sol, data){
  u <- sol$Rw/(0.5 + sol$Rw)
  objfun <- -sum(u)/data$nr
  objfun
}

# ---

Your objective: I have droppped the x1, x2, ... since they do not appear to
be needed. Also, you use an object 'nr', but never pass it. R will find it
if it is in the global environment, but I like it better to pass all objects
directly or make the dependence explicit by using environments. It remains
to run the algorithm.

# ---
# settings: see ?TAopt
algo <- list(x0 = x0, neighbour = neighbourU, nS = 2000L, nT = 10L,
             nD = 1000L, q = 0.20, printBar = FALSE, printDetail = FALSE)
system.time(res <- TAopt(Neu31, algo = algo,data = data))
res$OFvalue  # the obj fun value of the solution
res$xbest$w  # the best weights
sum(res$xbest$w)  # is budget satisfied?


# TA is a stochastic method: run restarts
allRes <- restartOpt(fun=TAopt, OF = Neu31, n = 20L, algo = algo, data =
data)
allF <- sapply(allRes, `[[`, "OFvalue")      # realised objective functions
allw <- sapply(allRes,`[[`, c("xbest","w"))  # weights



hope that helps,
Enrico
1 day later
#
Just to round out this thread on-list, I've attached a script for the
problem as described by Prof. Burkett that uses the Munger and Arrow
bounded utility function to construct a portfolio via PortfolioAnalytics
using both random portfolios and DEoptim as alternate optimization
engines.

I've modified the script to use the 'edhec' data series included with
PerformanceAnalytics for reproducibility.

As Enrico noted in an earlier post, this utility function tends towards
creating very concentrated portfolios in only a small number of assets.

Regards,

   - Brian
On Mon, 2011-08-01 at 22:57 -0400, John P. Burkett wrote:

  
    
3 days later
#
Brian G. Peterson wrote:
Thanks again for your ingenious script and illuminating comments.

The only contribution I can make to the discussion now is to provide 
some references on bounded utility functions.  Karl Menger (a 
mathematician and son of the Austrian economist Carl Menger) made a case 
for bounded utility functions in a 1934 article in German.  An English 
translation was published as "The Role of Uncertainty in Economics", ch. 
16 in Essays in Mathematical Economics in Honor of Oskar Morenstern, ed. 
Martin Shubik (Princeton University Press, 1967).  Other arguments for 
bounded utility functions have been made by Kenneth J. Arrow, Essays in 
the Theory of Risk-bearing (Markham, 1971) and by John W. Pratt, Howard 
Raiffa, and Robert Schlaifer, Introduction to Statistical Decision 
Theory (MIT Press, 1995), pp. 80-81.  The particular utility function I 
have been using is just one example of a bounded utility function, 
expressed in terms of gross return. To the best of my knowledge, none of 
the authorities cited above have recommended this particular function. 


Best regards,
John

  
    
#
Care to enlighten the list so someone else doesn't have the same problem
you had?  - Brian
On Mon, 2011-08-08 at 21:14 -0400, Yihao Lu aeolus_lu wrote:
3 days later
#
Sry about the threadjack. I got this working by using an older version of rJava however it looks the behavior has changed in the new version when you fetch a list of securities. 

The resultant data.frame has security1 followed by security 2..followed by security 3 rather than a column bind as in the previous versions (which didn't use rJava)

I must not have read the new documentation but is there a way to go back to the old RBloomberg which returned zoo objects and did the merge of the data than the append.

 Is there someplace where I can find the older win32 binary?.

Thanks & Best,
Krishna
On Aug 9, 2011, at 8:51 AM, Yihao Lu aeolus_lu wrote: