program define binphat

   * Quick-and-dirty program to compute predicted probabilities of
   * a binomial probit.  Outcomes are presumed to be {0,1,2} with
   * omitted category 0.  Value functions are:
   *
   *      V0 = 0
   *      V1 = XB1 + u1
   *      V2 = XB2 + u2
   *
   * This program computes P(Y=0) = P(V0 > max(V1,V2)),
   * P(Y=1) = P(V1 > max(V0,V2)), and P(Y=2) = P(V2 > max(V0,V1)).
   *
   * Usage:  binphat xb1 xb2 rho phat0 phat1 phat2
   *    xb1     [in] variable with XB1
   *    xb2     [in] variable with XB2
   *    rho     [in] scalar representing corr(u1,u2)
   *    phat0  [out] variable with predicted P(Y=0)
   *    phat1  [out] variable with predicted P(Y=1)
   *    phat2  [out] variable with predicted P(Y=2)
   *
   * Stan Panis, 31 March 2003

   version 7

   ****************************************
   *  Confirm that all arguments are OK.  *
   ****************************************
   local error = "Usage:  binphat xb1 xb2 rho phat0 phat1 phat2"
   if ("`6'"=="" | "`7'"!="") {
      display in red "`error'"
      exit 999
      }

   confirm variable `1'
   local xb1 "`1'"
   confirm variable `2'
   local xb2 "`2'"

   confirm number `3'
   local rho = `3'

   macro shift 3
   local varlist "required new min(3) max(3)"
   parse "`*'"
   parse "`varlist'", parse(" ")
   local phat0 "`1'"
   local phat1 "`2'"
   local phat2 "`3'"

   ***************************
   *  All arguments are OK.  *
   ***************************

   * Calculate predicted P(Y=0)
   qui replace `phat0' = binorm(-`xb1',-`xb2',`rho')

   * Calculate predicted P(Y=1)
   local newrho = (1-`rho')/sqrt(2-2*`rho')
   qui replace `phat1' = binorm(`xb1', (`xb1'-`xb2')/sqrt(2-2*`rho'), `newrho')

   * Calculate predicted P(Y=2)
   local newrho = (1-`rho')/sqrt(2-2*`rho')
   qui replace `phat2' = binorm(`xb2', (`xb2'-`xb1')/sqrt(2-2*`rho'), `newrho')

end