*! xtgraph.ado, Version 2.02
*! by PT Seed (paul.seed@kcl.ac.uk)

*! graphs of xt (cross-sectional time-series) data
*! Presents averages with error bars of a single variable, by time
*! data can be grouped.

* 14/4/2001 All major functions now seem to be working, ready for SUG
* Log & power transformations, model fitting, outputting data & graphs, 
 
* The only problem seems to be that the undeclared "show" does not 
* work with axis labels on the graph.

* -av(am)- now allowed as a synonym for -av(mean)-

* 19/6/2001 Add in error sd as an option following -model-

* version 2.02
* Include boxcox(#) option
* Correct a series of bugs:
* Half-bars facing the wrong way, 
* Missing values in string group variables,
* Spaces in temporary file names.

* To allow version 8.0 graphics


prog define xtgraph
version 8.0

preserve

 *****  parse the command ***************
	error 0
	syntax varlist(min=1 max=1 numeric) [if] [in]  , ///
[AVerage(string) bar(string) display ///
POWer(real 1) log(string) boxcox(real 1) model ///
i(string) t(string) GRoup(varlist max=1) Half ///
OFFset(real 0) show MINobs(real 1) level(real $S_level) ///
list savedat(string asis) by(string) listwise missing ///
baropts(string) lineopts(string) bartype(string) ///
L1title(string) Symbol(string) Connect(string) Weight(string) nograph v7 * ]

if "`display'" == ""  {
	local display = "*" 
}	

`display'"*****  ***** ***** ***** ***** ***** 

`display'"*****  Collect the I and T vars ***** 
	xt_iis `i'
	local ivar "`s(ivar)'"

	xt_tis `t'
	local tvar "`s(timevar)'"
`display'"*****  ***** ***** ***** ***** ***** 

`display'"*****  Handle the conditions *****

	local yvar "`varlist'"
	cap keep `if' `in'
	qui keep if `yvar' ~= . 
	cap keep if `group' ~= .
	cap keep if `group' ~= "." & `group' ~= ""

	tempvar n
	if `minobs' > 1 {
		qui egen `n' = count(`yvar'), by(`group' `by' `tvar') 
		cap keep if `n' >= `minobs'
		if _rc & "`1'" ~= ""  { 
			di in red "Invalid option minobs(" in ye"`minobs'" in red")"
			exit 198
		}
	drop `n'
	}

	if "`listwise'" ~= "" {
		qui egen `n' = count(`yvar'), by(`ivar') 
		summ `n' , meanonly
		keep if `n' == r(max)
	drop `n'
	}

	if "`missing'" == "" { 
		cap drop if `group' == . 
		cap drop if `group' == ""
		if "`by'" ~= "" {
			cap drop if `by' == .
			cap drop if `by' == ""
		}
	}


	qui egen `n' = count(`yvar'), by(`group' `by' `tvar')

	if "`saving'"  ~= "" {
		local graph = "" 
	}
***** ***** ***** ***** ***** ***** 

