OxMetrics - Modèle logit conditionnel

Question

J'essaie de développer un code OxMetrics pour un modèle logit multinomial (ou conditionnel) (comme décrit dans http://data.princeton.edu/wws509/notes/c6s2.html). Je suis nouveau à OxMetrics et j'ai encore du mal à comprendre comment fonctionne la routine d'optimisation (MaxBFGS).

Toute aide est plus que bienvenue!

// 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)); 
} 

Answer 1

Vous devriez lire le MaxBFGS officiel de l'aide (appuyez sur "F1" sur la fonction) ou here.

Cette fonction est déclarée comme:

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

• argument de func: fonction pour maximiser
• argument de avP: entrée -> matrice de valeurs de départ, sortie -> coefficients finals. Vous devriez donc d'abord stocker tous vos arguments qui doivent être optimisés dans une seule matrice (votre variable sv).
• argument amInvHess: mis à 0 pour une estimation simple.
• argument fNumDer: mis à 0: pour les dérivées premières analytiques et mis à 1 pour utiliser les dérivées premières numériques.

La fonction qui doit être optimisé doit ont le prototype suivant:

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

comme une simple introduction, vous pouvez jeter un oeil à estimer un modèle probit l'exemple suivant (OxMetrics7\ox\samples\maximize\probit1.ox) via MaxBFGS - il est documenté dans la documentation officielle "Introduction à Ox -Jurgen A. Doornik et 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: vous pouvez utiliser des variables globales pour vos X,Y,N,T,J.. les variables et améliorer votre code (voir probit2.ox, probit3.ox ...)

Answer 2

Merci beaucoup, votre réponse a été très utile (mon cerveau a du mal à passer de R à Ox) Je ne sais pas si c'est le bon endroit pour publier un Ox pour un logit mixte (ou des paramètres aléatoires, logit, ou densité de noyau), mais cela peut être utile aux autres. Bien sûr, les idées/pensées/astuces pour améliorer le code suivant sont les bienvenues!

DECLARE FORMAT DU DATASET

decl N = 433;  // Number of respondents 
decl R = 100;  // Number of draws 
decl T = 8; // Number of tasks per respondent (!!! same for all respondents) 
decl J = 2; // Number of options per task (!!! same for all tasks) 
decl Y = <0>;  // choice variable 
decl X = <1;2;3;4;5;6>; // 6 param, including both fixed effects + mean of random effects 
decl RC = <2;3;4;5;6>; // 5 param corresponding to SD of random effects 
decl H, mY, mX, mRC; 

FONCTION DE PROBABILITÉ MIXTES COMPUTE LOGIT

fn_MXL(const beta, const obj, const grad, const hess) 
{ 
decl i, seqprb, num, den, sllik; 
seqprb = zeros(N,R); 
for(i = 0; i < R; ++i) 
{ 
num = reshape(exp(mX * beta[0:5][] + mRC .* (H[N*T*J*i:(N*T*J*(i+1))-1][]) * beta[6:][]), N*T, J); // !!! Edit "beta" accordingly 
den = sumr(num); 
seqprb[][i] = prodr(reshape(sumr(num ./ den .* mY), N, T)); 
} 
sllik = sumc(log(meanr(seqprb))); 
obj[0] = double(sllik); 
return 1; 
} 

CODAGE MODELE/ESTIMATIONS

main() 
{ 
decl data, i, id1, id2, Indiv, prime, Halton, sv, ir, obj; 
// 1. Select variables 
data = loadmat("C:/Users/.../choice_data.xls"); 
mY = reshape(data[][Y], N*T, J); 
mX = data[][X]; 
mRC = data[][RC]; 
delete data; 
// 2. Generate halton draws 
id1 = id2 = zeros(N,R); // Create ID variables for both respondents and draws to re-order the draws 
for(i = 0; i < N; ++i) 
{ 
id1[i][] = range(1,R); // ID for draws 
id2[i][] = i + 1; // ID for respondents 
} 
Indiv = reshape(id1, N*R, 1) ~ reshape(id2, N*R, 1); 
Halton = zeros(N*R, sizeof(RC)); 
prime = <2;3;5;7;11;13;17;19;23;29;31;37;41;43;47;53;59;61;67;71;73;79;83;89;97;101;103;107;109;113>; // List of prime numbers 
for(i = 0; i < sizeof(RC); ++i) 
{ 
Halton[][i] = haltonsequence(prime[i], N*R); // prime number ; number of draws 
} 
Halton = Indiv ~ Halton; 
H = zeros(rows(Halton)*T*J,columns(Halton));  // Duplicate draws (T*J) times to match structure of "mX" (!!! possible to run out of memory) 
for(i = 0; i < (T*J); ++i) 
{ 
H[rows(Halton)*i:(rows(Halton)*(i+1)-1)][] = Halton; 
} 
H = sortbyc(H, <0;1>)[][2:]; // Draws are organised as following: Draw > Indiv (e.g. first 5 rows would correspond to 1st draw for first 5 indiv) 
H = quann(H); // Transform halton draws into normal draws 
// 3. Estimation 
sv = <0;0;0;0;0;0;0;0;0;0;0>;  // Starting values for X and RC 
MaxControl(-1, 1, 1); 
ir = MaxBFGS(fn_MXL, &sv, &obj, 0, TRUE); 
print("\n", MaxConvergenceMsg(ir), " using numerical derivatives", "\nFunction value = ", obj, "; parameters:", sv); 
}