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

[Jeff Pitblado (StataCorp)]

Thank you for providing a working example.

deriv() does not handle constraints. optimize() and moptimize() provide their own evaluators to deriv(), and it is these evaluators that handle the constraints for deriv().

In the following, I show the changes one might employ to get deriv() to compute the second order numerical derivatives. The changes involve applying the contraints to the predictors matrix and point estimates that are supplied to deriv(), then reversing the transformations on the resulting derivatives matrix.

My changes take advantage of the fact that probit posts the constraints implied by the factor variables to e(Cns). Command matcprob uses this matrix to construct the tranformation matrix T and offset vector a. For more information on this, in Stata type help makecns.

clear all

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

probit y i.x

* get the matrices that apply the linear constraints
matcproc T a C
matrix list T
matrix list a
matrix list C

mata:

// evaluator is not changed
void llprobit_v(c, RHS, LHS, lnf)
{
    lnf = ln(normal(RHS*c')) :* (LHS) :+ ln(normal(-RHS*c')) :* (1 :- LHS)
}

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

// Mata copy of linear constraint tranformation matrices
T = st_matrix("T")
a = st_matrix("a")

// values at which to compute the numerical derivatives
b0 = st_matrix("e(b)") 

// transform the data and values at which to compute the numerical derivatives
tXF = XF*T
tb0 = (b0 :- a)*T

D = deriv_init()
deriv_init_evaluator(D, &llprobit_v())
deriv_init_evaluatortype(D, "v")
// supply transformed data
deriv_init_argument(D, 1, tXF)  
deriv_init_argument(D, 2, Y)  

// supply transformed values
deriv_init_params(D, tb0)

// compute second order derivatives
td2 = deriv(D, 2)
// show
td2
// invert
itd2 = invsym(-td2)

// transform back to unconstrained space
V = T*itd2*T'
V

// show that above are equivalent to the VCE from -probit-
st_matrix("e(V)")

end