quadratic programming in Stata

[boriskaiser]

Hi everyone,

I want to find optimal weights to combine different forecasts by minimizing a quadratic form given some constraints. It's a quadratic programming problem of the form:

min (w' * S * w) subject to sum(w_i)=1 and 0<=w_i<=1 for all i ,

where w is a m-x-1-vector and w_i is an element of w. S is a m-x-m matrix.

I know that this could be solved in R using the quadprog package. However, my whole analysis is set up in Stata and I would prefer solving the problem in Stata. How can I do this?

Any pointers would be appreciated, thanks!

Boris

[Matthew J. Baker]

boriskaiser --

I'm not sure if it helps, but here is a little example that uses the optimize() function in Mata. I know one can do a lot better by specifying explicitly a quadratic problem, but hopefully this is a start. I think the idea is to work with a transformed version of the parameters so that they always lie in between 0 and 1, and sum to 1. Something like:

clear all
set seed 5150
mata:
S = 1 :+ rnormal(3,3,0,1)

void quadfun(todo, b, s1, s2, s3, v, g, H) {
    
    r1 = exp(b[1])/(1+exp(b[1]))
    r2 = (1-r1)*exp(b[2])/(1+exp(b[2]))
    r3 = 1 - r1 - (1 - r1)*r2
    
    w = r1 \ r2 \ r3

    S=s1, s2, s3
    
    v = w'S*w
}

Problem = optimize_init()
optimize_init_evaluator(Problem, &quadfun())
optimize_init_which(Problem, "min")
optimize_init_evaluatortype(Problem, "d0")
optimize_init_params(Problem, (.5, .5))
optimize_init_argument(Problem, 1, S[,1])
optimize_init_argument(Problem, 2, S[,2])
optimize_init_argument(Problem, 3, S[,3])
b = optimize(Problem)
    
r1 = exp(b[1])/(1+exp(b[1]))
r2 = (1-r1)*exp(b[2])/(1+exp(b[2]))
r3 = 1 - r1 - (1 - r1)*r2    
    
r1,r2,r3

end

The idea in my code is to let r1 be a fraction between zero one, r2 a fraction of the remaining fraction, and then r3 be whatever is left over. One problem that I have had in these situations is that sometimes the coefficients aren't that stable - they want to go to either zero or one. So, using an optimizer that does not require derivatives such as Nelder-Mead helps. In my example, things go through okay even though the second weight goes to zero, but that is something to watch out for.

Hope that helps!

Matt Baker

[boriskaiser]

Dear Matt

Thanks a lot for your post. That's actually a nice way of imposing the constraints. I hadn't thought of that.

The downside is perhaps that in my application m, the length of the vector w, is not known apriori and can be large (>100). In that case, it's difficult to implement the code you suggest.

I guess the only way is to implement a quadratic programming algorithm in Stata myself or to change software ...

[Christophe Kolodziejczyk]

I don't think the size of the w vector is an issue for optimize(). You have to adapt the first two lines of the evaluator in order to be able to work with more than two weights which sum to 1.

...
 w = exp(b)/(1:+quadsum(exp(b))

w = w,1-quadsum(w)

v = w*S*w'
...

The problem with this approach is that you cannot have solutions with either some of the weights eqal to 0 or one weight equal to 1. You would have to implement the algorithm yourself. I don't know if it is possible to adapt optimize() in a way that allows inequality constraints.