3 Nov 2005. Model estimation procedures that can be called by MAXLIKC

	==================================================

SORRY. THIS NEEDS SERIOUS WORK (it is woefully outdated).

	==================================================

These procedures are designed to be used with the GRBL maxlikc procedure.

The calling syntax for all these procedures is:

          {lnlik,agrad,dirmtx}=procname(b,narg)

where:
     B -- Kx1 coefficient vector.
  NARG -- flag dictating what to compute, where:
          0) Li,G and DirMtx=G'G
          1) Li, G and DirMtx=Hessian
          2) Hessian only
          3) Likelihood only
          4) G'G only
          5) G and Li only
          9) Return vector of likelihoods (1 row per observation) -- useful in conjunction with HESSP and GRADP
and
  lnlik -- scalar; the sum of loglikelihoods
  agrad  -- Kx1 vector; the sum of the score (gradient) matrix
 dirmtx -- KxK matrix; the hessian, or gradient outer product

Note that each procedure expects different globals (as noted below)


Also note that certain error conditions will throw an error message and stop execution


Todo: a slew of stuff.  





-------------------------------------------------------------------
MLE_PROBIT		 Probit Estimator.
-------------------------------------------------------------------

Required global matrices:

    XUSE  - matrix of independent variables
    YUSE - a vector of dependent variable (the response).
            YUSE should take values of 0 (NO) and 1 (YES).

Optional matrices:

      AUXUSE - auxillary variables  
      BAUX - fixed coefficients used with AUXUSE. If no AUXUSE vars, BAUX=0 is okay
      WUSE - weight (optional)

  Note: it is the caller's responsiblity to provide a 0/1 YUSE (and to
        check that YUSE is not all 0s or all 1s)


-------------------------------------------------------------------
MLE_LOGIT 		    LOGIT Estimator.
------------------------------------------------------------------- 

Required global matrices:

    XUSE  - NxK matrix of independent variables
    YUSE - Nx1vector of dependent variable (the response).
            YUSE should take values of 0 (NO) and 1 (YES).

Optional matrices:

    WUSE  - Nx1 matrix of weights. Or scalar 1, if no weighting
    AUXUSE - NxJ matrix of auxillary variables 
    BAUX - fixed coefficients used with AUXUSE. If no AUXUSE vars, BAUX=0 is okay


-------------------------------------------------------------------
MLE_MNL		 MNL, Non-nested (joint) estimator.
-------------------------------------------------------------------

Note that each observation consists of a sub-matrix with NALTS rows (NALTS= # of alternatives
available for this observation), with each row in the sub-matrix corresponding to an alternative.

NALTS need NOT be the same for each observation (the model can be unbalanced). 


Required global matrices:


    PANELS  -   Nx2 matrix of obseravtion ids.
 	        1 : first row for this observation
	        2 : last row for this observation

	    Thus 1+id[.,2]-id[.,1] is the number of alternatives per obseration

	    Note: totalts= sumc(1+id[.,3]-id[.,2])

    XUSE  - (N*totalts) x K matrix of independent variables 

	     The XUSE variables can contain "conditional" and "multinomial" variables
 
	    Conditional:
	        Values differ across alternatives. The coefficients on the XUSE variables
	        is the same for all alternatives.
  	        Conditional variables have an implicit normalization: one of the alternatives 
                can have its XUSE values set to 0, with the others values recentered 
	        around this "set to 0" value.     In other words, what drives the
 	        estimation is the differences between XUSE variables (across observations), not the
	        actual value.

	      Multinomial:
	        Values are the same across alternatives. The coefficients depend on the alternative.
	        A normalization is enforced, with the coefficient on one of  alternative dropped
		(implicitily set to 0).

	        Thus, if you include multinomial variables, they should have some "0" rows. For example,
		if there are 4 alternatives, for "multinomial" variable Z1,  observation i would look like:
			0      0     0
			Z1i    0     0
		        0    Z1i     0
			0      0   Z1i

		Note the first line is all 0s -- the first alternative is the one that is "normalized".

		Note that multinomial variables should only be used when the number of alternatives is balanced
	        (same set of alternatives for each observation).


    YUSE - (N*totalts) x 1 vector of dependent variable (the response).
	    Each obseration-specific "submatrix" of YUSE should be:
		A vector with NALTS-1 zeros, and 1 one.  The "one" value is the chosen alternative
	        A vector of percents; the sum adding to 1. This is a "shares" model
		A vector of counts (which can be fractional). This is an "aggregate" model.


Optional global matrices:

    WUSE  - (N*totalts) x 1 matrix of weights. 
             Or scalar 1, if no weighting (or don't specify)
	    Note that weights are analytically equivalent to using a counts type of YUSE. Thus, there is probably
	    little need to use weights in this MNL model.
            

    AUXUSE - (N*totalts) x J matrix of auxillary variables -- with fixed coefficient values

	    If AUXUSE variables aren't needed, set AUXUSE=0 (or don't specify).

	    Independent variables with fixed coefficients.  The AUX variables are
	    useful when using samples of the alternatives; such as inclusion of
	    a size factor when aggregating alternatives, a probability of selection
	    factor when using importance sampling of alternatives, or a sampling
	    proportion factor when choice-based samples are constructed
	    (see Ben-Akiva and Lerman, ch. 9).


    BAUX  - Coefficients used with AUXUSE. If AUXUSE is not specified, baux is ignored.


   BH_EXP  - controls whether to use expected value of BHHH matrix


