plot.efficientFrontier <- function(Dat, mufree=NULL, noSS=FALSE, maxMuP=NULL,
                                   plot.tagent=TRUE, xlim=NULL, ylim=NULL){
    ## plot of the efficient frontier including tagent portfolio
    ## Dat:         return data of the assets
    ## mufree:      input value of risk-free interest rate
    ## noSS:        no short selling if TRUE, short selling allowed if FALSE
    ## maxMuP:      maximum value of mu_P
    ## plot.tagent: plot tagent of tagent portfolio (=TRUE) or not (=FALSE)
    ## 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 deviation
    N <- length(riskH) # number of assets

    if(noSS){
        Amat <- cbind(rep(1,N), muH, diag(1,N))
                                        # set the constraints matrix
        ## no short selling --> target muP is between min and max of assets
        muP <- seq(min(muH) + 0.0001, max(muH) - 0.0001, length=300)
    } else {
        Amat <- cbind(rep(1,N), muH)  # set the constraints matrix
        if(is.null(maxMuP)) maxMuP <-  max(muH)*2
        muP <- seq(.01, maxMuP, 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
        if(noSS){
            bvec <- c(1,muP[i], rep(0,N))  # constraint vector
        } else {
            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
    if(is.null(xlim)) xlim <- range(c(0,riskP))
    if(is.null(ylim)) ylim <- 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
        if(plot.tagent){ # if =TRUE --> plot tagent
            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)
    if(is.null(mufree)){
        Value <- list(mvp=c(riskP[id2], muP[id2]), Wmvp=weights[id2,],
                  tp=NULL, Wtp=NULL, tpSharpe=NULL)
    } else{
        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, maxMu=NULL, noSS=FALSE, 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
    ## noSS:    TRUE = short selling not allowed , FALSE = short selling allowed
    ## nBoot:   number of bootstrap simulations
    ## setSeed: if not NULL, give a "seed number"
    ## 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=4)
    colnames(out) <- c("actual", "estimated", "muTP", "riskTP")
    ## mu.out <- matrix(1,nrow = nBoot, ncol=ncol(Dat))
    n <- nrow(Dat)
    N <- ncol(Dat)
    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)
        SigmaBoot <- estSigma(Rboot)
        if(noSS){
            muP <- seq(min(muBoot)+ 0.001, max(muBoot) - 0.001,
                       length=300)
            Amat <- cbind(rep(1,N), muBoot, diag(1,N))
        } else {
            MaxMu <- ifelse(missing(maxMu), 1.4, maxMu)
             muP <- seq(.01, maxMu, length=300)
            Amat <- cbind(rep(1,N),muBoot)
        }
        riskP <- muP
        weights <- matrix(0,nrow=300, ncol=N)
        for (i in 1:length(muP)) {
            if(noSS){
                bvec <- c(1,muP[i], rep(0,N))  # constraint vector
            } else {
                bvec <- c(1,muP[i])  # constraint vector
            }
            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,2:4] <- c(sharpe[id], muP[id], riskP[id])
        wTP <- weights[id,]
        sharpeActual <- (wTP%*%muTRUE - mufree) /sqrt(wTP %*%SigmaTRUE%*% wTP)
        out[iBoot,1] <- as.numeric(sharpeActual)
    }

return(as.data.frame(out))
}

## ======================================================================
BootTP.muFix <- function(Dat, mufree, nBoot=250, setSeed=NULL){
    ## Runs a bootstrap simulation of the tangency portfolio
    ## 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)
    muTRUE <- apply(Dat,2,mean)
    SigmaTRUE <- cov(Dat)
    estSigma <- function(x, mu=muTRUE){
        xx <- as.matrix(scale(x, center=mu, scale=FALSE))
        return((t(xx) %*% xx)/(nrow(x)-1))
    }

    ##
    out <- matrix(1, nrow=nBoot, ncol=4)
    colnames(out) <- c("actual", "estimated", "muTP", "riskTP")
    ##
    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 <- muTRUE
        ##
        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,2:4] <- c(sharpe[id], muP[id], riskP[id])
        wTP <- weights[id,]
        sharpeActual <- (wTP%*%muTRUE - mufree) /sqrt(wTP %*%SigmaTRUE%*% wTP)
        out[iBoot,1] <- as.numeric(sharpeActual)
    }

return(as.data.frame(out))
}
## =======================================================================
BootTP.sigFix <- function(Dat, mufree, maxMu=NULL, nBoot=250, setSeed=NULL){
    ## Runs a bootstrap simulation of the tangency portfolio
    ## 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)
    muTRUE <- apply(Dat,2,mean)
    SigmaTRUE <- cov(Dat)
    estSigma <- function(x, Sigma=SigmaTRUE) Sigma
    ##
    out <- matrix(1, nrow=nBoot, ncol=4)
    colnames(out) <- c("actual", "estimated", "muTP", "riskTP")
    ##
    n <- nrow(Dat)
    N <- ncol(Dat)
    MaxMu <- ifelse(missing(maxMu), 1.4, maxMu)
    muP <- seq(.01, maxMu, 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,2:4] <- c(sharpe[id], muP[id], riskP[id])
        wTP <- weights[id,]
        sharpeActual <- (wTP%*%muTRUE - mufree) /sqrt(wTP %*%SigmaTRUE%*% wTP)
        out[iBoot,1] <- as.numeric(sharpeActual)
    }

return(as.data.frame(out))
}
## =======================================================================

BootMVP <- function(Dat, nBoot=250, setSeed=NULL){
    ## Runs a bootstrap simulation of the minimum variance portfolio
    ## Dat:    return data of the assets
    ## Formula "Sigma^(-1) 1 /(1^T Sigma^(-1) 1)"  see, e.g., Zivot 2013
    ##
    n <- nrow(Dat)
    N <- ncol(Dat)
    one.vec <- rep(1, N)
    wMVP <- matrix(1, nrow=nBoot, ncol=N)
    out <- matrix(1, nrow=nBoot, ncol=2)
    colnames(out) <- c("riskMVP", "muMVP")

    if(!is.null(setSeed)) set.seed(setSeed)
    ##
    for (iBoot in (1:nBoot)){
        Rboot <- Dat[sample(1:n,n,replace=TRUE),]
        SigmaBoot <- cov(Rboot)
        SigmaBoot.inv <- solve(SigmaBoot)
        h.top.mat <- SigmaBoot.inv %*% one.vec
        h.bot.val <- as.numeric(t(one.vec)  %*%  h.top.mat)
        wMVP[iBoot,] <- h.top.mat / h.bot.val
        out[iBoot,] <- c( sqrt(wMVP[iBoot,]^T %*% SigmaBoot %*% wMVP[iBoot,]),
                         sum(apply(Rboot,2,mean) * wMVP[iBoot,]))
    }

    return(list(out=as.data.frame(out), weightsMVP=wMVP))
}