*! version 1.6.6 8/20/03 -- jf fix of bug where observed probabilities improperly computed for plot option
*! version 1.6.5 1/10/01

capture program drop prcounts
program define prcounts
    version 6.0
    tempname b alpha ai gai mu
    syntax newvarname [if] [in] [, Max(integer 9) Plot]

    if `max' > 99 | `max' < 0 {
        di in r "max() must be in range from 0 to 99"
        exit 198
    }

    local stem "`varlist'"
    local stem = substr("`stem'", 1, 28)
    cap version 7
        if _rc!=0 { local stem = substr("`stem'", 1, 4) }
    version 6.0
    local modelis "from `e(cmd)'"

*-> GENERATE PREDICTED RATE
    quietly predict `stem'rate `if' `in', n
    label variable `stem'rate "Predicted rate `modelis'"
    gen `mu' = `stem'rate

*-> GENERATE PREDICTED FREQUENCY OF ALWAYS ZERO
    if "`e(cmd)'"=="zip" | "`e(cmd)'"=="zinb" {
        quietly predict `stem'all0 `if' `in', p
        label variable `stem'all0 "Pr(always 0) `modelis'"
        * predict n, n computes: exp(x*beta)*(1-all0).
        * The last term needs to be removed.
        quietly replace `mu' = `mu'/(1-`stem'all0) `if' `in'
    }

*-> TAKE CARE OF ALPHA
    if "`e(cmd)'"=="nbreg" {
        sca `alpha' = e(alpha)
    }
    if "`e(cmd)'"=="zinb" {
        mat `b' = e(b)
        local temp = colsof(`b')
        sca `alpha' = `b'[1, `temp']
        sca `alpha' = exp(`alpha')
    }
    if "`e(cmd)'"=="nbreg" | "`e(cmd)'"=="zinb" {
        sca `ai' = 1/`alpha'
        sca `gai' = exp(lngamma(`ai'))
        if `gai'==. {
            di in r "problem with alpha prevents" /*
            */ " estimation of predicted probabilities."
            exit 198
        }
    }

*-> GENERATE PREDICTED PROBABILITIES
    local i 0
    while `i' <= `max' {
        local newvar "`stem'pr`i'"
        if "`e(cmd)'"=="poisson" {
            quietly gen `newvar' = /*
            */ ((exp(-`mu'))*(`mu'^`i'))/(round(exp(lnfact(`i'))), 1) /*
            */ `if' `in'
        }
        if "`e(cmd)'"=="nbreg" {
            quietly gen `newvar' = /*
            */ (exp(lngamma(`i'+`ai')) / /*
            */ (round(exp(lnfact(`i')),1) * exp(lngamma(`ai')))) /*
            */ * ((`ai'/(`ai'+`mu'))^`ai') * ((`mu'/(`ai'+`mu'))^`i') /*
            */ `if' `in'
        }
        if "`e(cmd)'"=="zip" {
            quietly gen `newvar' = (1-`stem'all0)*((exp(-`mu')) /*
            */ * (`mu'^`i'))/(round(exp(lnfact(`i'))), 1) `if' `in'
            if `i'==0 {
                quietly replace `newvar' = `newvar' + `stem'all0 `if' `in'
            }
        }
        if "`e(cmd)'"=="zinb" {
            quietly gen `newvar' = (1-`stem'all0)*(exp(lngamma(`i'+`ai')) /*
            */ / ( round(exp(lnfact(`i')),1) * exp(lngamma(`ai')) )  ) /*
            */ * ((`ai'/(`ai'+`mu'))^`ai') * ((`mu'/(`ai'+`mu'))^`i') /*
            */ `if' `in'
            if `i'==0 {
                quietly replace `newvar' = `newvar' + `stem'all0 `if' `in'
            }
        }
        label variable `newvar' "Pr(y=`i') `modelis'"
        local i = `i' + 1
    }

*-> GENERATE cuMULATIVE PROBABILITIES
    quietly gen `stem'cu0=`stem'pr0
    label variable `stem'cu0 "Pr(y=0) `modelis'"
    local i 1
    while `i' <= `max' {
        quietly egen `stem'cu`i' = rsum(`stem'pr0-`stem'pr`i') if `mu'~=.
        label variable `stem'cu`i' "Pr(y<=`i') `modelis'"
        local i = `i' + 1
    }

*-> GENERATE GREATER THAN VARIABLE
    quietly gen `stem'prgt = 1-`stem'cu`max' if `mu'~=.
    label variable `stem'prgt "Pr(y>`max') `modelis'"

*-> IF PLOT OPTION SPECIFIED
    if "`plot'"=="plot" {
        quietly gen `stem'val = .
        label variable `stem'val "Count"
        quietly gen `stem'obeq = .
        label variable `stem'obeq "Observed Pr(y=k) `modelis'"
        quietly gen `stem'preq = .
        label variable `stem'preq "Predicted Pr(y=k) `modelis'"
        quietly gen `stem'oble = .
        label variable `stem'oble "Observed Pr(y<=k) `modelis'"
        quietly gen `stem'prle = .
        label variable `stem'prle "Predicted Pr(y<=k) `modelis'"

        ** bug fix -- makes sure observed probabilities are computed on estimation sample
        if "`if'" == "" {
            local if "if e(sample)==1"
        }
        else {
            local if "`if' & e(sample)==1"
        }

        local i 0
        while `i' <= `max' {
            quietly {
                local obs = `i' + 1
                replace `stem'val = `i' in `obs'
                tempvar count1 count2
                gen `count1' = (`e(depvar)'==`i') `if' `in'
                sum `count1' `if' `in'
                replace `stem'obeq = r(mean) in `obs'
                gen `count2' = (`e(depvar)'<=`i') `if' `in'
                sum `count2' `if' `in'
                replace `stem'oble = r(mean) in `obs'
                sum `stem'pr`i' `if' `in'
                replace `stem'preq = r(mean) in `obs'
                sum `stem'cu`i' `if' `in'
                replace `stem'prle = r(mean) in `obs'
            }
            local i = `i' + 1
        }
    }
end