program t2t2innerv3
version 13
args lnf mu theta lsigma
tempvar sigma 
quietly gen double `sigma' = exp(`lsigma')
quietly replace `lnf' = ln((((1/2)*(1 + ((`mu') + sqrt((-1 + ($ml_y1+0.000000001) + exp(`theta')*($ml_y1+0.000000001))^2/(1 + (-1 + exp(`theta'))*($ml_y1+0.000000001))^2)/(sqrt(2)*sqrt(-((exp(`theta')*(-1 + ($ml_y1+0.000000001))*($ml_y1+0.000000001))/(1 + (-1 + exp(`theta'))*($ml_y1+0.000000001))^2))))/((`sigma')*sqrt(2 + ((`mu') + sqrt((-1 + ($ml_y1+0.000000001) + exp(`theta')*($ml_y1+0.000000001))^2/(1 + (-1 + exp(`theta'))*($ml_y1+0.000000001))^2)/(sqrt(2)*sqrt(-((exp(`theta')*(-1 + ($ml_y1+0.000000001))*($ml_y1+0.000000001))/(1 + (-1 + exp(`theta'))*($ml_y1+0.000000001))^2))))/(`sigma')^2)))) - ((1/2)*(1 + ((`mu') + sqrt((-1 + $ML_y1 + exp(`theta')*$ML_y1)^2/(1 + (-1 + exp(`theta'))*$ML_y1)^2)/(sqrt(2)*sqrt(-((exp(`theta')*(-1 + $ML_y1)*$ML_y1)/(1 + (-1 + exp(`theta'))*$ML_y1)^2))))/((`sigma')*sqrt(2 + ((`mu') + sqrt((-1 + $ML_y1 + exp(`theta')*$ML_y1)^2/(1 + (-1 + exp(`theta'))*$ML_y1)^2)/(sqrt(2)*sqrt(-((exp(`theta')*(-1 + $ML_y1)*$ML_y1)/(1 + (-1 + exp(`theta'))*$ML_y1)^2))))/(`sigma')^2)))))/0.000000001) /// 
   if (exp(`theta')*$ML_y1)/(1 + (-1 + exp(`theta'))*$ML_y1) > 0 & (exp(`theta')*$ML_y1)/(1 + (-1 + exp(`theta'))*$ML_y1) < 0.5
quietly replace `lnf' = ln((((1/2)*(1 + ((`mu') - sqrt((-1 + ($ml_y1+0.000000001) + exp(`theta')*($ml_y1+0.000000001))^2/(1 + (-1 + exp(`theta'))*($ml_y1+0.000000001))^2)/(sqrt(2)*sqrt(-((exp(`theta')*(-1 + ($ml_y1+0.000000001))*($ml_y1+0.000000001))/(1 + (-1 + exp(`theta'))*($ml_y1+0.000000001))^2))))/((`sigma')*sqrt(2 + ((`mu') - sqrt((-1 + ($ml_y1+0.000000001) + exp(`theta')*($ml_y1+0.000000001))^2/(1 + (-1 + exp(`theta'))*($ml_y1+0.000000001))^2)/(sqrt(2)*sqrt(-((exp(`theta')*(-1 + ($ml_y1+0.000000001))*($ml_y1+0.000000001))/(1 + (-1 + exp(`theta'))*($ml_y1+0.000000001))^2))))/(`sigma')^2)))) - ((1/2)*(1 + ((`mu') - sqrt((-1 + $ML_y1 + exp(`theta')*$ML_y1)^2/(1 + (-1 + exp(`theta'))*$ML_y1)^2)/(sqrt(2)*sqrt(-((exp(`theta')*(-1 + $ML_y1)*$ML_y1)/(1 + (-1 + exp(`theta'))*$ML_y1)^2))))/((`sigma')*sqrt(2 + ((`mu') - sqrt((-1 + $ML_y1 + exp(`theta')*$ML_y1)^2/(1 + (-1 + exp(`theta'))*$ML_y1)^2)/(sqrt(2)*sqrt(-((exp(`theta')*(-1 + $ML_y1)*$ML_y1)/(1 + (-1 + exp(`theta'))*$ML_y1)^2))))/(`sigma')^2)))))/0.000000001) /// 
   if (exp(`theta')*$ML_y1)/(1 + (-1 + exp(`theta'))*$ML_y1) >= 0.5 & (exp(`theta')*$ML_y1)/(1 + (-1 + exp(`theta'))*$ML_y1) < 1
quietly replace `lnf' = ln(`sigma'^2*exp(`theta')) if $ML_y1 ==0
quietly replace `lnf' = ln(`sigma'^2/(exp(`theta'))) if $ML_y1 ==1
end