rm(list=ls())
library(tseries)
library(zoo)
library(nloptr) 

startDate<-"2013-01-02"
endDate<-"2013-12-31"

max_risk<-1.0
risk_increment<-0.05

eval_f <- function(x, covar, rr, ii) {
  z<-0
  g<-0
  z<-NULL
  n<-nrow(covar)
  R<-sum(colMeans(rr)*x)
  s<-sum(x*colSums(covar*x))
  for(i in 1:n) {
    g<-0  
    for(j in 1:n) {
      g<-g+x[i]*covar[i,j]
    }
    g<-0.5*g-colMeans(rr)[i]
    z<-c(z, g)
  }
  return(list("objective"= 0.5*s-ii*R, "gradient"=z))
}


eval_g_eq <- function(x, covar, rr, ii) {
  n<-nrow(covar)
  s<-rep(1,n)
  return(list("constraints"= sum(x)-1, "jacobian"=s))
  return(sum(x)-1)
}

divers <- function (rtn, part.cov, max_risk, risk.increment) {
 
  rs<-rtn[,-1]
  n<-ncol(rs)
  n1<-n-1

  covariance <- cov(rs)
  opts <- list("algorithm"="NLOPT_LD_SLSQP","xtol_rel"=1.0e-8, "maxeval"=30000)

  x0<-rep(1/n, n)
  
  lb<-NULL
  ub<-NULL
  for(i in 1:n) {
    lb<-c(lb, 0)
    ub<-c(ub, 1)
  }
  # Calculate the number of loops
  loops <- max_risk/ risk.increment + 1
  loop <- 1
  
  # Initialize a matrix to contain allocation and statistics
  # This is not necessary, but speeds up processing and uses less memory
  eff <- matrix(nrow=loops, ncol=n+3)
  # Now I need to give the matrix column names
  colnames(eff) <- c(colnames(rs), "Std.Dev", "Exp.Return", "sharpe")
  
  # Loop through the quadratic program solver
  for (i in seq(from=0, to=max_risk, by=risk.increment)){
    #if(i>0) break
    sol <- nloptr(x0=x0, eval_f=eval_f, lb=lb, ub=ub, eval_g_eq=eval_g_eq, opts=opts, covar=covariance, rr=rs, ii=i)
    eff[loop,"Std.Dev"] <- sqrt(sum(sol$solution*colSums((covariance*sol$solution))))
    eff[loop,"Exp.Return"] <- as.numeric(sol$solution %*% colMeans(rs))
    eff[loop,"sharpe"] <- eff[loop,"Exp.Return"] / eff[loop,"Std.Dev"]
    eff[loop,1:n] <- sol$solution
    loop <- loop+1
  }
  
  return(as.data.frame(eff))
}

etf_list <- read.table("C:/R/SPYETF3.csv", header=TRUE, sep=",")
Netf<-nrow(etf_list)


hdr<-NULL
prc_list=get.hist.quote(instrument = etf_list[1,1], startDate, endDate, quote = "AdjClose",   provider = "yahoo" );
for(i in 1:Netf) {hdr<-c(hdr, sprintf("%s",etf_list[i,1]))}
for(i in 2:Netf) {
  prc <- get.hist.quote(instrument = etf_list[i,1], startDate, endDate, quote = "AdjClose",   provider = "yahoo" )
  prc_list<-cbind(prc_list, prc)
}
names(prc_list)<-hdr

Lt<-nrow(prc_list)
init_list<-NULL
for(i in 1:Lt) 
    init_list<-rbind(init_list, prc_list[i])



r_list<-diff(log(init_list[,1]))
for(i in 2:Netf)
  r_list<-cbind(r_list, diff(log(init_list[,i])))
names(r_list)<-hdr
eff <- divers(rtn=r_list, part.cov=PART.COV, max_risk=max_risk, risk.increment=risk_increment)

# Find the optimal portfolio
opt.sharpe <- eff[eff$sharpe==max(eff$sharpe),][1,]   
hdr2<-c(hdr, "std", "rtn", "sharpe")
hdr2<-hdr2[-1]
weights<-opt.sharpe

print(opt.sharpe)