*! version 2.0.1  30mar99  STB-49 sg97.1
*
* syntax: outreg [varlist] using filename [, SE Pvalue noLabel 
*                BDec(numlist int >=0 <12) Tdec(numlist int >=0 <12 max=1) 
*                noPAren BRacket noASter noCONs Xstats Onecol EForm 
*                CLAssic AT noNOBs noNI noR2 ADjr2 
*                RDec(numlist int >=0 <12 max=1) noNOTes COMma APpend REPLACE]
*

program define outreg
   version 6.0
   * write formatted regression output to file

   syntax [varlist(default=none)] using [, SE Pvalue noLabel BDec(numlist >=0 <=11) Tdec(numlist >=0 <=11) noPAren BRacket noASter noCONs Xstats Onecol EForm CLAssic AT(string) noNOBs noNI noR2 ADjr2 RDec(numlist >=0 <=11 max=1) noNOTes COMma APpend REPLACE]

   if "`se'"=="se" & "`pvalue'"=="pvalue" {
      di in red "cannot choose both se and pvalue options"
      exit 498
   }
   if "`bdec'"=="" {local bdec = 3}
   if "`tdec'"=="" {local tdec = 2}
   if "`rdec'"=="" {local rdec = `tdec'}

   tempname regN df_r rsq numi r2mat
   tempname varname coefcol b vc b_alone

   /* for svy commands with subpop(), N_sub is # of obs used for estimation */
   if e(N_sub) != . {scalar `regN'   = e(N_sub)}  
   else {scalar `regN'   = e(N)}
   scalar `df_r'   = e(df_r)
   scalar `numi'   = e(N_g)
   local   depvar  = e(depvar)
   local   robust  = e(vcetype)
   if "`robust'"=="." {local robust "none"}
   local   ivar    = e(ivar)
   mat     `b'     = e(b)
   mat     `vc'    = e(V)
   local   univar  = 1          /* univariate (single equation) regression */
   local   bcols   = colsof(`b')  /* cols of b */
   local   bocols  = `bcols'      /* cols of b only, w/o other stats */
   local   labdp   = ""           /* for use with mlogit (see) */
   capture local fracpol = (e(fp_cmd)=="fracpoly")
    
   /* logistic reports coeffients in exponentiated form (odds ratios) */
   if e(cmd)=="logistic" {local eform "eform"}
   if "`eform'"=="eform" {local cons "nocons"}
   
   /* multi-equation models */
   if e(cmd)=="mvreg" | e(cmd)=="sureg" | e(cmd)=="reg3" {
      if "`onecol'" != "onecol" {
         local univar  = 0      /* multivariate (multiple equation) regression */
         mat `r2mat'   = e(b)   /* get column labels */
         local neq     = e(k_eq)
         local eqlist  = e(eqnames)
         if e(cmd)=="mvreg" {
            local r2list  = e(r2)
            local depvar = "`eqlist'"
         }
         local eq = 1
         while `eq' <= `neq' {
            if e(cmd)=="mvreg" { 
               local r2str: word `eq' of `r2list'
               scalar `rsq' = real("`r2str'")
            }
            else {scalar `rsq' = e(r2_`eq')}
            mat `r2mat'[1,`eq'] = `rsq'
            local eq = `eq' + 1
         }
      }
      else {   /* if onecol */
         local r2 = "nor2"  
         scalar `rsq' = .
      }
   }     /* `rsq' after `r2list' to avoid type mismatch  */
   else if "`adjr2'"=="adjr2" {scalar `rsq' = e(r2_a)}
   else {scalar `rsq' = e(r2)}

   if (e(cmd)=="intreg" | e(cmd)=="svyintrg" | e(cmd)=="xtintreg"){ 
      local depvar : word 1 of `depvar'  /* 2 depvars listed */
   }

   if (e(cmd)=="arch" | e(cmd)=="arima") {
      local xcons = "xcons"  /* include "Constant" in extra stats */
   }

   /* nolabels for anova and fracpoly */
   else if (e(cmd)=="anova" | `fracpol' | e(cmd)=="nl") {  
      /* e(fp_cmd)!=. means fracpoly */
      local cons "nolabel"
   }

   /* dprobit transformation of b and V */
   else if e(cmd)=="dprobit" {
      if "`at'"!="" {local at = "at(`at')"}
      dprbcalc, `at' `classic' b(`b') vc(`vc') 
      local cons "nocons"
   }

   /* bivariate probit, heckman, heckprob */
   else if (e(cmd)=="biprobit" | e(cmd)=="heckman" | e(cmd)=="heckprob") {
      if ("`onecol'"!="onecol") {
         local univar  = 0      /* multivariate (multiple equation) regression */
         local neq     = e(k_dv)
         local eqlist  = e(depvar)
         if e(cmd)=="heckman" {
            local eqlist = "`eqlist'" + " select"
            local depvar = "`depvar'" + " select"
         }
      }
      else { /* if onecol */
         local depvar : word 1 of `depvar'
      }
   }

   /* multinomial regression */
   else if ((e(cmd)=="mlogit" | e(cmd)=="svymlog") & "`onecol'"!="onecol") {
      local univar = 0    /* multivariate (multinomial) regression */
      local labdp = "nolabdp"  /* no var label for mlogit categories */
      local neq   = e(k_cat)-1
      mat `vc' = e(V)
      local k = colsof(`vc')/`neq'
      local catlist : coleq(`vc')
      local depvar ""
      local cat = 1
      while `cat' <= `neq' {
         local catwrdn = (`cat'-1)*`k' + 1 
         local catname: word `catwrdn' of `catlist'
         local depvar = "`depvar' " + "`catname'"
         local cat = `cat' + 1
      }
      local eqlist  = "`depvar'"
   }

   /* when e(b) includes extra statistics (which don't have variable labels) */
   capture mat `b_alone' = `b'[1,"`depvar':"]
   if _rc==0 {local bocols = colsof(`b_alone')}
   else if (e(cmd)=="ologit" | e(cmd)=="oprobit") {
      local bocols = e(df_m)
      mat `b_alone' = `b'[1,1..`bocols']
   }
   else if (e(cmd)=="cnreg" | e(cmd)=="tobit") {
      local bocols = `bocols'-1 /* last element of e(b) is not est coef */
      mat `b_alone' = `b'[1,1..`bocols']
   }
   else if (e(cmd)=="intreg" | e(cmd)=="svyintrg"){
      mat `b_alone' = `b'[1,"model:"]
      local bocols = colsof(`b_alone')
   }
   if `univar' & `bcols'>`bocols' & "`xstats'"!="xstats" {
      mat `b' = `b_alone'
      mat `vc' = `vc'[1..`bocols',1..`bocols']
   }
   
   
   quietly {
      preserve
      
      if `univar' {
         coeftxt `varlist', `se' `pvalue' `label' bdec(`bdec') tdec(`tdec') `paren' `bracket' `aster' `cons' `eform' `nobs' `ni' `r2' `adjr2' rdec(`rdec') `notes' `append' regN(`regN') df_r(`df_r') rsq(`rsq') numi(`numi') ivar(`ivar') depvar(`depvar') robust(`robust') `labdp' borows(`bocols') `xcons' b(`b') vc(`vc') varname(`varname') coefcol(`coefcol')
   
         if "`append'"=="append" {
            appfile `using', varname(`varname') coefcol(`coefcol') 
            outsheet v* `coefcol' `using', nonames noquote `comma' replace
            drop v*
         }
         else {   
            noi outsheet `varname' `coefcol' `using', nonames noquote `comma' `replace'
         }
      }
      else {   /* multiple equation regression */
         tempname b_eq vc_eq
         local depvlst `depvar'
         local eq = 1
         while `eq' <= `neq' {
            restore, preserve  /* to access var labels after first equation */
            local eqname: word `eq' of `eqlist'
            local depvar: word `eq' of `depvlst'
            /* r2mat doesn't exist for mlogit => "capture" */
            capture scalar `rsq' = `r2mat'[1,`eq'] 
            mat `b_eq'  = `b'[.,"`eqname':"]
            mat `vc_eq' = `vc'["`eqname':","`eqname':"]
            local bocols = colsof(`b_eq')
   
            if `eq' == 1 & "`append'"~="append" {local apptmp ""}
            else {local apptmp "append"}
            coeftxt `varlist', `se' `pvalue' `label' bdec(`bdec') tdec(`tdec') `paren' `bracket' `aster' `cons' `eform' `nobs' `ni' `r2' `adjr2' rdec(`rdec') `notes' `apptmp' regN(`regN') df_r(`df_r') rsq(`rsq') numi(`numi') ivar(`ivar') depvar(`depvar') robust(`robust') `labdp' borows(`bocols') `xcons' b(`b_eq') vc(`vc_eq') varname(`varname') coefcol(`coefcol')
            if `eq' == 1 & "`append'"~="append" {
               noi outsheet `varname' `coefcol' `using', nonames noquote `comma' `replace'
            }
            else {
               appfile `using', varname(`varname') coefcol(`coefcol') 
               outsheet v* `coefcol' `using', nonames noquote `comma' replace
               drop v*
            }
            local eq = `eq' + 1
         }
      }
      
   }  /* for quietly */
end


program define coeftxt
   version 6.0
   * convert vectors b and vc into coefficient and t-statistics text + variable name text

   syntax [varlist(default=none)] [, SE Pvalue noLabel BDEC(numlist) TDEC(numlist) noPAren BRacket noASter noCONs EForm noNOBs noNI noR2 ADjr2 RDec(numlist) noNOTes APpend regN(string) df_r(string) rsq(string) numi(string) ivar(string) depvar(string) robust(string) noLabdp BOROWS(string) xcons b(string) vc(string) varname(string) coefcol(string)]

   tempvar beta st_err tstat astrix mrgrow 
   tempfile bcopy 
   tempname b_alone vc_alon b_xtra vc_xtra

   mat `b' = `b''
   mat `vc' = vecdiag(`vc')
   mat `vc' = `vc''
   local brows = rowsof(`b')

   if ("`label'"=="" & "`labdp'"=="" & "`depvar'"!="" ) {
      local depname : var label `depvar'
      }
   if "`depname'"=="" {local depname "`depvar'"}

   if (`numi'~=. & "`ni'"~="noni") {
      if ("`label'"=="" & "`ivar'"~=".") {local iname : var label `ivar'}
      if "`iname'"=="" {local iname "`ivar'"}
      if "`iname'"=="." {local iname "groups"}
   }

   /* in case varlist is specified:  */
   mat `b_alone' = `b'[1..`borows',1] /* use to avoid _cons in xtra stats */
   if `brows'>`borows' {  /* allow for xtra stats */
      local borows1 = `borows'+1
      mat `b_xtra' = `b'[`borows1'...,1...]
      mat `vc_xtra' = `vc'[`borows1'...,1...]
   }
   if "`varlist'"~="" {
      tempname arow testnew newb newvc
        /* add the constant unless "nocons" is chosen */
      if "`cons'"~="nocons" {local varlist = "`varlist'" + " _cons"}
      local vname : word 1 of `varlist'
      local i=1
      while "`vname'"~="" {
         local j = rownumb(`b_alone',"`vname'")
         if `j'!=. {
            matrix `arow'  = `b'[`j',1...]  /* "..." needed to get rownames */
            matrix `newb'  = nullmat(`newb')\ `arow'
            matrix `arow'  = `vc'[`j',1...]
            matrix `newvc' = nullmat(`newvc')\ `arow'
         }
         else if (e(cmd)!="sureg" & e(cmd)!="reg3" & "`vname'"~="_cons") {
            di in red "`vname' not found in regression coefficients"
            exit 111
         }
         local i = `i'+1
         local vname : word `i' of `varlist'
      }
      mat `b_alone' = `newb'
      if `brows'>`borows' {  /* allow for xtra stats */
         mat `newb'  = `newb'\ `b_xtra'
         mat `newvc' = `newvc'\ `vc_xtra'
      }
      mat `b' = `newb'
      mat `vc' = `newvc'
   }
   else if "`cons'"=="nocons" {   /* delete the constant if "nocons" is chosen */
      local j_1 = rownumb(`b_alone',"_cons")-1
      if `j_1'!=. {  /* in case there is no constant in b */
         mat `b_alone' = `b_alone'[1..`j_1',1...]
         mat `vc_alon' = `vc'[1..`j_1',1...]
         if `brows'==`borows' {  
            mat `b' = `b_alone'
            mat `vc' = `vc_alon'
         }
         else {  /* allow for xtra stats */
            mat `b' = `b_alone' \ `b_xtra'
            mat `vc' = `vc_alon' \ `vc_xtra'
         }
      }
   }
   local borows = rowsof(`b_alone')
   local brows = rowsof(`b')  /* reset brows */

   gen double `beta' = matrix(`b'[_n, 1]) in 1/`brows'
   keep if `beta'~=.
   gen double `st_err' = matrix(`vc'[_n, 1]) in 1/`brows'
   replace `st_err' = sqrt(`st_err')
   gen double `tstat' = abs(`beta'/`st_err')
   if "`pvalue'"=="pvalue" {
      if `df_r'==. {replace `tstat' = 2*(1-normprob(`tstat'))}
      else{replace `tstat' = tprob(`df_r', `tstat')}
   }
   if "`eform'"=="eform" {
      replace `beta' = exp(`beta')
      replace `st_err' = `beta'*`st_err'
   }

   gen str31 `varname' = ""
   local bnames : rownames(`b_alone')
   local bname : word 1 of `bnames'
   local i 1
   while "`bname'"~="" {
      if "`bname'"~="_cons" {
         if "`label'"!="nolabel" {
            local dotloc = index("`bname'", ".") 
            if `dotloc'!=0 {    /* deal w/time series prefixes */
               local tspref = substr("`bname'",1,`dotloc')
               local bname = substr("`bname'",`dotloc'+1,.)
               local blabel : var label `bname'
               local blabel = "`tspref' " + "`blabel'"
            }
            else {local blabel : var label `bname'}
         }
         else {local blabel ""}
      }
      else { 
         local blabel "Constant" 
      }
      if "`blabel'"=="" { 
         local blabel "`bname'" 
      }
      replace `varname' = "`blabel'" if `i'==_n
      local i = `i'+1
      local bname : word `i' of `bnames'
   }
   if `brows'>`borows' {  /* allow for xtra stats */
      local borows1 = `borows'+1
      mat `b_xtra' = `b'[`borows1'...,1...]
      local bnames : rownames(`b_xtra')
      local beqs : roweq(`b_xtra')
      local bname : word 1 of `bnames'
      local beq : word 1 of `beqs'
      local i 1
      while "`bname'"~="" {
         if "`bname'"~="_cons" {
            if "`beq'"=="_" {local blabel "`bname'"}
            else {local blabel "`beq':`bname'"}
         }
         else if "`xcons'"=="xcons" {local blabel "`beq':Constant"}
         else {local blabel "`beq'"}
         replace `varname' = "`blabel'" if `i'+`borows'==_n
         local i = `i'+1
         local bname : word `i' of `bnames'
         local beq : word `i' of `beqs'
      }
   }

     /* get rid of original data since labels accessed */
   keep `varname' `beta' `tstat' `st_err'  
   
   gen `mrgrow' = 2*_n + 1 + ("`append'"=="append")
   if "`aster'"~="noaster" {
      if "`pvalue'"=="pvalue" {
         gen str2 `astrix' = "*" if (`tstat'<=0.05 & `tstat'!=.)
         replace `astrix' = "**" if (`tstat'<=0.01 & `tstat'!=.)
      }
      else {
         if `df_r'==. {
            scalar Tcrit1 = invnorm(0.995)
            scalar Tcrit5 = invnorm(0.975)
         }
         else {
            scalar Tcrit1 = invt(`df_r',0.99)
            scalar Tcrit5 = invt(`df_r',0.95)
         }
         gen str2 `astrix' = "*" if (`tstat'>=Tcrit5 & `tstat'!=.)
         replace `astrix' = "**" if (`tstat'>=Tcrit1 & `tstat'!=.)
      }
   }
   else {
      gen str2 `astrix' = ""
   }      
   if "`paren'"!="noparen" {
      if "`bracket'"=="bracket" {
         local lparen "["
         local rparen "]"
      }
      else {
         local lparen "("
         local rparen ")"
      }
   }
   else {
      local lparen ""
      local rparen ""
   }
   local bdeccnt : word count `bdec'   
   if "`se'"~="se" {gen str12 `coefcol' = string(`tstat',"%12.`tdec'f")}  
   else if `bdeccnt'==1 {
      gen str12 `coefcol' = string(`st_err',"%12.`bdec'fc")
   }
   else {
      gen str12 `coefcol' = ""
      local abdec : word 1 of `bdec'
      local i 1
      while `i'<=_N {
         replace `coefcol' = string(`st_err',"%12.`abdec'fc") if `i'==_n
         local i = `i'+1
         local prvbdec = `abdec'
         local abdec : word `i' of `bdec'
         if "`abdec'"=="" {local abdec = `prvbdec'}
      }
   }
   replace `coefcol' = "`lparen'" + `coefcol' + "`rparen'" + `astrix'
   sort `mrgrow'
   save `bcopy', replace
   replace `mrgrow' = `mrgrow'-1
   if `bdeccnt'==1 {replace `coefcol' = string(`beta',"%12.`bdec'fc")}
   else {
      local abdec : word 1 of `bdec'
      local i 1
      while `i'<=_N {
         replace `coefcol' = string(`beta',"%12.`abdec'fc") if `i'==_n
         local i = `i'+1
         local prvbdec = `abdec'
         local abdec : word `i' of `bdec'
         if "`abdec'"=="" {local abdec = `prvbdec'}
      }
   }
   sort `mrgrow'
   merge `mrgrow' using `bcopy'
   replace `varname' = " " if _merge==2
   drop `beta' `st_err' `tstat' `astrix' _merge
   sort `mrgrow'
   save `bcopy', replace

   drop `varname' `coefcol'
   local coefrow = 2*`brows'+1+("`append'"=="append")
   local totrows = `coefrow' + ("`nobs'"!="nonobs") + /*
    */(`numi'!=.) + ("`r2'"!="nor2"&`rsq'!=.&`df_r'!=.) + /*
    */("`notes'"!="nonotes"&"`append'"!="append")*(1+("`aster'"!="noaster"))

   local minobs = max(1000,`brows')
   set obs `minobs'  /* obs must be set to >= 100 (depending on matsize?) */
   keep if _n <= `totrows'
   replace `mrgrow' = _n
   sort `mrgrow'
   merge `mrgrow' using `bcopy'
   sort `mrgrow'
   drop _merge
   
   replace `coefcol' = "`depname'" if _n==1 & ("`append'"~="append")
   replace `coefcol' = "`depname'" if _n==2 & ("`append'"=="append")
   if "`nobs'"~="nonobs" {
      local coefrow = `coefrow'+1
      replace `varname' = "Observations" if _n==`coefrow'
      replace `coefcol' = string(`regN') if _n==`coefrow'
   }
   if (`numi'~=. & "`ni'"~="noni") {
      local coefrow = `coefrow'+1
      replace `varname' = "Number of `iname'" if _n==`coefrow'
      replace `coefcol' = string(`numi') if _n==`coefrow'
   }
   if "`r2'"~="nor2" & `rsq'~=. & `df_r'!=. { /* if df_r=., not true r2 */
      local coefrow = `coefrow'+1
      replace `coefcol' = string(`rsq',"%12.`rdec'f") if _n==`coefrow'
      replace `varname' = "R-squared" if _n==`coefrow'
      if "`adjr2'"=="adjr2" {
         replace `varname' = "Adjusted " + `varname' if _n==`coefrow'
      }
   }
   if "`notes'"~="nonotes" & "`append'"~="append" {
      local coefrow = `coefrow'+1
      if "`bracket'"=="bracket" {local par_bra "brackets"}
      else {local par_bra "parentheses"}
      if "`se'"~="se" {
         if `df_r'~=. {local t_or_z "t"}
         else {local t_or_z "z"}
         if "`robust'"=="none" {
            replace `varname' = "Absolute value of `t_or_z'-statistics in `par_bra'" if _n==`coefrow'
         }
         else{
            replace `varname' = "`robust' `t_or_z'-statistics in `par_bra'" if _n==`coefrow'
         }
      }
      else {
         if "`robust'"=="none" {
            replace `varname' = "Standard errors in `par_bra'" if _n==`coefrow'
         }
         else{
            replace `varname' = "`robust' standard errors in `par_bra'" if _n==`coefrow'
         }
      }
      if "`aster'"~="noaster" {
         local coefrow = `coefrow'+1
         replace `varname' = "* significant at 5% level; ** significant at 1% level" if _n==`coefrow'
      }
    }
