*! version 2.2.5 SRH 7 Sept 2011
program define gllarob, eclass
    version 6.0

    syntax [,CLuster(varname) DOts First Macs SCorefil(string) temp noROb]

    tempname Vr b V
    matrix `V' = e(V)
    matrix `b' = e(b)

/* first: have not computed robust se before */
/* macs: macros are already there */
/* rob: do not compute/report robust ses */

*disp in re "first = `first' and mac= `macs' "

/* sort out depvar */
    local depv "`e(depvar)'"    
    global ML_y1 "`depv'"

    
/* sort out scorefil */

    if "`scorefil'"~=""{
        if "`e(scorefil)'"~=""{
            disp in re "There is already a score file `e(scorefil)', option ignored"
            local score
        }
        else{
            * check if file can be opened
            capture postfile junk a using "`scorefil'"
            if _rc~=0{
                disp in re "cannot open `scorefil'.dta for writing"
                exit 198
            }
            postclose junk
            local score "`scorefil'"
        }
    
    }
/* sort out first and calc */
    if "`first'"~=""{
        local first = 1
        local calc = 1
        
    }
    else{
        local first = 0
    }

    local robu = 1
    if "`rob'"~=""{
        local robu = 0
    }
    if `first' == 0 {
        global HG_const = e(const)
        matrix `V' = e(Vs)
        local robclus  "`e(robclus)'"
        local calc = 0
        * disp "`robclus' == `cluster' ?"
        capture matrix `Vr' = e(Vr)
            if _rc>0|"`robclus'"~="`cluster'"{
            local calc = 1
        }
        if "`score'"~=""{
            local calc = 1
        }
    }
    * disp "first = "  `first' " and calc = " `calc'
    *if `first'==0& `calc' {
    if "`macs'"==""& `calc' {
        if "`cluster'"~=""{
            qui count if `cluster'==.&e(sample)
            if r(N)>0{
                disp in re "`cluster' has missing values in the estimation sample"
                exit(198)
            }
        }
        preserve
        qui keep if e(sample)
/* set all global macros needed by gllam_ll */
        setmacs 
/* sort out temporary variables */
/* sort out weight */
        local weight "`e(weight)'"
        global HG_weigh "`e(weight)'"
        local pweight "`e(pweight)'"
*10/29/06
        local numlv: word count `e(clus)'
        tempvar wt
        quietly gen double `wt'=1
*End 10/29/06
        local i = 1
        while `i'<= `numlv'{ /* 10/29/06 not $HG_tplv{ because wrong for init option */
            tempvar wt`i'
            qui gen double `wt`i''=1
            global HG_wt`i' "`wt`i''"
            capture confirm variable `weight'`i'
            if _rc==0 {
                qui replace `wt`i'' = `wt`i'' * `weight'`i'
            }
            capture confirm variable `pweight'`i'
            if _rc==0{
                qui replace `wt`i'' = `wt`i'' * `pweight'`i'
            }
*10/29/06
            qui replace `wt' = `wt'*`wt`i''
            local i = `i'+1
        }
*10/29/06
        if `e(init)'{
            qui replace `wt1' = `wt'
        }

/* sort out level 1 clus variable */
        local clus `e(clus)'
        global HG_clus `clus'
        tempvar id
        if $HG_exp~=1&$HG_comp==0{
            gen long `id'=_n
            if $HG_tplv == 1{
                global HG_clus "`id'"
            }
            else{
                tokenize "`clus'"
                local l= $HG_tplv
                local `l' "`id'"
                global HG_clus "`*'"
            }
        }
        * disp in re "HG_tplv = " $HG_tplv " and HG_clus: $HG_clus" 
/* sort out denom */
        local denom "`e(denom)'"
        if "`denom'"~=""{
            capture confirm variable `denom'
            if _rc>0{
                tempvar den
                qui gen byte `den'=1
                global HG_denom "`den'"
            }
            else{
                global HG_denom `denom'
            }
        }

/* sort out HG_ind */
        capture confirm variable $HG_ind
        if _rc>0{
            tempname junk
            global HG_ind "`junk'"
            gen byte $HG_ind=1
        }

