capture program drop switchoprobit_work program define switchoprobit_work, eclass local gradient global neqs = 5+2*$ncut forvalues i = 1/$neqs { local gradient "`gradient' g`i'" } args todo b lnfj `gradient' //H tempvar eta1 eta2_0 eta2_1 tempname r1 r0 forvalues i = 1/$ncut { tempname c1`i' tempname c0`i' } mleval `eta1' = `b', eq(1) mleval `eta2_0' = `b', eq(2) mleval `eta2_1' = `b', eq(3) mleval `r0' = `b', eq(4) scalar mleval `r1' = `b', eq(5) scalar //di "ncut = $ncut" local r 0 forvalues i = 0/1 { forvalues j = 1/$ncut { local `++r' local k = 5+`r' mleval `c`i'`j'' = `b', eq(`k') scalar } } local touse $touse local rho0 = tanh(`r0') local rho1 = tanh(`r1') tempvar eta2i rhoi qui g double `eta2i'=`eta2_0' if $ML_y1==0 qui replace `eta2i'=`eta2_1' if $ML_y1==1 qui g double `rhoi' = `rho0' if $ML_y1==0 qui replace `rhoi' = `rho1' if $ML_y1==1 qui levelsof($ML_y2) global nchoices: word count `r(levels)' //local imat0 //local imat1 forvalues i = 1/$nchoices { //local j = `i'+$nchoices tempvar m0`i' m1`i' ll0`i' ll1`i' //ll`j' qui g double `m0`i'' = $ML_y1==0 & ($ML_y2==`i') qui g double `m1`i'' = $ML_y1==1 & ($ML_y3==`i') qui g double `ll0`i'' = 0 qui g double `ll1`i'' = 0 /*qui g double `ll`j'' = 0*/ local imat0 "`imat0' `m0`i''" local imat1 "`imat1' `m1`i''" } local imat `imat0' `imat1' local neginf = minfloat() local posinf = maxfloat() tempname cutpts0 cutpts1 cutptst mat `cutpts0' = J(1,$ncut,0) mat `cutpts1' = J(1,$ncut,0) forvalues i = 1/$ncut { mat `cutpts0'[1,`i'] = `c0`i'' mat `cutpts1'[1,`i'] = `c1`i'' } mat `cutpts0' = `neginf', `cutpts0', `posinf' mat `cutpts1' = `neginf', `cutpts1', `posinf' mat `cutptst' = `neginf', 0, `posinf' tempname fv_tmp qui g double `fv_tmp'=0 forv d = 0/1 { local e = `d'+2 local f = `d'+1 local cut_h1 = `cutptst'[1,`e'] local cut_l1 = `cutptst'[1,`f'] forv i = 1/$nchoices { local j = `i'+$nchoices local k = `i'+1 local cut_h2 = `cutpts`d''[1,`k'] local cut_l2 = `cutpts`d''[1,`i'] qui replace `ll`d'`i'' = `m`d'`i''*(((binorm(`cut_h1'-`eta1',`cut_h2'-`eta2_`d'',`rho`d''))- /// (binorm(`cut_l1'-`eta1',`cut_h2'-`eta2_`d'',`rho`d''))-(binorm(`cut_h1'-`eta1',`cut_l2'-`eta2_`d'',`rho`d''))+ /// (binorm(`cut_l1'-`eta1',`cut_l2'-`eta2_`d'',`rho`d'')))) qui replace `fv_tmp' = `fv_tmp'+`ll`d'`i'' } } qui replace `lnfj' = ln(`fv_tmp') if (`todo'==0) exit local gradfun forv i = 1/$neqs { tempvar g`i'_tmp qui g double `g`i'_tmp'=0 if `touse' local gradfun "`gradfun' `g`i'_tmp'" } local r_init_0 = tanh(`r0') local r_init_1 = tanh(`r1') tempname treat qui g `treat' = $ML_y1 mata: mata_gradfun_switchoprobit("`eta1'", "`eta2_0'", "`eta2_1'", "`cutptst'", "`cutpts0'", /// "`cutpts1'", "`r_init_0'", "`r_init_1'", "`imat'", "`lnfj'", "`touse'", "`gradfun'") local gfun forvalues i = 1/$neqs { qui replace `g`i''=`g`i'_tmp' if `touse' local gfun "`gfun' `g`i''" } if (`todo'==1) exit end mata: void mata_gradfun_switchoprobit(string scalar xb1, string scalar xb2_0, string scalar xb2_1, /// string scalar cutptst, string scalar cutpts0, string scalar cutpts1, string scalar r0, /// string scalar r1, string scalar indmat, string scalar llike, string scalar touse, /// string scalar gfunout) { n = st_nobs() real matrix gradfun, imat real colvector eta1, eta2_0, eta2_1, lnlike real scalar rho1, rho0 eta1=st_data(.,xb1,touse) eta2_0=st_data(.,xb2_0,touse) eta2_1=st_data(.,xb2_1,touse) imat=st_data(.,tokens(indmat),touse) lnlike=st_data(.,llike,touse) cg = cols(tokens(gfunout)) rho0 = atanh(strtoreal(r0)) rho1 = atanh(strtoreal(r1)) //printf("rho0= %3.2f \n", rho0) //printf("rho1= %3.2f \n", rho1) cutpts_0 = st_matrix(cutpts0) cutpts_1 = st_matrix(cutpts1) cutpts_t = st_matrix(cutptst) fvi = exp(lnlike) ci = cols(imat) cc = ci/2 imat0 = imat[|1,1\n,cc|] imat1 = imat[|1,cc+1\n,ci|] //printf("ci=%2.0f \n", ci) choices = ci/2 //printf("choices = %2.0f \n", choices) neqs = 5+2*(choices-1) gradfun = J(n,neqs,0) gradfun_10 = J(n,cc,0) gradfun_11 = J(n,cc,0) gradfun_2 = J(n,cc,0) gradfun_3 = J(n,cc,0) delta0 = sqrt(1-rho0^2) delta1 = sqrt(1-rho1^2) /*derivatives wrt to beta1, beta2_0, and beta2_1*/ cut_ht = cutpts_t[1,1] cut_lt = cutpts_t[1,2] Atp1 = cut_ht:-eta1 At = cut_lt:-eta1 for(j=1;j<=choices;j++){ cut_h0 = cutpts_0[1,j+1] cut_l0 = cutpts_0[1,j] A0p1 = cut_h0:-eta2_0 A0 = cut_l0:-eta2_0 gradfun_10[.,j]=(normalden(Atp1):*(normal((A0p1:-rho0:*Atp1):/(delta0)))):- /// (normalden(At):*(normal((A0p1:-rho0:*At):/(delta0)))):- /// (normalden(Atp1):*(normal((A0:-rho0:*Atp1):/(delta0)))):+ /// (normalden(At):*(normal((A0:-rho0:*At):/(delta0)))) gradfun_2[.,j]=normalden(A0p1):*(normal((Atp1:-rho0:*A0p1):/(delta0))):- /// normalden(A0):*(normal((Atp1:-rho0:*A0):/(delta0))):- /// normalden(A0p1):*(normal((At:-rho0:*A0p1):/(delta0))):+ /// normalden(A0):*(normal((At:-rho0:*A0):/(delta0))) } cut_ht = cutpts_t[1,2] cut_lt = cutpts_t[1,3] Atp1 = cut_ht:-eta1 At = cut_lt:-eta1 for(j=1;j<=choices;j++){ cut_h1 = cutpts_1[1,j+1] cut_l1 = cutpts_1[j] A1p1 = cut_h1:-eta2_1 A1 = cut_l1:-eta2_1 gradfun_11[.,j]=(normalden(Atp1):*(normal((A1p1:-rho1:*Atp1):/(delta1)))):- /// (normalden(At):*(normal((A1p1:-rho1:*At):/(delta1)))):- /// (normalden(Atp1):*(normal((A1:-rho1:*Atp1):/(delta1)))):+ /// (normalden(At):*(normal((A1:-rho1:*At):/(delta1)))) gradfun_3[.,j]=normalden(A1p1):*(normal((Atp1:-rho1:*A1p1):/(delta1))):- /// normalden(A1):*(normal((Atp1:-rho1:*A1):/(delta1))):- /// normalden(A1p1):*(normal((At:-rho1:*A1p1):/(delta1))):+ /// normalden(A1):*(normal((At:-rho1:*A1):/(delta1))) } gradfun[.,1]=((rowsum(imat1:*gradfun_11):/fvi):+(rowsum(imat0:*gradfun_10):/fvi)) gradfun[.,2]=(rowsum(imat0:*gradfun_2):/fvi) gradfun[.,3]=(rowsum(imat1:*gradfun_3):/fvi) /*derivatives wrt r's, rho's */ gradfun_4 = J(n,cc,0) gradfun_5 = J(n,cc,0) gamma0 = (1-rho0^2) gamma1 = (1-rho1^2) nu0 = 1/sqrt(gamma0) nu1 = 1/sqrt(gamma1) zeta0 = ((2*pi())*sqrt(1-rho0^2))^-1 zeta1 = ((2*pi())*sqrt(1-rho1^2))^-1 cut_h1 = cutpts_t[1,2] cut_l1 = cutpts_t[1,1] A1p1 = cut_h1:-eta1 A1 = cut_l1:-eta1 for(j=1;j<=choices;j++){ cut_h2 = cutpts_0[1,j+1] cut_l2 = cutpts_0[1,j] A2p1 = cut_h2:-eta2_0 A2 = cut_l2:-eta2_0 gradfun_4[.,j]=(exp((-.5):*(((A1p1:^2):+(A2p1:^2):-(rho0:*A1p1:*A2p1:*2)):/gamma0)):- /// exp((-.5):*(((A1:^2):+(A2p1:^2):-(rho0:*A1:*A2p1:*2)):/gamma0)):- /// exp((-.5):*(((A1p1:^2):+(A2:^2):-(rho0:*A1p1:*A2:*2)):/gamma0)):+ /// exp((-.5):*(((A1:^2):+(A2:^2):-(rho0:*A1:*A2:*2)):/gamma0))):*zeta0 } rr0 = tanh(rho0) drho_dr0 = 1/(cosh(rr0)^2) // derivative of rho wrt r gradscalar = drho_dr0 //zeta gradfun[.,4] = gradscalar:*((rowsum(imat0:*gradfun_4)):/fvi) cut_h1 = cutpts_t[1,3] cut_l1 = cutpts_t[1,2] A1p1 = cut_h1:-eta1 A1 = cut_l1:-eta1 for(j=1;j<=choices;j++){ cut_h2 = cutpts_1[1,j+1] cut_l2 = cutpts_1[1,j] A2p1 = cut_h2:-eta2_1 A2 = cut_l2:-eta2_1 gradfun_5[.,j]=(exp((-.5):*(((A1p1:^2):+(A2p1:^2):-(rho1:*A1p1:*A2p1:*2)):/gamma1)):- /// exp((-.5):*(((A1:^2):+(A2p1:^2):-(rho1:*A1:*A2p1:*2)):/gamma1)):- /// exp((-.5):*(((A1p1:^2):+(A2:^2):-(rho1:*A1p1:*A2:*2)):/gamma1)):+ /// exp((-.5):*(((A1:^2):+(A2:^2):-(rho1:*A1:*A2:*2)):/gamma1))):*zeta1 } rr1 = tanh(rho1) //printf("rho1=%4.3f \n", rr1) drho_dr1 = 1/(cosh(rr1)^2) // derivative of rho wrt r gradscalar = drho_dr1 gradfun[.,5] = gradscalar:*((rowsum(imat1:*gradfun_5)):/fvi) /*derivatives wrt cutpoints*/ ncuts = choices-1 gradfun_6_tmp = J(n,cc,0) gradfun_6 = J(n,ncuts,.) for(i=1;i<=ncuts;i++){ cut_h1 = cutpts_t[1,2] cut_l1 = cutpts_t[1,1] A1p1 = cut_h1:-eta1 A1 = cut_l1:-eta1 for(k=1;k<=choices;k++){ delta_kj1 = (i==k) delta_kj = (i==k-1) cut_h2 = cutpts_0[1,k+1] cut_l2 = cutpts_0[1,k] A2p1 = cut_h2:-eta2_0 A2 = cut_l2:-eta2_0 gradfun_6_tmp[.,k]=((normal((A1p1:-rho0:*A2p1):/delta0):-normal((A1:-rho0:*A2p1):/delta0)):*normalden(A2p1):*delta_kj1):- /// ((normal((A1p1:-rho0:*A2):/delta0):-normal((A1:-rho0:*A2):/delta0)):*normalden(A2):*delta_kj) } gradfun_6[.,i] = (rowsum(imat0:*gradfun_6_tmp)):/fvi } gradfun_7_tmp = J(n,cc,0) gradfun_7 = J(n,ncuts,.) for(i=1;i<=ncuts;i++){ cut_h1 = cutpts_t[1,3] cut_l1 = cutpts_t[1,2] A1p1 = cut_h1:-eta1 A1 = cut_l1:-eta1 for(k=1;k<=choices;k++){ delta_kj1 = (i==k) delta_kj = (i==k-1) cut_h2 = cutpts_1[1,k+1] cut_l2 = cutpts_1[1,k] A2p1 = cut_h2:-eta2_1 A2 = cut_l2:-eta2_1 gradfun_7_tmp[.,k]=((normal((A1p1:-rho1:*A2p1):/delta1):-normal((A1:-rho1:*A2p1):/delta1)):*normalden(A2p1):*delta_kj1):- /// ((normal((A1p1:-rho1:*A2):/delta1):-normal((A1:-rho1:*A2):/delta1)):*normalden(A2):*delta_kj) } gradfun_7[.,i] = (rowsum(imat1:*gradfun_7_tmp)):/fvi } gradfun_cuts = gradfun_6, gradfun_7 //gradfun_cuts[|1,1\3,cols(gradfun_cuts)|] gradfun[|1,6\n,cols(gradfun)|]=gradfun_cuts st_store(.,(tokens(gfunout)),gradfun) } end