/* fig2_2.do

   Figure 2.2 (from Exercise 2.3)

   programs required
       binscrn1.ado  (v 1.0.4 or later)
       bvnellip.ado  (v 1.1.0 or later)

   data required: est1.dta  (CASS exercise stress test data)

   i)   elliptical 95% conf region for (DLR+,DLR-) based on
         Var(ln(DLR+),ln(DLR-)) per Simel (?)
          (is Cov also per Simel??, or just var ?)

   ii)  transformation of joint confidence region for (FPR,TPR)
             TPR/FPR     => DLR+
         (1-TPR)/(1-FPR) => DLR-

   iii) the rectangle

           TPR_l   TPR_u       1-TPR_u  1-TPR_l
         ( ----- , ----- ) x ( -------, ------- )
           FPR_u   FPR_l       1-FPR_l  1-FPR_u


   last update:  02 May 2003

*/
version 7
set more off
cap semt_profile
cap log close
pause on

use ${semt_data}est1,clear
preserve

tempname sigm_pn cov Var Lndlr

qui binscrn1 y1 d

* vector for ln(dlr) and var-cov matrix Var[ln(dlr)]

sca `cov' = -(1/r(nd) + 1/r(ndb))
matrix `Var' = (r(selndlrp)^2, `cov'\ `cov', r(selndlrn)^2)
matrix `Lndlr' = (ln(r(dlrpos))\ln(r(dlrneg)))

local df = r(nd) + r(ndb) - 2

matrix l `Var'

matr l `Lndlr'

di "df: " `df'
* pause

bvnellip, pl(lndlrpci lndlrnci) cov(`Var') beta(`Lndlr') df(`df') sort(pltordr1) replace

gen dlrpci = exp(lndlrpci)
gen dlrnci = exp(lndlrnci)

keep if pltordr1~=.

* pause

gen byte expn = _n==1
expand 2 if expn
gen dlrneg = .
gen dlrpos = .
bys expn : replace dlrpos = exp(`Lndlr'[1,1]) if expn==1 & _n==1
bys expn : replace dlrpci = exp(`Lndlr'[1,1]) if expn==1 & _n==1
bys expn : replace dlrneg = exp(`Lndlr'[2,1]) if expn==1 & _n==1
bys expn : replace dlrnci = .                 if expn==1 & _n==1

sort pltordr1

#delimit ;
gr dlrneg dlrnci dlrpci, s(p.) c(.l)
   yscale(.2,.4) ylabel(.2,.3,.4)
   xscale(2,4) xlab(2(.5)4)
;
#delimit cr
pause q <return> to continue . . .

keep pltordr1 dlrnci dlrpci dlrneg dlrpos
order pltordr1 dlrnci dlrpci dlrneg dlrpos
sort pltordr1
tempfile f1
save `f1'

replace dlrpci = . if dlrnci==.
save ${semt_log}fig2_2i,replace   /* save file for conversion to SigmaPlot */
***********************
* ii)

restore, preserve

qui binscrn1 y1 d

local nx = 100

sca tp_lb = r(tpj_lb)
sca tp_ub = r(tpj_ub)
sca fp_lb = r(fpj_lb)
sca fp_ub = r(fpj_ub)

drop _all

set obs 4
gen byte intv = _n
expand `nx'
bys intv: gen pltordr = _n
sort intv pltordr

gen     tpr = tp_ub if intv==1
replace tpr = tp_lb if intv==3
by intv: replace tpr = tp_ub - ((_n-1) * ( (tp_ub-tp_lb)/`nx') ) if intv==2
by intv: replace tpr = tp_lb + ((_n-1) * ( (tp_ub-tp_lb)/`nx') ) if intv==4

gen     fpr = fp_ub if intv==2
replace fpr = fp_lb if intv==4
by intv: replace fpr = fp_lb + ((_n-1) * ( (fp_ub-fp_lb)/`nx') ) if intv==1
by intv: replace fpr = fp_ub - ((_n-1) * ( (fp_ub-fp_lb)/`nx') ) if intv==3

expand 2 if _n==1
replace intv=4 if _n==1
replace pltordr = `nx' + 1 if _n==1
sort intv pltordr
gen pltordr2 = _n

gen dlrpci = tpr/fpr
gen dlrnci2 = (1-tpr)/(1-fpr)

keep pltordr2 dlr*
order pltordr2 dlrnci2 dlrpci
save ${semt_log}fig2_2ii,replace   /* save file for conversion to SigmaPlot */
tempfile f2
save `f2'
drop _all  /* the needed scalars remain */
************************************************
* iii)

set obs 4
gen int pltordr3 = _n
gen     dlrnci3 = (1-tp_ub)/(1-fp_lb) if _n==1 | _n==2
replace dlrnci3 = (1-tp_lb)/(1-fp_ub) if _n==3 | _n==4
gen     dlrpci  = (tp_lb)/(fp_ub) if _n==2 | _n==3
replace dlrpci  = (tp_ub)/(fp_lb) if _n==1 | _n==4
expand 2 if pltordr3==1
sort pltordr3
replace pltordr3 = 5 if _n==1
sort pltordr3

#delimit ;
gr dlrnci3 dlrpci, s(.) c(l)
   yscale(.2,.4) ylabel(.2,.3,.4)
   xscale(2,4) xlab(2(.5)4)
;
#delimit cr
pause q <return> to continue . . .


************************************************
append using `f1'
append using `f2'

la var dlrneg

sort pltordr3 pltordr2 pltordr1

#delimit ;
gr dlrneg dlrnci dlrnci2 dlrnci3 dlrpci, s(p...) c(.lll)
   yscale(.2,.4) ylabel(.2,.3,.4)
   l1("DLR-") b1("      DLR+") b2(" ")
   t1(95% joint confidence regions)
   key2(pen(3) c(l) "BVN (ln(DLR+),ln(DLR-))")
   key3(pen(4) c(l) "fr. joint reg. (FPR,TPR)")
   key4(pen(5) c(l) "from univar. CI's TPR, FPR")
   xscale(2,4) xlab(2(.5)4) saving(${semt_log}fig2_2,replace)
;
#delimit cr

gprefs set custom1 pen5_color 3 121 255
gprefs set custom1 pen1_thick 3
gprefs set custom1 pen2_thick 4
gprefs set custom1 pen3_thick 4
gprefs set custom1 pen4_thick 4
gprefs set custom1 pen5_thick 4

translator set gph2wmf scheme custom1
translate ${semt_log}fig2_2.gph fig2_2stcolor.wmf, /* translator(gph2wmf) */ replace