******************************************************************************** *! "pr_pred_cub", v.19, Cerulli, 02oct2020 ******************************************************************************** prog pr_pred_cub , eclass syntax varlist(min=1) [if] [in] [fweight pweight] , prob(name) local y `varlist' tempvar p M R S ******************************************************************************** marksample touse ******************************************************************************** cub informat [`weight'`exp'] if `touse' , pi() xi() vce(oim) // esimates "cub00" tempvar theta1 predict `theta1' , equation(pi_beta) tempvar theta2 predict `theta2' , equation(xi_gamma) local m=e(M) // number of categories ******************************************************************************** quietly generate double `p' = 1/(1+exp(-`theta1')) local c = exp(lnfactorial(`m'-1)) mat cmb = J(`m',1,.) forv i=1/`m' { sca d = (exp(lnfactorial(`i'-1))*exp(lnfactorial(`m'-`i'))) mat cmb[`i',1] = `c'/d } qui gen double `M' = cmb[`y',1] quietly generate double `R' = ((exp(-`theta2'))^(`y'-1))/((1+exp(-`theta2'))^(`m'-1)) quietly generate double `S' = 1/`m' cap gen `prob' = (`p'*(`M'*`R'-`S')+`S') ******************************************************************************** di as text in red "" di as text in red "{hline}" di as text in red "{bf:******************************************************************************}" di as text in red "{bf:********************* PROBABILITY BY CATEGORY ********************************}" di as text in red "{bf:******************************************************************************}" di as text in red "{hline}" tab `y' `prob' di as text in red "{hline}" ******************************************************************************** end ********************************************************************************