> 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: [EMAIL PROTECTED] 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." ______________________________________________ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html