`display'"*****  Check the syntax ***** 


	if "`log'" ~= "" { 
		cap confirm number `log' 
		if _rc == 7 {
			di in red "Option log() must contain a number." 
			exit 198
		}
	}

	if "`log'" ~= "" & "`power'" ~= "1" {
		di in red "Options log() and power() canot be used together	
		exit 198
	}
	if "`log'" ~= "" & "`boxcox'" ~= "1" {
		di in red "Options log() and boxcox() canot be used together	
		exit 198
	}
	if "`boxcox'" ~= "1" & "`power'" ~= "1" {
		di in red "Options boxcox() and power() canot be used together	
		exit 198
	}

	if "`power'" == "0" | "`boxcox'" == "0" {
		local log = 0
		local power = 1
		local boxcox = 1
	}


	if "`average'" == "" |  "`average'" == "am" {
		local average = "mean" 
	}
	if ("`average'" == "gm" | "`average'" == "hm" ) & "`power'" ~= "1" {
		di in red "Option power() should not be used with average(`average')"
		exit 198
	}
	if ("`average'" == "gm" | "`average'" == "hm" ) & "`boxcox'" ~= "1" {
		di in red "Option boxcox() should not be used with average(`average')"
		exit 198
	}

	if "`average'" == "gm" | "`boxcox'" == "0" | "`power'" == "0" {
		local average "mean" 
		local power = 1
		local boxcox = 1
		local log = 0
	}

	if "`average'" == "hm" {
		local average "mean" 
		local power = -1
	}

	if "`bar'" == "" {
		local bar "ci" 
	}

	if "`average'" == "median" {
		if "`bar'" ~= "ci" & "`bar'" ~= "no" & "`bar'"~= "iqr" & "`bar'" ~= "rr" { 
			di in red "Option bar(`bar') invalid with average(`average')."
			exit 198
		}
	}

	if "`average'" == "mean" {
		if "`bar'" ~= "ci" & "`bar'" ~= "no" & "`bar'"~= "se" /*
		*/ & "`bar'"~= "sd" & "`bar'" ~= "rr" & "`bar'" ~= "error" { 
			di in red "Option bar(`bar') invalid with average(`average')."
			exit 198
		}
	}


	if "`group'" == "" {
		tempvar group
		qui gen `group' = 1
	}

	cap assert ~(( "`bar'" == "no") & ("`half'" == "half" ))
	if _rc { 
		di in gr "Note: There is no point in asking for half-bars if no bars are requested." 
		local half
	}


	if "`bartype'" == "" local bartype rcap

`display'"*****  Symbol ***** 

	tempvar gpno 
	sort `group'
	qui by `group' : gen `gpno' = _n == 1
	qui replace `gpno' = sum(`gpno')
	sort `gpno'
	local ngrp = `gpno'[_N]

	while length("`symbol'") < `ngrp' & `ngrp' ~= . { 
		local symbol "`symbol'OSTodp" 
	}
	local symbol = substr("`symbol'",1,`ngrp') 
	local symbol "`symbol'ii"


`display'"*****  Connect ***** 

	while length("`connect'") < `ngrp' { 
		local connect "`connect'l" 
	}
	local connect "`connect'II"

`display'"*****  Half ***** 

	if "`half'" ~= "" & `ngrp' ~= 2 {
		di in red "Half-bars are not appropriate with more than two groups."
		di in red "Perhaps you should use the offset()option."
		exit 198
	}

`display'"*****  Collect the Y and T labels ***** 

	local y_vrlab : var label `varlist'
	if "`y_vrlab'" == "" { 
		local y_vrlab = "`varlist'" 
	}
	local y_vllab : val label `varlist'

	local t_vrlab : var label `tvar'
	if "`t_vrlab'" == "" { 
		local t_vrlab = "`tvar'" 
	}
	local t_vllab : val label `tvar'
	

***** ***** ***** ***** ***** ***** 

`display'"*****  Carry out the transformations *****

	if "`power'" ~= "1" & "`model'" == "" { 
		qui replace `yvar' = `yvar'^`power' 
	}

	else if "`boxcox'" ~= "1" & "`model'" == "" { 
		qui replace `yvar' = (`yvar'^`boxcox'-1)/`boxcox' 
	}

	else if "`log'" ~= "" & "`model'" == "" { 
		cap assert `yvar' - `log' > 0 
		if _rc { 
			local ln_sign = "-" 
			cap assert -`yvar' - `log' > 0 
			if _rc {
				di in red "Non-positive values of `yvar' - `log' encountered." 
				exit _rc
			}
		}
		qui replace `yvar' = log(`ln_sign'1*`yvar'-`log')
	}
	
******************************************

`display'"*****  Calculate the averages and bars ***********

	tempvar av 
	if "`bar'" ~= "no" { 
		tempvar lb ub 
	}

`display'"*****  ***** ***** means ***** ***** ***** 
	if "`average'" == "mean" { 
		if "`model'" == "" { 
			qui egen `av' = `average'(`yvar'), by(`group' `by' `tvar') 
		}
		else { 
			cap predict `av' 
			if _rc { 
				di in red "The last model fitted cannot be detected."
				exit _rc
			}
		}

`display'"*****  ***** sd & rr ***** ***** 
		if "`bar'" == "sd" | "`bar'" == "rr" {
			tempvar sd 
			if "`model'" == "" { 
				qui egen `sd' = sd(`yvar'), by(`group' `by' `tvar') 
			}
			else { 
				cap predict `sd' , stdf 
				if _rc { 
					di in red "The last model fitted does not permit estimates of stdf,
					di in red "needed for standard deviation bars"
					exit _rc
				}
			}
			if "`bar'" == "sd" {
				qui gen `lb' = `av' - `sd'
				qui gen `ub' = `av' + `sd'
			}
			if "`bar'" == "rr" {
				qui gen `lb' = `av' - invnorm(.5+`level'/200)*`sd'
				qui gen `ub' = `av' + invnorm(.5+`level'/200)*`sd'
			}
		}

`display'"*****  ***** se & ci ***** ***** 

		if "`bar'" == "ci"  | "`bar'" == "se" {
 			tempvar se
			if "`model'" == "" {
				qui egen `se' = sd(`yvar'), by(`group' `by' `tvar') 
				qui replace `se' = `se'/`n'^.5
			}
			else { 
				qui predict `se', stdp  
				if _rc { 
					di in red "The last model fitted does not permit estimates of stdp,
					di in red "(the standard deviation of the forecast) needed for standard error bars"
					exit _rc
				}
			}

			if "`bar'" == "ci" {
				qui gen `lb' = `av' - (invttail(`n'-1,.5-`level'/200)*`se')
				qui gen `ub' = `av' + (invttail(`n'-1,.5-`level'/200)*`se')
			}
			if "`bar'" == "se"  {
				qui gen `lb' = `av' - `se'
				qui gen `ub' = `av' + `se'
			}
		}

