*! version 1.0.0  09May2016
program define mynlprobit1, eclass sortpreserve
    version 14.1

    syntax varlist [if] [in] [,  noCONStant ]
    marksample touse

    gettoken depvar indeps : varlist

    tempname b V N rank

    mata: mywork("`depvar'", "`indeps'", "`touse'", "`constant'", ///
       "`b'", "`V'", "`N'", "`rank'" )

    if "`constant'" == "" {
        local indeps "`indeps' _cons"
    }
    matrix colnames `b' = `indeps'
    matrix colnames `V' = `indeps'
    matrix rownames `V' = `indeps'

    ereturn post `b' `V', esample(`touse')
    ereturn scalar N       = `N'
    ereturn scalar rank    = `rank'
    ereturn local  cmd     "mynlexp"

    ereturn display

end

mata:

void MYNLProbit(real scalar todo, real vector b,  ///
	real vector y, real matrix X,	 	  ///
	val, grad, hess)
{
	real vector  r, f, xb
	real matrix  df

	xb  = X*b'
	f   = normal(xb)
	r   = y - f
	val = -(r:^2)
	df  = normalden(xb):*X

	if (todo>=1) {
		grad = r:*df
	}
	if (todo==2) {
		hess = -1*quadcross(df, df)
	}
}

void mywork( string scalar depvar,  string scalar indeps,
             string scalar touse,   string scalar constant,
             string scalar bname,   string scalar Vname,
             string scalar nname,   string scalar rname)
{

    real vector y, b
    real matrix X, V
    real scalar n, p, ssr
    transmorphic S

    y = st_data(., depvar, touse)
    n = rows(y)
    X = st_data(., indeps, touse)
    if (constant == "") {
        X = X,J(n, 1, 1)
    }
    p = cols(X)

    S = optimize_init()
    optimize_init_argument(S, 1, y)
    optimize_init_argument(S, 2, X)
    optimize_init_evaluator(S, &MYNLProbit())
    optimize_init_params(S, J(1, p, .01))
    optimize_init_evaluatortype(S, "gf2")
    optimize_init_conv_vtol(S, 1e-10)
    b   = optimize(S)
    V   = invsym(-1*optimize_result_Hessian(S))
    ssr = (-1/(n-p))*optimize_result_value(S)
    V   = ssr*V

    st_matrix(bname, b)
    st_matrix(Vname, V)
    st_numscalar(nname, n)
    st_numscalar(rname, p)
}
end