plot.efficientFrontier <- function(Dat, mufree=NULL){
    ## plot of the efficientfrontier
    ## Dat:    return data of the assets
    ## mufree: input value of risk-free interest rate
    ## Impelementation is based on R code of David Ruppert.
    require(quadprog)
    muH <- apply(Dat,2,mean)
    SigmaH <- cov(Dat)
    riskH <- sqrt(diag(SigmaH)) # risk = standard devaiation
    N <- length(riskH) # number of assets

    Amat <- cbind(rep(1,N), muH)  # set the constraints matrix
    muP <- seq(.01, max(muH)*1.4, length=300) # set of 300 possible target
                                      # values for the expect portfolio return
    riskP <- muP # set up storage for standard deviations of portfolio returns
    weights <- matrix(0, nrow=300, ncol=N) # storage for portfolio weights

    for (i in 1:length(muP)){# find the optimal portfolios for each target
                             # expected return
        bvec <- c(1,muP[i])  # constraint vector
        result <- solve.QP(Dmat=2*SigmaH, dvec=rep(0,N), Amat=Amat,
                           bvec=bvec, meq=2)
        riskP[i] <- sqrt(result$value)
        weights[i,] <- result$solution
    }
     ##
    ## Display of the result
    xlim <- if(is.null(mufree)) range(riskP) else range(c(0,riskP))
    ylim <- if(is.null(mufree)) range(muP) else range(c(mufree,muP))
    ## plot of the efficient frontier (and inefficient frontier)
    plot(riskP, muP, type="l", lty=3, las=1, xlim=xlim, ylim=ylim)
    if(! is.null(mufree))
        points(0, mufree, cex=4, pch="*", col="blue")  # show risk-free asset
    sharpe <- (muP-mufree)/riskP # compute Sharpe ratios
    id <- which(sharpe == max(sharpe)) # Find maximum Sharpe ratio
    xlop <- c(0,riskP[id]*1.2) # x range of line of optimal portfolios
    lines(xlop, mufree+xlop*(muP[id]-mufree)/riskP[id], lwd=2, col="gray")
           # show line of optimal portfolios
    points(riskP[id],muP[id], cex=2, pch=16, col="blue") # show tangency portfolio
    id2 <- which(riskP == min(riskP)) # find the minimum variance portfolio
    points(riskP[id2], muP[id2], cex=2, pch="+", col="red")
                                        # show minimum variance portfolio
    id3 <- (muP > muP[id2])
    lines(riskP[id3], muP[id3], type="l", lwd=2, col="magenta")
                                        # mark the efficient frontier
    text(riskH, muH, labels=colnames(Dat), cex=0.7)
    Value <- list(mvp=c(riskP[id2], muP[id2]), Wmvp=weights[id2,],
                  tp=c(riskP[id],muP[id]), Wtp=weights[id,],
                  tpSharpe=sharpe[id])
    invisible(Value)
}
## -------------------------------------------------------------------------
h.covByFacA <- function(x, NoFAC){
    sDcov <- diag(sqrt(diag(cov(x))))
    h.fa <- factanal(x, factors=NoFAC, rotation="none")
    h.l <- loadings(h.fa)
    CorA <- h.l %*% t(h.l) + diag(h.fa$uniquenesses)
    sDcov %*% CorA %*% sDcov
}

## -------------------------------------------------------------------------

BootTP <- function(Dat, mufree, nBoot=250, setSeed=NULL, NoFAC=NULL){
    ## Runs a bootstrap simulation of the tangency portfolio
    ## Dat:    return data of the assets
    ## mufree: input value of risk-free interest rate
    ## NoFAC:  an integer > 0. If not NULL, a "factor analysis" approximation
    ##         with NoFAC components is used to estimate the covariance matrix.
    ## Impelementation is based on R code of David Ruppert.
    ##
    require(quadprog)
    muTRUE <- apply(Dat,2,mean)
    SigmaTRUE <- cov(Dat)
    estSigma <- {if(is.null(NoFAC)) cov
                 else {function(x) h.covByFacA(x, NoFAC=NoFAC)}}
    ##
    out <- matrix(1, nrow=nBoot, ncol=2)
    mu.out <- matrix(1,nrow = nBoot, ncol=ncol(Dat))
    n <- nrow(Dat)
    N <- ncol(Dat)
    muP <- seq(.01, max(muTRUE)*1.4, length=300)

    if(!is.null(setSeed)) set.seed(setSeed)
    ##
    for (iBoot in (1:nBoot)){
        Rboot <- Dat[sample(1:n,n,replace=TRUE),]
        muBoot <- apply(Rboot,2,mean)
        mu.out[iBoot,] <- muBoot
        SigmaBoot <- estSigma(Rboot)
        Amat <- cbind(rep(1,N),muBoot)
        riskP <- muP
        weights = matrix(0,nrow=300, ncol=N)
        for (i in 1:length(muP)) {
            bvec <- c(1,muP[i])
            result <- solve.QP(Dmat=2*SigmaBoot, dvec=rep(0,N), Amat=Amat,
                               bvec=bvec, meq=2)
            riskP[i] <- sqrt(result$value)
            weights[i,] <- result$solution
        }
        sharpe <- (muP-mufree)/riskP
        id <- (sharpe == max(sharpe))
        out[iBoot,1] <- sharpe[id]
        wTP <- weights[id,]
        sharpeActual <- (wTP%*%muTRUE - mufree) /sqrt(wTP %*%SigmaTRUE%*% wTP)
        out[iBoot,2] <- as.numeric(sharpeActual)
    }

return(data.frame(actual=out[,2], estimated=out[,1]))
}