@ This program replicates the Monte Carlo Experiment design of section    @
@ 6.1 in Golan, Judge, & Perloff JASA, 1996 article.                      @
@                                                                         @
library  maxlik,co;
#include maxlik.ext;
#include co.ext;

coset;
maxset;
begprog     = hsec;

__title     = "GME Estimation";

@                                                                           @
@       Set Simulation Parameters                                           @
@                                                                           @
iterate = 2;
n       = 1000;
@                                                                           @
@       Defining the Error Support Vector                                   @
@                                                                           @
NV      = 9;
bhat    = zeros(iterate,16);
retcode = 0;
eps     = 2.5*(10^(-16));
VMAX    = 1/(alpha*SQRT(N));
vstep   = vmax*2/(nv-1);
jj      = 1;
v       = zeros(nv,1);
DO      until jj > nv;
        v[jj,1] = -vmax + vstep*(jj-1);
        jj = jj+1;
endo    ;
data    = zeros(n,3);
@                                                                           @
@       Set Multinomial Parameters                                          @
@                                                                           @
stvall  = ones(8,1);         @ Starting Values for ML MNL  @
stvalg  = ones(8,1);         @ Starting Values for GME     @
dataset = "logitdta";
b1      = {0, 0, 0, 0};
b2      = {-1, 1, 2, -1};
b3      = {1, -1, -2, 1};
x       = ones(n,1)~rndn(n,3);
z0      = x[.,1];
z0      = z0|z0|z0;
z1      = x[.,2];
z1      = z1|z1|z1;
z2      = x[.,3];
z2      = z2|z2|z2;
z3      = x[.,4];
z3      = z3|z3|z3;
yak1    = x*b1;
yak2    = x*b2;
yak3    = x*b3;

@                                                                           @
@       Define Logit Likelihood Procedure                                   @
@                                                                           @
proc llogit(bl,dta);
    local   xb,llik,y,x,bhat,xb0,xb1,xb2;
    y       = dta[.,1];
    x       = dta[.,2:5];
    xb0     = x*zeros(4,1);
    xb1     = x*bl[1:4];
    xb2     = x*bl[5:8];
    xb       = xb0[1:n,.]|xb1[n+1:2*n,.]|xb2[2*n+1:3*n,.];
    llik    = -ln(1+exp(xb0)+exp(xb1)+exp(xb2))+y.*xb;
    retp    (llik);
endp;

@                                                                           @
@       Define GME Procedure                                                @
@                                                                           @
proc fctgme(bgme);
    local   omega,lambda1,lambda2,pnorm,M,psi1,psi2,J,I,xx,y,psi,psi3,
             lambda0;
    y       = data[.,1];
    XX      = data[1:n,2:5];
    M       = rows(v);
    J       = 3;
    I       = eye(J);
    lambda0 = zeros(4,1);
    lambda1 = bgme[1:4];
    lambda2 = bgme[5:8];
    omega   = (ln(exp(-XX*lambda1) +exp(-XX*lambda2)
                 +exp(-XX*lambda0)))'*ones(n,1);
    psi1    = ln(exp(-XX*lambda1*v'*ones(M,1)+eps))'*ones(n,1);
    psi2    = ln(exp(-XX*lambda2*v'*ones(M,1)+eps))'*ones(n,1);
    psi3    = (ln(M*ones(rows(y),1)))'*ones(rows(y),1);
    psi     = psi1+psi2 +psi3;
    pnorm   = lambda0|lambda1|lambda2;
    retp    (y'*(I.*.XX)*pnorm + omega + psi);
endp;

@                                                                           @
@       Outer Loop on Simulations (indexed by ct)                           @
@                                                                           @
ct      = 1;
totit   = 0;
do      while ct < (iterate+1);

    @                                                                       @
    @   Construct Data for current interation                               @
    @
                                                                            @
    r1    = rndu(n,3);
    r2    = rndu(n,3);
    sigma = 1;
    nu1   = -ln(-ln(r1[.,1]))~-ln(-ln(r1[.,2]))~-ln(-ln(r1[.,3]));
    nu2   = -ln(-ln(r2[.,1]))~-ln(-ln(r2[.,2]))~-ln(-ln(r2[.,3]));
    err1   = sigma*(nu1[.,1] - nu2[.,1]);
    err2   = sigma*(nu1[.,2] - nu2[.,2]);
    err3   = sigma*(nu1[.,3] - nu2[.,3]);
    wtp1   = yak1 + err1;
    wtp2   = yak2 + err2;
    wtp3   = yak3 + err3;
    Yes1    = (wtp1 .> wtp2).*(wtp1.>wtp3);
    Yes2    = (wtp2 .> wtp1).*(wtp2.>wtp3);
    Yes3    = 1-Yes1-Yes2;
    YES     = Yes1|Yes2|Yes3;
    Z       = x|x|x;
    let dep = YES;
    let ind = z0 z1 z2 z3;
    vnames  = { YES z0 z1 z2 z3 };
    vars    = dep|ind;
    data    = YES~Z;
    call    saved(data,dataset,vnames);

    @                                                                       @
    @   Run Logit and GME for current interation                            @
    @                                                                       @
    {bl,logl,g,vc,retlgt} = maxprt(maxlik(dataset,vars,&llogit,stvall));
    {bgme,f,g,retgme}     = co(&fctgme,stvalg);

    @                                                                       @
    @   Store results from current interation if convergence achieved       @
    @                                                                       @
    if  ((retlgt<1) and (retgme<1));
        bhat[ct,1:8]  = bl';
        bhat[ct,9:16]  = -bgme';
        ct          = ct+1;
    endif;
   totit=totit+1;
 endo;

@                                                                           @
@       Compute average simulation results                                  @
@                                                                           @
avgbs   = meanc(bhat)';
sdevb   = stdc(bhat)' ;
output file = c:\gauss\entropy\cv\golan.txt reset;
print   "logistic-distributed errors.  ";
print   "AVERAGE Bs for Logit";
print   avgbs[1:8];
print   "AVERAGE Bs for GME";
print   avgbs[9:16];
print   "Average Sbs in the same order";
print   sdevb[1:8];
print   sdevb[9:16];
print   "total number of iterations";
print   totit;
progtime = (hsec - begprog)/100;
print   "Program Elapse Time in Seconds";
print   progtime;
print "error support vectors";
print v';
end;