R help - Aug 2004 - linear constraint optim with bounds/reparametrization

Hello All,

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
                          ---
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.

Is there an alternative to constrOptim that can handle such strict
(equality) linear constraints?

2) Another option of course would be to reparametrize the problem as
follows. I will illustrate with an example:

I have parameters:
> p
      [,1]
 [1,]  0.8
 [2,]  0.2
 [3,]  0.2
 [4,]  0.8
 [5,]  0.6
 [6,]  0.1
 [7,]  0.3
 [8,]  0.1
 [9,]  0.3
[10,]  0.6
[11,]  0.5
[12,]  0.5

and the following constraints (all these constraint amount to certain
probabilities summing to one, ie the first two parameters sum to one, the
second pair of parameters  sum to one etc):
> A
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
[1,]    1    1    0    0    0    0    0    0    0     0     0     0
[2,]    0    0    1    1    0    0    0    0    0     0     0     0
[3,]    0    0    0    0    1    1    1    0    0     0     0     0
[4,]    0    0    0    0    0    0    0    1    1     1     0     0
[5,]    0    0    0    0    0    0    0    0    0     0     1     1

and the bounds on  the constraints are:
> ci
      [,1]
 [1,]    1
 [2,]    1
 [3,]    1
 [4,]    1
 [5,]    1

The bounds on the parameters are all zero and one. In order to reparametrize
I could use z=ginv(b)%*%(p-ginv(A)%*%cc) with b=Null(t(A)) and optimize z
instead of p using p=ginv(a)%*%ci + b%*%z

My question is however what the bounds on z are?

In the above example z is:
> z
            [,1]
[1,] -0.23333333
[2,] -0.03333333
[3,] -0.57974349
[4,] -0.37974349
[5,] -0.07974349
[6,] -0.18856181
[7,] -0.18856181

which conforms to the constraints, so these values can be used as an intial
estimate. If I knew the bounds on z I could use optim with
method="L-BFGS".

I hope this problem is sufficiently clear to elicit response, if not please
let me know.

ingmar

ps: to reproduce above example use the following:

p=matrix(c(0.8, 0.2, 0.20, 0.8, 0.6, 0.1, 0.3, 0.1, 0.3, 0.6, 0.5,
0.5),nc=1)

A = matrix(c(1, 1, 0, 0, 0, 0, 0, 0, 0,  0,  0,  0
, 0, 0, 1, 1, 0, 0, 0, 0, 0,  0,  0,  0
, 0, 0, 0, 0, 1, 1, 1, 0, 0,  0,  0,  0
, 0, 0, 0, 0, 0, 0, 0, 1, 1,  1,  0,  0
, 0, 0, 0, 0, 0, 0, 0, 0, 0,  0,  1,  1),nc=12,byrow=T)

ci=matrix(rep(1,5),nc=1)

b=Null(t(A))

z=ginv(b)%*%(p-ginv(A)%*%cc)

pp=ginv(a)%*%ci + b%*%z

> 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."

On Mon, 9 Aug 2004, Kahra Hannu wrote:
> >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.
>
Ick. You do not want to use constrOptim for equality constraints.
constrOptim is a log-barrier interior-point method, meaning that it adds
a multiple of log(A%*%theta-c) to the objective function. This is a really
bad idea as a way of faking equality constraints.

Use Lagrange multipliers and optim.

	-thomas

On 8/9/04 4:52 PM, "Thomas Lumley" <tlumley at u.washington.edu>
wrote:
> On Mon, 9 Aug 2004, Kahra Hannu wrote:
> 
>>> 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.
>> 
> 
> Ick. You do not want to use constrOptim for equality constraints.
> constrOptim is a log-barrier interior-point method, meaning that it adds
> a multiple of log(A%*%theta-c) to the objective function. This is a really
> bad idea as a way of faking equality constraints.
> 
> Use Lagrange multipliers and optim.
Is there a package that does all that for me? Or is there example code that
does something similar?

ingmar

>From Spencer Graves:
>However, for an equality constraint, I've had good luck by with an
objective function that adds something like the
>following to my objective function:
constraintViolationPenalty*(A%*%theta-c)^2, where
"constraintViolationPenalty" is
>passed via "..." in a call to optim.
I applied Spencer's suggestion to a set of eight different constrained
portfolio optimization problems. It seems to give a usable practice to solve the
portfolio problem, when the QP optimizer is not applicable. After all, practical
portfolio management is more an art than a science.
>I may first run optim with a modest value for constraintViolationPenalty
then restart it with the output of the
>initial run as starting values and with a larger value for
constraintViolationPenalty.
I wrote a loop that starts with a small value for the penalty and stops when the
change of the function value, when increasing the penalty, is less than epsilon.
I found that epsilon = 1e-06 provides a reasonable accuracy with respect to
computational time.

Spencer, many thanks for your suggestion.

Hannu Kahra