Numerical derivatives with deriv() when using factor variables

[Harald Tauchmann]

Dear all,
We are writing a program that involves calculating second derivatives numerically using mata's deriv() function. We want to allow factor variables to enter the function for which the second derivatives are calculated. Creating a mata matrix from stata factor variables means that this matrix will contain columns of zeros. In this case, deriv() aborts with the error message "could not calculate numerical derivatives -- flat or discontinuous region encountered", which makes perfect sense, since the function needs to be flat with respect to the coefficients of the omitted (base level) variables. See the code below for an illustration of this problem using a super-simple probit example..

clear
set seed 152507
set obs 99
gen byte x = rbinomial(1,0.5)
gen byte y = rbinomial(1,normal(-0.5+x))
gen byte one = 1

mata : XF = st_data(.,"i.x one")
mata : X = st_data(.,"x one")
mata : Y = st_data(.,"y")

mata
: void llprobit_v(c, RHS, LHS, lnf)
{
 lnf = ln(normal(RHS*c')) :* (LHS) :+ ln(normal(-RHS*c')) :* (1 :- LHS)
 }
end

probit y i.x

mata
 D = deriv_init()
 deriv_init_evaluator(D, &llprobit_v())
 deriv_init_evaluatortype(D, "v")
 deriv_init_params(D, st_matrix("e(b)"))
 deriv_init_argument(D, 1, XF)  
 deriv_init_argument(D, 2, Y)  
 deriv(D, 2)
end

However, one can make deriv() do what one want it to do, i.e. set the (second) derivatives with respect to the coefficients of the omitted variables to zero, by specifying deriv_init_search(D, "off"), i.e. preventing deriv() from searching for optimal delta values. In the "playing around" example above, this works well, i.e. the derivatives of the other coefficients are the same as those obtained by not specifying x as a factor variable and not specifying deriv_init_search(D, "off").

mata
 D = deriv_init()
 deriv_init_evaluator(D, &llprobit_v())
 deriv_init_evaluatortype(D, "v")
 deriv_init_search(D, "off")
 deriv_init_params(D, st_matrix("e(b)"))
 deriv_init_argument(D, 1, XF)  
 deriv_init_argument(D, 2, Y)  
 deriv(D, 2)
end

However – and not surprisingly - in more complex applications, specifying deriv_init_search(D, "off") seems to cause deriv() to calculate some rather 'strange' values for the derivatives.
So my question is:is there another way to get deriv()to set the derivatives with respect to the coefficients of omitted variable coefficients to zero?
(This question is obviously related to the problem of how to deal with factor variables when using optimize(). However, the solution of 'telling' optimize() by specifying optimize_init_constraints()that the respective coefficients are constrained to zero cannot, as far as I got it, be applied to deriv()).

Thank you all for your suggestions,
Harald