/* sort out HG_lvolo (from gllapred) */
        if $HG_nolog>0{
            tempname junk
            global HG_lvolo "`junk'"
            qui gen byte $HG_lvolo = 0
            local no = 1
            if "$HG_lv"==""{
                local olog = M_olog[1,`no']
                qui replace $HG_lvolo = 1       
            }
            else{
                while `no'<=$HG_nolog{
                    local olog = M_olog[1,`no']
                    qui replace $HG_lvolo = 1 if $HG_lv == `olog'
                    local no = `no' + 1
                }
            }
        }
        
/** (from gllapred) **/  
        if $HG_mlog>0{
            if $HG_nolog==0{
                tempname junk
                global HG_lvolo "`junk'"
                qui gen byte $HG_lvolo = 0
            }
            if "$HG_lv"~=""{ /* more than one link */
                qui replace $HG_lvolo = 1 if $HG_lv == $HG_mlog
            }
            else{
                qui replace $HG_lvolo = 1
            }
        }



/* sort out constraints */
        if $HG_const {
            * disp "constraints used"
            matrix `b' = `b'*M_T
            matrix `V' = M_T'*`V'*M_T
            *matrix junk = M_T*`V'*M_T'
            * matrix list junk
        }

        tempname junk
        capture matrix `junk' = inv(`V')
        if _rc>0{
            disp in re "parameter covariance matrix not invertible"
            exit(198)
        }  

/* set up HG_MU, HG_SD and HG_Cij */
        local rf = 1
        while `rf'<=$HG_tprf{
            tempname junk
            global HG_MU`rf' "`junk'"
            tempname junk
            global HG_SD`rf' "`junk'"
            gen double ${HG_MU`rf'}=0
            gen double ${HG_SD`rf'}=1
            local rf2 = `rf' + 1
            while `rf2' < $HG_tprf {
                tempname junk
                global HG_C`rf2'`rf' "`junk'"
                gen double ${HG_C`rf2'`rf'}=0
                local rf2 = `rf2' + 1
            }
            local rf = `rf' + 1
        }
    } 

/* endif set-up macros and variables */

/* compute scores and/or robust standard errors */
    if `calc'{ 
        *set trace on
        tempvar tpwt
        gen int `tpwt' = 1
        global HG_tpwt "`tpwt'"
        local weight "$HG_weigh"
* 10/20/06
        if $HG_init{
            local numlv: word count `e(clus)'
            local i = 1
            while `i'<=`numlv'{
                capture confirm variable `weight'`i'
                if _rc==0{
                    qui replace `tpwt' = `tpwt'*`weight'`i'
                }
                local i = `i' + 1
            }
        }
        else{
            local l = $HG_tplv
            capture confirm variable `weight'`l'
            if _rc==0 {
                /* there are frequency weights at the top level */
                qui replace `tpwt' = `tpwt' * `weight'`l'
            }
        }
        *summ `tpwt'
        * disp "before comprob, HG_const = " $HG_const " and first = " `first'
        noi cap noi comprob3 "`b'" "`V'" "`cluster'" "`score'" "`dots'" `robu'
        if _rc>0{
            disp in re "something went wrong in comprob3"
            delmacs
            exit 198
        }
        matrix `Vr' = `V'
    }

/* post results */
    *disp in re "temp = `temp' score = `score' "
    if "`score'"~=""&"`temp'"==""{
        est local scorefil "`score'"
    }
    if `robu'{
        * disp "after comprob, HG_const = " $HG_const " and first = " `first'
        if $HG_const&`first' == 0{
            matrix `V' = `Vr'
            est matrix Vr  `V'
            tempname M_T
            matrix `M_T'=e(T)
            matrix `Vr' = `M_T'*`Vr'*`M_T''
            * matrix list `Vr'
        }
        else{
            matrix `V' = `Vr'
            est matrix Vr  `V'
            * matrix list e(Vr)
        }
        estimates repost V = `Vr'
        est local robclus "`cluster'"
    }
    if `first' == 0&`calc'{
        delmacs
    }

end


program define comprob3
version 6.0
args coeffs var cluster score dots rob
*set trace on
    *disp in re "in comprob3: dots=`dots', score = `score', rob=`rob' "
    if "`cluster'"~=""{
        local cluster cluster(`cluster')
    }
    *disp "cluster is `cluster'"
    *matrix list `coeffs'

    tempvar where last
    local l = $HG_tplv
    local weight "$HG_tpwt"

    if $HG_tplv>1{
        local top: word 1 of $HG_clus
    }
    else{
        local k: word count $HG_clus
        local top: word `k' of $HG_clus
    }
    sort `top' $HG_ind
    qui by `top': gen byte `last' = _n==_N
    qui by `top': gen int `where' = _n==1
    qui replace `where' = sum(`where')
    local n = colsof(`coeffs')

/* compute scores */
if "`e(scorefil)'"==""{
    *disp "Computing scores"
    tempname b1 lnf 
    tempname S g fm fp fm0 S0 h Sgoal0 Sgoal1 Sming
    tempvar llp llm lls  done dfvar ok goal0 goal1 mingoal

/* sort out adapt */
    local adapt = 0
    if $HG_adapt==1{
        local adapt = 1
    }
    if `adapt'{
        tempname llast lnf
        global HG_adapt=1
        noi gllam_ll 1 `coeffs' "`lnf'" "junk" "junk" 1
        disp in gre "Non-adaptive log-likelihood: " in ye `lnf'
        scalar `llast' = 0
        local i = 1
                qui `noi' disp in gr "Updating posterior means and variances"
                qui `noi' disp in gr "log-likelihood: "
        while abs((`llast'-`lnf')/`lnf')>1e-8&`i'<60{
                scalar `llast' = `lnf'
                qui gllam_ll 1 "`coeffs'" "`lnf'" "junk" "junk" 1
                disp in ye %10.4f `lnf' " " _c
                        if mod(`i',6)==0 {disp " " }
                local i = `i' + 1
        }
    }

    qui gen double `llm'=.
    qui gen double `llp'=.
    qui gen double `lls'=.
    gllam_ll 1 `coeffs' "`lnf'" "junk" "junk" 3 "`lls'"
    * with mlogit link, all but last lls are missing:
    sort `top' `last'
    qui by `top': replace `lls' = `lls'[_N]

    if abs((`lnf'-`e(ll)')/`e(ll)')>1e-5{
        disp in re "can't get correct log-likelihood: " `lnf' " should be " `e(ll)'
        exit 198
    }

    *disp in re "got the right likelihood: " `e(ll)'
    
    qui count if `last'>0
    local N = r(N)

    
    local l 1e-8
    local u 1e-7
    local m 1e-10
    local epsf 1e-3
    local j = 1
    * disp "N = " `N'
    
    * CHANGED TO VARIABLES         
    gen double `goal0' = (abs(`lls')+`l')*`l'
    gen double `goal1' = (abs(`lls')+`u')*`u'
    gen double `mingoal' = (abs(`lls')+`m')*`m'

    qui gen byte `done' = 0
    qui gen byte `ok' = 0
    qui gen double `dfvar' = 0

    local ds
    local ss

    local i = 1
    while `i'<=`n'{ /* loop over parameters */
        * if "`dots'"~=""{ noi disp in gr "." _c}
*JUNK
        * disp in re " "
        * disp in re "parameter " `i'
        scalar `S' = 1 /* ML_d0_S[1,`i'] */
        scalar `h' = (abs(`coeffs'[1,`i'])+`epsf')*`epsf'
        matrix `b1' = `coeffs'

/* find zero derivatives */
       
        matrix `b1'[1,`i'] = `coeffs'[1,`i']+500*`h'*`S'
*JUNK
*        noi matrix list `b1'
        gllam_ll 1 `b1' "`lnf'" "junk" "junk" 3 "`llp'"
        
        matrix `b1'[1,`i'] = `coeffs'[1,`i']-500*`h'*`S'
*JUNK
*        noi matrix list `b1'
        gllam_ll 1 `b1' "`lnf'" "junk" "junk" 3 "`llm'"
        
        qui replace `done' = abs(`llp' - `llm')<`goal0' /* derivative is zero */
        sort `top' `last'
        qui by `top': replace `done' = `done'[_N]
        