`display'"*****  ***** error ***** ***** 
		if "`bar'" == "error" {
			local error = e(rmse)
			if "`error'" == "." { 
				di in red "The last model fitted does not give a root-mean square error."
				exit 198
			}
			qui gen `lb' = `av' - `error'
			qui gen `ub' = `av' + `error'
		}
	}
	***** ***** ***** ***** ***** ***** ***** 

`display'"*****  ***** ***** medians ***** ***** ***** 
	if "`average'" == "median" & "`bar'" == "ci" {
		qui gen `av' = .
		qui gen `lb' = .
		qui gen `ub' = .
		sort `group' `by' `tvar'
		tempvar gp_t
		qui by `group' `by' `tvar': gen `gp_t' = _n == 1
		qui replace `gp_t' = sum(`gp_t')
		local n_gp_t = `gp_t'[_N]
		local i = 1
		while `i' <= `n_gp_t' {
			qui centile `yvar' if `gp_t' == `i', level(`level')
			qui replace `av' = r(c_1) if `gp_t' == `i'
			qui replace `lb' = r(lb_1) if `gp_t' == `i'
			qui replace `ub' = r(ub_1) if `gp_t' == `i'
			local i = `i' + 1
		}
	}

	if "`average'" == "median" & "`bar'" == "iqr" {
		qui egen `av' = pctile(`yvar') , by(`group' `by' `tvar')
		qui egen `lb' = pctile(`yvar') , by(`group' `by' `tvar') p(25)
		qui egen `ub' = pctile(`yvar') , by(`group' `by' `tvar') p(75)
		}

	if "`average'" == "median" & "`bar'" == "rr" {
		local l_cent = 50-`level'/2
		local u_cent = 50+`level'/2
		qui egen `av' = pctile(`yvar') , by(`group' `by' `tvar')
		qui egen `lb' = pctile(`yvar') , by(`group' `by' `tvar') p(`l_cent')
		qui egen `ub' = pctile(`yvar') , by(`group' `by' `tvar') p(`u_cent')
	}
	if "`average'" == "median" & "`bar'" == "no" {
		qui egen `av' = pctile(`yvar') , by(`group' `by' `tvar')
	}

