*
*  MiscProb
*  Assorted probit/logit
*
all 400
*
*   Generate multinomial logit data
*
set x = %rangamma(10.0)
set i1 =  .3 + .07*x
set i2 = -.5 - .06*x
set norm = 1+exp(i1)+exp(i2)
set y = (u=%uniform(0.0,norm)),%if(u<exp(i1),1,%if(u<exp(i1)+exp(i2),2,0))
*
*   Estimate multinomial logit
*
nonlin b11 b12 b21 b22
compute b11=b12=b21=b22=0.0
frml index1 = b11 + b12*x
frml index2 = b21 + b22*x
frml multilgt = (f1=index1), (f2=index2), $
       %if(y==2,f2,%if(y==1,f1,0.0))-log(1.0+exp(f1)+exp(f2))
maximize(method=bfgs) multilgt
*
*   Generate sequential probit data
*
set x1 = %rangamma(10.0)
set x2 = %ran(20.0)
set i13 = -.5 + .01*x1 + %ran(1.0)
set i12 = 2.0 + .02*x2 + %ran(1.0)
set i34 =     - .04*x1 + .03*x2 + %ran(1.0)
set y = %if(i13>=0.0,%if(i12>=0.0,1,2),%if(i34>=0.0,3,4))
*
*   Estimate sequential probit model
*
nonlin a0 a1 a2 b0 b1 b2 c0 c1 c2
frml indexa = a0 + a1*x1 + a2*x2
frml indexb = b0 + b1*x1 + b2*x2
frml indexc = c0 + c1*x1 + c2*x2
frml logl = pa=%cdf(indexa), pb=%cdf(indexb), pc=%cdf(indexc), $
   %if(y==1,log(pa*pb), $
   %if(y==2,log(pa*(1-pb)), $
   %if(y==3,log((1-pa)*pc), $
          log((1-pa)*(1-pc)) )))
maximize(method=bfgs) logl
*
*   Generate ordered probit data
*
set x1 = %rangamma(10.0)
set x2 = %ran(20.0)
set io = .03*x1+.05*x2+%ran(1.0)
set y = 1+(io>-1.0)+(io>0.0)+(io>1.0)
*
*   Estimate ordered probit model
*
nonlin b1 b2 a1 a2 a3
compute b1=b2=0.0,a1=-2.0,a2=0.0,a3=2.0
frml index = b1*x1+b2*x2
frml logl = ix=index(t),$
   %if(y==1,log(%cdf(a1-ix)),$
   %if(y==2,log(%cdf(a2-ix)-%cdf(a1-ix)),$
   %if(y==3,log(%cdf(a3-ix)-%cdf(a2-ix)),$
                     log(%cdf(ix-a3)))))
maximize(method=bfgs,iters=500) logl