end  /* coeftxt */


program define appfile
   version 6.0
   * append regression results to pre-existing file
   
   syntax using/, varname(string) coefcol(string)

   tempname vorder vartmp varsml
   tempfile tmpf1
   gen str80 `vartmp' = `varname'
   replace `vartmp' = "0" if _n==1
   replace `vartmp' = `vartmp'[_n-1]+"1" if trim(`vartmp')==""
   keep `varname' `coefcol' `vartmp'
   sort `vartmp'
   save `tmpf1', replace
   if index("`using'", ".")==0 {local using = "`using'"+".out"}
   insheet using `using', nonames clear
   gen str80 `vartmp' = v1
   replace `vartmp' = "01" if _n==1 & v2~="(1)"
   replace `vartmp' = "0" if _n==1 & v2=="(1)"
   replace `vartmp' = `vartmp'[_n-1]+"1" if `vartmp'==""
   gen `vorder' = _n+(1*(v2[1]~="(1)"))
   gen str10 `varsml' = `vartmp'
   count if (`varsml'=="R-squared" | `varsml'=="Adjusted R"  /*
       */| `varsml'=="Observatio" | `varsml'=="Number of " )
   local oldstat = r(N)
   count if (`varsml'=="Constant")
   local oldcon = r(N)
   drop `varsml'
   sort `vartmp'

   merge `vartmp' using `tmpf1'
   gen str10 `varsml' = `vartmp'
   count if (`varsml'=="Constant" | `varsml'=="Constant1" /*
       */| `varsml'=="Observatio" | `varsml'=="R-squared"  /*
       */| `varsml'=="Adjusted R" | `varsml'=="Robust t-s" /*
       */| `varsml'=="Robust z-s" | `varsml'=="Absolute v" /*
       */| `varsml'=="Robust sta" | `varsml'=="Standard e" /*
       */| `varsml'=="* signific" | `varsml'=="Number of " /*
       */| `varsml'=="Newey-West" )
   local noncoef = r(N)
   count if (`varsml'=="R-squared" | `varsml'=="Adjusted R"  /*
       */| `varsml'=="Observatio" | `varsml'=="Number of " )
   local newstat = r(N) - `oldstat'
   count if (`varsml'=="Constant")
   local newcon = r(N) - `oldcon'
   count if _merge==2
   local newcoef = r(N) - 2*`newcon' - `newstat'
   count if v2=="(1)"
   local endcoef = _N + 1*(r(N)==0) - `noncoef' - `newcoef'
   replace `vorder' = `vorder'+_N-`endcoef'+2*`newcon' if `vorder'>`endcoef' & `vorder'~=.
   replace `vorder' = _n if `vorder'==.
   if `newcon' {  /* when constant in appended reg, but not prev reg */
      replace `vorder' = _N+1 if `vartmp'=="Constant" 
      replace `vorder' = _N+2 if `vartmp'=="Constant1"
   } 
   if `newstat' {  /* when R2, N, or Ni in appended reg, but not prev reg */
                   /* put after Constant & Constant1 (if any) */
      replace `vorder' = _N + 2*(`oldcon'|`newcon') + 0.2 /*
         */ if `varsml'=="Observatio"
      replace `vorder' = _N + 2*(`oldcon'|`newcon') + 0.5 /*
         */ if `varsml'=="Number of "
      replace `vorder' = _N + 2*(`oldcon'|`newcon') + 0.7 /*
         */ if (`varsml'=="R-squared" | `varsml'=="Adjusted R")
   } 
   replace `vorder' = 1 if `vartmp'=="0"
   sort `vorder'
   replace v1 = `varname' if v1=="" & `varname'~=""
   drop `vorder' `varname' `vartmp' `varsml' _merge
   describe
   local numcol = r(k)-1
   if (`numcol'==2) {
      replace v2 = "(1)" if _n==1
      replace `coefcol' = "(2)" if _n==1
   }
   else {
      replace `coefcol' = "(" + string(`numcol') + ")" if _n==1
   }