* Check that the model correctly predicts one value only per group.
	if "`model'" ~= "" {
		cap bys `group' `by' `tvar' (`av') : assert `av' == `av'[1] | `av' == .
		if _rc {
			di in red "The most recent model does not producer a single average for each group"
			di "as defined by `group' `by' `tvar'."
			exit _rc
		}
		cap bys `group' `by' `tvar' (`lb') : assert `lb' == `lb'[1] | `lb' == .
		if _rc & "`bar'" ~= "no "{
			di in red "The most recent model does not producer a single lower bound for each group"
			di "as defined by `group' `by' `tvar'."
			exit _rc
		}
		cap bys `group' `by' `tvar' (`ub') : assert `ub' == `ub'[1] | `ub' == .
		if _rc & "`bar'" ~= "no "{
			di in red "The most recent model does not producer a single upper bound for each group"
			di "as defined by `group' `by' `tvar'."
			exit _rc
		}


	}
	
	***** ***** ***** ***** ***** ***** 


`display'"*****  Carry out the back transformations *****
	if "`power'" ~= "1" & "`power'" ~= "0" { 
		qui replace `av' = `av'^(1/`power')
		if "`bar'" ~= "no" & `power' >=0 {

			qui replace `lb' = `lb'^(1/`power')
			qui replace `ub' = `ub'^(1/`power')
		}
		if "`bar'" ~= "no" & `power' <0 {
			tempvar lb_temp
			qui gen `lb_temp' = `ub'^(1/`power')
			qui replace `ub' = `lb'^(1/`power')
			qui replace `lb' = `lb_temp' 
		}
	}

	if "`boxcox'" ~= "1" & "`boxcox'" ~= "0" { 
		qui replace `av' = (`av'*`boxcox'+1)^(1/`boxcox')
		if "`bar'" ~= "no" & `boxcox' >=0 {

			qui replace `lb' = (`lb'*`boxcox'+1)^(1/`boxcox')
			qui replace `ub' = (`ub'*`boxcox'+1)^(1/`boxcox')
		}
		if "`bar'" ~= "no" & `boxcox' <0 {
			tempvar lb_temp
			qui gen `lb_temp' = (`ub'*`boxcox'+1)^(1/`boxcox')
			qui replace `ub' = (`lb'*`boxcox'+1)^(1/`boxcox')
			qui replace `lb' = `lb_temp' 
		}
	}

	if "`log'" ~= "" { 
		qui replace `av' = `ln_sign'1*(exp(`av')+`log')
		if "`bar'" ~= "no" {
			qui replace `lb' = `ln_sign'1*(exp(`lb')+`log')
			qui replace `ub' = `ln_sign'1*(exp(`ub')+`log')
		}
	}


***** ***** ***** ***** ***** 

`display'"*****  Offset the tvar *****

	tempvar t2
	if `offset' ~= 0 {
		qui tab `gpno', 
		qui gen `t2' = `tvar' + (`gpno'-(r(r)+1)/2)*`offset'
	}
	else {
		qui gen `t2' = `tvar' 
	}
	label var `t2' "`t_vrlab'"
	label val `t2' `t_vllab'
	
***** ***** ***** ***** ***** 

`display'"*****  Drop the extra variables & values *****

	sort `group' `by' `tvar'
	tempvar touse
	qui by `group' `by' `tvar' : gen `touse' = _n == 1 
	qui keep if `touse'

 	keep `ivar' `tvar' `n' `group' `gpno' `by' `yvar' `av' `avlist' `lb' `ub' `t2'
	if `group'[_N] == `group'[1] { 
		drop `group' 
		local group 	
	}
***** ***** ***** ***** ***** 

`display'"*****  Rename the variables ***** 
	cap rename `lb' lb
		if _rc == 0 { 
		local lb lb  
	}
	cap rename `ub' ub
		if _rc == 0 { 
			local ub ub  
		}
	cap rename `t2' t2
		if _rc == 0 { 
			local t2 t2 
		}
	cap rename `av' `average'
		if _rc == 0 { 
			local av `average'
		}
	cap noi rename `n' n
		if _rc == 0 { 
			local n n
		}
cap mac li _t2
cap desc `t2'
***** ***** ***** ***** ***** 

`display' "*****  Label the groups *****
	cap assert "`group'" ~= ""
	if _rc {
		local avlist "`av'"
		label var `avlist' "`y_vrlab'"
		local c "l"
	}

	else { 
		tempvar grp_num
		cap confirm numeric var `group'
		if _rc == 0 { 
			local grp_num "`group'" 
 		}
		else { 
			encode `group' , gen(`grp_num') 
		}
* grp_num is the numeric form of group.
		
		local i = 1
		while `i' <= `ngrp' {
			qui summ `grp_num' if `gpno' == `i', meanonly
			local gr_val = r(mean)
			local label  : label (`grp_num') `gr_val'
* First try to make name for average from label
			local name = lower(substr("`label'",1,8))
			while index("`name'", " ") {
				local l_name = substr("`name'", 1, index("`name'", " ")-1)
				local r_name = substr("`name'", index("`name'", " ")+1, .)
				local name "`l_name'_`r_name'"
			}
			cap confirm new var `name'
			while _rc { 
* Then try adding underscores
				local j = `j' + 1
				local name = "_" + substr("`name'",1,7)
				cap confirm new var `name'
				if `j' == 8 {
* Then try numbered groups
					local name "Group_`gr_val'"
				}
				if `j' == 10 {
* Then give up
					di in red "Problem giving meaningful names to all groups."
					di "No solution."
					exit 198
				}
			*End of [not] new var `name'
			}
`display' "*****  Rename the group average for group `i' *****
			qui gen `name' = `av' if `gpno' == `i'
			label var `name' "`label'"
			local avlist  "`avlist' `name'"
			local c "`c'l"
			local i = `i' + 1
		* End of consideration of name for group `i'
		}		
	* End of consideration of group names for multiple groups
	}
