options crt;
name nest31 'General 2-Level Nested Logit using ML Procedure' ;

? general 3-level nested logit program structure
? also handles special cases, like 2-level or standard mnl
? by restrictions on inclusive value coefficients

? 	par (Paul Ruud) 11/29/94
? 	based on a 2-level example by bhh (nlogit, 11/87, 2/93)
?	considerably cleaned up by clint 11/30/94
?	more comments by par 12/1/94
?	more automation by par 12/5/94 
?		(1) frml's for Lt, Bt., Tt..
?		(2) param's for tau., theta..

?============================================================

? using the nomenclature of Axel Boersch-Supan, L stands for Limbs, 
? B for Branches, and T for Twigs

? this particular version has the following structure
? for twelve alternatives

?                   |
?           +-------+-------+            (3 Limbs)
?           |       |       |
?           +   +---+---+   +-+          (1, 3, 2 Branches)
?           |   |   |   |   | | 
?         +-+ +-+ +-+ +-+-+ + +-+        (2, 2, 2, 3, 1, 2 Twigs)
?         | | | | | | | | | | | |
?         1 2 3 4 5 6 7 8 9 a b c

? to modify this code for a particular structure, one must
? 	(1) change the list statements below to specify no. of
?		options at each node
?	(2) impose normalizations
?		(a) tau's are coefficients for Limbs; tau1 for Limb 1,
?			tau2 for Limb 2, and so on
?		(b) theta's are coefficients for Branches; theta11 for
?			Limb 1 on Branch 1, theta12 for Limb 2 on Branch
?			1, and so on
?	(4) specify the deteriminstic utility functions xb... and create the
?		explanatory variables with the appropriate numerical
?		suffixes: e.g., DIST123 for characteristic DIST of
?		alternative on Limb1, Branch 2, Twig 3
?	(5) create a dependent variable called y containing the index of
?		the chosen alternative; be sure to define this variable
?		before the section of code that uses it (above the hist
?		commmand)

? The following is an alternative schematic for the tree pictured
? above.  This tree shows how the "subscripting" in the program works.
? It is necessary to adopt this subscripting in order to interpret
? coefficients and to correctly specify the list statements immediately
? below.

?    L      B       T       Alternative
? \_____1
?        \______11______111	1
?                 \_____112	2
?
? \_____2
?        \______21______211	3
?                 \_____212	4
?
?        \______22______221	5
?                 \_____222	6
?
?	 \______23______231	7
?                 \_____232	8
?                 \_____233	9
?
? \_____3
?        \______31______311	a
?
?        \______32______321	b
?                 \_____322	c

?                                   In this tree there are:
?  3 Limbs
?           1 Branch
?                   2 Twigs
?           3 Branches
?                   2 Twigs
?                   2 Twigs
?                   3 Twigs
?           2 Branches
?                   1 Twig
?                   2 Twigs


?                                  These translate into the 
?                                  the following list statements:
list listL 1-3 ;
 list list1B 1 ;                   ? *See NOTE
  list list11T 1-2 ;
 list list2B 1-3 ;
  list list21T 1-2 ;
  list list22T 1-2 ;
  list list23T 1-3 ;
 list list3B 1-2 ;
  list list31T 1 ;                 ? *See NOTE
  list list32T 1-2 ;

? *NOTE: Wherever there is a single branch or twig, the list
?        has a solitary "1".  The coefficients at these points
?        must be held constant in estimation.  Below you will
?        find "const tau1" and "const theta31" statements for
?        the two cases illustrated here.

?            DO NOT ALTER ANYTHING BELOW THIS LINE
?============================================================

const null 0 ;
frml Lt null ;
frml eq_llf -log(Lt) ;
dot listL ;
	frml Lo Lt + L. ; 
	eqsub Lo Lt ;
	frml Lt Lo ;
	eqsub Lt Lo ;

	frml llfe eq_llf + ll. ;
	eqsub llfe eq_llf ;
	frml eq_llf llfe ;
	eqsub eq_llf llfe ;

	frml Bt null ;
	dot list.B ;
		frml Bo Bt + B.. ;
		eqsub Bo Bt ;
		frml Bt Bo ;
		eqsub Bt Bo ;

		frml llfe eq_llf + ll.. ;
		eqsub llfe eq_llf ;
		frml eq_llf llfe ;
		eqsub eq_llf llfe ;

		frml Tt null ;
		dot list..T ;
			frml To Tt + T... ;
			eqsub To Tt ;
			frml Tt To ;
			eqsub Tt To ;

			frml llfe eq_llf + ll... ;
			eqsub llfe eq_llf ;
			frml eq_llf llfe ;
			eqsub eq_llf llfe ;
		enddot ;
		frml Tt.. Tt ;
		eqsub Tt.. Tt ;

		param theta.. 1 ;
	enddot ;
	frml Bt. Bt ;
	eqsub Bt. Bt ;

	param tau. 1 ;
