*
* Example 21.7.8 from pp 729-732
*
open data tablef21-2[1].txt
data(format=prn,org=columns) 1 840 mode ttme invc invt gc hinc psize
*
set indiv  = (t+3)/4
set which  = %clock(t,4)
set chosen = mode
*
set dair   = which==1
set dtrain = which==2
set dbus   = which==3
set airhinc = dair*hinc
*
* Conditional logit model
*
ddv(type=conditional,indiv=indiv,modes=which) chosen
# dair dtrain dbus gc ttme airhinc
compute covfull=%xx,betafull=%beta
*
* Conditional logit restricted to those not choosing air, to test
* for IIA. To make the Hausman test easier, the dair and airhinc
* are included in the regressors. They'll estimate zero with zero
* standard errors, since they'll have no effect on the likelihood
* for the restricted sample.
*
ddv(type=conditional,indiv=indiv,modes=which,smpl=which<>1) chosen
# dair dtrain dbus gc ttme airhinc
*
test(all,vector=betafull,covmat=%xx-covfull,title="Hausman Test for IIA Among Ground Modes")
*
* Nested logit
*
* Standard utility index
*
frml(reg,vector=ba) fa
# dair dbus dtrain gc ttme
*
* Additional shift terms for group choice
*
frml(reg,vector=bg) fg
# airhinc
*
* Set group categories
*
compute choices=4,groups=2
dec rect grouplook(groups,choices)
input grouplook
 1 0 0 0
 0 1 1 1
*
* frmllist(1) is for the index of the choice given the branch
* frmllist(2) is for the index of the branch. For a given individual, this should
* produce the same value for all choices within a branch.
*
dec vect[frml] frmllist(2)
frml frmllist(1) = fa
frml frmllist(2) = fg
*
* tau is the adjustment parameter for the groups
*
dec vect tau(groups)
*
* This evaluates the nested logit probability for the individual whose data begins
* at "time". As with the conditional logit, the data points for time through
* time+choices-1 give the attributes of the choices. This is written to look for
* the choice itself in the series called "chosen"
*
function SequentialLogit time
type real SequentialLogit
type integer time
local integer i j
declare vect fullprobs relprobs ones
local real ez totprobs
*
dim fullprobs(choices) relprobs(choices) ones(choices)
compute fullprobs=%const(0.0)
*
*  Loop over groups on the outside
*
compute total=0.0
do i=1,groups
   compute within=0.0
   compute relprobs=%const(-1.0)
   *
   * Loop over the choices that belong to group "i". Compute their relative probabilities.
   *
   do j=1,choices
      if grouplook(i,j) {
         compute ez=exp(frmllist(1)(time+j-1)+frmllist(2)(time+j-1))
         compute relprobs(j)=ez
         compute within=within+ez
      }
   end do j
   *
   * Divide by within to convert to relative probabilities, and multiply by within**tau(i)
   * to weight the branch. In the end, fullprobs will have to be divided by total to get
   * probabilities.
   *
   ewise fullprobs(j)=%if(relprobs(j)==-1,fullprobs(j),relprobs(j)*within**(tau(i)-1))
   compute total=total+within**tau(i)
end do i
*
*  Return the probability of the chosen value
*
do j=1,choices
   if chosen(time+j-1) {
      compute SequentialLogit=fullprobs(j)/total
      return
   }
end do j
end
*
frml NLogit = log(SequentialLogit(t))
*
nonlin ba bg tau
compute ba=%const(0.0),bg=%const(0.0),tau=%const(1.0)
maximize(title="Nested Logit Model",pmethod=simplex,piters=10,method=bfgs,smpl=which==1) NLogit