***** ***** ***** ***** ***** ***** 

`display'"*****  Label the Y axis ***** 
	if "`l1title'" == "" { 
		local l1title "`y_vrlab'"
	}
	local ytitle "ytitle(`l1title')"
	local l1title "l1(`l1title')"

`display'"*****  List the data *****

	if "`list'" ~= "" {
		if  "`average'" == "mean" {
			if "`log'" == "0" { 
				local av_type "Geometric mean" 
			}
			else if `power' == -1 { 
				local av_type "Harmonic mean" 
			}
			else if "`power'" ~= "1" { 
				local av_type = "Transformed mean [x^`power']" 
			}
			else if "`boxcox'" ~= "1" { 
				local av_type = "Box-Cox transformed mean [x^`power']" 
			}
			else if "`log'" ~= "" { 
				local av_type = "Transformend mean [log(x -`log')]" 
			} 
			else { 
				local av_type "Arithmetic mean" 
			}
		}
		else local av_type = "`average'"
	
		di _n "`av_type' of `varlist' " _c

		if "`bar'" == "ci" |  "`bar'" == "rr" { 
			local l = "`level'% " 
		}
		if "`bar'" ~= "no" {
			local up_bar = upper("`bar'")
			di "with bars based on `l'`up_bar'"
		}
		if "`group'" ~= "" & "`by'" == "" { 
			di "by `group' and `tvar'." 
		}
		else if "`group'" ~= "" & "`by'" ~= "" { 
			di "by `group', `by' and `tvar'." 
		}
		else if "`group'" == "" & "`by'" ~= "" { 
			di "by `by' and `tvar'." 
		}
		else di "by `tvar'."
		if "`group'" ~= "" | "`by'" ~= "" { 
			sort `group' `by' 
			by `group' `by' : list `tvar' `n' `av' `lb' `ub' , noobs 
		}
		else  { 
			list `tvar' `n' `av' `lb' `ub' , noobs 
		}
	}
***** ***** ***** ***** ***** 

`display'"*****  Weight the points ***** 


	if "`weight'" ~= "" { 
		if "`weight'" === "n" { 
			local wtvar "`n'"
		}
		else if "`weight'" == "bar" { 
			tempvar wtvar
			qui gen wtvar = `ub' - `lb'
		}
		else local wtvar = `weight'
		local weight "[fw=`wtvar']"
	}
	

`display'"*****  Save the data ***** 
	cap mac li _savedat
	if _rc == 0 {
		tempfile temp
		qui save "`temp'"
		keep `group' `by' `tvar' `n' `av' `lb' `ub' `wtvar'
		cap noi save `savedat'
		if _rc {	
			di in ye "Make sure that you have put quotation marks round the temporary file name." 
		}
		qui use "`temp'", replace
	}
***** ***** ***** ***** ***** 

