#include <oxstd.h>
#import <modelbase>

/*------------------ Probit : Modelbase -----------------*/
class Probit : Modelbase
{
    Probit();                                // constructor
    fProbit(const vP, const adFunc, const avScore, const amHess);
    InitPar();                     // initialize parameters
    DoEstimation(vPar);                // do the maximation
    Covar();          // compute variance-covariance matrix
};
/*----------------- END Probit : Modelbase --------------*/

Probit::Probit()
{
    Modelbase();                   // initialize base class
}

Probit::fProbit(const vP, const adFunc, const avScore, const amHess)
{
    decl prob = probn(m_mX * vP);
    decl tail = 1 - prob;

    adFunc[0] = double(
        meanc(m_mY .* log(prob) + (1-m_mY) .* log(tail)));

    if (avScore)
    {
        decl weight = (m_mY - prob) .* densn(m_mX * vP)
            ./ (prob .* tail);
        avScore[0] = meanc(weight .* m_mX)';
    }

return 1;                            // 1 indicates success
}

Probit::InitPar()
{
	Modelbase::InitPar();   // first call Modelbase version
	SetParCount(m_cX);
	SetPar(zeros(m_cX, 1));              // start from zero
    return TRUE;
}
Probit::DoEstimation(vStart)
{
    m_iResult = MaxBFGS(fProbit, &vStart, &m_dLogLik, 0, FALSE);
       // m_iResult has return code, vStart will have new params
    m_dLogLik *= m_cT;             // change to unscaled log-lik
    return {vStart, "BFGS", FALSE};       // three return values
}
Probit::Covar()
{
    decl mhess;
    if (Num2Derivative(fProbit, m_vPar, &mhess))
    {
        m_mCovar = -invert(mhess) / m_cT;
    }
}   
main()
{
    decl probitobj = new Probit();        // create object
                             // load data file into object
    probitobj.LoadIn7("data/finney.in7");
    probitobj.Deterministic(FALSE);     // create constant
                                    // Formulate the model
    probitobj.Select(Y_VAR, { "vaso",0,0 } );
    probitobj.Select(X_VAR, { "Constant",0,0,
        "Lrate",0,0, "Lvolume",0,0 } );
    probitobj.SetSelSample(-1, 1, -1, 1);   // full sample

    probitobj.Estimate();                      // estimate

    delete probitobj;
}
