w ̂=argmin w'Σ w
s. t.  w'I = 1           #weights sum up to 1
       w'μ=ρ             #target expected return
       w≥0               #non-negativity(short-sale) constraint

optimization<-function(returns) {
  p <- ncol(x)                    #number of assets
  n <- nrow(x)                    #number of observations
  x <- matrix(data$return,n,assets)
  mean <- colMeans(na.rm=FALSE,x)
  M <- as.integer(10)             #nuber of ports on the eff.front.
  S <- cov(x)                     #covariance matrix
  Rmax<- 0.01                     #max monthly return value
  Dmat   <- solve(S)              #inverse of covariance matrix
  u <- rep(1,p)                   #vector of ones

  These codes are for the Lagrange solutions

  a<- matrix(rep(0,4), nrow=2)
  a[1,1] <- t(u)%*% Dmat %*%u
  a[1,2] <- t(mean)%*%Dmat%*%u
  a[2,1] <- a[1,2]
  a[2,2] <- t(mean)%*%Dmat%*%mean
  d <- a[1,1]*a[2,2]-a[1,2]*a[1,2]
  f <- (Dmat%*%(a[2,2]*u-a[1,2]*mean))/d
  g <- (Dmat%*%(-a[1,2]*u+a[1,1]*mean))/d
  r <- seq(0, Rmax, length=M)
  w <- matrix((rep(0, p*M)), nrow=p)

for(i in 1:M) { w[,i] = f+r[i]*g                    #portfolio weights
      if (w[,i] <0) {w[,i]=0} else {w[,i]=w[,i]}
    }

while (w> 0)
  { for(i in 1:M) { w[,i] = f+r[i]*g }
  print(w)
  }