end  /* appfile */


program define dprbcalc
   version 6.0
   * calculate probability change for dprobit results 
   * adapted from dprobit.ado

   syntax , [at(string) classic] b(string) vc(string)

   tempname mns Vhat v xb fxb c

   if "`at'"!="" { 
      if colsof(`b')==colsof(matrix(`at')) { 
         if `at'[1,colsof(`b')] != 1 { 
            di in red "`at':  last (_cons) element not 1"
            exit 598
         }
         mat `mns' = `at'
      }
      else if colsof(`b')-1 == colsof(matrix(`at')) { 
         mat `mns' = nullmat(`at'),1
      }
      else { 
         di in red "`at':  conformability error" 
         exit 503
      }
   }
   else {mat `mns' = get(mns)}

   if "`classic'"!="" {
      local cons "_cons"
   }
   _evlist  /* returns s(varlist) containing names of e(b), _cons removed */
   tokenize `e(depvar)' `s(varlist)' `cons'  /* puts names into # macros */
   sret clear

   Vhatdiag `mns' `b' `vc' -> `Vhat'

   mat `xb' = `mns'*`b''
   scalar `xb' = `xb'[1,1]
   scalar `fxb' = normd(`xb')

   local lhs "`1'"
   mac shift

   local i 1
   while `"``i''"'!="" {
      local isdum : word `i' of `e(dummy)'
      if `isdum' & `"`classic'"'=="" {
         scalar `c' = /*
         */ normprob(`xb'+(1-`mns'[1,`i'])*`b'[1,`i']) - /*
         */ normprob(`xb'-(`mns'[1,`i'])*`b'[1,`i']) 
         vdummy `i' `mns' `b' `vc' -> `v'
      }
      else {
         scalar `c' = `b'[1,`i']*`fxb'
         scalar `v' = `Vhat'[1,`i']
      }

      mat `b'[1,`i'] = `c'
      mat `vc'[`i',`i'] = `v'  /* only diagonals of vc are correct */

      local i=`i'+1
   }
end


program define Vhatdiag /* mns b V -> res */
   args xbar b V ARROW res

   tempname tmp dgdb

   local k = colsof(`V')

   mat `tmp' = `xbar'*`b'' + e(offbar)
   mat `dgdb' = normd(`tmp'[1,1]) * ( I(`k') - `tmp'[1,1]*`b''*`xbar' )
   mat `res' = vecdiag(`dgdb'*`V'*`dgdb'')

end


program define vdummy /* i mns b V -> res */
   args i mns b V ARROW res

   tempname x xb fxb Delp Vhat
   mat `x' = `mns'
   mat `x'[1,`i']=1
   mat `xb' = `x'*`b'' + e(offbar)
   scalar `fxb' = normd(`xb'[1,1])
   mat `Delp' = `fxb' * `x'

   mat `x'[1,`i']=0
   mat `xb' = `x'*`b'' + e(offbar)
   scalar `fxb' = normd(`xb'[1,1])
   mat `xb' = `fxb'*`x' 
   mat `Delp' = `Delp' - `xb'

   mat `Vhat' = `Delp'*`V'*`Delp''
   scalar `res' = `Vhat'[1,1]
end