*JUNK
*        disp in re "non-zero derivatives for following clusters (showing lls llm llp and goal0) step = " 1000*`h'*`S'
*        noi list `top' `lls' `llm' `llp' `goal0' if `last'==1&`done'==0
        
        tempvar d`i'
        qui gen double `d`i''= 0 if `done' == 1&`last'==1
        local ds "`ds' `d`i''"

        if "`score'"~=""{
            tempvar s`i'
            qui gen double `s`i''= 0 if `done' == 1&`last'==1
            local ss "`ss' `s`i''"
        }

        qui count if `done'==0&`last'==1
        local num = r(N)
    
        while `num'>0{
            **disp "`num' clusters left to do"
            scalar `S' = 1 /* ML_d0_S[1,`i'] */
            if "`dots'"~=""{ noi disp in ye "." _c}
            sort `done' `where'
            local nxt = `where'[1]
            scalar `lnf' = `lls'[1]
  
            scalar `Sgoal0' = (abs(scalar(`lnf'))+`l')*`l'
            scalar `Sgoal1' = (abs(scalar(`lnf'))+`u')*`u'
            scalar `Sming' = (abs(scalar(`lnf'))+`m')*`m'
            
            preserve

            ***disp in re " "
*JUNK
*            disp in re "sorting out cluster `nxt': `top' = " `top'[1]
            qui keep if `where' == `nxt'

