[R-SIG-Finance] Static Portfolio Optimization

[Thomas Etheber]

Brian G. Peterson wrote:
> Jesse Velez wrote:
>> Is there any function or example in R or Rmetrics of static portfolio
>> optimization, where I have a vector of expected returns for N assets
>> and a
>> expected covariance matrix of said N assets all at a fixed time (say
>> generated from a MFM risk and return model).
>>
>> fPortfolio, Portfolio, portfolio.optim appear to all require time
>> series of
>> returns to generate the expected return and historical covariance
>> matrix for
>> use in creating weights.
>>
>> Ideally, I hope to find an example that allows easily allows Long/Short
>> weights to make the portfolio  market neutral (i.e. Summation of Weights
>> =0).
>>   
> All the implementations of Markowitz style mean/variance optmization
> use quadprog in R.
>
> Plenty of information on the list archives from before all these
> packages existed about using quadprog for optimization.
>
> Regards,
>
>  - Brian
>

Hi there,

I also had the problem with fixed parameter inputs some time ago.
Implementing methods to perform this tasks would certainly be a nice
improvement of the library (as would be some help/error messages if the
covariance matrix is not positive semidefinite).  
Although Brian's comment is helpful as usual, using basic quadprog
sounds like reinventing the wheel, but might nevertheless be needed to
solve your second task of a market-neutral portfolio.

In order to use prespecified estimates as inputs I helped myself with
overwriting some of the methods. It's not a nice solution, but it worked
for me. You will find the methods attached below.
I didn't check the code again, but I think it should work. Please note,
some other methods of Rmetrics and fPortfolio might rely on the
timeseries objects and might not work properly.

Hth
Thomas

>>>
require(MBESS)
require(fPortfolio)
rm(list=ls())
spec <- portfolioSpec()
constraints <- NULL

portfolioData <- function (data, spec = portfolioSpec())
{  
    ans = NULL
    if(class(data) == "timeSeries") {
       data = sort(data)
       nAssets = dim(data)[2]
       statistics = portfolioStatistics(data, spec)
       tailRisk = spec at model$tailRisk
       ans <- new("fPFOLIODATA", data = list(series = data, nAssets =
nAssets),
           statistics = statistics, tailRisk = tailRisk)
    }
    if(class(data) == "list") {
      statistics = list(mu = data$mu, Sigma = data$Sigma )
      attr(statistics, "estimator") = spec at model$estimator
      ans <- new("fPFOLIODATA", data = list( nAssets = length(data$mu)
), statistics = statistics, tailRisk = list() )
    }
    ans
}


####################################################################################

.efficientConstrainedMVPortfolio <- function (data, spec, constraints)
{
    if (!inherits(data, "fPFOLIODATA"))
        data = portfolioData(data, spec)
    mu = getMu(data)
    Sigma = getSigma(data)
    nAssets = getNumberOfAssets(data)
    targetAlpha = getTargetAlpha(spec)
    solver = getSolver(spec)
    stopifnot(solver == "quadprog" | solver == "Rdonlp2")
    if (solver == "quadprog") {
        portfolio = solveRQuadprog(data, spec, constraints)
    }
    else if (solver == "Rdonlp2") {
        portfolio = solveRDonlp2(data, spec, constraints)
    }
    weights = portfolio$weights
    attr(weights, "status") <- portfolio$status
    names(weights) = names(mu)
    targetReturn = matrix(as.numeric(mu %*% weights), nrow = 1)
    colnames(targetReturn) <- getEstimator(spec)[1]
    covTargetRisk = sqrt(as.numeric(weights %*% Sigma %*% weights))
#   x = getSeries(data)@Data %*% weights
#   VaR = quantile(x, targetAlpha, type = 1)
#   CVaR = VaR - 0.5 * mean(((VaR - x) + abs(VaR - x)))/targetAlpha
#   targetRisk = matrix(c(covTargetRisk, CVaR, VaR), nrow = 1)
#   colnames(targetRisk) <- c("cov", paste(c("CVaR.", "VaR."),
#   targetAlpha * 100, "%", sep = ""))
    targetRisk = matrix(c(covTargetRisk), nrow = 1)
    ## is needed to use the plotting functions....
    targetRisk = matrix(c(covTargetRisk, covTargetRisk ), nrow = 1)
    colnames(targetRisk) <- c( "cov", "dummy" )
    new("fPORTFOLIO", call = match.call(), data = list(data = data),
        spec = list(spec = spec), constraints = as.character(constraints),
        portfolio = list(weights = weights, targetReturn = targetReturn,
            targetRisk = targetRisk, targetAlpha = targetAlpha,
            status = portfolio$status), title = paste("Constrained MV
Portfolio - Solver:",
            solver), description = .description())
}

