*! version 1.0.0 2013-10-06

//  tpr - VanHeuvelen & Long (tvanheuv@indiana.edu)
//  tpr_pop_design01: logistic; homoscedasticity; 4 x's

program define tpr_pop_design01

    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: homoscedastic logistic error"
	global error_dist "logistic homoscedastic"
    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
    gen xb = $pop_b0 + $pop_b1*x1 + $pop_b2*x2 + $pop_b3*x3 + $pop_b4*x4
    label var xb "xb: xb"

    * ystar: latent variable
    gen ystar = xb + e
    label var ystar "ystar: xb + e"
    * ycat: categorized outcome
    local catsNm1 = $y_catsN - 1
    local icat = 1
    foreach iprob in $y_dist {
        local pct`icat' = `iprob'*100
        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: $y_catsN-category outcome"

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

    noi tabulate ycat, matcell(ycat)
    forvalues icat = 1/$y_catsN {
        global ncat`icat' = ycat[`icat',1]
        local gnm ncat`icat'
        global y_est`icat'prob = $`gnm'/r(N)
    }
}

end
exit