[R-SIG-Finance] Quadratic programming solve.QP's Lagrangians

[Ton Jean-Claude]

Hi I am currently playing a little with Markowitz's portfolio optimization problem and from my computations I get that if we start from the following problem :

$$    \underset{w}{\arg \max} \quad w'r - \frac{\gamma}{2} w' Q w \quad \text{s.t}  \sum_i w_i=1, \quad w_i \geq 0\quad
$$

Then $\gamma$ should be of the form :
$$\gamma = \frac{w'(r-\beta)-\alpha}{w' \Sigma w},$$

where $\alpha$ is the Lagrangian multiplier due to the sum of weights being equal to one and $\beta$ the Lagrangian multipliers due to the no short selling condition.
Using `solve.QP` from the package `quadprog` I have coded the mean-variance formulation. And the function `ComputeRho` which you may find attached below computes the $\gamma$ (Sorry for the bad notations)


My problem is the following
depending on the r.a (risk aversion that is $\gamma$ in my case), `ComputeRho` gets me the correct answer or not.

    rm(list=ls())
        ###################################### begin : packages      ####################################
        require(tseries)
        require(quantmod)
        require(quadprog)
        require(MASS)
        require(PerformanceAnalytics)
        require(corpcor)
        require(quantmod)
        require(parallel)
        require(ggplot2)
        #require(mvtnorm)
        ###################################### end : packages  ######################################
        LogRet <- function(x,lag=5){
                #require(quantmod)
                #return(na.omit(apply(x,2,function(x) Delt(x,k=lag))) )
                require(zoo)
                return(na.omit(diff(log(x),lag=lag)))
        }
        ##########################################################################
        # Input :        sol = solution from solve.qp
        #                lr = logreturns
        # Output :      rho
        ComputeRho<-function(sol,lr){
                lr<-as.matrix(lr)
                alpha<-sol$lagrange_mults[1] # Lagrangian due to normalizing constraint
                beta<-sol$lagrange_mults[2:length(sol$lagrange_mults)] # Lagrangians due to no short sell constraint

                weights<-sol$weights
                cov.mat<-cov(lr)
                # unconstrained form
                uncons <-(t(weights)%*%colMeans(lr))/(t(weights)%*%cov.mat%*%weights)
                # due to Lagrangian multipliers
                cons1 <- (t(weights)%*%beta)/(t(weights)%*%cov.mat%*%weights)
                cons2 <- alpha/(t(weights)%*%cov.mat%*%weights)
                cons<-cons1+cons2

                rho<-uncons+cons # lagrangians change sign...
                rho
        }


        ############################################################
        # Input :       data = usually logreturns
        #               r.aversion = risk aversion coefficient usually between 1 to 10
        # Output :      er = expected return
        #               sd = standard deviation
        #               weights = weights
        #               Lagrangians = Lagrange multipliers

        MeanVariance<-function(data,r.aversion=10,method="MeanVariance")
        {
                # compute probability beliefs mean and risk ( here sample covariance)
                er<-colMeans(data)
                cov.mat<-cov(data)
                # get length of data
                n_asset <- length(er)

                # quad problem
                Dmat <- r.aversion*as.matrix(cov.mat)
                dvec <- er
                meq <-  1

                bvec <- c(1, rep(0,n_asset))
                Amat <- matrix( c(rep(1,n_asset), diag(nrow=n_asset)), nrow=n_asset)

                sol<-solve.QP(Dmat, dvec,Amat, bvec=bvec, meq=meq)

                # fetch returns from solve.QP
                weights<-sol$solution
                exp.ret <- t(er)%*%weights
                std.dev <- sqrt(weights %*% cov.mat %*% weights)

                # return results
                ret <- list(er = as.vector(exp.ret),
                            sd = as.vector(std.dev),
                            weights = weights,
                            lagrange_mults=sol$Lagrangian,
                            weights2=sol$unconstrained.solution)

        }
        # load asset prices
        data <- read.table("~/data.csv", header=TRUE, quote="\"")
        # risk aversion
        r.a<- 17 # for values below 18 r.a is different from rho
        lr<-LogRet(as.matrix(data))

        sol<-MeanVariance(lr,r.aversion = r.a)
        rho<-ComputeRho(sol,lr)
        rho

        r.a2<- 18 # for values below 18 r.a is different from rho
        sol2<-MeanVariance(lr,r.aversion = r.a2)
        rho2<-ComputeRho(sol2,lr)
        rho2


I was able to cern the problem to the sign of the Lagranians (here `cons1` and `cons2` in `ComputeRho` ) which switch signs depending on the `r.a.` I choose.

Does someone knows how the Lagrangians are computed explicitly ? Since It seems that this is done in FORTRAN.

Best,
Jc

here is the data set
https://www.dropbox.com/s/ynlnmdgfm9w4rq0/data.csv?dl=0

	[[alternative HTML version deleted]]