program t2t2outerv3
version 13
args lnf mu theta lsigma
tempvar sigma 
quietly gen double `sigma' = exp(`lsigma')
quietly replace `lnf' = ln(((exp(`theta')*((`mu') - (`sigma')*sqrt(2 + ((`mu') + sqrt((1 - 2*($ML_y1 + 0.000001))^2)/(sqrt(2)*sqrt((1 - ($ML_y1 + 0.000001))*($ML_y1 + 0.000001))))^2/(`sigma')^2) + sqrt((1 - 2*($ML_y1 + 0.000001))^2)/(sqrt(2)*sqrt((1 - ($ML_y1 + 0.000001))*($ML_y1 + 0.000001)))))/((-1 + exp(`theta'))*(`mu') - (-1 + exp(`theta'))*(-(sqrt((1 - 2*($ML_y1 + 0.000001))^2)/(sqrt(2)*sqrt((1 - ($ML_y1 + 0.000001))*($ML_y1 + 0.000001))))) - (1 + exp(`theta'))*(`sigma')*sqrt(((`mu')^2 + 2*(`sigma')^2 - 2*(`mu')*(-(sqrt((1 - 2*($ML_y1 + 0.000001))^2)/(sqrt(2)*sqrt((1 - ($ML_y1 + 0.000001))*($ML_y1 + 0.000001))))) + (sqrt((1 - 2*($ML_y1 + 0.000001))^2)/(sqrt(2)*sqrt((1 - ($ML_y1 + 0.000001))*($ML_y1 + 0.000001))))^2)/(`sigma')^2)) - (exp(`theta')*((`mu') - (`sigma')*sqrt(2 + ((`mu') + sqrt((1 - 2*$ML_y1 )^2)/(sqrt(2)*sqrt((1 - $ML_y1 )*$ML_y1 )))^2/(`sigma')^2) + sqrt((1 - 2*$ML_y1 )^2)/(sqrt(2)*sqrt((1 - $ML_y1 )*$ML_y1 ))))/((-1 + exp(`theta'))*(`mu') - (-1 + exp(`theta'))*(-(sqrt((1 - 2*$ML_y1 )^2)/(sqrt(2)*sqrt((1 - $ML_y1 )*$ML_y1 )))) - (1 + exp(`theta'))*(`sigma')*sqrt(((`mu')^2 + 2*(`sigma')^2 - 2*(`mu')*(-(sqrt((1 - 2*$ML_y1 )^2)/(sqrt(2)*sqrt((1 - $ML_y1 )*$ML_y1 )))) + (sqrt((1 - 2*$ML_y1 )^2)/(sqrt(2)*sqrt((1 - $ML_y1 )*$ML_y1 )))^2)/(`sigma')^2)))/0.000001) /// 
   if $ML_y1 > 0 & $ML_y1 < 0.5
quietly replace `lnf' = ln(((exp(`theta')*((`mu') - (`sigma')*sqrt(2 + ((`mu') - sqrt((1 - 2*($ML_y1 + 0.000001))^2)/(sqrt(2)*sqrt((1 - ($ML_y1 + 0.000001))*($ML_y1 + 0.000001))))^2/(`sigma')^2) - sqrt((1 - 2*($ML_y1 + 0.000001))^2)/(sqrt(2)*sqrt((1 - ($ML_y1 + 0.000001))*($ML_y1 + 0.000001)))))/((-1 + exp(`theta'))*(`mu') - (-1 + exp(`theta'))*(sqrt((1 - 2*($ML_y1 + 0.000001))^2)/(sqrt(2)*sqrt((1 - ($ML_y1 + 0.000001))*($ML_y1 + 0.000001)))) - (1 + exp(`theta'))*(`sigma')*sqrt(((`mu')^2 + 2*(`sigma')^2 - 2*(`mu')*(sqrt((1 - 2*($ML_y1 + 0.000001))^2)/(sqrt(2)*sqrt((1 - ($ML_y1 + 0.000001))*($ML_y1 + 0.000001)))) + (sqrt((1 - 2*($ML_y1 + 0.000001))^2)/(sqrt(2)*sqrt((1 - ($ML_y1 + 0.000001))*($ML_y1 + 0.000001))))^2)/(`sigma')^2)) - (exp(`theta')*((`mu') - (`sigma')*sqrt(2 + ((`mu') - sqrt((1 - 2*$ML_y1 )^2)/(sqrt(2)*sqrt((1 - $ML_y1 )*$ML_y1 )))^2/(`sigma')^2) - sqrt((1 - 2*$ML_y1 )^2)/(sqrt(2)*sqrt((1 - $ML_y1 )*$ML_y1 ))))/((-1 + exp(`theta'))*(`mu') - (-1 + exp(`theta'))*(sqrt((1 - 2*$ML_y1 )^2)/(sqrt(2)*sqrt((1 - $ML_y1 )*$ML_y1 ))) - (1 + exp(`theta'))*(`sigma')*sqrt(((`mu')^2 + 2*(`sigma')^2 - 2*(`mu')*(sqrt((1 - 2*$ML_y1 )^2)/(sqrt(2)*sqrt((1 - $ML_y1 )*$ML_y1 ))) + (sqrt((1 - 2*$ML_y1 )^2)/(sqrt(2)*sqrt((1 - $ML_y1 )*$ML_y1 )))^2)/(`sigma')^2)))/0.000001) /// 
   if $ML_y1 >= 0.5 & $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