PRJMultinomial - predicted probs for a multinomial logit

[TomDoan]

This computes the predicted probabilities and marginal effects for the choices in a multinomial logit model given an input set of explanatory variables.

Code: Select all

*
* @PRJMultinomial( options )
* computes predicted probabilities and marginal effects for choices in a
* multinomial logit. This should be only be used after doing DDV with
* TYPE=MULTINOMIAL.
*
* Options:
*  XVECTOR=vector of test values for the explanatory variables. Make sure you include
*    1 for the CONSTANT.
*  LABELS=VECTOR[STRINGS] with labels for the choices. These should correspond to the
*    numeric codings from smallest to largest.
*  TITLE=title for output
*
* Revision Schedule:
*  07/2008 Written by Tom Doan, Estima.
*
procedure prjmultinomial
option vector        xvector
option vect[strings] labels
option string        title
*
local real    p
local vect    ez pz dg m
local rect    betamat
local integer i j row
*
if %rows(%betasys)==0.or.%clock(%rows(%betasys),%nreg)<>%nreg {
   disp "@PRJMULTINOMIAL must follow DDV(TYPE=MULTINOMIAL)"
   return
}
*
if .not.%defined(xvector) {
   disp "Syntax: @PRJMULTINOMIAL(XVECTOR=xvector,other options)"
   return
}
*
* Reformat %betasys into an %NREG x (choices - 1) array
*
compute betamat=%vectorect(%betasys,%nreg)
*
* Compute the probabilities of the choices (other than 1)
*
compute ez=%exp(tr(xvector)*betamat)
compute pz=ez/(1+%sum(ez))
compute dg=betamat*pz
report(action=define)
compute row=0
if %defined(title)
   report(atrow=(row=row+1),atcol=1,span) title

do i=1,%rows(ez)+1
   *
   * Handle the first choice by taking the remainders
   *
   if i==1 {
      compute p=1-%sum(pz)
      compute m=-p*dg
   }
   else {
      compute p=pz(i-1)
      compute m=p*(%xcol(betamat,i-1)-dg)
   }
   compute row=row+1
   if %defined(labels)
      report(atrow=row,atcol=1) "Choice: "+labels(i)
   else
      report(atrow=row,atcol=1) "Choice "+i
   report(atrow=(row=row+1),atcol=1) "Predicted P" p
   do j=1,%nreg
      if %eqnreglabels(0)(j)=="Constant"
         next
      report(atrow=(row=row+1),atcol=1) "Marginal "+%eqnreglabels(0)(j) m(j)
   end do j
   compute row=row+1
end do i
report(action=show)
end