`display'"*****  Handle the half-bars ***** 

	if "`half'" ~= "" { 
		tempvar lbs ubs
		qui gen str1 `lbs' = ""
		qui gen str1 `ubs' = ""
		sort `by' `tvar' `av' `group' 
		qui by `by' `tvar' : replace `lb' = `av' if `group' ==  `group'[_N]
		qui by `by' `tvar' : replace `ub' = `av' if `group' ==  `group'[1]

}
***** ***** ***** ***** ***** 


`display'"*****  Sort the data *****
	sort `by' `group' `t2' 

***** ***** ***** ***** ***** 

`display'"*****  Draw the graph *****
* set trace on

	if "`v7'" ~= "" {
		if "`by'" ~= "" { 
				local by "by(`by')" 
		}
		if "`show'" ~= "" {
			cap noi di in wh "version 7: graph `avlist' `lb' `ub' `t2' `weight', connect(`connect') symbol(`symbol') `xlab' `ylab' `by' `options' `saving' "
		}

		if "`graph'" == "" {
			version 7: graph `avlist' `lb' `ub' `t2' `weight', connect(`connect') symbol(`symbol') `xlab' `ylab' `l1title' `by' `saving' `options'
		}
	}

	else  {

		if "`show'" ~= "" {
                	 di "twoway line `avlist' `t2', `lineopts' `ytitle' || `bartype' `lb' `ub' `t2' , `baropts'
		}
		if "`graph'" == "" {
                	 twoway line `avlist' `t2', `lineopts' `ytitle'  || `bartype' `lb' `ub' `t2' , `baropts'
		}
	}


*************************************


end xtgraph

*************************************
*************************************
*************************************


* version 2.2.0  01jul1999
program define centile, rclass
	version 6
	syntax [varlist] [if] [in] [, CCi /*
		*/ Centile(numlist >=0 <=100) /*
		*/ Level(real $S_level) Meansd Normal display]
	
if "`display'" == ""  {
	local display = "*" 
}	
	

	tempvar touse notuse
	mark `touse' `if' `in'
	qui gen byte `notuse' = .  /* will be reset for each variable */

	if "`centile'"=="" { 
		local centile 50 
	}
/*
	Parse `centile' and count the requested centiles, placing them
	into c1, c2, ....
*/
	local nc 0
	tokenize "`centile'"
	while "`1'" != "" {
		local nc = `nc' + 1
		local c`nc' `1'
		local cents "`cents' `1'" /* to return in r() */
		mac shift
	}

	local tl1 "    Obs "
	local ttl "Percentile"
	if "`meansd'"=="" { 
		if "`normal'"=="" { 
			if "`cci'"~="" {
				di in gr _n _col(52) "-- Binomial Exact --"
			}
			else {
				di in gr _n _col(52) "-- Binom. Interp. --"
			}
		}
		else {
			di in gr _n _col(32) /*
			*/ "-- Normal, based on observed centiles --"
		}
	}
	else {
		di in gr _n _col(32) "-- Normal, based on mean and std. dev.--"
	}
	#delimit ;
	di in gr 
	"Variable | `tl1' `ttl'      Centile        [`level'% Conf. Interval]" 
	_n "---------+" _dup(61) "-" ;
	#delimit cr

	local anymark 0
	local alpha2 = (100-`level')/200
	local zalpha2 = -invnorm(`alpha2')
/*
	Run through varlist
*/
	tokenize `varlist'
	local vl
	while "`1'" != "" {
		capt conf str var `1'
		if _rc { 
			local vl "`vl' `1'"
		}
		mac shift 
	}

	local nobs 0			/* in case loop aborted */
	tokenize `vl'
	while "`1'" ~= "" {
		local yvar "`1'"

		/* notuse: 0 = use, 1 = not use -- sort will put use first */
		qui replace `notuse' = ~`touse'
		qui replace `notuse' = 1 if `yvar'==.
		sort `notuse' `yvar'

		qui sum `yvar' if ~`notuse'
		local nobs = r(N)
		local mean = r(mean)
		local sd = sqrt(r(Var))
/*
	Formatting
*/
		local skip = 8 - length("`yvar'")
		local fmt : format `yvar'
		if substr("`fmt'",-1,1)=="f" { 
			local ofmt="%9."+substr("`fmt'",-2,2)
		}
		else if substr("`fmt'",-2,2)=="fc" {
			local ofmt = "%9." + substr("`fmt'",-3,3)
		}
		else	local ofmt "%9.0g"
