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

decl g_mY;                                  // global data
decl g_mX;                                  // global data


///////////////////////////////////////////////////////////////////
// Possible improvements:
// 1) Use log(tailn(g_mX * vP)) instead of log(1-prob)
//    in fProbit. This is slower, but a bit more accurate.
// 2) Use analytical first derivatives. This is a bit more robust.
//    Add numerical computation of standard errors.
// 3) Use analytical second derivatives and Newton instead
//    of Quasi-Newton (BFGS). Because the log-likelihood is concave,
//    we don't really need the robustness of BFGS.
// 4) Encapsulate in a class to avoid global variables.
// 5) General starting value routine, etc. etc.
//
// probit2.ox implements (2)
// probit3.ox implements (4)
// probit4.ox implements (4) in a more general way
///////////////////////////////////////////////////////////////////

fProbit(const vP, const adFunc, const avScore,
    const amHessian)
{
    decl prob = probn(g_mX * vP);   // vP is column vector

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

return 1;                           // 1 indicates success
}

main()
{
    decl vp, dfunc, ir;

    print("Probit example 1, run on ", date(), ".\n\n");

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

    g_mY = mx[][0];       // dependent variable: 0,1 dummy
    g_mX = 1 ~ mx[][3:4]; // regressors: 1, Lrate, Lvolume
    delete mx;

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

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

    print("\n", MaxConvergenceMsg(ir),
        " using numerical derivatives",
        "\nFunction value = ", dfunc * rows(g_mY),
		"; parameters:", vp);
}
