[R] linear constraint optim with bounds/reparametrization

[Kahra Hannu]

> from Ingmar Visser:

>I would like to optimize a (log-)likelihood function subject to a number of
>linear constraints between parameters. These constraints are equality
>constraints of the form A%*%theta=c, ie (1,1) %*% 0.8,0.2)^t = 1 meaning
>that these parameters should sum to one. Moreover, there are bounds on the
>individual parameters, in most cases that I am considering parameters are
>bound between zero and one because they are probabilities.
>My problems/questions are the following:

>1) constrOptim does not work in this case because it only fits inequality
>constraints, ie A%*%theta > =  c
                          ---
I was struggling with the same problem a few weeks ago in the portfolio optimization context. You can impose equality constraints by using inequality constraints >= and <= simultaneously. See the example bellow. 

>As a result when providing starting values constrOptim exits with an error
>saying that the initial point is not feasible, which it isn't because it is
>not in the interior of the constraint space.

In constrOptim the feasible region is defined by ui%*%theta-ci >=0. My first attempt to obtain feasible starting values was based on solving for theta from ui%*%theta = ci. Some of the items in the right hand side of the feasible region are, however, very often very small negative numbers, and hence, theta is not feasible. Next, I started to increase ci by epsilon ("a slack variable") and checked if the result was feasible. The code is in the example bellow. If ui is not a square matrix, theta is obtained by simulation. It is helpfull to know the upper and lower limits of the theta vector. In my case they are often [0,1]. Usually only 2-3 simulations are required.

In the example, the weights (parameters) have limits [0,1] such that their sum is limited to unity, as in your case. See, how Amat and b0 are defined.

V <- matrix(c(0.12,0,0,0,0.21,0,0,0,0.28),3,3)              # variances
Cor <- matrix(c(1,0.25,0.25,0.25,1,0.45,0.05,0.45,1),3,3)   # correlations
sigma <- t(V)%*%Cor%*%V                                     # covariance matrix
ER <- c(0.07,0.12,0.18)                                     # expected returns
Dmat <- sigma                                               # adopted from solve.QP
dvec <- rep(0,3)                                            #    "      "     "
k <- 3                                                      # number of assets
reslow <- rep(0,k)                                          # lower limits for portfolio weights
reshigh <- rep(1,k)                                         # upper limits for portfolio weights
pm <- 0.10                                                  # target return                      
rfree <- 0.05                                               # risk-free return
risk.aversion <- 2.5                                        # risk aversion parameter
####### RISKLESS = FALSE; SHORTS = TRUE ; CONSTRAINTS = TRUE ########################################
a1 <- rep(1, k)                
a2 <- as.vector(ER)+rfree  
a3 <- matrix(0,k,k)
diag(a3) <- 1    
Amat <- t(rbind(a1, a2, a3, -a3))
b0 <- c(1,pm,reslow, -reshigh)

objectf.mean <- function(x)                                                 # object function
{
return(sum(t(x)%*%ER-1/2*risk.aversion*t(x)%*%sigma%*%x)-(t(x)%*%ER-pm)) 
}

# Getting starting values for constrOptim

        ui <- t(Amat)                         
        ci <- b0
        if (dim(ui)[1] == dim(ui)[2])     
        {
        test <- F
        cieps <- ci
            while(test==F) 
            {
            theta <- qr.solve(ui,cieps)
            foo <- (ui%*%theta-ci)              # check if the initial values are in the feasible area
            test <- all(foo > 0)
            if(test==T) initial <- theta        # initial values for portfolio weights
            cieps <- cieps+0.1
            }
         } 
         if (dim(ui)[1] != dim(ui)[2])          # if Amat is not square, simulate the starting values        
         {
        test <- F
        i <- 1
        while(test==F)
            {
            theta <- runif(k, min = 0, max = 1)
            foo <- (ui%*%theta-ci)
            test <- all(foo > 0)                # and check that theta is feasible
            if (test == T) initial <- theta
            i <- i+1
            }    
        }

initial

res <- constrOptim(initial, objectf.mean, NULL, ui=t(Amat), ci=b0, control = list(fnscale=-1))
res$par                                     # portfolio weights (parameters)                    
y <- t(res$par)%*%ER                        # return on the portfolio
y



I hope this helps.

Hannu Kahra 
Progetti Speciali 
Monte Paschi Asset Management SGR S.p.A. 
Via San Vittore, 37 - 20123 Milano, Italia 
Tel.: +39 02 43828 754 
Mobile: +39 333 876 1558 
Fax: +39 02 43828 247 
E-mail: kahra at mpsgr.it 
Web: www.mpsam.it 

"Le informazioni trasmesse sono da intendersi per esclusivo uso del destinatario e possono contenere informazioni e materiale confidenziale e privilegiato. Qualsiasi correzione, inoltro e divulgazione in qualsiasi forma e modo anche a tenore generale è del tutto proibita. Se avete ricevuto per errore il presente messaggio, cortesemente contattate subito il mittente e cancellate da qualsiasi supporto il messaggio e gli allegati a voi giunti. Tutti gli usi illegali saranno perseguiti penalmente e civilmente da Monte Paschi Asset Management SGR S.p.A."

"The information transmitted are intended only for the person or entity to wich it is addressed and may contain confidential and/or privileged material. Any review, retransmission, dissemination, taking of any action in reliance upon, or other general use are strictly prohibited. If you received this in error, please contact immediately the sender and delete the material from any computer. All the illegal uses will be persecuted by Monte Paschi Asset Management SGR S.p.A."