/*
	Calc required centile(s)
*/
		local j 1
		local s 7
		while `j' <= `nc' {
			local mark ""
			local cj "c`j'"
			local quant = ``cj''/100
			if "`meansd'" ~= "" & (`nobs' > 0) {
			/*
				Normal distribution (parametric estimates)
			*/
				local z = invnorm(`quant')
				local centil = `mean'+`z'*`sd'
				local se = `sd'*sqrt(1/`nobs'+(`z')^2/(2*(`nobs'-1)))
				local cLOWER = `centil'-`zalpha2'*`se'
				local cUPPER = `centil'+`zalpha2'*`se'
			}
			else if `nobs' > 0 {
/*
	If `normal' is not set, calc centile and exact nonparametric confidence
	limits (for example, see Gardner and Altman 1989, pp 72-74),
	interpolating when not at ends of distribution.  An iterative process
	is required for each limit.  As starting values, find ranks for lower
	and upper limits using a normal approximation.

	If `normal' is set, use normal distribution formula for variance, 
	(10.29) in Kendall & Stuart (1969) p 237.

	`alpha2' is for example .025 for a 95% CI.
*/
				local frac1 = (`nobs'+1)*`quant'
`display'"Central fraction and rank = " `quant',`frac1'
				local i1 = int(`frac1')
				local frac1 = `frac1'-`i1'
				if `i1' >= `nobs' {
					local centil = `yvar'[`nobs']
				}
				else if `i1' < 1 {
					local centil = `yvar'[1]
				}
				else {
					local centil = `yvar'[`i1']+ /*
					*/ `frac1'*(`yvar'[`i1'+1]-`yvar'[`i1'])
				}
				if "`normal'" == "" {
					local nq = `nobs'*`quant'
					local z = sqrt(`nq'*(1-`quant'))*`zalpha2'
					local rzLOW = int(.5+`nq'-`z')
					local rzHIGH = 1+int(.5+`nq'+`z')
`display'"lower and upper approx ranks = " `rzLOW',`rzHIGH'
					local r1 `rzHIGH'
					if `r1' > `nobs'+1 { 
						local r1 = `nobs'+1
					}
					local r0 = `r1'-1
					local p0 = Binomial(`nobs',`r0',`quant')
					local p1 = Binomial(`nobs',`r1',`quant')
					local done 0
					while ~`done' {
						if `p0' > `alpha2' {
							if `p1' <= `alpha2' {
								local done 1
							}
							else {
								local r0 = `r1'
								local p0 = `p1'
								local r1 = `r1'+1
								local p1 = Binomial(`nobs',`r1',`quant')
							}
						}
						else if `p0' == `alpha2' {
							local r1 = `r0'
							local p1 = `p0'
							local done 1
						}
						else {
							local r1 = `r0'
							local p1 = `p0'
							local r0 = `r0'-1
							local p0 = Binomial(`nobs',`r0',`quant')
						}
					}
`display'"Upper r0, p0, r1, p1, interp =" `r0',`p0',`r1',`p1',(`p0'-`alpha2')/(`p0'-`p1')
/*
	Upper limit. Interpolate between r0 and r1, r1 being
	conservative. Note that p0>=p1 (upper tail areas).
*/
					if `r0' >= `nobs' {
						local cUPPER = `yvar'[`nobs']
						local mark "*"
						local anymark 1
					}
					else if `r0' < 1 {
						local cUPPER = `yvar'[1]
						local mark "*"
						local anymark 1
					}
					else {
						if "`cci'" == "" {
							local cUPPER = `yvar'[`r0']		/*
							*/	+((`p0'-`alpha2')/(`p0'-`p1'))	/*
							*/	*(`yvar'[`r1']-`yvar'[`r0'])
						}
						else {
							local cUPPER = `yvar'[`r1']
						}
					}
					local r1 `rzLOW'
					if `r1' < 0 { 
						local r1 0 
					}
					local r0 = `r1'-1
					local p0 = 1-Binomial(`nobs',`r0'+1,`quant')
					local p1 = 1-Binomial(`nobs',`r1'+1,`quant')
					local done 0
					while ~`done' {
						if `p1' > `alpha2' {
							if `p0' <= `alpha2' {
								local done 1
							}
							else {
								local r1 = `r0'
								local p1 = `p0'
								local r0 = `r0'-1
								local p0 = 1-Binomial(`nobs',`r0'+1,`quant')
							}
						}
						else if `p0' == `alpha2' {
							local r0 = `r1'
							local p0 = `p1'
							local done 1
						}
						else {
							local r0 = `r1'
							local p0 = `p1'
							local r1 = `r1'+1
							local p1 = 1-Binomial(`nobs',`r1'+1,`quant')
						}
					}
/*
	Interpolate between r1 and r1+1, r1 being conservative.
	Note that p0<=p1 (both are lower tail areas).
*/
					if `r1' >= `nobs' {
						local cLOWER = `yvar'[`nobs']
						local mark "*"
						local anymark 1
					}
					else if `r1' < 1 {
						local cLOWER = `yvar'[1]
						local mark "*"
						local anymark 1
					}
					else {
						if "`cci'" == "" {
							local cLOWER = `yvar'[`r1'] /*
							*/ +((`alpha2'-`p0')/(`p1'-`p0')) /*
							*/ *(`yvar'[`r1'+1]-`yvar'[`r1'])
						}
						else {
							local cLOWER = `yvar'[`r1']
						}
					}
`display'"Lower r0, p0, r1, p1, interp =" `r0',`p0',`r1',`p1',(`alpha2'-`p0')/(`p1'-`p0')
`display'"Var = `yvar', centile `cj'=" `centil' ", CI = " `cLOWER',`cUPPER'
				}
				else {
/*
	Normal distribution, observed centiles
*/
					local dens = exp(-0.5*((`centil'-`mean')/`sd')^2)/(`sd'*sqrt(2*_pi))
					local se = sqrt(`quant'*(1-`quant')/`nobs')/`dens'
					local cLOWER = `centil'-`zalpha2'*`se'
					local cUPPER = `centil'+`zalpha2'*`se'
				}
			}

			if (`j' == 1) & (`nobs' > 0)  {
				di in gr _skip(`skip') "`yvar' |" _col(10) /*
				*/ in yel %8.0f `nobs'		/*
		 		*/ _col(23) %7.0g ``cj''	/*
		 		*/ _col(35) `ofmt' `centil'	/*
		 		*/ _col(51) `ofmt' `cLOWER'	/*
		 		*/ _col(63) `ofmt' `cUPPER' in gr "`mark'"
			}
			else if `nobs' > 0 {
				di in gr "         |"		/*
		 		*/ in yel _col(23) %7.0g ``cj''	/*
		 		*/ _col(35) `ofmt' `centil'	/*
		 		*/ _col(51) `ofmt' `cLOWER'	/*
		 		*/ _col(63) `ofmt' `cUPPER' in gr "`mark'"
			}
			else {  /* `nobs' == 0 */
				if (`j' == 1) {
					di in gr _skip(`skip') "`yvar' |" /*
					*/ _col(10) in yel %8.0f `nobs'
				}
				local centil .
				local cLOWER .
				local cUPPER .
			}
/*
	Store centile in S_# macros, starting at 7.
	Store centiles also in r(c_#) where # starts at 1.
	Store lower and upper bounds on centiles in r(lb_#) and r(ub_#)
*/
			local tmp = `s' - 6
			ret scalar c_`tmp' = `centil'    /* save in r() */
			global S_`s' `centil'      /* also save in S_# */
			/* save confidence limits also in r(), but not S_# */
			ret scalar lb_`tmp' = `cLOWER'
			ret scalar ub_`tmp' = `cUPPER'
			local j = `j'+1
			local s = `s'+1
		}
		mac shift
	}
	if "`anymark'" == "1" {
		di in gr _n /*
*/ "`mark' Lower (upper) confidence limit held at minimum (maximum) of sample"
	}
/*
	Store quantities at final point; S_4 is duplicate of S_[`nc'+6].
*/
	ret scalar N = `nobs'
	ret scalar n_cent = `nc'
	ret local centiles `cents'

	/* double save in S_# */
	global S_1 `nobs'
	global S_2 `nc'
	global S_3 ``cj''
	global S_4 `centil'
	global S_5 `cLOWER'
	global S_6 `cUPPER'
end 


exit