*! 1.3.1 Anders Alexandersson 20030211
program define ellip7, rclass sortpreserve
	version 7.0
	if _caller() <= 6 {
		ellip6 `0'
		exit
	}
	syntax varlist(min=1 max=2 numeric) [if] [in] [, COEFS MEANS  /*
	*/ Level(int $S_level) NPoints(integer 400) Pool(int -1)      /*
	*/ Constant(str) CONNect(str) SYmbol(str) Generate(str)       /*
	*/ Add(varlist min=2 max=2 numeric) B2title(str) L1title(str) /*
	*/ T1title(str) T2title(str) YFormat XFormat YFormat2(str) /*
	*/ XFormat2(str) SAving(str) EVR(real 1) NOGRaph REPLACE /*
        */ BY(varname) Need(str) *]
	* parse generate(str)
	if "`generate'" ~= "" {
		local wc : word count `generate'
		if `wc' ~= 2 {
			di as err "exactly TWO variables required in generate()"
			exit 198
		}
		tokenize `generate'
		args new1 new2
		if "`replace'" ~= "" {
			capture drop `new1' `new2'
		}
		confirm new var `new1'
		confirm new var `new2'
		marksample touse, novarlist	/* include missing obs */
	}
        else {
		marksample touse		/* exclude missing obs */
	}
	* parse constant(str)
	if "`constant'" ~= "" {
		local wc : word count `constant'
		if `wc' >= 3 {
			di as err "max TWO arguments allowed in constant()"
			exit 198
		}
		tokenize "`constant'"
		args stat p
		confirm names `stat'
		if "`p'" ~= "" {
			confirm number `p'
		}
	}

	* set common defaults
	if "`nograph'" ~= "" {
		local graph ""
	}
	else {
		local graph "graph"
	}
	local symbol = cond("`symbol'"=="", "......", "`symbol'")
	local symbol "symbol(`symbol')"
	local connect = cond("`connect'"=="", "llllll", "`connect'")
	local connect "connect(`connect')"
	if "`saving'" ~= "" {
		local saving "saving(`saving')"
	}
	set textsize 125
	if `"`t1title'"' == `""' {
		local t1title `" "'
	}
	if `"`t2title'"' == `""' {
		local t2title `" "'
	}

	* check for errors
	if "`generate'" ~= "" & "`by'" ~= "" {
		di as err "may not combine options generate() and by()"
		exit 198
	}
	if `level' < 10 | `level' > 99 {
		di as err `"level(`level') invalid"'
		di as text "level() must be an integer between 10 and 99"
		exit 198
	}
	if `npoints' < 20 | `npoints' > 9999 {
		di as err `"npoints(`npoints') invalid"'
		di as text "npoints() must be an integer between 20 and 9999"
		exit 198
	}
	if `evr' < 0 | `evr' > 10^36 {
		di as err `"evr(`evr') invalid"'
		di as text "evr() must be a floating-point number between " /*
		*/ "0 and 10^36"
		exit 198
	}
	if "`means'" ~= "" & "`coefs'" ~= "" {
		di as err "may not specify both means and coefs"
		exit 198
	}
	if "`yformat'" ~= "" & `"`yformat2'"' ~= "" {
		di as err "may not specify both yformat and yformat()"
		exit 198
	}
	if `"`yformat2'"' ~= "" {
		capt local tmp : display `yformat2' 1
		if _rc {
			di as err `"invalid %fmt in yformat(`yformat2')"'
			exit 120
		}
	}
	if "`xformat'" ~= "" & `"`xformat2'"' ~= "" {
		di as err "may not specify both xformat and xformat()"
		exit 198
	}
	if `"`xformat2'"' ~= "" {
		capt local tmp : display `xformat2' 1
		if _rc {
			di as err `"invalid %fmt in xformat(`xformat2')"'
			exit 120
		}
	}
	if "`constant'" == "sd" & `level' ~= $S_level {
		di as err "may not specify both constant(sd) and level()"
		exit 198
	}
	if "`constant'" == "sq" & `level' ~= $S_level {
		di as err "may not specify both constant(sq) and level()"
		exit 198
	}
	if "`pool'" ~= "-1" {
		ChkPool "`pool'" "`add'" "`if'" "`in'" "`graph'"
	}

	* define macros, and tokenize
	tempvar t new1 new2 yy xx id tmpy tmpx iwt
	tempname n f F r a b c xbar ybar s_x s_y x_a y_a A B m1 m2 /*
	*/ var_y var_x cov_xy C D lambda1 lambda2
	tempfile new tmp
	local old_N = _N
	if "`generate'" ~= "" {
		tokenize `generate'
		args new1 new2
	}
	if "`add'" ~= "" {
		tokenize `add'
		args old1 old2
	}
	if "`constant'" ~= "" {
		tokenize `constant'
		args stat p
	}
	tokenize "`varlist'"
	args y x

	* set centering-specific defaults
	if "`coefs'" == "coefs" {
		if "`e(cmd)'" ~= "regress" {
			exit 198
			* exit 301
		}
		if "`x'" == "" {
			local x "_cons"
		}
		else {
			local x "`2'"
			if "`e(depvar)'" == "`x'" | "`e(depvar)'" == "`y'" {
				di as err "`e(depvar)' is not independent variable"
				exit 111
			}
		}
		if "`l1title'" == "" {
			local l1title "Estimated `y'"
		}
		if "`b2title'" == "" {
			local b2title "Estimated `x'"
		}
	}
	else {
		local ylab : variable label `y'
		local xlab : variable label `x'
		if "`l1title'" == "" {
			if "`ylab'" == "" {
				local l1title "`y'"
			}
			else {
				local l1title "`ylab'"
			}
		}
		if "`b2title'" == "" {
			if "`xlab'" == "" {
				local b2title "`x'"
			}
			else {
				local b2title "`xlab'"
			}
		}
	}

	if "`by'" == "" {
		* calculate centering-specific basic ingredients
		if "`coefs'" == "coefs" {
			local ybar "_b[`y']"
			local s_y "_se[`y']"
			local var_y "_se[`y']^2"
			local xbar "_b[`x']"
			local s_x "_se[`x']"
			local var_x "_se[`x']^2"
			qui cor if `touse', _coef
			scalar `n' = e(N)
		}
		else {
			qui summarize `y' if `touse'
			scalar `ybar' = r(mean)
			scalar `s_y' = r(sd)
			scalar `var_y' = r(Var)
			qui summarize `x' if `touse'
			scalar `xbar' = r(mean)
			scalar `s_x' = r(sd)
			scalar `var_x' = r(Var)
			qui corr `y' `x' if `touse'
			scalar `n' = r(N)
		}
		scalar `r' = r(rho)
		scalar `cov_xy' = `r' * `s_x' * `s_y'

		* calculate c in subroutine, and then amplitudes of arcs x and y
		DoCons "`p'" `n' `level' "`stat'" "`coefs'"
		scalar `c' = r(c)
		scalar `x_a' = `s_x' * sqrt(`c')
		scalar `y_a' = `s_y' * sqrt(`c')

		* bring back data, calculate ellipse, and return results
		preserve
		qui range `t' 0 2*_pi `npoints'				/* change data */
		qui gen `new2' = `xbar' + `x_a' * cos(`t')
 		qui gen `new1' = `ybar' + `y_a' * cos(`t' + acos(`r'))
		qui summarize `new1'
		ret scalar min_y = r(min)	/* min of y */
		ret scalar max_y = r(max)	/* max of y */
		qui summarize `new2'
		ret scalar min_x = r(min)	/* min of x */
		ret scalar max_x = r(max)	/* max of x */

		* create macros for readabaility: ylist, xlist
		if "`generate'" == "" {
			if "`add'" == "" {
				local ylist "`new1'"
				local xlist "`new2'"
			}
			else {
				stack `new1' `new2' `old1' `old2', into(`yy' `xx') clear
				qui gen `new1' = `yy' if _stack==1  /* (obs 001-400) */
				qui gen `old1' = `yy' if _stack==2  /* (obs 401-800) */
				keep `new1' `old1' `xx'             /*    >= 800 obs */
				local ylist "`new1' `old1'"
				local xlist "`xx'"
			}
		}
		else {							/* gen() */
			if "`add'" == "" {
				local ylist "`new1'"
				local xlist "`new2'"
			}
			else {						/* add() */
				qui save "`new'" /* save dataset to allow damage */
				stack `new1' `new2' `old1' `old2', into(`yy' `xx') clear
				qui gen `new1' = `yy' if _stack==1 /* (obs 001-400) */
				qui gen `old1' = `yy' if _stack==2 /* (obs 401-800) */
				keep `new1' `old1' `xx'            /*    >= 800 obs */
				if "`pool'" == "-1" {
					local ylist "`new1' `old1'"
					local xlist "`xx'"
				}
				else {					/* pool() */
					qui gen long `id' = _n
					sort `id'
					local ylist "`new1' `old1'"
					local xlist "`xx'"
					qui save "`tmp'"		/* save tmp dataset */
				}
			}
		}

		* format ylist and xlist
		if `"`yformat2'"' ~= "" {
			format `ylist' `yformat2'
		}
		else format `ylist' %9.0g /* includes fy */
		if `"`xformat2'"' ~= "" {
			format `xlist' `xformat2'
		}
		else format `xlist' %9.0g /* includes fx */

		* graph ellipse
		if "`generate'" ~= "" & "`add'" ~= "" & "`pool'" ~= "-1" {
			* retrieve tmp dataset, and create fractionally pooled datasets
			use "`tmp'", clear
			* get rhs-variables (excluding _cons)
			tempname b
			mat `b' = e(b)
			local rhs : colnames(`b')
			local rhs : subinstr local rhs "_cons" "", word
			forvalues i = 0/`pool' {
				restore, preserve /* N >= 400 */
				qui gen double `iwt' = 1 if `touse'==1
				qui replace `iwt' = [`i']/`pool' if `touse'==0

				* bugfix in v1.3.1: `rhs' = iv1 iv2 ...
				qui reg `e(depvar)' `rhs' [iw=`iwt'], /*
				*/ level(`level') noheader
				qui gen `tmpy' = `ybar'
				qui gen `tmpx' = `xbar'
				qui keep `iwt' `tmpy' `tmpx'
				qui keep in l/l
				tempfile wt`i'
				qui save "`wt`i''"
			}

			* append into one file which contains the locus curve
			use "`wt0'", clear
			forvalues i = 1/`pool' {
				append using "`wt`i''"
			}

			* save locus curve, and merge with tmp dataset
			qui gen long `id' = _n
			sort `id'
			qui save "`wt`pool''", replace

			* return matrix for locus curve, as in -corrgram-
			tempname IWT YVAR XVAR
			mkmat `iwt', matrix(`IWT')
			mkmat `tmpy', matrix(`YVAR')
			mkmat `tmpx', matrix(`XVAR')
			mat colnames `IWT' = iwt
			mat colnames `YVAR' = "`y'"
			mat colnames `XVAR' = "`x'"
			ret matrix IWT `IWT'
			ret matrix YVAR `YVAR'
			ret matrix XVAR `XVAR'
			merge `id' using "`tmp'"

			* use the -gr20- STB-34 package within gph
			gph open, `saving'
			* pen() option in graph is invalid in gphdt -> r(198)
			gra `ylist' `xlist', `connect' `symbol' /*
			*/ t1("`t1title'") t2("`t2title'") /*
			*/ b2("`b2title'") l1("`l1title'") `options'
			gphsave
			local end = `pool' + 1
			local y_b = `tmpy' in 1/1			/*  b(y;_) */
			local y_bp = `tmpy' in `end'/`end'		/* bp(y;_) */
			local x_b = `tmpx' in 1/1			/*  b(_;x) */
			local x_bp = `tmpx' in `end'/`end'		/* bp(_;x) */
			gphdt text `y_b' `x_b' 0 1 b			/* b text */
			gphdt text `y_bp' `x_bp' 0 1 bp			/* bp text */
			qui gphdt vline `tmpy' `tmpx' in 1/`end'	/* var line */
			qui gphdt vpoint `tmpy' `tmpx' in 1/`end'	/* var pts */
			gph close
			qui keep if `iwt' ~= .
			rename `iwt' iweight
			rename `tmpy' `y'
			di						/* display table */
			di as text "The Fractionally Pooled Estimates"
			di as text "(i.e., the dots in the locus curve):"
			if "`x'" == "_cons" {
				rename `tmpx' constant			/* _cons is system variable */
				tabdisp iweight, cell(`y' constant) center
			}
			else {
				rename `tmpx' `x'
				tabdisp iweight, cell(`y' `x') center
			}
		}
		else {  /* NOT pool() */
			if "`graph'" ~= "" {
				gra `ylist' `xlist', by(`by') `connect' `symbol'  /*
				*/ t1("`t1title'") t2("`t2title'")  /*
				*/ b2("`b2title'") l1("`l1title'") `saving' `options'
			}
		}
		* make sure that dataset is not damaged
		if "`generate'" ~= "" & "`add'" ~= "" {
			use "`new'", clear
		}

		* all went well, keep any generated vars
		if "`generate'" ~= "" {
			if _N > `old_N' {
				di as text "obs was `old_N', now " _N
			}
			restore, not
		}
	}
	else { /* by() */
		preserve
		tempvar Y id
		tempname byname
		local by "`by'"
        	tempvar group
        	sort `touse' `by'
        	qui by `touse' `by': gen byte `group' = _n == 1 if `touse'
        	qui replace `group' = sum(`group')
        	local max = `group'[_N]
		qui range `t' 0 2*_pi `npoints'			/* 400 obs */
		qui tab `by', matrow(`byname')
		svmat `byname', name(`id') /* 0,1,2,... */

		* create yvars and xvars (ellipse variables for each group)
		tempvar subuse
		qui gen byte `subuse' = .
		* get rhs-variables (excluding _cons)
		tempname b
		mat `b' = e(b)
		local rhs : colnames(`b')
		local rhs : subinstr local rhs "_cons" "", word

		forvalues i = 1/`max' {
			if "`coefs'" == "coefs" {	/* bugfix in v 1.3.1 */
				* -reg ...by()- remembers only last estimates
				qui replace `subuse' = `touse' & `group' == `i'
				qui reg `e(depvar)' `rhs' if `subuse'
				local ybar "_b[`y']"
				local s_y "_se[`y']"
				local var_y "_se[`y']^2"
				local xbar "_b[`x']"
				local s_x "_se[`x']"
				local var_x "_se[`x']^2"
				qui cor if `touse', _coef
				scalar `n' = e(N)
			}
			else {
				qui summarize `y' if `group' == `i'
				scalar `ybar' = r(mean)
				scalar `s_y' = r(sd)
				scalar `var_y' = r(Var)
				qui summarize `x' if `group' == `i'
				scalar `xbar' = r(mean)
				scalar `s_x' = r(sd)
				scalar `var_x' = r(Var)
				qui corr `y' `x' if `group' == `i'
				scalar `n' = r(N)
			}

			scalar `r' = r(rho)
			scalar `cov_xy' = `r' * `s_x' * `s_y'
			* calculate c in subroutine, and then amplitudes of arcs x and y
			DoCons "`p'" `n' `level' "`stat'" "`coefs'"
			scalar `c' = r(c)
			scalar `x_a' = `s_x' * sqrt(`c')
			scalar `y_a' = `s_y' * sqrt(`c')
			tempvar ynew`i' xnew`i'
			gen `ynew`i'' = `ybar' + `y_a' * cos(`t' + acos(`r'))
			gen `xnew`i'' = `xbar' + `x_a' * cos(`t')

			* return min_x and min_y?
			qui summarize `ynew`i''
			ret scalar min_y = r(min)	/* min of y */
			ret scalar max_y = r(max)	/* max of y */
			qui summarize `xnew`i''
			ret scalar min_x = r(min)	/* min of x */
			ret scalar max_x = r(max)	/* max of x */
			local stargs "`stargs' `ynew`i'' `xnew`i'' `id' `by'"
		}


		* default titles: if no varlab, use varnam
		local bylab : variable label `by'
		if "`title'" == "" {
			if "`bylab'" == "" {
				local title "`by'"
			}
			else {
				local title "`bylab'"
			}
		}

		* stack does not allow into (_cons )
		if "`x'" == "_cons" {
			tempvar z
			gen `z' = 1
			stack `stargs', into (`Y' `z' `id' `by') wide clear
			local x "`z'"
		}
		else {
			stack `stargs', into (`Y' `x' `id' `by') wide clear
		}
		forvalues i = 1/`max' {
			tempvar sy_`i'
			qui gen `sy_`i'' = `Y' if _stack==`i'
			local ylist "`ylist' `sy_`i''"
			qui replace _stack = `id'[`i'] if _stack==`i'
		}
		local xlist "`x'"
		drop `by'
		rename _stack `by'
		sort `by'
		gra `ylist' `xlist', by(`by') `connect' `symbol' /*
		*/ l1("`l1title'") t1("`t2title'") /*
		*/ b2("`b2title'") b1(Graphs by "`title'") `saving' `options'
	}
	* calculate scalars to be returned
	scalar `A' = `var_y' - `evr' * `var_x'
	scalar `B' = `cov_xy'
	scalar `m1' = (`A' + sqrt((`A')^2 + 4*`evr' * (`B')^2)) / (2 * `B')
	scalar `C' = `var_y' + `var_x'
	scalar `D' = ((`C'^2) - 4*(`var_y' * `var_x' - `cov_xy'^2))^.5
	scalar `lambda1' = (`C' + `D') / 2
	scalar `lambda2' = `C' - `lambda1'  /* or: (`C' - `D') / 2 */

	* save results as scalars in r()
	ret scalar a = (1 / `lambda1')^-.5 * sqrt(`c')	/* l of semi-major */
	ret scalar b = (1 / `lambda2')^-.5 * sqrt(`c')	/* l of semi-minor */
	ret scalar c       = `c'		/* boundary constant */
	ret scalar r       = `r'		/* correlation coefficient */
	ret scalar m1      = `m1'		/* slope of major axis */
	ret scalar m2      = - 1 / `m1'		/* slope of minor axis */
	ret scalar lambda1 = `lambda1'		/* first eigenvalue of inverse V */
	ret scalar lambda2 = `lambda2'		/* second eigenvalue of inverse V */
	ret scalar mean_y  = `ybar'		/* mean of y */
	ret scalar s_y     = `s_y'		/* standard error of y */
	ret scalar mean_x  = `xbar'		/* mean of x */
	ret scalar s_x     = `s_x'		/* standard error of x */
end

program define ChkPool, rclass
        args pool add if in graph
	capture which gphdt
	if _rc ~= 0 {
		di in re "pool() must have gphdt.ado installed"
		exit _rc
	}
	capture which gphsave
	if _rc ~= 0 {
		di in re "pool() must have gphsave.ado installed"
		exit _rc
	}
	if "`add'" == "" {
		di as err "pool() must be used with add()"
		exit 198
	}
	if "`if'" == "" & "`in'" == "" {
		di as err "must use pool() with [if exp] or [in range]"
		exit 198
	}
	if "`graph'" == "" {
		di as err "may not specify both nograph and pool()"
		exit 198
	}
	if `pool' < 1 | `pool' > 100 {
		di as err `"pool(`pool') invalid"'
		di as text "pool() must be an integer between 1 and 100"
		exit 198
	}
	capture confirm new variable _merge
	if _rc ~= 0 {
		di as err "drop _merge before running ellip with pool()"
		exit _rc
	}
end

program define DoCons, rclass	/* returns r(c) */
	args p n level stat coefs
	tempname Fadj F c Chi T2
	local means_d `"("`stat'" == "" & "`coefs'" ~= "coefs")"'
	local coefs_d `"("`stat'" == "" & "`coefs'" == "coefs")"'
	* calculate p
	if "`stat'" == "sd" | (`means_d') {
		if "`p'" == "" { local p = "4" }
		else { local p = (`p')^2 }
	}
	else if "`stat'" == "sq" {
		if "`p'" == "" { local p = "4" }
		else { local p = `p' }
	}
	else {
		if "`p'" == "" { local p = "2" }
		else { local p = "`p'" }
	}

	* calculate c
	if "`stat'" == "fadj" {
		scalar `Fadj' = invfprob(2,(`n'-`p'),(100-`level')/100)
		ret scalar c = (2 * `Fadj')
	}
	else if "`stat'" == "f" | `coefs_d' {
		scalar `F' = invfprob(`p',(`n'-`p'),(100-`level')/100)
		ret scalar c = (2 * `F')
	}
	else if "`stat'" == "chisq" {
		ret scalar c = invchi2tail(2,(100-`level')/100) /* v7 */
	}
	else if "`stat'" == "chisqn" {
		scalar `Chi' = invchi2tail(2,(100-`level')/100) /* v7 */
		ret scalar c = (`Chi' / `n')
	}
	else if "`stat'" == "tsq" | "`stat'" == "hotel" {
		scalar `F' = invfprob(`p',(`n'-`p'),(100-`level')/100)
		scalar `T2' = `p' * (`n'-1) / (`n'-`p') * `F'
		ret scalar c = `T2'
	}
	else if "`stat'" == "tsqn" | "`stat'" == "hoteln" {
		scalar `F' = invfprob(`p',(`n'-`p'),(100-`level')/100)
		scalar `T2' = `p' * (`n'-1) / (`n'-`p') * `F'
		ret scalar c = (`T2' / `n')
	}
	else if "`stat'" == "ptsq" | "`stat'" == "photel" {
		scalar `F' = invfprob(`p',(`n'-`p'),(100-`level')/100)
		scalar `T2' = `p' * (`n'-1) / (`n'-`p') * `F'
		ret scalar c = (`T2' * (`n' + 1) / `n')
	}
	else if "`stat'" == "sd" | "`stat'" == "sq" | (`means_d') {
		ret scalar c = `p'
	}
end
exit

HISTORY
1.3.1 Anders Alexandersson 10feb2003			SSC Archive
       - Bug fixes of ellip7:
       1. coefs by()- would give the results for -means by()-;
         this is fixed in ellip7 (N/A for ellip6 and ellip5)
       2. -coefs by()- failed to return results; this is fixed
       3. -pool()- would unly use 2 independent variables in reg;
         this is fixed in ellip7 but not in ellip6 (no bug in ellip5)
       - replaces version 1.3.0 with same requirements
1.3.0 Anders Alexandersson 16jan2003			(outdated)
       - new nograph and c(sq) option; fixed by() option for means
       - requires Stata 7, gphdt for pool
1.2.0 Anders Alexandersson 16jan2003                    SSC Archive
       - fixes several minor bugs such as some r() in 1.1.0
       - is tested to reproduce 20+ published confidence ellipses
       - allows more options in c()
       - uses more modular code and calls -graph- less times
       - replaces version 1.1.0 with same requirements
1.1.0 Anders Alexandersson 12nov2001                    (outdated)
       - allows centering around means and more types of ellipses
       - requires Stata 6, gphdt for pool()
       - version 1.0.0 renamed into ellip_5.ado
1.0.0 Anders Alexandersson 7aug1998                     STB-46: gr32
       - graphs confidence ellipses around regression coefficients
       - requires Stata 5, gphdt for pool(), parsoptp