#
# Written by:
# -- 
# John L. Weatherwax                2009-04-21
# 
# email: wax@alum.mit.edu
# 
# Please send comments and especially bug reports to the
# above email address.
#
# Computes the Efficient Frontier.
#
#-----

library(quadprog)

efficient_frontier = function(returns, muP, mu_free=0.0, w_lower_limit=-Inf, w_upper_limit=+Inf){
  #
  # w_lower_limit and w_upper_limit (if finite) are bounds on the porfolio weights such that
  # 
  # w_lower_limit <= w <= w_upper_limit
  #
  # Note: if solve.QP gives errors about convergence one can try a more restricted range for the values of muP:
  #
  # muP = seq(min(mean_vect), max(mean_vect), length.out=500)
  #
  #--
  n_stocks = dim(returns)[2]
  
  # Extract individual equity (mean return, std return) values: 
  mean_vect = apply(returns, 2, mean)
  cov_mat = cov(returns)
  sd_vect = sqrt(diag(cov_mat))
  
  # Portfolio weight constraints:
  #
  #     first condition enforce \sum w_j = 1
  #     second condition enforce \sum mu_j w_j = mu_target
  #             
  Amat = cbind(rep(1,n_stocks),mean_vect)
  bvec = c(1,NaN)
  if( is.finite(w_lower_limit) ){
    Amat = cbind(Amat, diag(1, nrow=n_stocks))
    bvec = c(bvec, w_lower_limit*rep(1,n_stocks))
  }
  if( is.finite(w_upper_limit) ){
    Amat = cbind(Amat, -diag(1, nrow=n_stocks))
    bvec = c(bvec, -w_upper_limit*rep(1,n_stocks))
  }

  # storage for results: 
  sdP = muP
  weights = matrix(0, nrow=length(muP), ncol=n_stocks)

  # find the optimal portfolios for each target expected return:
  #
  for( i in 1:length(muP) ){
    bvec[2] = +muP[i] # enforce portfolio mean constraint 
    result = solve.QP(Dmat=2*cov_mat, dvec=rep(0,n_stocks), Amat=Amat, bvec=bvec, meq=2)
    sdP[i] = sqrt(result$value)
    weights[i,] = result$solution 
  }

  # Find maximum Sharpe portfolio: 
  sharpe = ( muP - mu_free ) / sdP 
  ind_ms = which.max(sharpe)

  # Find minimum variance portfolio: 
  ind_mv = which.min(sdP)
  
  list( muP=muP, sdP=sdP, weights=weights, sharpe=sharpe, max_sharpe_ind=ind_ms, min_variance_ind=ind_mv ) 
}