*! version 1.0.0  31Jan2016
program define mypoisson1, eclass sortpreserve
    version 14.1

    syntax varlist(numeric ts min=2) [if] [in] [, noCONStant ]
    marksample touse

    gettoken depvar indepvars : varlist

    _rmcoll `indepvars', `constant' forcedrop
    local indepvars  "`r(varlist)'"

    tempname b V N rank

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

    if "`constant'" == "" {
        local cnames "`indepvars' _cons"
    }
    else {
        local cnames "`indepvars'"
    }
    matrix colnames `b' = `cnames'
    matrix colnames `V' = `cnames'
    matrix rownames `V' = `cnames'

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

    ereturn display

end

mata:

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

    real vector y, b, mo_v, cv
    real matrix X, V, Cm
    real scalar n, p, rank, ko

    y = st_data(., depvar, touse)
    n = rows(y)
    X = st_data(., indepvars, 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, &plleval2())
    optimize_init_params(S, J(1, p, .01))

    b    = optimize(S)
    V    = optimize_result_V_oim(S)
    rank = p - diag0cnt(invsym(V))

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

void plleval2(real scalar todo, real vector b,     ///
              real vector y,    real matrix X,     ///
              val, grad, hess)
{
    real vector  xb

    xb = X*b'
   val = sum(-exp(xb) + y:*xb - lnfactorial(y))
}

end