enddot;
eqsub eq_llf Lt ;


dot listL ;
	frml ll. y.*(log(L.) - log(Bt.));
	frml L.  exp(w.*tau.) ;
	frml w. log(Bt.) ;
	eqsub eq_llf ll. L. w. Bt. ;
	frml tau tau. ;
	dot list.B ;
		frml ll.. y..*(log(B..) - log(Tt..)) ;
		frml B.. exp(v..*theta../tau) ;
		frml v.. log(Tt..) ;
		eqsub eq_llf ll.. B.. v.. Tt.. ;
		frml theta theta.. ;
		dot list..T ;
			frml ll... y...*log(T...) ;
			frml T... exp(xb.../theta) ;
			frml xb... x1...*b1 + x2...*b2 ; ? utility index
			eqsub eq_llf ll... T... xb... ;
		enddot ;
		eqsub eq_llf theta ;
	enddot ;
	eqsub eq_llf tau ;
enddot ;

? print eq_llf ; ? this will print out the LLF in all its glory

?============================================================
?            DO NOT ALTER ANYTHING ABOVE THIS LINE

? generate a test data set:  users should place their code
? here for their own data set and application

?	(1) create a variable y, which contains the index of
?		the chosen alternative
?	(2) create the variables that hold the characteristics
?		of the alternatives

? we will just create:
?	normally distributed characteristics
?	error components according to structure of the tree
?	coefficients equal magnitude and opposite sign

random(seedin=498724) ;
smpl 1 1000 ;
genr vmax = -1e30 ;
set alt = 0 ; genr y = 0 ;
dot listL ;
 random u ; ? Limb error component
 dot list.B ;
  random uu ; ? Branch error component
  dot list..T ;
	set alt = alt + 1;
	random x1... ;
        random x2... ;
	random uuu ; ? Twig error component
        genr v... = 3*(x1... - x2...) + u + 2*uu + 2*uuu ;
	genr test = vmax>=v... ;
	genr vmax = vmax*test + v...*(1-test) ;
	genr y = y*test + alt*(1-test) ;
  enddot ;
 enddot ;
enddot ;

delete u u uuu ; compress ;

hist (discrete) y ;

?            DO NOT ALTER ANYTHING BELOW THIS LINE
?============================================================

? generate the node switches:  y contains index of chosen
? alternative (in this example, 1,2,3,...,12); this section
? generates variables that describe the path through the tree
? structure to the chosen alternative

set alt = 0 ;
dot listL ;
 genr u = 0 ;
 dot list.B ;
  genr uu = 0 ;
  dot list..T ;
	set alt = alt + 1;
	genr y... = (y=alt) ;
	genr uu = uu + y... ;
  enddot ;
  genr y.. = uu ;
  genr u = u + uu ;
 enddot ;
 genr y. = u ;
enddot ;

?============================================================
?            DO NOT ALTER ANYTHING ABOVE THIS LINE

? try it out
? 	starting value strategy: begin with standard mnl 
?	estimates and then work your way down through the
?	tree in Gauss-Seidel fashion.  One may need several
?       cycles of Gauss-Seidel, and one may need to truncate
?	the number of iterations on the theta and tau
?	estimation steps.

const tau1 1 ;
const theta31 1 ;

supres @coef ; ? lower output volume

do i = 1 to 1 ; ? increase the number of cycles, if desired
	param b1 b2 ;
	const tau2 tau3 ;
	const theta11 theta21 theta22 theta23 theta32 ; 
	ml (silent) eq_llf ;

	param theta11 theta21 theta22 theta23 theta32 ;
	const b1 b2 ;
	ml (silent,maxit=4) eq_llf ;

	param tau2 tau3 ;
	const theta11 theta21 theta22 theta23 theta32 ;
	ml (silent,maxit=4) eq_llf ;
enddo ;

nosupres @coef ;

? fully unconstrained optimization:

param b1 b2 ;
param theta11 theta21 theta22 theta23 theta32 ;
param tau2 tau3 ;
ml eq_llf ;