What is causing this Mata error: 3200 conformability error?

[Abdul Adam]

Dear all,
I am trying to translate a SAS code into Stata using Mata. Part of the code implements a probit model where the gradient and Hessian are calculated by hand. I am using Stata 14.2, and the following is my code used on the auto dataset. I know I shouldn't be hardcoding a command that is always available in Stata, but this is part of a bigger code that I want to change from SAS to Stata, and I am new to Mata.
"
Executing the code ends with "3200 conformability error". Help appreciated.

sysuse auto, clear
probit foreign price mpg weight

*De…fine y and x
generate cons = 1
local y foreign
local xlist price mpg weight cons

mata
st_view(y=., ., "`y'") // read in stata data to y and X
st_view(X=., ., tokens("`xlist'"))
b = J(cols(X),1,0) // compute starting values
n = rows(X)
iter = 1   // initialize number of iterations
cha = 1   // initialize stopping criterion
do   {

 //real vector  r, f, xb
 //real matrix  df

  xb  = X*b
  f   = normal(xb)
  d = normalden(xb)
  r1  = y - f
  r2  = 1 - f
  df  = normalden(xb):*X'
  grad_numer = df:*r1
  grad_denom = f:*r2
  grad = grad_numer/grad_denom // kx1 gradient vector
  hes  = (d:^2)*cross(X',*X)/grad_denom // negative of the kxk hessian matrix
  bold = b
  b = bold + cholinv(hes)*(grad)
  cha = (bold-b)'(bold-b)/(b'b)
  iter = iter + 1
} while (cha > 1e-16)  // end of iteration loops
end

The formulas for the gradient and Hessian are based on Wooldridge (2002, p. 393).
Wooldridge, J.M. (2002). Econometric Analysis of Cross Section and Panel Data. The MIT Press

[Scott Merryman]

I would suggest trying the code without the do while loop while you are still debugging the code.

Matrix f is 74x1

r2 = 1 - f

is the source of the first conformability error.

[Abdul Adam]

Dear Scott,
Thank you for your answer. After solving the issue you pointed out, I have another conformability issue: see the "hes" line: Although I know the Hessian must be nxn, clearly, I am missing something, as I get error: <istmt>: 3200 conformability error. Here is the new code:

sysuse auto, clear
probit foreign price mpg weight

*De…fine y and x
generate cons = 1
local y foreign
local xlist price mpg weight cons

mata
st_view(y=., ., "`y'") // read in stata data to y and X
st_view(X=., ., tokens("`xlist'"))
X = X, J(rows(X), 1, 1)
b = J(cols(X),1,0)
n = rows(X)
  xb  = X*b
  f   = normal(xb)
  d = normalden(xb)
  r1  = y :- normal(xb)
  r2  = 1 :- f
  d2  = normalden(xb):*X
  grad_numer = d2:*r1
  grad_denom = f:*r2
  grad = grad_numer:/grad_denom // kx1 gradient vector
  hes  = (d:^2):*cross(X,X):/grad_denom // negative of the kxk hessian matrix
  bold = b
  b = bold + cholinv(hes)*(grad)
end

Also after the

b = bold + cholinv(hes)*(grad)
I get:
_cholinv(): 3205 square matrix required

[Abdul Adam]

Just an update:
I checked and derived the gradient formula by hand, and it seems there is no problem with the formula; however, when I try to code in Mata, I still get the conformability error. I can see that I should not be multiplying matrices of different sizes; But I wonder how I can make this Mata code replicate the results from Stata's probit command:

---

clear all
sysuse auto, clear
probit foreign price mpg weight

*De…fine y and x
generate cons = 1
local y foreign
local xlist price mpg weight cons

mata
st_view(y=., ., "`y'") // read in stata data to y and X
st_view(X=., ., tokens("`xlist'"))
b = J(cols(X),1,0)
  xb  = X*b
  grad = normalden(xb)*X'(y :- normal(xb)):/normal(xb)*(1 - normal(xb)) // should return kx1 gradient vector
  grad 
  hes  = (normalden(xb):^2)*cross(X,X):/normal(xb)*(1 - normal(xb)) // should return negative of the kxk hessian matrix
  hes
end

[Christophe Kolodziejczyk]

several comments on your code

1. you include the constant two times. Once in your local xlist and second by adding a column of ones in X. X becomes then singular.
2. Your gradient should be kx1 (k=4) and it is 74x1. The formula should be

grad = quadcolsum(grad_numer:/grad_denom )

3. Now the Hessian. It is not conformable. X'X is 4x4 and grad_denom is 74x1. I took a look at the formulas. In your cross-product each element of X'X is weighted by the vector (d^2/f*r2). So I think it should be

hes  = cross(X, d:^2:/(f:*r2),X)

4. In your last statement you shoud transpose the gradient

 b = bold + cholinv(hes)*(grad)'

5. It's better (numerically speaking) to code

