capture program drop _all program define treatoprobit_work, eclass local gradient global neqs = 3+$ncut forvalues i = 1/$neqs { local gradient "`gradient' g`i'" } args todo b lnfj `gradient' H tempvar eta1 eta2 tempname r forvalues i = 1/$ncut { tempname c`i' } mleval `eta1' = `b', eq(1) mleval `eta2' = `b', eq(2) mleval `r' = `b', eq(3) scalar forvalues i = 1/$ncut { local j = 3+`i' mleval `c`i'' = `b', eq(`j') scalar } local rho = (tanh(`r')) 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' ll`i' ll`j' qui g double `m0`i'' = $ML_y1==0 & $ML_y2==`i' qui g double `m1`i'' = $ML_y1==1 & $ML_y2==`i' qui g double `ll`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 cutpts0r cutptst mat `cutpts0' = J(1,$ncut,0) forvalues i = 1/$ncut { mat `cutpts0'[1,`i'] = `c`i'' } mat `cutpts0r' = `neginf', `cutpts0', `posinf' mat `cutptst' = `neginf', 0, `posinf' tempname fv_tmp qui g double `fv_tmp'=0 local count 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 `++count' local k = `i'+1 local cut_h2 = `cutpts0r'[1,`k'] local cut_l2 = `cutpts0r'[1,`i'] qui replace `ll`count'' = `m`d'`i''*(((binorm(`cut_h1'-`eta1',`cut_h2'-`eta2',`rho'))- /// (binorm(`cut_l1'-`eta1',`cut_h2'-`eta2',`rho'))-(binorm(`cut_h1'-`eta1',`cut_l2'-`eta2',`rho'))+ /// (binorm(`cut_l1'-`eta1',`cut_l2'-`eta2',`rho')))) qui replace `fv_tmp' = `fv_tmp'+`ll`count'' } } qui replace `lnfj' = ln(`fv_tmp') if (`todo'==0 | `lnfj'==.) exit local gradfun forv i = 1/$neqs { tempvar g`i'_tmp qui g double `g`i'_tmp'=0 local gradfun "`gradfun' `g`i'_tmp'" } local r_init = `rho' mata: mata_gradfun_sobiprobit("`eta1'", "`eta2'", "`cutpts0r'", "`cutptst'", /// "`r_init'", "`imat'", "`lnfj'","`gradfun'") local gfun forvalues i = 1/$neqs { qui replace `g`i''=`g`i'_tmp' local gfun "`gfun' `g`i''" } if (`todo'==1) exit local hess forv i = 1/3 { forv j = `i'/$neqs { tempvar g`i'`j'_tmp qui g double `g`i'`j'_tmp'=0 local hess "`hess' `g`i'`j'_tmp'" } } forvalues i = 4/$neqs{ tempvar g`i'`i'_tmp qui g double `g`i'`i'_tmp'=0 local hess "`hess' `g`i'`i'_tmp'" } global offdiag = $neqs-1 forvalues i = 4/$offdiag { local j = `i'+1 tempvar g`i'`j'_tmp qui g double `g`i'`j'_tmp'=0 local hess "`hess' `g`i'`j'_tmp'" } mata: mata_hessfun_sobiprobit("`eta1'", "`eta2'", "`cutpts0r'", "`cutptst'", /// "`r_init'", "`imat'", "`lnfj'", "`touse'", "`gfun'", "`hess'") forv i = 1/$neqs { forv j= `i'/$neqs { tempname g`i'`j' } } mlmatsum `lnfj' `g11' = `g11_tmp', eq(1) forv i = 2/$neqs { mlmatsum `lnfj' `g1`i'' = `g1`i'_tmp', eq(1,`i') } mlmatsum `lnfj' `g22' = `g22_tmp', eq(2) forv i = 2/$neqs { mlmatsum `lnfj' `g2`i'' = `g2`i'_tmp', eq(2,`i') } mlmatsum `lnfj' `g33' = `g33_tmp', eq(3) forv i = 3/$neqs { mlmatsum `lnfj' `g3`i'' = `g3`i'_tmp', eq(3,`i') } forv i = 4/$neqs { forv j = 4/$neqs { if `i'==`j' { mlmatsum `lnfj' `g`i'`j'' = `g`i'`j'_tmp', eq(`i') } if `i'==`j'-1 { mlmatsum `lnfj' `g`i'`j'' = `g`i'`j'_tmp', eq(`i',`j') } else { continue } } } *form hessian matrix* tempname H1 H2 H3 A B C d O mat `A' = `g11', `g12', `g13'\ `g12'', `g22', `g23'\ `g13'', `g23'', `g33' mat `H1' = `g14' mat `H2' = `g24' mat `H3' = `g34' forvalues i = 5/$neqs { mat `H1' = nullmat(`H1'), `g1`i'' mat `H2' = nullmat(`H2'), `g2`i'' mat `H3' = nullmat(`H3'), `g3`i'' } mat `B' = nullmat(`H1') \ `H2' mat `B' = nullmat(`B') \ `H3' mat `d' = `g44' forv i = 5/$neqs { mat `d' = nullmat(`d'), `g`i'`i'' } mat `C' = diag(`d') global cutM1 = $ncut-1 mat `O' = J(rowsof(`C'),colsof(`C'),0) forv i = 1/$cutM1 { forv j = 2/$ncut { if `i'==`j'-1 { local m = `i'+3 local n = `j'+3 mat `O'[`j',`i']=`g`m'`n'' mat `O'[`i',`j']=`g`m'`n'' } else { continue } } } mat `C' = `C'+`O' *hessian mat `H' = `A', `B' \ `B'', `C' end mata: void mata_gradfun_sobiprobit(string scalar xb1, string scalar xb2, string scalar cutpts0r, /// string scalar cutptst, string scalar r, string scalar indmat, string scalar llike, /// string scalar gfunout) { n = st_nobs() real matrix gradfun, imat real colvector eta1, eta2, lnlike real scalar rho st_view(eta1,.,xb1) st_view(eta2,.,xb2) st_view(imat,.,tokens(indmat)) st_view(lnlike,.,llike) cg = cols(tokens(gfunout)) rho =(strtoreal(r)) cutpts0r = st_matrix(cutpts0r) cutptst = st_matrix(cutptst) fvi = exp(lnlike) //lnfvisum = sum(lnlike) //printf("Log Like = %9.2f\n", lnfvisum) ci = cols(imat) //printf("ci=%2.0f \n", ci) choices = ci/2 //printf("rho = %3.2f \n", rho) neqs = 3+(choices-1) gradfun = J(n,neqs,0) gradfun_1 = J(n,ci,0) gradfun_2 = J(n,ci,0) delta = sqrt(1-rho^2) //printf("delta = %3.2f \n", delta) /*derivatives wrt to beta1 and beta2*/ for(i=1;i<=2;i++) { cut_h1 = cutptst[1,i+1] cut_l1 = cutptst[1,i] A1p1 = cut_h1:-eta1 A1 = cut_l1:-eta1 for(j=1;j<=choices;j++){ k = ((i-1)*choices)+ j cut_h2 = cutpts0r[1,j+1] cut_l2 = cutpts0r[1,j] A2p1 = cut_h2:-eta2 A2 = cut_l2:-eta2 gradfun_1[.,k]=(normalden(A1p1):*(normal((A2p1:-rho:*A1p1):/(delta)))):- /// (normalden(A1):*(normal((A2p1:-rho:*A1):/(delta)))):- /// (normalden(A1p1):*(normal((A2:-rho:*A1p1):/(delta)))):+ /// (normalden(A1):*(normal((A2:-rho:*A1):/(delta)))) gradfun_2[.,k]=normalden(A2p1):*(normal((A1p1:-rho:*A2p1):/(delta))):- /// normalden(A2):*(normal((A1p1:-rho:*A2):/(delta))):- /// normalden(A2p1):*(normal((A1:-rho:*A2p1):/(delta))):+ /// normalden(A2):*(normal((A1:-rho:*A2):/(delta))) } } gradfun[.,1]=(-1):*rowsum(imat:*gradfun_1):/fvi gradfun[.,2]=(-1):*rowsum(imat:*gradfun_2):/fvi /*derivatives wrt r, rho*/ gradfun_3 = J(n,ci,0) gamma = (1-rho^2) nu = 1/sqrt(gamma) zeta = ((2*pi())*sqrt(1-rho^2))^-1 for(i=1;i<=2;i++) { cut_h1 = cutptst[1,i+1] cut_l1 = cutptst[1,i] A1p1 = cut_h1:-eta1 A1 = cut_l1:-eta1 for(j=1;j<=choices;j++){ k = ((i-1)*choices)+ j cut_h2 = cutpts0r[1,j+1] cut_l2 = cutpts0r[1,j] A2p1 = cut_h2:-eta2 A2 = cut_l2:-eta2 gradfun_3[.,k]=(exp((-.5):*(((A1p1:^2):+(A2p1:^2):-(rho:*A1p1:*A2p1:*2)):/gamma)):- /// exp((-.5):*(((A1:^2):+(A2p1:^2):-(rho:*A1:*A2p1:*2)):/gamma)):- /// exp((-.5):*(((A1p1:^2):+(A2:^2):-(rho:*A1p1:*A2:*2)):/gamma)):+ /// exp((-.5):*(((A1:^2):+(A2:^2):-(rho:*A1:*A2:*2)):/gamma))) //:*zeta } } rr = atanh(rho) //strtoreal(r) drho_dr = 1/(cosh(rr)^2) // derivative of rho wrt r gradscalar = drho_dr*zeta gradfun[.,3] = gradscalar:*((rowsum(imat:*gradfun_3)):/fvi) ncuts = choices-1 gradfun_4_tmp = J(n,ci,0) gradfun_4 = J(n,ncuts,.) for(i=1;i<=ncuts;i++){ for(j=1;j<=2;j++){ cut_h1 = cutptst[1,j+1] cut_l1 = cutptst[1,j] A1p1 = cut_h1:-eta1 A1 = cut_l1:-eta1 for(k=1;k<=choices;k++){ q = ((j-1)*choices)+ k delta_kj1 = (i==k) delta_kj = (i==k-1) cut_h2 = cutpts0r[1,k+1] cut_l2 = cutpts0r[1,k] A2p1 = cut_h2:-eta2 A2 = cut_l2:-eta2 gradfun_4_tmp[.,q]=((normal((A1p1:-rho:*A2p1):/delta):-normal((A1:-rho:*A2p1):/delta)):*normalden(A2p1):*delta_kj1):- /// ((normal((A1p1:-rho:*A2):/delta):-normal((A1:-rho:*A2):/delta)):*normalden(A2):*delta_kj) } } gradfun_4[.,i] = (rowsum(imat:*gradfun_4_tmp)):/fvi } gradfun[|1,4\n,cols(gradfun)|]=gradfun_4 st_store(.,(tokens(gfunout)),gradfun) } void mata_hessfun_sobiprobit(string scalar xb1, string scalar xb2, string scalar cutpts0r, /// string scalar cutptst, string scalar r, string scalar indmat, string scalar llike, /// string scalar touse, string scalar gfunin, string scalar hess) { n = st_nobs() real matrix gradfun, imat real colvector eta1, eta2, lnlike real scalar rho //touse = st_data(.,tokens(touse)) eta1=st_data(.,xb1,touse) eta2=st_data(.,xb2,touse) imat=st_data(.,tokens(indmat),touse) lnlike=st_data(.,llike,touse) cg = cols(tokens(gfunin)) gradfunc = st_data(., tokens(gfunin),touse) z = cols(gradfunc) /*gradients*/ g1 = gradfunc[.,1] g2 = gradfunc[.,2] g3 = gradfunc[.,3] g4 = gradfunc[|1,4\n,z|] /*matrix to store hessian*/ hess_c = cols(tokens(hess)) hessfun = J(n,hess_c,0) rho = strtoreal(r) cutpts0r = st_matrix(cutpts0r) cutptst = st_matrix(cutptst) fvi = exp(lnlike) fvi2 = fvi:^2 //lnfvisum = sum(lnlike) //printf("Log Like = %9.2f\n", lnfvisum) ci = cols(imat) choices = ci/2 neqs = 3+(choices-1) gradfun = J(n,neqs,0) gamma = 1-rho^2 delta = (sqrt(1-rho^2))^-1 /*g11*/ zeta = ((2*pi())*sqrt(1-rho^2))^-1 g11 = J(n,1,0) gh11 = J(n,ci,0) gh11_1 = J(n,ci,0) gh11_2 = J(n,ci,0) gh11_3 = J(n,ci,0) gh11_4 = J(n,ci,0) for(i=1;i<=2;i++) { cut_h1 = cutptst[1,i+1] cut_l1 = cutptst[1,i] A1p1 = (cut_h1:-eta1) A1 = (cut_l1:-eta1) for(j=1;j<=choices;j++){ k = ((i-1)*choices)+ j cut_h2 = cutpts0r[1,j+1] cut_l2 = cutpts0r[1,j] A2p1 = cut_h2:-eta2 A2 = cut_l2:-eta2 gh11_1[.,k]=zeta:*(exp((-.5):*(((A1p1:^2):+(A2p1:^2):-(rho:*A1p1:*A2p1:*2)):/gamma))):*(-rho):+ /// normal((A2p1:-rho:*A1p1):*delta):*normalden(A1p1):*(-A1p1) gh11_2[.,k]=zeta:*(exp((-.5):*(((A1:^2):+(A2p1:^2):-(rho:*A1:*A2p1:*2)):/gamma))):*(-rho):+ /// normal((A2p1:-rho:*A1):*delta):*normalden(A1):*(-A1) gh11_3[.,k]=zeta:*(exp((-.5):*(((A1p1:^2):+(A2:^2):-(rho:*A1p1:*A2:*2)):/gamma))):*(-rho):+ /// normal((A2:-rho:*A1p1):*delta):*normalden(A1p1):*(-A1p1) gh11_4[.,k]=zeta:*(exp((-.5):*(((A1:^2):+(A2:^2):-(rho:*A1:*A2:*2)):/gamma))):*(-rho):+ /// normal((A2:-rho:*A1):*delta):*normalden(A1):*(-A1) gh11[.,k] = gh11_1[.,k]:-gh11_2[.,k]:-gh11_3[.,k]:+gh11_4[.,k] } } g11 = ((rowsum(imat:*gh11)):/fvi - (g1:^2)) /*g12*/ zeta = ((2*pi())*sqrt(1-rho^2))^-1 g12 = J(n,1,0) gh12 = J(n,ci,0) gh12_1 = J(n,ci,0) gh12_2 = J(n,ci,0) gh12_3 = J(n,ci,0) gh12_4 = J(n,ci,0) for(i=1;i<=2;i++) { cut_h1 = cutptst[1,i+1] cut_l1 = cutptst[1,i] A1p1 = cut_h1:-eta1 A1 = cut_l1:-eta1 for(j=1;j<=choices;j++){ k = ((i-1)*choices)+ j cut_h2 = cutpts0r[1,j+1] cut_l2 = cutpts0r[1,j] A2p1 = cut_h2:-eta2 A2 = cut_l2:-eta2 /*gh12_1[.,k]=normalden(A1p1):*normalden((A2p1:-rho:*A1p1):*delta) gh12_2[.,k]=normalden(A1):*normalden((A2p1:-rho:*A1):*delta) gh12_3[.,k]=normalden(A1p1):*normalden((A2:-rho:*A1p1):*delta) gh12_4[.,k]=normalden(A1):*normalden((A2:-rho:*A1):*delta)*/ gh12_1[.,k]=exp((-.5):*(((A1p1:^2):+(A2p1:^2):-(rho:*A1p1:*A2p1:*2)):/gamma)) gh12_2[.,k]=exp((-.5):*(((A1:^2):+(A2p1:^2):-(rho:*A1:*A2p1:*2)):/gamma)) gh12_3[.,k]=exp((-.5):*(((A1p1:^2):+(A2:^2):-(rho:*A1p1:*A2:*2)):/gamma)) gh12_4[.,k]=exp((-.5):*(((A1:^2):+(A2:^2):-(rho:*A1:*A2:*2)):/gamma)) gh12[.,k] = (gh12_1[.,k]:-gh12_2[.,k]:-gh12_3[.,k]:+gh12_4[.,k]):*zeta } } g12 = ((rowsum(imat:*gh12)):/fvi):-((g1:*g2)) sigma = (1-rho^2)^(-3/2) rr = atanh(rho) drho_dr = 1/((cosh(rr))^2) // derivative of rho wrt r /*g13*/ g13 = J(n,1,0) gh13 = J(n,ci,0) gh13_1 = J(n,ci,0) gh13_2 = J(n,ci,0) gh13_3 = J(n,ci,0) gh13_4 = J(n,ci,0) for(i=1;i<=2;i++) { cut_h1 = cutptst[1,i+1] cut_l1 = cutptst[1,i] A1p1 = cut_h1:-eta1 A1 = cut_l1:-eta1 for(j=1;j<=choices;j++){ k = ((i-1)*choices)+ j cut_h2 = cutpts0r[1,j+1] cut_l2 = cutpts0r[1,j] A2p1 = cut_h2:-eta2 A2 = cut_l2:-eta2 gh13_1[.,k]=-normalden(A1p1):*normalden((A2p1:-rho:*A1p1):*delta):*((rho:*A2p1:-A1p1):*sigma) gh13_2[.,k]=-normalden(A1):*normalden((A2p1:-rho:*A1):*delta):*((rho:*A2p1:-A1):*sigma) gh13_3[.,k]=-normalden(A1p1):*normalden((A2:-rho:*A1p1):*delta):*((rho:*A2:-A1p1):*sigma) gh13_4[.,k]=-normalden(A1):*normalden((A2:-rho:*A1):*delta):*((rho:*A2:-A1):*sigma) gh13[.,k] = (gh13_1[.,k]:-gh13_2[.,k]:-gh13_3[.,k]:+gh13_4[.,k] ):*drho_dr } } g13 = ((rowsum(imat:*gh13)):/fvi):-((g1:*g3)) ncuts = choices-1 /*g14*/ g14 = J(n,ncuts,.) g14_tmp = J(n,ncuts,0) for(i=1;i<=ncuts;i++){ gh14 = J(n,ci,0) gh14_1 = J(n,ci,0) gh14_2 = J(n,ci,0) gh14_3 = J(n,ci,0) gh14_4 = J(n,ci,0) for(j=1;j<=2;j++){ cut_h1 = cutptst[1,j+1] cut_l1 = cutptst[1,j] //printf("cut_l1=%4.3f \n",cut_l1) A1p1 = cut_h1:-eta1 A1 = cut_l1:-eta1 for(k=1;k<=choices;k++){ q = ((j-1)*choices)+ k delta_kj1 = (i==k) delta_kj = (i==k-1) cut_h2 = cutpts0r[1,k+1] cut_l2 = cutpts0r[1,k] A2p1 = cut_h2:-eta2 A2 = cut_l2:-eta2 gh14_1[.,q]=((normalden(A1p1):*normalden((A2p1:-rho:*A1p1):*delta)):*delta):*delta_kj1 gh14_2[.,q]=((normalden(A1):*normalden((A2p1:-rho:*A1):*delta)):*delta):*delta_kj1 gh14_3[.,q]=((normalden(A1p1):*normalden((A2:-rho:*A1p1):*delta)):*delta):*delta_kj gh14_4[.,q]=((normalden(A1):*normalden((A2:-rho:*A1):*delta)):*delta):*delta_kj gh14[.,q] = (-gh14_1[.,q]:+gh14_2[.,q]:+gh14_3[.,q]:-gh14_4[.,q]) //minus sign for zeros } } g14_tmp[.,i] = (rowsum(imat:*gh14)):/fvi } g14= g14_tmp:-(g1:*g4) /*g22*/ g22 = J(n,1,0) gh22 = J(n,ci,0) gh22_1 = J(n,ci,0) gh22_2 = J(n,ci,0) gh22_3 = J(n,ci,0) gh22_4 = J(n,ci,0) for(i=1;i<=2;i++) { cut_h1 = cutptst[1,i+1] cut_l1 = cutptst[1,i] A1p1 = cut_h1:-eta1 A1 = cut_l1:-eta1 for(j=1;j<=choices;j++){ k = ((i-1)*choices)+ j cut_h2 = cutpts0r[1,j+1] cut_l2 = cutpts0r[1,j] A2p1 = cut_h2:-eta2 A2 = cut_l2:-eta2 gh22_1[.,k]=zeta:*(exp((-.5):*(((A1p1:^2):+(A2p1:^2):-(rho:*A1p1:*A2p1:*2)):/gamma))):*(-rho):+ /// normal((A1p1:-rho:*A2p1):*delta):*normalden(A2p1):*(-A2p1) gh22_2[.,k]=zeta:*(exp((-.5):*(((A1:^2):+(A2p1:^2):-(rho:*A1:*A2p1:*2)):/gamma))):*(-rho):+ /// normal((A1:-rho:*A2p1):*delta):*normalden(A2p1):*(-A2p1) gh22_3[.,k]=zeta:*(exp((-.5):*(((A1p1:^2):+(A2:^2):-(rho:*A1p1:*A2:*2)):/gamma))):*(-rho):+ /// normal((A1p1:-rho:*A2):*delta):*normalden(A2):*(-A2) gh22_4[.,k]=zeta:*(exp((-.5):*(((A1:^2):+(A2:^2):-(rho:*A1:*A2:*2)):/gamma))):*(-rho):+ /// normal((A1:-rho:*A2):*delta):*normalden(A2):*(-A2) gh22[.,k] = gh22_1[.,k]:-gh22_2[.,k]:-gh22_3[.,k]:+gh22_4[.,k] } } g22 = ((rowsum(imat:*gh22)):/fvi - (g2:^2)) /*g23*/ //printf("drho_dr = %3.2f \n", drho_dr) sigma = (1-rho^2)^-1 g23 = J(n,1,0) gh23 = J(n,ci,0) gh23_1 = J(n,ci,0) gh23_2 = J(n,ci,0) gh23_3 = J(n,ci,0) gh23_4 = J(n,ci,0) for(i=1;i<=2;i++) { cut_h1 = cutptst[1,i+1] cut_l1 = cutptst[1,i] A1p1 = cut_h1:-eta1 A1 = cut_l1:-eta1 for(j=1;j<=choices;j++){ k = ((i-1)*choices)+ j cut_h2 = cutpts0r[1,j+1] cut_l2 = cutpts0r[1,j] A2p1 = cut_h2:-eta2 A2 = cut_l2:-eta2 gh23_1[.,k]=-zeta:*(exp((-.5):*(((A1p1:^2):+(A2p1:^2):-(rho:*A1p1:*A2p1:*2)):/gamma))):*((rho:*A1p1:-A2p1):*sigma) gh23_2[.,k]=-zeta:*(exp((-.5):*(((A1:^2):+(A2p1:^2):-(rho:*A1:*A2p1:*2)):/gamma))):*((rho:*A1:-A2p1):*sigma) gh23_3[.,k]=-zeta:*(exp((-.5):*(((A1p1:^2):+(A2:^2):-(rho:*A1p1:*A2:*2)):/gamma))):*((rho:*A1p1:-A2):*sigma) gh23_4[.,k]=-zeta:*(exp((-.5):*(((A1:^2):+(A2:^2):-(rho:*A1:*A2:*2)):/gamma))):*((rho:*A1:-A2):*sigma) gh23[.,k] = (gh23_1[.,k]:-gh23_2[.,k]:-gh23_3[.,k]:+gh23_4[.,k]):*drho_dr } } g23 = ((rowsum(imat:*gh23)):/fvi):-((g2:*g3)) /*g24*/ g24 = J(n,ncuts,.) for(i=1;i<=ncuts;i++){ gh24 = J(n,ci,0) gh24_1 = J(n,ci,0) gh24_2 = J(n,ci,0) gh24_3 = J(n,ci,0) gh24_4 = J(n,ci,0) for(j=1;j<=2;j++){ cut_h1 = cutptst[1,j+1] cut_l1 = cutptst[1,j] A1p1 = cut_h1:-eta1 A1 = cut_l1:-eta1 for(k=1;k<=choices;k++){ q = ((j-1)*choices)+ k delta_kj1 = (i==k) delta_kj = (i==k-1) cut_h2 = cutpts0r[1,k+1] cut_l2 = cutpts0r[1,k] A2p1 = cut_h2:-eta2 A2 = cut_l2:-eta2 gh24_1[.,q]=((normalden(A1p1):*normalden((A2p1:-rho:*A1p1):*delta):*-rho*delta):+ /// (normal((A1p1:-rho:*A2p1):*delta):*normalden(A2p1):*(-A2p1))):*delta_kj1 gh24_2[.,q]=((normalden(A1):*normalden((A2p1:-rho:*A1):*delta):*-rho*delta):+ /// (normal((A1:-rho:*A2p1):*delta):*normalden(A2p1):*(-A2p1))):*delta_kj1 gh24_3[.,q]=((normalden(A1p1):*normalden((A2:-rho:*A1p1):*delta):*-rho*delta)+ /// (normal((A1p1:-rho:*A2):*delta):*normalden(A2):*(-A2))):*delta_kj gh24_4[.,q]=((normalden(A1):*normalden((A2:-rho:*A1):*delta):*-rho*delta):+ /// (normal((A1:-rho:*A2):*delta):*normalden(A2):*(-A2))):*delta_kj gh24[.,q] = -gh24_1[.,q]:+gh24_2[.,q]:+gh24_3[.,q]:-gh24_4[.,q] } } g24[.,i] = ((rowsum(imat:*gh24)):/fvi) } g24 = g24:-((g2:*g4)) /*g33*/ rr=atanh(rho) zeta = ((2*pi())*sqrt(1-rho^2))^-1 zetainv = zeta^-1 d2rho_d2r = -2*sinh(rr)/cosh(rr)^3 //-2*tanh(rr)/cosh(rr)^3 //8*exp(2*`r')*(1-exp(2*`r'))/(1+exp(2*`r'))^3 d2zeta_d2r = rho/(2*pi()*((1-rho^2)^(3/2))) //((-2*pi()*rho)/((1-rho^2)^(3/2))) drho_dr = 1/(cosh(rr)^2) gradscalar = drho_dr*zeta d2cons_d2r = zeta*d2rho_d2r+drho_dr^2*d2zeta_d2r gamma = (1-rho^2) g3k = (g3:*(zetainv):*(drho_dr^-1)) g33 = J(n,1,0) gh33 = J(n,ci,0) gh33_1 = J(n,ci,0) gh33_2 = J(n,ci,0) gh33_3 = J(n,ci,0) gh33_4 = J(n,ci,0) for(i=1;i<=2;i++) { cut_h1 = cutptst[1,i+1] cut_l1 = cutptst[1,i] A1p1 = cut_h1:-eta1 A1 = cut_l1:-eta1 for(j=1;j<=choices;j++){ k = ((i-1)*choices)+ j cut_h2 = cutpts0r[1,j+1] cut_l2 = cutpts0r[1,j] A2p1 = cut_h2:-eta2 A2 = cut_l2:-eta2 gh33_1[.,k]=(exp((-.5):*(((A1p1:^2):+(A2p1:^2):-(rho:*A1p1:*A2p1:*2)):/gamma)):* /// (-.5:*((((-2:*A1p1:*A2p1):*gamma):+2:*rho:*((A1p1:^2)+(A2p1:^2)-2:*rho:*(A1p1:*A2p1))):/(gamma:^2)))) gh33_2[.,k]=(exp((-.5):*(((A1:^2):+(A2p1:^2):-(rho:*A1:*A2p1:*2)):/gamma)):* /// (-.5:*((((-2:*A1:*A2p1):*gamma):+2:*rho:*((A1:^2)+(A2p1:^2)-2:*rho:*(A1:*A2p1))):/(gamma:^2)))) gh33_3[.,k]=(exp((-.5):*(((A1p1:^2):+(A2:^2):-(rho:*A1p1:*A2:*2)):/gamma)):* /// (-.5:*((((-2:*A1p1:*A2):*gamma):+2:*rho:*((A1p1:^2)+(A2:^2)-2:*rho:*(A1p1:*A2))):/(gamma:^2)))) gh33_4[.,k]=(exp((-.5):*(((A1:^2):+(A2:^2):-(rho:*A1:*A2:*2)):/gamma)):* /// (-.5:*((((-2:*A1:*A2):*gamma):+2:*rho:*((A1:^2)+(A2:^2)-2:*rho:*(A1:*A2))):/(gamma:^2)))) gh33[.,k] = (gh33_1[.,k]:-gh33_2[.,k]:-gh33_3[.,k]:+gh33_4[.,k]):*drho_dr*gradscalar } } g33 = ((((rowsum(imat:*gh33)):/fvi):-((g3:^2)))):+g3k:*d2cons_d2r /*g34*/ // printf("g34") g34 = J(n,ncuts,.) for(i=1;i<=ncuts;i++){ gh34 = J(n,ci,0) gh34_1 = J(n,ci,0) gh34_2 = J(n,ci,0) gh34_3 = J(n,ci,0) gh34_4 = J(n,ci,0) for(j=1;j<=2;j++){ cut_h1 = cutptst[1,j+1] cut_l1 = cutptst[1,j] A1p1 = cut_h1:-eta1 A1 = cut_l1:-eta1 for(k=1;k<=choices;k++){ q = ((j-1)*choices)+ k delta_kj1 = (i==k) delta_kj = (i==k-1) cut_h2 = cutpts0r[1,k+1] cut_l2 = cutpts0r[1,k] A2p1 = cut_h2:-eta2 A2 = cut_l2:-eta2 gh34_1[.,q]=(exp((-.5):*(((A1p1:^2):+(A2p1:^2):-(rho:*A1p1:*A2p1:*2)):/gamma)):* /// (-.5:*(((2:*A2p1):-(2:*rho:*A1p1)):/gamma))):*delta_kj1 gh34_2[.,q]=(exp((-.5):*(((A1:^2):+(A2p1:^2):-(rho:*A1:*A2p1:*2)):/gamma)):* /// (-.5:*(((2:*A2p1):-(2:*rho:*A1)):/gamma))):*delta_kj1 gh34_3[.,q]=(exp((-.5):*(((A1p1:^2):+(A2:^2):-(rho:*A1p1:*A2:*2)):/gamma)):* /// (-.5:*(((2:*A2):-(2:*rho:*A1p1)):/gamma))):*delta_kj gh34_4[.,q]=(exp((-.5):*(((A1:^2):+(A2:^2):-(rho:*A1:*A2:*2)):/gamma)):* /// (-.5:*(((2:*A2):-(2:*rho:*A1)):/gamma))):*delta_kj gh34[.,q] = (gh34_1[.,q]:-gh34_2[.,q]:-gh34_3[.,q]:+gh34_4[.,q]):*(zeta*drho_dr) } } g34[.,i] = ((rowsum(imat:*gh34)):/fvi) } g34 = g34:-((g3:*g4)) /*g44*/ //printf("g44") g44_diag = J(n,ncuts,.) for(i=1;i<=ncuts;i++){ gh44 = J(n,ci,0) gh44_1 = J(n,ci,0) gh44_2 = J(n,ci,0) gh44_3 = J(n,ci,0) gh44_4 = J(n,ci,0) for(j=1;j<=2;j++){ cut_h1 = cutptst[1,j+1] cut_l1 = cutptst[1,j] A1p1 = cut_h1:-eta1 A1 = cut_l1:-eta1 for(k=1;k<=choices;k++){ q = ((j-1)*choices)+ k delta_kj1 = (i==k) delta_kj = (i==k-1) cut_h2 = cutpts0r[1,k+1] cut_l2 = cutpts0r[1,k] A2p1 = cut_h2:-eta2 A2 = cut_l2:-eta2 gh44_1[.,q]=(normalden(A2p1):*normalden((A1p1:-rho:*A2p1):*delta):*(-rho):+ /// normal((A1p1:-rho:*A2p1):*delta):*normalden(A2p1):*(-A2p1)):*delta_kj1 gh44_2[.,q]=(normalden(A2p1):*normalden((A1:-rho:*A2p1):*delta):*(-rho):+ /// normal((A1:-rho:*A2p1):*delta):*normalden(A2p1):*(-A2p1)):*delta_kj1 gh44_3[.,q]=(normalden(A2):*normalden((A1p1:-rho:*A2):*delta):*(-rho):+ /// normal((A1p1:-rho:*A2):*delta):*normalden(A2):*(-A2)):*delta_kj gh44_4[.,q]=(normalden(A2):*normalden((A1:-rho:*A2):*delta):*(-rho):+ /// normal((A1:-rho:*A2):*delta):*normalden(A2):*(-A2)):*delta_kj gh44[.,q] = gh44_1[.,q]:-gh44_2[.,q]:-gh44_3[.,q]:+gh44_4[.,q] } } g44_diag[.,i] = ((rowsum(imat:*gh44)):/fvi) } g44_diag = g44_diag:-((g4:^2)) /*g44 off diagonal*/ zeros = J(n,1,0) offdiag=ncuts-1 g44_od = J(n,offdiag,.) for(i=1;i<=offdiag;i++){ t=i+1 gh44 = J(n,ci,0) gh44_1 = J(n,ci,0) gh44_2 = J(n,ci,0) gh44_3 = J(n,ci,0) gh44_4 = J(n,ci,0) for(j=1;j<=2;j++){ cut_h1 = cutptst[1,j+1] cut_l1 = cutptst[1,j] A1p1 = cut_h1:-eta1 A1 = cut_l1:-eta1 for(k=1;k<=choices;k++){ q = ((j-1)*choices)+ k delta_kj1 = (i==k) delta_kj = (i==k-1) cut_h2 = cutpts0r[1,k+1] cut_l2 = cutpts0r[1,k] A2p1 = cut_h2:-eta2 A2 = cut_l2:-eta2 gh44_1[.,q]=(normalden(A2p1):*normalden((A1p1:-rho:*A2p1):*delta):*(-rho*delta):+ /// normal((A1p1:-rho:*A2p1):*delta):*normalden(A2p1):*(-A2p1)):*delta_kj gh44_2[.,q]=(normalden(A2p1):*normalden((A1:-rho:*A2p1):*delta):*(-rho*delta):+ /// normal((A1:-rho:*A2p1):*delta):*normalden(A2p1):*(-A2p1)):*delta_kj gh44_3[.,q]=(normalden(A2):*normalden((A1p1:-rho:*A2):*delta):*(-rho*delta):+ /// normal((A1p1:-rho:*A2):*delta):*normalden(A2):*(-A2)):*delta_kj1 gh44_4[.,q]=(normalden(A2):*normalden((A1:-rho:*A2):*delta):*(-rho*delta):+ /// normal((A1:-rho:*A2):*delta):*normalden(A2):*(-A2)):*delta_kj1 gh44[.,q] = gh44_1[.,q]:-gh44_2[.,q]:-gh44_3[.,q]:+gh44_4[.,q] //- + + - } } g44_od[.,i] = ((rowsum(imat:*gh44)):/fvi) g44_od[.,i] = zeros:-(g4[.,i]:*g4[.,t]) } g44 = g44_diag, g44_od zeros = J(n,1,0) hessfun = g11, g12, g13, g14, g22, g23, g24, g33, g34, g44 st_store(.,tokens(hess),hessfun) } end