-------------------------------------------------------------------
MLE_MNL2_A	Nested MNL Estimator. First stage estimator:
-------------------------------------------------------------------

Note that each observation consists of a sub-matrix with NALTS rows (NALTS= # of alternatives
available for this observation), with each row in the sub-matrix corresponding to an alternative.
Furthemore, each alternative beonds to on of NNESTS "nests".

NALTS and NNESTS need NOT be the same for each observation (the model can be unbalanced). 

For balanced designs, then number of nests and alternatives must be the same.
Furthermore, the number of alternatives per nest must be the same (within and
across observations)


Required global matrices:


    PANELS  -   Nx2 matrix of obseravtion ids.
 	        1 : first row for this observation
	        2 : last row for this observation

	    Thus 1+id[.,2]-id[.,1] is the number of alternatives per obseration

	    Note: totalts= sumc(1+id[.,3]-id[.,2])

   NEST_CT: Nx2 matrix of nest ionfo
     Column 
         1 :   first row, in NEST_IDS, of info on this observation's nests
	 2 :   # of nests for this observation
 
    NEST_IDS: N*B x 1 matrix of nest ideenetivers
       NEST_IDS is structured into multi-row blocks, with each block corresponding to a
       an observation (with a row in NEST_CT).
            The first several values in a block are the # of observations
            for each nest. In particular, 
		i) for the ith obseravation, lookup the value of the first column 
	           of the ith row of NEST_CT. Call this MM1.
		   Then, lookup the value of the 2nd column (of the ith row). Call this
		   MM2
	       ii) The MM2 to MM2+MM1-1 rows of NEST_IDS identify the number of alternatives
		   in the MM1 nests identified for this observation
	      iii) From MM2+MM1 are the identifiers for the nests. These identifiers
		   are pointers into the rows  of XUSE and YUSE.
        Example,N_NESTS=3

  	   NESTCT1
	   NESTCT2
	   NESTCT3
           nest1_1 pointer
		..
	   nest1_nestct1  pointer
           nest2_1 pointer
		..
	   nest2_nestct2  pointer
           nest3_1 pointer
		..
	   nest3_nestct3  pointer

      For example,if the nest1_2 pointer equals 37, then
	  "the 2nd alternative in nest 1 can be found in row 37 of the data"
      ... that is, the 37th row of the X, Y, etc matrices.

      Example:
        if observation 4 has 3 nests, with 2,3, and 2 alternatives,
        and the last row of obseration 3 information is as row 40 of nest_ids
        and the nest ids for these 7 alternatives are:
	   12 12 14 51 14 51 14   (two alts for nest id=12, 3 for nest id=14, two for nest id=14)
        and the first alternative of observation 4 is at row 31 of the dataset

        Nest_ct[4,.] = 41~3

        nest_ids[41]=2		# alts in nest1
        nest_ids[42]=3	        # alts in nest2
        nest_ids[43]=2          # alts in nest3
        nest_ids[44]=31     - row 31 of data is  in nest 1
        nest_ids[45]=32     - row 32 of data is in nest 1
        nest_ids[46]=33     - row 33 of data is in nest 2
        nest_ids[47]=35     - row 35 of data is in nest 2
        nest_ids[48]=37     - row 37 of data is in nest 2
        nest_ids[49]=34     - row 34 of data is in nest 3
        nest_ids[50]=36     - row 36 of data is in nest 3

      NOtes:
          * NEST_CT[5,1] will be 51.
          * The order of the "pointers", within a subset, doesn't matter.
	    For example:
                nest_ids[46]=33 
                nest_ids[47]=35 
                nest_ids[48]=37 
	    and
                nest_ids[46]=37 
                nest_ids[47]=33 
                nest_ids[48]=35 
	    are equivalent.



    XUSE  - (N*totalts) x K matrix of independent variables 
	    
            The XUSE variables an contain "conditional" and "multinomial" variables
	    "totalts" is the number of alternatives per observation (or the "average", if 
	    not balanced).

	     See "MNL, 1 stage, Estimator" for the details

    YUSE - (N*totalts) x 1 vector of dependent variable (the response).
	    Each obseration-specific "submatrix" of YUSE should be:
		A vector with NALTS-1 zeros, and 1 one.  The "one" value is the chosen alternative
	        A vector of percents; the sum adding to 1. This is a "shares" model
		A vector of counts (which can be fractional). This is an "aggregate" model.


