Markowitz model / portfolio optimization using local search in R

Question

I am taking baby steps to use metaheuristics for solving constrained optimization problems. I am trying to solve basic Markowitz Mean-Variance optimization model (given below) using NMOFpackage in R.

  Min 
   lambda * [sum{i=1 to N}sum{j = 1 to N}w_i*w_i*Sigma_ij] - (1-lambda) * [sum{i=1 to N}(w_i*mu_i)]
 subject to 
    sum{i=1 to N}{w_i} = 1
    0 <= w_i <= 1; i = 1,...,N

where, lambda takes values between 0 and 1, N is number of assets.

Following is my code (Based on Book: Numerical Methods and Optimization in Finance):

    library(NMOF)

na <- dim(fundData)[2L]
ns <- dim(fundData)[1L]
Sigma <- cov(fundData)
winf <- 0.0
wsup <- 1.0
m <- colMeans(fundData)


resample <- function(x,...) x[sample.int(length(x),...)]


data <- list(R = t(fundData), 
             m = m,
             na = dim(fundData)[2L],
             ns = dim(fundData)[1L], 
             Sigma = Sigma, 
             eps = 0.5/100,
             winf = winf, 
             wsup = wsup,
             nFP = 100)

 w0 <- runif(data$na); w0 <- w0/sum(w0)

 OF <- function(w,data){    
   wmu <- crossprod(w,m)
   res <- crossprod(w, data$Sigma)
   res <- tcrossprod(w,res)
   result <- res - wmu
  }


 neighbour <- function(w, data){
   toSell <- w > data$winf
   toBuy <- w < data$wsup
   i <- resample(which(toSell), size = 1L)
   j <- resample(which(toBuy), size = 1L)
   eps <- runif(1) * data$eps
   eps <- min(w[i] - data$winf, data$wsup - w[j], eps)
   w[i] <- w[i] - eps
   w[j] <- w[j] + eps
   w
 }


 algo <- list(x0 = w0, neighbour = neighbour, nS = 5000L)
 system.time(sol1 <- LSopt(OF, algo, data))

I am not sure how to include lambda in the objective function (OF). The above code does not include lambda in OF. I tried using for loop but it resulted in following error:

OF <- function(w,data){

  lambdaSeq <- seq(.001,0.999, length = data$nFP)
  for(lambda in lambdaSeq){
  wmu <- crossprod(w,m)
  res <- crossprod(w, data$Sigma)
  res <- tcrossprod(w,res)
  result <- lambda*res - (1-lambda)*wmu 
  }
}

Error:

 Local Search.
    Initial solution:   
      |                                                                                          |   0%
    Error in if (xnF <= xcF) { : argument is of length zero
    Timing stopped at: 0.01 0 0.03

It would be nice if someone could help me in this regard.

P.S: I am also aware that this can be solved using quadratic programming. This is just an initiation to include other constraints.

Answer 1

If I understand correctly, you want to replicate the mean--variance efficient frontier by Local Search? Then you need to run a Local Search for every value of lambda that you want to include in the frontier.

The following example should help you get going. I start by attaching the package and setting up the list data.

require("NMOF")
data <- list(m = colMeans(fundData), ## expected returns
             Sigma = cov(fundData),  ## expected var of returns
             na = dim(fundData)[2L], ## number of assets
             eps = 0.2/100,          ## stepsize for LS
             winf = 0,               ## minimum weight
             wsup = 1,               ## maximum weight
             lambda = 1)

Next I compute a benchmark for the minimum-variance case (i.e. lambda equals one).

## benchmark: the QP solution
## ==> this will only work with a recent version of NMOF,
## which you can get by saying:
##   install.packages('NMOF', type = 'source',
##                    repos = c('http://enricoschumann.net/R',
##                              getOption('repos')))
##
require("quadprog")
sol <- NMOF:::minvar(data$Sigma, 0, 1)

Objective function and neighbourhood function. I have slightly simplified both functions (for clarity; using crossprod in the objective function would probably be more efficient).

OF <- function(w, data){
    data$lambda * (w %*% data$Sigma %*% w) -
    (1 - data$lambda) * sum(w * data$m)
}

neighbour <- function(w, data){
    toSell <- which(w > data$winf)
    toBuy  <- which(w < data$wsup)
    i <- toSell[sample.int(length(toSell), size = 1L)]
    j <- toBuy[sample.int(length(toBuy), size = 1L)]
    eps <- runif(1) * data$eps
    eps <- min(w[i] - data$winf, data$wsup - w[j], eps)
    w[i] <- w[i] - eps
    w[j] <- w[j] + eps
    w
}

Now we can run Local Search. Since it is a fairly large dataset (200 assets), you will need a relatively large number of steps to reproduce the QP solution.

w0 <- runif(data$na) ## a random initial solution
w0 <- w0/sum(w0)
algo <- list(x0 = w0, neighbour = neighbour, nS = 50000L)
sol1 <- LSopt(OF, algo, data)

You can compare the weights you get from Local Search with the QP solution.

par(mfrow = c(3,1), mar = c(2,4,1,1), las = 1)
barplot(sol, main = "QP solution")
barplot(sol1$xbest, main = "LS solution")
barplot(sol - sol1$xbest,
        ylim = c(-0.001,0.001)) ## +/-0.1%

Finally, if you want to compute the whole frontier, you need to rerun this code for different levels of data$lambda.