####################################################################################

.minvarianceConstrainedMVPortfolio <- function (data, spec, constraints)
{
    if (!inherits(data, "fPFOLIODATA"))
        data = portfolioData(data, spec)
    mu = getMu(data)
    Sigma = getSigma(data)
    nAssets = getNumberOfAssets(data)
    targetAlpha = getTargetAlpha(spec)
    .minVariancePortfolioFun = function(x, data, spec, constraints) {
        spec at portfolio$targetReturn = x
        ans = .efficientConstrainedMVPortfolio(data = data, spec = spec,
            constraints = constraints)
        f = getTargetRisk(ans)[1]
        attr(f, "targetReturn") <- getTargetReturn(ans)
        attr(f, "targetRisk") <- getTargetRisk(ans)[1]
        attr(f, "weights") <- getWeights(ans)
        f
    }
    minVar = optimize(.minVariancePortfolioFun, interval = range(mu),
        data = data, spec = spec, constraints = constraints,
        tol = .Machine$double.eps^0.5)
    weights = attr(minVar$objective, "weights")
    names(weights) = names(mu)
    targetReturn = spec at portfolio$targetReturn =
as.numeric(attr(minVar$objective,
        "targetReturn"))
    targetReturn = matrix(targetReturn, nrow = 1)
    colnames(targetReturn) <- spec at model$estimator[1]
    covTargetRisk = as.numeric(attr(minVar$objective, "targetRisk"))
    # x = getSeries(data)@Data %*% weights
    # VaR = quantile(x, targetAlpha, type = 1)
    # CVaR = VaR - 0.5 * mean(((VaR - x) + abs(VaR - x)))/targetAlpha
    #targetRisk = matrix(c(covTargetRisk, CVaR, VaR), nrow = 1)
    #colnames(targetRisk) <- c("cov", paste(c("CVaR.", "VaR."),
    targetRisk = matrix(c(covTargetRisk), nrow = 1)
    ## is needed to use the plotting functions....
    targetRisk = matrix(c(covTargetRisk, covTargetRisk ), nrow = 1)
    colnames(targetRisk) <- c( "cov", "dummy" )
    new("fPORTFOLIO", call = match.call(), data = list(data = data),
        spec = list(spec = spec), constraints = as.character(constraints),
        portfolio = list(weights = weights, targetReturn = targetReturn,
            targetRisk = targetRisk, targetAlpha = targetAlpha,
            status = 0), title = "Minimum Variance Portfolio",
        description = .description())
}

