Imposing inequality constraints in mata optimize

[George D]

Hi all,

I posted this question under another thread and am moving it to a new thread to keep things organized. I am trying to impose a constraint using optimize_init_constraints(S,Cc) for the following optimization:

void mvopt(todo, p, fml, g, H) {
external m
external v
fml= 253*(p*m) - sqrt(253)*sqrt(p*v*p')
}

and I was wondering if it is possible to impose an inequality constraint. I know how to impose constraints such as P₁ + P₂ + P₃ =1 by using optimize_init_constraints(S,(1,1,1,1)) but I don't know how to do something like P₁ >= 0.

Thank you.
George

[Sergiy Radyakin]

George, if you want to minimize F(x) subject to P1>=0, then you can minimize F(x)+(P1<0)*10^200 unconditionally. See 'Penalty method'.

Best, Sergiy

[Matthew J. Baker]

George --

Another option to Sergiy's method - which is sometimes the only way one can impose really difficult constraints - might be to just build the inequality constraint into your code via variable transform, as one often does when estimating a standard deviation. In your example:

void mvopt(todo, p, fml, g, H) {
external m
external v
v1=exp(v)
fml=253*(p*m)-sqrt(253)*sqrt(p*v1*p')
}

The only problem is that if you are doing this as part of an estimation problem, you will have to transform standard errors appropriately. You can also impose an inequality constraint between two values using similar ideas:

void mvopt(todo, p, fml, g, H) {
external m
external v
external A
external B
v1=1/(1+exp(v))*A+exp(v)/(1+exp(v))*B
fml=253*(p*m)-sqrt(253)*sqrt(p*v1*p')
}

I know I've seen this in a Stata reference book or manual somewhere - perhaps in Gould, Pitblado, and Poi's Maximum Likelihood with Stata?

Best,

Matt

[George D]

Sergiy, can you explain why the penalty function is (P1<0)*10^200 ?

Matt, can you explain how you know to transform V? I am getting results that don't really make sense if I do that. Do I have to transform my results in p back into another form after the optimization?

Thanks guys.

[Matthew J. Baker]

George --

Sorry I should have said something about that. I think the usual way of reverse-transforming a transformed coefficient is the delta method. The Wikipedia link pretty much sums up the idea, which is to use a linear approximation of the transforming function around the estimated parameter values. That is, if h(p) is the parameter function, and the estimated parameters are P, then one uses h(P)+(p-P)J(P) as an approximation with J(P) the Jacobian. Then, the variance-covariance matrix in terms of the sought-after "original" parameter(s) is J(P)VJ(P)'.

As an example, in the following linear-regression model, I use p2=ln(s^2) instead of just using s directly for reasons of numerical stability:

clear all
sysuse auto
mata: 
function lregeval(M,todo,b,crit,s,H)
{
real colvector p1, p2
real colvector y1
    p1=moptimize_util_xb(M,b,1)
    p2=moptimize_util_xb(M,b,2)
    y1=moptimize_util_depvar(M,1)
    crit=-(y1:-p1)'*(y1:-p1)/(2*exp(p2))- ///
        rows(y1)/2*p2
}
M=moptimize_init()
moptimize_init_evaluator(M,&lregeval())
moptimize_init_evaluatortype(M,"d0")
moptimize_init_depvar(M,1,"mpg")
moptimize_init_eq_indepvars(M,1,"price weight displacement")
moptimize_init_eq_indepvars(M,2,"")
moptimize(M)
end

After estimation, I want to express results in terms of s, not p2, so I observe that s=exp(p2/2), and compute the Jacobian, which in this case is a diagonal matrix with ones on the diagonal for the regular parameters, and only has a change for the standard deviation entry:

mata:
bOrig=moptimize_result_coefs(M)
bNew=bOrig
bNew[5]=exp(bOrig[5]/2)
VOrig=moptimize_result_V(M)
VNew =VOrig
J=diag(1\1\1\1\ 1/2*exp(bOrig[5]/2))
VNew=J*VOrig*J'
end

Then, bNew holds the "actual" estimated coefficients, with s instead of ln(s^2), and VNew holds the variance matrix of the "actual" estimated parameters.

I hope I got that right, and that the explanation helps!

Matt