*! version 1  18Aug2000 by Duolao Wang
program define koopmani, rclass
	version 6
	
	gettoken a 0 : 0, parse(" ,")
	gettoken b 0 : 0, parse(" ,")
	gettoken c 0 : 0, parse(" ,")
	gettoken d 0 : 0, parse(" ,")

	confirm integer number `a'
	confirm integer number `b'
	confirm integer number `c'
	confirm integer number `d'

	if `a'<0 | `b'<0 | `c'<0 | `d'<0 { 
		di in red "negative numbers invalid"
		exit 498
	}

	syntax [,Level(real 95)]

	if `level'<.01 | `level'>99.99 { 
		local level 95
	}
	di

	scalar x=`a'
	scalar m=`b'
	scalar y=`c'
	scalar n=`d'
	scalar cinv = invchi(1,1-`level'/100)

/* determine interval via Koopman's Chi-squared method */

	scalar   p1=x/m
	scalar   p2=y/n
	scalar   T = p1/p2

/* lower bound - first handle extreme values */

	if x == 0 & y == 0 {scalar chlower = 0}
	else {
	
		if x == m {scalar xtemp = x - 0.5}
		else {scalar xtemp = x}

		if y == n {scalar ytemp = y - 0.5}
 		else {scalar ytemp = y}
	}

	scalar  bisect1 = 0.0001
	scalar  bisect2 = T
	scalar  bisect = (bisect1 + bisect2)/2
	scalar  tol = 1

	local i 1

	while `i'<=50 & abs(tol) > 0.0001 {

		scalar p1hat = (bisect*(m+ytemp) + x + n -  /*
 		*/ sqrt(  (bisect*(m+ytemp)+xtemp+n)^2 - (4*bisect*(m+n)*(xtemp+ytemp))  ) ) /(2*(m+n))
 		scalar p2hat = p1hat/bisect

		if x == m & y == n & T <=1 {scalar  Utheta0 = m*(1-bisect)/bisect}
    		else { 
		scalar Utheta0 =(((xtemp-m*p1hat)^2) / (m*p1hat*(1-p1hat))) + (((ytemp-n*p2hat)^2) / (n*p2hat*(1-p2hat)))
      	} 

		if Utheta0 < cinv {scalar bisect2 = bisect}
		else {scalar bisect1 = bisect}

		scalar bisect = (bisect1+bisect2)/2
		scalar tol = Utheta0 - cinv

		scalar chlower = bisect
		local i=`i'+1
	}

/* upper bound - first handle extreme values */

	if x == 0 & y == 0 { scalar chupper = 99999}
	else {
		if x == m {scalar xtemp = x - 0.5}	
		else {scalar xtemp = x}

 		if y == n {scalar ytemp = y - 0.5}
		else {scalar ytemp = y}
	}


	scalar bisect1 = 999
	scalar bisect2 = T
	scalar bisect = (bisect1 + bisect2)/2
	scalar tol = 1


	local i=1

	while `i'<=100 & abs(tol) > 0.0001 {
		scalar p1hat = (bisect*(m+ytemp) + xtemp + n - /*
	           */ sqrt(  (bisect*(m+ytemp)+xtemp+n)^2 - (4*bisect*(m+n)*(xtemp+ytemp))  ) ) / (2*(m+n))
		scalar p2hat = p1hat/bisect
	
		if x == m & y == n & T >=1 {scalar Utheta0 = n*(bisect-1)}
		else {
		scalar Utheta0 = (((xtemp-m*p1hat)^2) / (m*p1hat*(1-p1hat))) + (((ytemp-n*p2hat)^2) / (n*p2hat*(1-p2hat)))
		}

		if Utheta0 < cinv {scalar bisect2 = bisect}
		else { scalar bisect1 = bisect}

		scalar bisect = (bisect1+bisect2)/2
		scalar tol = Utheta0 - cinv
		scalar chupper = bisect
		local i=`i'+1
	}	 


	scalar  PE = 100 * (1 - T)
	scalar  PEchil = 100 * (1-chupper)
	scalar  PEchiu = 100 * (1-chlower)

*display the results

*	wcrc4fld x m-x y n-y `"`col'"' "" 
	wcrc4fld x m-x y n-y "    Event    " " " /*
		*/ Yes No Group1 Group2 /*
		*/ "" yes
	di in gr _col(18) "|" _col(43) "|"
	di in gr _col(18) /*
	*/ "|      Point estimate    |  [`level'% Conf. Interval]" _n /*
	*/_col(1) _dup(17)"-" _col(18) "|" _dup(24) "-" "+" _dup(22) "-"  
	
	di in gr "      Odds Ratio | " in ye /* 
		*/ _col(27) %9.0g T in gr _col(43) "| " /*
		*/ in ye %9.0g  chlower "   " %9.0g chupper 

	di in gr _col(1) _dup(65) "-" 

*save results

	return scalar  x      = x
	return scalar  y      = y
	return scalar  m      = m
	return scalar  n      = n
 	return scalar  theta      = T
	return scalar  theta_l    = chlower 
	return scalar  theta_u    = chupper
 	return scalar  level  = `level'

end

*! version 1  18Aug2000
program define wcrc4fld
	version 6
	args a b c d stub1 stub2 col1 col2 row1 row2 notot addpro
	if `"`stub1'"'!="" { 
		local header yes
		local s1 : piece 1 24 of `"`stub1'"'
		di _col(18) in gr `"|            `s1'"' _col(43) "|" _c
		local s1 : piece 2 24 of `"`stub1'"'
		local j 2
		while `"`s1'"' != "" {
			di _n _col(18) in gr `"| `s1'"' _col(43) "|" _c
			local j = `j' + 1
			local s1 : piece `j' 24 of `"`stub1'"'
		}
	}
	local s1 = 9-length(`"`col1'"')
	local s2 = 9-length(`"`col2'"')
	if `"`addpro'"'!="" { 
		if `"`header'"'=="yes" {
			di in gr _skip(12) "Proportion"
		}
		else	di in gr _col(57) "Proportion"
		local total "Total         Yes"
		local dd 22
		local cc "_c"
	}
	else {
		if `"`header'"'=="yes" { di }
		local total "Total"
		local dd 10
	}
	di in gr `"`stub2'"' _col(18) "| " /*
		*/ _skip(`s1') `"`col1'"' "   " /*
		*/ _skip(`s2') `"`col2'  |     `total'"' _n /*
		*/ _dup(17) "-" "+" _dup(24) "-" "+" _dup(`dd') "-"
	local row1=substr(`"`row1'"',1,16)
	local s1 = 16-length(`"`row1'"')
	di in gr _skip(`s1') `"`row1' | "' /*
		*/ in ye %9.0g `a' "   " %9.0g `b' /*
		*/ in gr "  | " in ye %9.0g `a'+`b' `cc'
	if `"`addpro'"'!="" { 
		di %12.4f `a'/(`a'+`b')
	}
	local row2=substr(`"`row2'"',1,16)
	local s1 = 16-length(`"`row2'"')
	di in gr _skip(`s1') `"`row2' | "' /*
		*/ in ye %9.0g `c' "   " %9.0g `d' /*
		*/ in gr "  | " in ye %9.0g `c'+`d' `cc'
	if `"`addpro'"'!="" { 
		di %12.4f `c'/(`c'+`d')
	}
	di in gr /*
		*/ _dup(17) "-" "+" _dup(24) "-" "+" _dup(`dd') "-"
	if `"`notot'"'!="" { exit } 
	di in gr _col(12) "Total | " in ye /*
		*/ %9.0g `a'+`c' "   " %9.0g `b'+`d' /*
		*/ in gr "  | " in ye %9.0g `a'+`b'+`c'+`d' `cc'
	if `"`addpro'"'!="" { 
		di %12.4f (`a'+`c')/(`a'+`b'+`c'+`d')
	}
end