*! 1.1.13 MLB 06Sep2010 *! 1.1.8 MLB 08Apr2010 *! 1.1.0 MLB 04Jun2009 program define seqlogit_uh_lf version 9.2 /*count number of equations*/ forvalues i = 1/$S_eqs { local xb "`xb' xb`i'" } args lnf `xb' tokenize $S_treelevels local N_levs : word count $S_treelevels tempvar touse2 qui gen byte `touse2' = 0 local start = 0 local end = $S_drawstart forvalues i = 1/`N_levs' { qui replace `touse2' = $ML_y1 == ``i'' & $ML_samp == 1 local start = `start' + `end' qui count if `touse2' local end = `start' + r(N) mata comp_ll(``i'', `start') } end mata: void comp_ll(real scalar lev, real scalar start) { // get xb eqs = strtoreal(st_global("S_eqs")) xbnames = J(1,eqs,"") for(i=1; i<= eqs; i++) { xbnames[1,i] = st_local("xb" + strofreal(i)) } xb = . st_view(xb,.,xbnames,st_local("touse2")) // get y y = . st_view(y,.,st_global("ML_y1"),st_local("touse2")) // get sd_var if (st_global("S_sd_var") != "") { sd_var = . st_view(sd_var,.,st_global("S_sd_var"),st_local("touse2")) } // get lnf lnf = . st_view(lnf,.,st_local("lnf"),st_local("touse2")) // get other globals t = tokens(st_global("S_t" + strofreal(lev))) c = strtoreal(tokens(st_global("S_c" + strofreal(lev)))) draws = strtoreal(st_global("S_draws")) eqstr = J(length(t),1,"") Nchoice = J(length(t),1,.) for(i=1; i <= length(t); i++) { Nchoice[i,1] = strtoreal(st_global("S_Nchoice" + t[i])) eqstr[i] = st_global("S_eqstr" + t[i]) } sigma = strtoreal(tokens(st_global("S_sigma" ))) deltasd = strtoreal(tokens(st_global("S_deltasd"))) rho = . rho = strtoreal( st_global("S_rho" ) ) pr = . pr = strtoreal(tokens(st_global("S_pr" ))) mpnts = strtoreal(tokens(st_global("S_mpnts" ))) m_mn = strtoreal(tokens(st_global("S_m_mn" ))) pr_mn = strtoreal(tokens(st_global("S_pr_mn" ))) sd_mn = strtoreal(tokens(st_global("S_sd_mn" ))) // create e if (st_global("S_pr")!= "") { u = halton(draws*rows(xb), 1, start) e = J(draws*rows(xb),1,0) cumpr = 0 for (i = 1; i <= length(pr); i++){ prev_cumpr = cumpr cumpr = cumpr + pr[i] e = e :+ ( (u:=prev_cumpr) ) :* mpnts[i] } } else if (st_global("S_pr_mn") != "" ) { u = halton(draws*rows(xb), 2, start) e = J(draws*rows(xb),1,0) cumpr = 0 for (i = 1; i <= length(pr_mn); i++){ prev_cumpr = cumpr cumpr = cumpr + pr_mn[i] e = e :+ ( (u[.,1]:=prev_cumpr) ) :* ( ( invnormal(u[.,2]) :* sd_mn[i] ) :+ m_mn[i] ) } } else if (st_global("S_uniform") != "") { e = (halton(draws*rows(xb), 1, start) :- .5) :/ sqrt(1/12) } else { e = invnormal(halton(draws*rows(xb), 1, start)) } if (rho == .) { } else { x = . st_view(x,.,st_global("S_ofinterest"),st_local("touse2")) e = colshape(e,draws) sdx = strtoreal(st_global("S_sdx")) mx = strtoreal(st_global("S_mx")) e = (rho:*(x :- mx):/sdx :+ sqrt(1-rho^2):*e) e = colshape(e,1) } // compute probabilities w = J(rows(e),1,1) p = J(rows(xb),1,1) for(i = 1; i <= length(t); i++) { denom = J(rows(xb),draws,1) teqstr = strtoreal(tokens(eqstr[i])) for(j = 1; j <= (Nchoice[i,1]-1); j++) { if (st_global("S_deltasd") == "") { t_sigma = sigma[1,teqstr[j]] } else { t_sigma = sigma[1,teqstr[j]] :+ (deltasd[1,teqstr[j]] :* sd_var) } denom = denom :+ exp( xb[.,teqstr[j]]:+ (t_sigma:*colshape(e,draws)) ) } if (c[i]==0) { temp = colshape(1:/denom, 1) } else { if (st_global("S_deltasd") == "") { t_sigma = sigma[1,teqstr[c[i]]] } else { t_sigma = sigma[1,teqstr[c[i]]] :+ (deltasd[1,teqstr[c[i]]] :* sd_var) } temp = colshape(exp(xb[.,teqstr[c[i]]]:+ (t_sigma:*colshape(e,draws))) :/ denom, 1) } p = p:* ( rowsum(colshape(w:*temp, draws)):/rowsum(colshape(w,draws)) ) w = w:*temp } // compute and replace the log-likelihood lnf[.,.] = ln(p) } end