Optional global matrices:

    WUSE  - (N*totalts) x 1 matrix of weights. 
             Or scalar 1, if no weighting (or don't specify)
	    Note that weights are analytically equivalent to using a counts type of YUSE. Thus, there is probably
	    little need to use weights in this MNL model.
            

    AUXUSE - (N*totalts) x J matrix of auxillary variables -- with fixed coefficient values

	    If AUXUSE variables aren't needed, set AUXUSE=0 (or don't specify).

	    Independent variables with fixed coefficients.  The AUX variables are
	    useful when using samples of the alternatives; such as inclusion of
	    a size factor when aggregating alternatives, a probability of selection
	    factor when using importance sampling of alternatives, or a sampling
	    proportion factor when choice-based samples are constructed
	    (see Ben-Akiva and Lerman, ch. 9).


    BAUX  - Coefficients used with AUXUSE. If AUXUSE is not specified, baux is ignored.


   BH_EXP  - controls whether to use expected value of BHHH matrix




-------------------------------------------------------------------
MLE_MNL2_I		Nested MNL Estimator. Compute inclusive values
-------------------------------------------------------------------

Note that each observation consists of a sub-matrix with NALTS rows (NALTS= # of alternatives
available for this observation), with each row in the sub-matrix corresponding to an alternative.
Furthemore, each alternative beonds to one of NNESTS "nests".
(see  the  "Nested MNL Estimator. First stage estimator" for the details)


B1 should be the output of a first-stage MNL estimator (say, MLE_MNL_a)

Returns a two  vector N * anests x 1 vectors
  Inclus: the inclusive values 
  Choices: # times this nest chosen (sumc of choices within the nest)
where anests is the (average) number of nests per observation.

Required global matrices:

  PANELS  -   Nx2 matrix of obseravtion ids.
 	        1 : first row for this observation
	        2 : last row for this observation

	    Thus 1+id[.,2]-id[.,1] is the number of alternatives per obseration

	    Note: totalts= sumc(1+id[.,3]-id[.,2])

   NEST_CT: Nx2 matrix of nest info
     Column 
         1 :   first row, in NEST_IDS, of info on this observation's nests
	 2 :   # of nests for this observation
 
   NEST_IDS: N*B x 1 matrix of nest identifiers
       NEST_IDS is structured into multi-row blocks, with each block corresponding to a
       an observation (with a row in NEST_CT).
       B is the (average) number of nests + number of alternatives per observation.
       (see  the  "Nested MNL Estimator. First stage estimator" for the details)

    XUSE,YUSE,AUXUSE,BAUX,WUSE
        First stage variables
        (see  the  "Nested MNL Estimator. First stage estimator" for the details)

-------------------------------------------------------------------
MLE_MNL2_B		Nested MNL Estimator. 2nd stage estimator:
-------------------------------------------------------------------


Note that each observation consists of a sub-matrix with NALTS rows (NALTS= # of alternatives
available for this observation), with each row in the sub-matrix corresponding to an alternative.
Furthemore, each alternative beonds to one of NNESTS "nests".
(see  the  "Nested MNL Estimator. First stage estimator" for the details)


Required global matrices:


  PANELS  -   Nx2 matrix of obseravtion ids.
 	        1 : first row for this observation
	        2 : last row for this observation

	    Thus 1+id[.,2]-id[.,1] is the number of alternatives per obseration

	    Note: totalts= sumc(1+id[.,3]-id[.,2])

   NEST_CT: Nx2 matrix of nest info
     Column 
         1 :   first row, in NEST_IDS, of info on this observation's nests
	 2 :   # of nests for this observation
 
   NEST_IDS: N*B x 1 matrix of nest identifiers
       NEST_IDS is structured into multi-row blocks, with each block corresponding to a
       an observation (with a row in NEST_CT).
       B is the (average) number of nests + number of alternatives per observation.
       (see  the  "Nested MNL Estimator. First stage estimator" for the details)


    YUSE2  -- N*nests x 1 vector of choice variables (indicating whether, or how many times, this
	     nest was chosen)
	     You can use MLE_MNL_I to create this "nest choice" vector

    XUSE2 -- N*nests x K2 matrix of independent variables
             Nests is the number of nests per observation (or the average number in unbalanced designs)
	
	     XUSE2 MUST CONTAIN THE NEST SPECIFIC INCLUSIVE VALUES.

	     You can use MLE_MNL_I to create a vector of inclusive values.

	     XUSE2 may also contain conditional and multinomial variables


 Optional:

    AUXUSE2 -- N*nests x  J2 matrix of auxillary variables
    BAUX2   -- J2 x 1 vector of fixed coefficient values


/***********************/
 Compute corrected nested mnl cv matrix
  This is a score matrix: is based on expected value of gradient.
  THus, actual choice is not used.

MLE_MNL2_F
 bb1: first stage beta
 bb2: second stage beta
 rho: coefficient on inclusive value (rho will also be one of the bb2 parameters).




/*******************/
Nested MNL fiml and liml ala Small and Brownstone.
Note that iarg is ignored (for now) -- impose a BHHH method.


mle_mnl2_fiml(b1,iarg);


