rm(list=ls())
library(tseries)
library(zoo)
library(quadprog) #Needed for solve.QP

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

SHORT<-'no'
max_risk<-1.0
risk_increment<-0.05
max_alloc<-NULL #0.25

divers <- function (rtn, part.cov, short, risk.premium.up, risk.increment, max.allocation){

  rs<-rtn[,-1]
  n<-ncol(rs)
  n1<-n-1
  covariance <- cov(rs)
    
  # Create initial Amat and bvec assuming only equality constraint
  # (short-selling is allowed).
  Amat <- matrix (1, nrow=n)
  bvec <- 1
  meq <- 1
  
  # Then modify the Amat and bvec if short-selling is prohibited
  if(short=="no"){
    Amat <- cbind(1, diag(n))
    bvec <- c(bvec, rep(0, n))
  }
  
  # And modify Amat and bvec if a max allocation (concentration) is specified
  if(!is.null(max.allocation)){
    if(max.allocation > 1 | max.allocation <0){
      stop("max.allocation must be greater than 0 and less than 1")
    }
    if(max.allocation * n < 1){
      stop("Need to set max.allocation higher; not enough assets to add to 1")
    }
    Amat <- cbind(Amat, -diag(n))
    bvec <- c(bvec, rep(-max.allocation, n))
  }
  
  # Calculate the number of loops
  loops <- risk.premium.up / 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=risk.premium.up, by=risk.increment)){
    dvec <- colMeans(rs) * i # This moves the solution along the EF
    sol <- solve.QP(covariance, dvec=dvec, Amat=Amat, bvec=bvec, meq=meq)
    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, short=SHORT, risk.premium.up=max_risk, 
                    risk.increment=risk_increment, max.allocation=max_alloc)

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