//  tpr - VanHeuvelen & Long (tvanheuv@indiana.edu) | 2013-09-29
//  tpr_pop_geny: create y outcome

program define tpr_pop_geny_form1

    version 13

    capture drop e
    capture drop xb
    capture drop ystar
    capture drop ycat

qui {

    * logistic error
    gen double e = -log(1/runiform() - 1)
    label var e "e: logistic error"
    sum e
    scalar sdstep1 = r(sd)
    scalar sdtheory = sqrt((c(pi)^2)/3)
    scalar sddiv = sdstep1 / sdtheory
    qui replace e = (e - r(mean)) / sddiv

    * xb: use rescaled betas to be equal as std coefficients
    foreach rhsnm of varlist $rhsnms {
        noi sum `rhsnm'
        local sd`rhsnm' = r(sd)
    }
    gen xb = $b0 + ($b1/`sdx1')*x1 + ($b2)*x2 ///
            + ($b3/`sdx3')*x3 + ($b4/`sdx4')*x4
    label var xb "xb: xb"

    * ystar: latent variable
    gen ystar = xb + e
    label var ystar "ystar: xb + e"

    * ycat: categorized outcome
    local catsNm1 = $catsN - 1
    local icat = 1
    foreach c in $y_dist {
        local ipct : word `icat' of $y_dist
        local pct`icat' = `ipct'
        local ++icat
    }
    forvalues c = `icat'/5 { // dummy out categories not used
        local pct`c' = 0
    }
    local pct_yeq1 = `pct1'
    local pct_yeq2 = `pct1' + `pct2'
    local pct_yeq3 = `pct1' + `pct2' + `pct3'
    local pct_yeq4 = `pct1' + `pct2' + `pct3' + `pct4'
    local pct_yeq5 = `pct1' + `pct2' + `pct3' + `pct4' + `pct5'

    forvalues icat = 1/`catsNm1' {
        centile ystar, centile(`pct_yeq`icat'')
        local cut`icat' = r(c_1)
    }
    gen ycat = .
    replace ycat = 0 if ystar <= `cut1'
    forvalues icat = 1/`catsNm1' {
        replace ycat = `icat' if ystar > `cut`icat''
    }
    label var ycat "ycat: $catsN-category outcome"

    gen id = _n
    label var id "population observation number"

    noi tabulate ycat, matcell(ycat)
    forvalues icat = 1/$catsN {
        global ncat`icat' = ycat[`icat',1]
        global pct`icat' = $ncat`icat'/r(N)
    }
}

end
exit