*JUNK            
*            if `nxt'==241{disp in re "lnf = " `lnf' " goal0 = " `Sgoal0' " goal1 = " `Sgoal1'}
            * disp in re "got here!"
            GetStep `coeffs' `h' `S' `i' `Sgoal0' `Sgoal1' `Sming' `lnf' /* "`coeffs'" */
*JUNK
*            if `nxt'==241{ disp in re "after GetStep: S = " `S'  " h = " `h'}

            restore
            matrix `b1'[1,`i'] = `coeffs'[1,`i']-`h'*`S'
            gllam_ll 1 `b1' "`lnf'" "junk" "junk" 3 "`llm'"
            matrix `b1'[1,`i'] = `coeffs'[1,`i']+`h'*`S'
            gllam_ll 1 `b1' "`lnf'" "junk" "junk" 3 "`llp'"
            qui replace `dfvar' = abs(`lls'-`llm')
            qui replace `ok' = `goal0'<`dfvar'&`dfvar'<`goal1'
*JUNK
*            if `nxt'==241{
*                disp in re "goal0 = " `Sgoal0' " goal1 = " `Sgoal1'
*                noi list `goal0' `goal1' if `where'==`nxt'
*                disp in re "llm llp lls"
*                noi list `llm' `llp' `lls' if `where'==`nxt'
*            }
            * first derivatives
            qui replace `d`i'' = (`llp' - `llm')/(2*`h'*`S'*`weight') if `ok' == 1&`last'==1
            if "`score'"~=""{
                * second derivatives
                qui replace `s`i'' = (`llp' + `llm' -2*`lls' )/(`h'*`S'*`h'*`S'*`weight') /*
                */ if `ok' == 1&`last'==1
            }
            qui replace `done' = 1 if `ok'==1
            sort `top' `last'
            qui by `top': replace `done' = `done'[_N]
            
*JUNK?
            * if program didn't crash, derivative for `nxt' is valid:
            qui replace `done' = 1 if `where'==`nxt'
/*
            qui summ `done' if `where' == `nxt', meanonly
            if abs(r(mean)-1)>1e-3{
                disp in re "Problem with derivative for parameter `i', cluster `nxt'" 
            }
*/
            qui count if `done'==0&`last'==1
            local num = r(N)
        }
        *summ `d`i''
        local i = `i'+1
    }
*set trace on
    if "`score'"~=""{
        preserve
        format `ds' `ss' %16.11g
        qui outfile `top' `ds' `ss' if `last'==1 using "`score'.dta", replace
        local tp `top'
        if $HG_tplv==1{
            local tp __idno
        }
        qui infile `tp' double (d1-d`n' s1-s`n') using "`score'.dta", clear
        qui sort `tp'
        qui save "`score'", replace
        restore
    }
} /*endif*/
/* compute robust ses */
    if `rob'{
        local es
        local dds
        local i = 1
        while `i'<=`n'{
            local dds "`dds' d`i'"
            local es "`es' e`i'"
            local i = `i' + 1
        }
        * disp " "
        local eqs: coleq(`var')
        local rown: rownames(`var')
        local coln: colnames(`var')
        * matrix list `var'
        matrix coleq `var' = `es'
        matrix roweq `var' = `es'
        matrix colnames `var' = _cons
        matrix rownames `var' = _cons

        if "`e(scorefil)'"~=""{
            preserve
            local tp `top'
            if $HG_tplv == 1{
                rename $HG_clus __idno
                local tp __idno
            }
            sort `tp' `last'
            merge `tp' using `e(scorefil)'
            noi _robust `dds' [fweight=`weight'] if `last'>0, `cluster' variance(`var')
            restore
        }
        else{
            * matrix list `var'
            *disp "before _robust, _rc = " _rc
            *corr `ds', cov
            * disp "_robust `ds' [fweight=`weight'] if `last'>0, `cluster' variance(`var')"
            noi _robust `ds' [fweight=`weight'] if `last'>0, `cluster' variance(`var')
            *disp "after _robust, _rc = " _rc
            * matrix list `var'
        }

        * disp "setting equation names"
        matrix coleq `var' = `eqs'
        matrix roweq `var' = `eqs'
        matrix colnames `var' = `coln'
        matrix rownames `var' = `rown'
        * matrix list `var'
    }
    exit(0)
end

program define GetStep
    version 6.0

    args a h S i goal0 goal1 mingoal lnf /* coeffs */

    ***disp in re "In GetStep: h = " `h' " S = " `S' " i = " `i' " goal0 = " `goal0' 
    ***disp in re "goal0 = " `goal1' " mingoal = " `mingoal' " lnf = " `lnf' 
    tempname b1 S0 fm0 fp fm coeffs
    matrix `coeffs' = `a'
    matrix `b1' = `coeffs'
    * matrix list `coeffs'
    ***disp in re " "

/***** from here, stolen from ml_adj */
    matrix `b1'[1,`i'] = `coeffs'[1,`i']-`h'*`S'
    gllam_ll 1 `b1' "`fm'" 
    ***disp in re "in iteration 0, lnf = " `lnf' " and fm = " `fm' 

/* Save initial values of S and fm. */

    scalar `S0' = `S'
    scalar `fm0' = `fm'

    if `fm'==.{
        * disp in re "calling MisStep now"
        * disp in re "step before: " `S'
        MisStep `coeffs' `h' `S' `i' `lnf'
        * disp in re "step before: " `S'
        exit
    }

/* Compute df.  We want goal0 <= df <= goal1. */

    local df = abs(scalar(`lnf')-`fm')

    local Sold1 0
    local dfold1 0
    local iter 1

    local itmax = 20

    while (`df'<`goal0' | `df'>`goal1') & `iter'<=`itmax' {

        GetS `mingoal' `goal0' `goal1' `S' `df' /* interpolate ...
        */ `Sold1' `dfold1' `Sold2' `dfold2'
        local Sold2 `Sold1'
        local dfold2 `dfold1'
        local Sold1 = `S'
        local dfold1 `df'

        scalar `S' = r(S)

        matrix `b1'[1,`i'] = `coeffs'[1,`i']-`h'*`S'
*JUNK
*        noi disp in re "changing paramter by " `h'*`S' " to " `b1'[1,`i']
        gllam_ll 1 `b1' "`fm'" 
*JUNK
*        disp in re "in iteration `iter', S= " `S' " fm = " `fm'

        if `fm'==. {
            * disp in re "calling MisStep now"
            * disp in re "step before: " `S'
            MisStep `coeffs' `h' `S' `i' `lnf'
            * disp in re "step after: " `S'
            exit
        }

        local df = abs(scalar(`lnf')-`fm')
        local iter = `iter' + 1
    }

    if `df'<`goal0' | `df'>`goal1' { /* did not meet goal */
        scalar `S' = `S0'    /* go back to initial values */
        scalar `fm' = `fm0'  /* guaranteed to be nonmissing */
    }

    matrix `b1'[1,`i'] = `coeffs'[1,`i']+`h'*`S'
    gllam_ll 1 `b1' "`fp'" 

    if `fp'==. {
        * disp in re "calling MisStep now"
        * disp in re "step before: " `S'
        MisStep `coeffs' `h' `S' `i' `lnf'
        * disp in re "step after: " `S'
        exit
    }

    if `df'<`goal0' | `df'>`goal1' { /* did not meet goal; we redo
                    stepsize adjustment looking at
                    both sides; starting values are
                    guaranteed to be nonmissing
                 */
*JUNK
*        disp in re  "calling TwoStep now, lnf = " `lnf' 
*        disp in re "lnf = " `lnf' " fm = " `fm' " df  = " `df' 
        * disp in re "step before: " `S'
        TwoStep `fp' `fm' `coeffs' `h' `S' `i' `lnf'
        * disp in re "step after: " `S'

    }
/***** up to here, stolen from ml_adj */
end

program define GetS, rclass
* stolen from ml_adj
    args mingoal goal0 goal1 S df Sold1 dfold1 Sold2 dfold2

    if `df' < `mingoal' {
    /* di "GetS: below mingoal, doubling S --> 2*S" */  /* diag */
        return scalar S = 2*`S'
        exit
    }

/* Interpolate to get f(newS)=mgoal.

   When `Sold2' and `dfold2' are empty (on the first iteration), we do
   linear interpolation of f(S)=df, f(0)=0.

   Thereafter, we do quadratic interpolation with the current and previous
   two positions.
*/
    tempname newS
    local mgoal = (`goal0' + `goal1')/2

    Intpol `newS' `mgoal' `S' `df' `Sold1' `dfold1' `Sold2' `dfold2'

    if `newS'==. | `newS'<=0  | (`df'>`goal1' & `newS'>`S') /*
    */                        | (`df'<`goal0' & `newS'<`S') {

        return scalar S = `S'*cond(`df'<`goal0',2,.5)
    }
    else    return scalar S = `newS'
end


program define Intpol
* stolen from ml_adj
    args y x y0 x0 y1 x1 y2 x2

    if "`y2'"=="" { local linear 1 }
    else if `y2'==. | `x2'==. { local linear 1 }

    if "`linear'"!="" {
        scalar `y' = ((`y1')-(`y0'))*((`x')-(`x0'))/((`x1')-(`x0')) /*
        */           + (`y0')
        exit
    }

    scalar `y' = /*
*/   (`y0')*((`x')-(`x1'))*((`x')-(`x2'))/(((`x0')-(`x1'))*((`x0')-(`x2'))) /*
*/ + (`y1')*((`x')-(`x0'))*((`x')-(`x2'))/(((`x1')-(`x0'))*((`x1')-(`x2'))) /*
*/ + (`y2')*((`x')-(`x0'))*((`x')-(`x1'))/(((`x2')-(`x0'))*((`x2')-(`x1')))
end




program define MisStep  /* This routine is called if missing values were
               encountered in GetStep.
            */

    /* di "in MisStep!"  */                 /* diag */
    *args h S caller i fpout fmout x0
    args coeffs h S i lnf
    *macro shift 7
    *local list "`*'"

    local itmax 50

    tempname fm fp b1
    scalar `fm' = .
    scalar `fp' = .
    local iter 1
    while (`fm'==. | `fp'==.) & `iter'<=`itmax' {
        scalar `S' = `S'/2

        matrix `b1' = `coeffs'
        matrix `b1'[1,`i'] = `coeffs'[1,`i']-`h'*`S'
        gllam_ll 1 `b1' "`fm'" 

        if `fm'!=. {
            matrix `b1'[1,`i'] = `coeffs'[1,`i']+`h'*`S'
            gllam_ll 1 `b1' "`fp'" 
        }

        local iter = `iter' + 1
    }

    if `fm'==. | `fp'==. {
        di as err "could not calculate numerical derivatives" _n /*
        */ "discontinuous region with missing values encountered"
        exit 430
    }

    TwoStep `fp' `fm' `coeffs' `h' `S' `i' `lnf'
end

program define TwoStep  /* This routine is called if

                (1) goal was not reached, or

                (2) missing values were encountered
                    and MisStep then found nonmissing
                    values.

               Note: Input is guaranteed to be nonmissing
                     on both sides.
            */

    /* di "in two-step"  */                 /* diag */
    *args fp fm h S caller i fpout fmout x0
    args fp fm coeffs h S i lnf
    *macro shift 9
    *local list "`*'"

    tempname bestS b1

    local ep0   1e-8
    local ep1   1e-7
    local epmin 1e-12
    local itmax 40

    local goal0   = (abs(scalar(`lnf'))+`ep0')*`ep0'
    local goal1   = (abs(scalar(`lnf'))+`ep1')*`ep1'
    local mingoal = (abs(scalar(`lnf'))+`epmin')*`epmin'

    local df = (abs(scalar(`lnf')-`fp')+abs(scalar(`lnf')-`fm'))/2
    local bestdf `df'
    scalar `bestS' = `S'
    local Sold1 0
    local dfold1 0
    local iter 1

    while (`df'<`goal0' | `df'>`goal1') & `iter'<=`itmax' {

*di "TwoStep   iter = `iter'   df = " %12.4e `df' "   S = "  %12.3e `S'

        GetS `mingoal' `goal0' `goal1' `S' `df' /* interpolate ...
        */ `Sold1' `dfold1' `Sold2' `dfold2'

        local Sold2 `Sold1'
        local dfold2 `dfold1'
        local Sold1 = `S'
        local dfold1 `df'

        scalar `S' = r(S)

        matrix `b1' = `coeffs'
        matrix `b1'[1,`i'] = `coeffs'[1,`i']-`h'*`S'
        gllam_ll 1 `b1' "`fm'" 

        *Lik`caller' -`h'*`S' `i' `fm' `fmout' `x0' `list'

        if `fm'!=. {
            matrix `b1'[1,`i'] = `coeffs'[1,`i']+`h'*`S'
            gllam_ll 1 `b1' "`fp'" 
            * Lik`caller' `h'*`S' `i' `fp' `fpout' `x0' `list'
        }
        if `fm'==. | `fp'==. {
            if `bestdf' >= `mingoal' { /* go with best value */
                scalar `S' = `bestS'

            matrix `b1' = `coeffs'
            matrix `b1'[1,`i'] = `coeffs'[1,`i']-`h'*`S'
            gllam_ll 1 `b1' "`fm'" 
            matrix `b1'[1,`i'] = `coeffs'[1,`i']+`h'*`S'
            gllam_ll 1 `b1' "`fp'" 

                di as txt /*
                */ "numerical derivatives are approximate" /*
                */ _n "nearby values are missing"
                exit
            }

            di as err /*
            */ "could not calculate numerical derivatives" /*
            */ _n "missing values encountered"
            exit 430
        }

        local df = (abs(scalar(`lnf')-`fp')+abs(scalar(`lnf')-`fm'))/2

        if `df'>1.1*`bestdf' | (`df'>=0.9*`bestdf' & `S'<`bestS') {
            local bestdf `df'
            scalar `bestS' = `S'
        }

        local iter = `iter' + 1
    }

*JUNK
*disp in re "TwoStep   df = " %12.4e `df' "   S = "  %12.3e `S' " goal0 = " `goal0'  " goal1 = " `goal1' " lnf = " `lnf'

    if `df'<`goal0' | `df'>`goal1' { /* did not reach goal */

disp in re "TwoStep: did not reach goal"

        if `bestdf' >= `mingoal' { /* go with best value */
            scalar `S' = `bestS'

            matrix `b1' = `coeffs'
            matrix `b1'[1,`i'] = `coeffs'[1,`i']-`h'*`S'
            gllam_ll 1 `b1' "`fm'" 
            matrix `b1'[1,`i'] = `coeffs'[1,`i']+`h'*`S'
            gllam_ll 1 `b1' "`fp'" 

            di as txt "numerical derivatives are approximate" /*
            */ _n "flat or discontinuous region encountered"
        }
        else {
            di as err "could not calculate numerical derivatives" /*
            */ _n "flat or discontinuous region encountered"
            exit 430
        }
    }
end



program define setmacs
version 6.0
/* may not work for higher level models yet */
/* tplv */
    global HG_tplv = e(tplv)

/* link and family-related macros */
    global HG_famil "`e(famil)'"
    global HG_link "`e(link)'"
    global HG_linko "`e(linko)'"
    global HG_nolog = `e(nolog)'
    global HG_ethr = `e(ethr)'
    global HG_lv "`e(lv)'"
    global HG_fv "`e(fv)'"
    global HG_oth = e(oth)
    global HG_mlog = e(mlog)
    global HG_smlog = `e(smlog)'
    capture matrix M_resp=e(mresp)
    capture matrix M_respm=e(mrespm)
    capture matrix M_frld=e(frld)
    capture matrix M_olog=e(olog)
    capture matrix M_oth=e(moth)
    capture matrix ML_d0_S=e(do_S)
    global HG_exp = e(exp)
    global HG_expf = e(expf)
    global HG_ind = "`e(ind)'"
    global HG_init = e(init)
    global HG_lev1 = e(lev1)
    global HG_comp = e(comp)
    capture local coall "`e(coall)'"
    if $HG_comp~=0{
        local i = 1
        while `i'<=$HG_comp{
            local k: word `i' of `coall'
            global HG_co`i' `k'
            local i = `i' + 1
        }
    }

/* prior related global macros */
    global HP_prior = e(prior)
    if $HP_prior == 1{
            global HP_sprd = 1
            global HP_invga = e(invga)
            global HP_invwi = e(invwi)
            global HP_foldt = e(foldt)
            global HP_logno = e(logno)
            global HP_gamma = e(gamma)
            global HP_corre = e(corre)
            global HP_boxco = e(boxco)
            global HP_spect = e(spect)
            global HP_wisha = e(wisha)
            if $HP_invga==1{
                global shape = e(shape)
                global rate = e(rate)
            }
            if $HP_invwi==1{
                global df = e(df)
                matrix scale = e(scale)
            } 
            if $HP_foldt==1{
                global df = e(df)
                global scale = e(scale)
                global location = e(location)
            }
            if $HP_logno==1{
                global meanlog = e(meanlog)
                global sdlog = mean(sdlog)
            }
            if $HP_gamma==1{
                global HP_scale = e(scale)
                global HP_var = e(var)
                global HP_shape = e(shape)
            }
            if $HP_corre==1{
                global alpha = e(alpha)
                global beta = e(beta)
            }
            if $HP_boxco==1{
                global scale = e(scale)
                global lambda = e(lambda)
            }
            if $HP_spect==1{
                global alpha = e(alpha)
                global beta = e(beta)
            }
            if $HP_wisha==1{
                global wisha = e(wisha)
                global df = e(df)
                matrix scale =e(scale)
            }
            matrix M_nu = e(nu)
    }


/* set all other global macros */
    global HG_nats = `e(nats)'
    global HG_noC = `e(noC)'
    global HG_noC1 = `e(noC1)'
    global HG_adapt = `e(adapt)'
    global HG_tplv = e(tplv)
    global HG_tpff = `e(tpff)'
    global HG_tpi = `e(tpi)'
    global HG_free = e(free)
    global HG_mult = e(mult)
    global HG_lzpr lzprobg
    global HG_zip zipg
    if $HG_mult{
        global HG_lzpr lzprobm
    }
    else if $HG_free{
        global HG_lzpr lzprobf
        global HG_zip zipf
        matrix M_np=e(mnp)
    }
    global HG_cip = e(cip)
    global which = 4
    global HG_off "`e(offset)'"
    global HG_error = 0
    global HG_cor = `e(cor)'
    global HG_bmat = e(bmat)
    global HG_tprf = e(tprf)
    global HG_const = e(const)
    if $HG_const==1{
        matrix M_T = e(T)
        matrix M_a = e(a)
    }
    global HG_ngeqs = e(ngeqs)
    global HG_inter = e(inter)
    global HG_dots = 0
    matrix M_nbrf = e(nbrf)
    matrix M_nrfc = e(nrfc)
    matrix M_ip =  J(1,$HG_tprf+2,1)
    matrix M_nffc =  e(nffc)
    local tprf = $HG_tprf
    if `tprf'<2 { local tprf = 2 }
    matrix M_znow =J(1,`tprf'-1,1)
    matrix M_nip = e(nip)
    capture matrix M_ngeqs = e(mngeqs)
    capture matrix M_b=e(mb)
    *matrix M_chol = e(chol)
    local l = M_nrfc[1,1] + 1  /* loop */
    local k = M_nrfc[2,1] + 1  /* r. eff. */
    if $HG_tplv>1{
    while `l'<=M_nrfc[1,2]{
        while `k'<=M_nrfc[2,2]{
            * disp "loop " `l' " random effect " `k'
            local w = M_nip[2,`k']

            /* same loc and prob as before? */
            local found = 0
            local ii=M_nrfc[2,1] + 1
            while `ii'<`k'{
                if `w'==M_nip[2,`ii']{
                    local found = 1
                }
                local ii = `ii'+1
            }


            capture matrix M_zps`w' =e(zps`w')

            * disp "M_zlc`w'"
            matrix M_zlc`w'=e(zlc`w')
            local k = `k' + 1
        }
        local l = `l' + 1
    }
    }
end

program define delmacs
    version 6.0
/* deletes all global macros and matrices*/
    tempname var
    if "$HG_tplv"==""{
        * macros already gone
        exit
    }
    local nrfold = M_nrfc[2,1]
    local lev = 2
    while (`lev'<=$HG_tplv){
        local i2 = M_nrfc[2,`lev']
        local i1 = `nrfold'+1
        local i = `i1'
        local nrfold = M_nrfc[2,`lev']
        while `i' <= `i2'{
            local n = M_nip[2,`i']
            if `i' <= M_nrfc[1,`lev']{
                capture matrix drop M_zps`n'
            }
            capture matrix drop M_zlc`n'
            local i = `i' + 1
        }
        local lev = `lev' + 1
    }
    if $HG_free==0{
        capture matrix drop M_chol
    }
    matrix drop M_nrfc
    matrix drop M_nffc
    matrix drop M_nbrf
    matrix drop M_ip
    capture matrix drop M_b
    capture matrix drop M_resp
    capture matrix drop M_respm
    capture matrix drop M_frld
    matrix drop M_nip
    matrix drop M_znow
    capture matrix drop M_ngeqs
    capture matrix drop CHmat

    /* globals defined in gllam_ll */
    local i=1
    while (`i'<=$HG_tpff){
        global HG_xb`i'
        local i= `i'+1
    }
    local i = 1
    while (`i'<=$HG_tprf){
        global HG_s`i'
        local i= `i'+1
    }
    local i = 1
    while (`i'<=$HG_tplv){
        global HG_wt`i'
        local i = `i' + 1
    }
    global HG_nats
    global HG_noC
    global HG_noC1
    global HG_noC2
    global HG_adapt
    global HG_fixe
    global HG_lev1
    global HG_bmat
    global HG_tplv 
    global HG_tprf
    global HG_tpi
    global HG_tpff
    global HG_clus
    global HG_weigh
    global which   
    global HG_gauss 
    global HG_free
    global HG_famil 
    global HG_link 
    global HG_linko
    global HG_nolog
    global HG_lvolo
    global HG_oth
    global HG_mlog
    global HG_exp
    global HG_expf
    global HG_lv
    global HG_fv
    global HG_nump
    global HG_eqs
    global HG_obs
    global HG_off
    global HG_denom
    global HG_cor
    global HG_s1
    global HG_init
    global HG_ind
    global HG_const
    global HG_dots
    global HG_inter
    global HG_ngeqs
    global HG_tpwt
    global HG_ethr
    global HG_mult
    global HG_lzpr
    global HG_zip
    global HG_cip
    global HG_comp
    global HG_pwt
    global HG_befB
    global HG_smlog
    global HG_cn
    capture drop macro HG_co*
end