r2=normal(-xb)

instead of

r2=1:-f

6. In your original code try to use

cha =  quadcross(bold-b,bold-b)/quadcross(b,b)

instead. It works the way you have implemented it, but use cross() to compute cross-products whenever you can.
7. Finally, I'll recommend to use optimize() or better moptimize() instead of this algorithm. It's safer and a better implementation of the Newton-Raphson algorithm, although it is a good exercise to go through these different steps .

Changing steps 1-4 makes the code work and you should obtain the same parameters as Stata's probit.

[Abdul Adam]

Dear Christophe,
Thank you very much for your detailed answers. The code is now working. May I ask you (and others) two follow-up questions: 1) regarding standard errors, and 2) a follow-up question on the use of use optimize() or better moptimize() which you suggested?

/****************Begin**************/
sysuse auto, clear

*De…fine y and x
generate cons = 1
local y foreign
local xlist price mpg weight cons

mata
st_view(y=., ., "`y'") // read in stata data to y and X
st_view(X=., ., tokens("`xlist'"))
b = J(cols(X),1,0)
iter = 1   // initialize number of iterations
cha = 1   // initialize stopping criterion
do   {
  xb  = X*b
  f   = normal(xb)
  d = normalden(xb)
  r1  = y :- normal(xb)
  r2  = normal(-xb)
  d2  = normalden(xb):*X
  grad_numer = d2:*r1
  grad_denom = f:*r2
  grad = quadcolsum(grad_numer:/grad_denom) // kx1 gradient vector
  grad
  hes  = cross(X, d:^2:/(f:*r2),X) // negative of the kxk hessian matrix
  hes
  bold = b
  b = bold + cholinv(hes)*(grad)'
  cha =  quadcross(bold-b,bold-b)/quadcross(b,b)
    iter = iter + 1
} while (cha > 1e-16)  // end of iteration loops

/*Compute the variance-covariance matrix of b */
vb = cholinv(hes)
iter   // num iterations
cha   // stopping criterion
st_matrix("b",b') // pass results from Mata to Stata
st_matrix("V",vb) // pass results from Mata to Stata
st_matrix("hes",hes) // pass results from Mata to Stata
st_matrix("grad",grad) // pass results from Mata to Stata

end

/*save results for later use*/
matrix list b
matrix list grad
matrix list hes
matsave b, p("matrices") dropall replace
matsave grad, p("matrices") dropall replace
matsave hes, p("matrices") dropall replace

*Present results, nicely formatted using Stata command ereturn
matrix colnames b = `xlist'
matrix colnames V = `xlist'
matrix rownames V = `xlist'
ereturn post b V
ereturn display

//use probit for comparison
sysuse auto, clear
probit foreign price mpg weight

/****************End**************/

My questions:
1) I observe slight differences between the standard errors by the probit command and those I get from Mata calculations? Any idea for this slight discrepancy?

Second question: the reason I want to use Mata is because; I want to save betas, gradient, and hessians for later use in another model, hence my use of the matsave command above. Can I do the same in optimize() & moptimise()?

Once again thank you very much for your help.

[Christophe Kolodziejczyk]

The computation of the Hessian you use is the expected value of the Hessian (information matrix). According to Wooldridge the computation of the Hessian is more complicated hence the use of the expectation. There might still be some discrepancies because the algorithms are not completely the sames.

You can certainly retrieve the paramater vector, the gradient and the Hessian from moptimize(). See below for an example. See help moptimize() (and optimize()) for more details. Notice that if you use the observed information matrix (oim) with moptimize() you will get closer results to the ones obtained with your algorithm.

clear mata
mata:

// void myprobit(M, todo, b, lnf, g, H)
void myprobit(M, b, lnf)
{
        
        real colvector y, xb

        xb = moptimize_util_xb(M,b,1)
        y  = moptimize_util_depvar(M,1)
        
        // lnf = moptimize_util_sum(M,y:*ln(normal(-xb))+(1:-y):*ln(normal(xb)))
        lnf = y:*ln(normal(-xb))+(1:-y):*ln(normal(xb))

}


M = moptimize_init()
moptimize_init_evaluator(M,&myprobit())
moptimize_init_evaluatortype(M,"lf")
moptimize_init_depvar(M, 1,"foreign")
moptimize_init_eq_indepvars(M, 1,"price mpg weight")
moptimize_init_technique(M,"nr")
// moptimize_init_eq_coefs(M,1,J(1,4,0))
moptimize_init_vcetype(M,"oim")
moptimize(M)
moptimize_result_display(M)

b = moptimize_result_coefs(M)
g = moptimize_result_gradient(M) 
H = moptimize_result_Hessian(M)
b
g
H
end

probit foreign price mpg weight, vce(oim)

[Abdul Adam]

Dear @Christophe Kolodziejczyk,
Thank you very much. I will definitely have a look at help moptimize() (and optimize()). For now, I have no more questions. Thanks.