#include <oxstd.h>
#import <database>
#import <maximize>


/*-------------------------- Probit : Database -----------------------------*/
class Probit : Database
{
    decl m_mY;                                /* dependent variable [cT][1] */
    decl m_mX;                          /* regressor data vector [cT][m_cX] */

    Probit();                                                /* constructor */
    Estimate();                                        /* do the estimation */
    fProbit(const vP, const adFunc, const avScore, const amHessian);  /* LL */
};
/*------------------------ END Probit : Database ---------------------------*/

/*-------------------------- Probit : Database -----------------------------*/
Probit::Probit()
{
    this.Database();            // intialize base class
    println("Probit example 3, object created on ", date());
}

Probit::fProbit(const vP, const adFunc, const avScore,
    const amHessian)
{
    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::Estimate()
{
    decl  ir, vp, dfunc, mhess, cn;

    decl mx = loadmat("data/finney.in7");

    m_mY = mx[][0];       // dependent variable: 0,1 dummy
    m_mX = 1 ~ mx[][3:4]; // regressors: 1, Lrate, Lvolume
    cn = rows(m_mY);
    delete mx;

    vp = <-0.465; 0.842; 1.439>;        // starting values

    MaxControl(-1, 1);             // print each iteration
    ir = MaxBFGS(fProbit, &vp, &dfunc, 0, FALSE);

    print("\n", MaxConvergenceMsg(ir),
        " using analytical derivatives",
        "\nLog-likelihood = ", dfunc * cn,
        "; parameters:", vp);

                  // if converged: compute standard errors
    if (ir == MAX_CONV || ir == MAX_WEAK_CONV)
    {
        if (Num2Derivative(fProbit, vp, &mhess))
        {
            mhess = -invert(mhess) / cn;
            print("standard errors:", sqrt(diagonal(mhess)'));
        }
    }
}

main()
{
    decl probitobj;

    // create an object of class Probit
    probitobj = new Probit();

    // load the data, estimate the model
    probitobj.Estimate();

    // done with object
    delete probitobj;
}