show.fPORTFOLIO <- function (object)
{
    cat("\nTitle:\n ")
    cat(getTitle(object), "\n")
    cat("\nCall:\n ")
    print.default(getCall(object))
    cat("\nPortfolio Weight(s):\n")
    weights = round(getWeights(object), digits = 4)
    if (length(weights) == 1) {
        cat(" ", weights, "\n")
    }
    else {
        print.table(weights)
    }
    cat("\nRiskBudget(s):\n")
    riskBudgets = round(getCovRiskBudgets(object), digits = 4)
    if (length(riskBudgets) == 1) {
        cat(" ", riskBudgets, "\n")
    }
    else {
        print.table(riskBudgets)
    }
    if (FALSE) {
        if (!is.na(getTailRiskBudgets(object))) {
            cat("\nRiskBudget(s):\n")
            riskBudgets = round(getTailRiskBudgets(object), digits = 4)
            if (length(riskBudgets) == 1) {
                cat(" ", riskBudgets, "\n")
            }
            else {
                print.table(riskBudgets)
            }
        }
    }
    targetReturn = object at portfolio$targetReturn
    targetRisk = object at portfolio$targetRisk
    spec = getSpec(object)
    cat("\nTarget Risk(s) and Return(s):\n")
    if (is.null(dim(targetReturn))) {
        targetReturn = matrix(targetReturn, nrow = 1)
        colnames(targetReturn) = getEstimator(spec)[1]
    }
    if (is.null(dim(targetRisk))) {
        targetRisk = matrix(targetRisk, nrow = length(targetRisk) )
        colnames(targetRisk) = getEstimator(spec)[2]
    }
    target = cbind(targetReturn, targetRisk)
    colnames(target) = c(colnames(targetReturn), colnames(targetRisk) )
    if (nrow(target) == 1) {
        print(target[1, ])
    }
    else {
        print(target)
    }
    cat("\nDescription:\n ")
    cat(getDescription(object), "\n")
    invisible(object)
}

setMethod("show", "fPORTFOLIO", show.fPORTFOLIO)

####################################################################################

.portfolioConstrainedMVFrontier <- function (data, spec, constraints)
{
    if (!inherits(data, "fPFOLIODATA"))
        data = portfolioData(data, spec)
    mu = getMu(data)
    Sigma = getSigma(data)
    nAssets = getNumberOfAssets(data)
    targetAlpha = getTargetAlpha(spec)
    nFrontierPoints = getNFrontierPoints(spec)
    targetReturn = targetRisk = targetWeights = error = NULL
    Spec = spec
    solver = spec at solver$solver
    Spec at portfolio$weights = rep(1/nAssets, nAssets)
    k = 0
    solverType = spec at solver$solver
    status = NULL
    for (nTargetReturn in seq(min(mu), max(mu), length = nFrontierPoints)) {
        k = k + 1
        setTargetReturn(Spec) <- nTargetReturn
        nextPortfolio = .efficientConstrainedMVPortfolio(data = data,
            spec = Spec, constraints = constraints)
        Spec at portfolio$weights = nextPortfolio at portfolio$weights
        targetReturn = rbind(targetReturn,
nextPortfolio at portfolio$targetReturn)
        targetRisk = rbind(targetRisk, nextPortfolio at portfolio$targetRisk)
        nextWeights = nextPortfolio at portfolio$weights
        names(nextWeights) = names(mu)
        status = c(status, nextPortfolio at portfolio$status)
        targetWeights = rbind(targetWeights, t(nextWeights))
    }
    Index = (1:length(status))[status == 0]
    weights = targetWeights
    colnames(weights) = names(mu)
    weights = weights[Index, ]
    DIM = dim(targetReturn)[2]
    targetReturn = targetReturn[Index, ]
    targetReturn = matrix(targetReturn, ncol = DIM)
    colnames(targetReturn) = getEstimator(spec)[1]
    targetRisk = targetRisk[Index, ]
    new("fPORTFOLIO", call = match.call(), data = list(data = data),
        spec = list(spec = spec), constraints = as.character(constraints),
        portfolio = list(weights = weights, targetReturn = targetReturn,
            targetRisk = targetRisk, targetAlpha = targetAlpha,
            status = status), title = "Constrained MV Frontier",
        description = .description())
}

####################################################################################

# You should be able to specify the data in this form:
mu <- c( 0.1, 0.08, 0.065)
sigma <- c( 0.18, 0.12, 0.09 )

correlationMatrix <- rbind( c( 1, 0.8, 0.9 ),
                                              c( 0.8, 1, 0.75),
                                              c( 0.9, 0.75, 1) )

covarianceMatrix <- cor2cov(correlationMatrix, sigma )


data = list( mu = mu, Sigma = covarianceMatrix )

# And then do the optimisation
frontier <- portfolioFrontier(data, spec = spec, constraints )