OxMetrics - Conditional logit model

I am trying to develop a OxMetrics code for a multinomial (or conditional) logit model (as described in http://data.princeton.edu/wws509/notes/c6s2.html). I am new to OxMetrics and still struggle to understand how the optimisation routine (MaxBFGS) works.

Any help is more than welcome!

// FUNCTION FOR MULTINOMIAL LOGIT MODELLING
//X = Independent variables (e.g. product features)
//Y = Dependent variable (i.e. discrete choices (0/1))
//N = Total number of individuals (N=500)
//T = Number of observations (choice tasks) per respondent (T=20)
//J = Number of products per task (J=2 => choice between two options)
//sv => starting values
//llik => model log-likelihood

LOGIT(const sv, const llik, X, Y, N, T, J)
{
decl num, den, prbj, prbi, maty;
num = reshape(exp(X*sv[0:6]), N*T, J);
den = sumr(num);
prbj = num ./ den;
maty = reshape(Y .== 1, N*T, J);
prbi = sumr(prbj .* maty);
llik[0] = -sumc(log(prbi));
return llik[0];
}

main()
{
decl data, N, T, J, X, Y, sv, dFunc, fit;
data = loadmat("C:/Users/.../Data/data.csv");
X = data[][33] ~ data[][5:10];
Y = data[][12];
N = 500;
T = 20;
J = 2;
sv = <0;0;0;0;0;0;0>;
println ("\nEstimating using MaxBFGS");
fit = MaxBFGS(LOGIT, X=X, Y=Y, N=N, T=T, J=J, &sv, &dFunc, 0, TRUE);
println (MaxConvergenceMsg(fit), " at parameters ", sv', "with function value ", double(dFunc));
}

Solution

You should read the official MaxBFGS help (press "F1" on the function ) or here.

This function is declared as :

MaxBFGS(const func, const avP, const adFunc, const amInvHess, const fNumDer);

argument func : function to maximize
argument avP : input--> matrix of starting values, output -> final coefficients. So you should first store all your arguments that has to be optimized in a single matrix (your sv variable) .
argument amInvHess : set to 0 to simple estimation.
argument fNumDer : set to 0:for analytical first derivatives and set to 1 to use numerical first derivatives.

The function that has to be optimized must have the following prototype :

LOGIT(const vP, const adFunc, const avScore, const amHessian)

As a simple introduction you can have a look to the the following example (OxMetrics7\ox\samples\maximize\probit1.ox) to estimate a probit model via MaxBFGS - it is documented in the official documentation "Introduction to Ox -Jurgen A. Doornik and Marius Ooms -2006").

// file : ...\OxMetrics7\ox\samples\maximize\probit1.ox
#include <oxstd.oxh>
#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);
}

Ps: you can use global variables for your X,Y,N,T,J.. variables and then improve your code (see probit2.ox, probit3.ox...)