*! version 1.0.2  19apr1996
program define svy_x
	version 4.0
	local cmd "`1'"
	macro shift
	if substr("`1'",1,1)=="," | "`*'"=="" { 
		if "$S_E_cmd"!="svy`cmd'" { 
			error 301
		}
		svy_disp `*'  /* display results */
		exit
	}
	capture noisily {
		tempvar doit w nh subvar
		tempname b V Vdeff Vmeff npop nsub osub error

		svy_pars `cmd' `doit' `w' `nh' `subvar' `npop' `nsub' `osub' `*'

		local  options  "$S_VYopt"
		global S_VYb    "`b'"
		global S_VYv    "`V'"
		global S_VYdeff "`Vdeff'"
		global S_VYmeff "`Vmeff'"
		global S_VYerr  "`error'"

		if "$S_VYfpc"!="" {
			tempname Vdeft
			global S_VYdeft "`Vdeft'"
		}

		if "$S_VYcomp"!="" { /* complete cases */

			if "`cmd'"=="mean" { local cmd "ratio" }

			All`cmd' `doit' `nh' `subvar' $S_VYvl

			Labelmat $S_VYb $S_VYv $S_VYosub $S_VYnsub /*
			*/	 $S_VYdeff $S_VYdeft $S_VYmeff $S_VYerr

			local cmd = upper(substr("$S_VYcmd",1,1)) + /*
			*/	    substr("$S_VYcmd",2,.)
			
			local dof = $S_VYnpsu - $S_VYnstr

			matrix post `b' `V', /*
			*/	dof(`dof') obs($S_VYnobs) depn(`cmd')
		}
		else { /* available; do one variable at a time */
			tempname mstr mpsu
			global S_VYmstr "`mstr'"
			global S_VYmpsu "`mpsu'"

			OneByOne `doit' `subvar' $S_VYvl

			Labelmat $S_VYb $S_VYv $S_VYosub $S_VYnsub /*
			*/	 $S_VYdeff $S_VYdeft $S_VYmeff     /* 
			*/	 $S_VYerr $S_VYmstr $S_VYmpsu 
			
			discard  /* clear previous estimates */
		}

		svy_pars save  /* save macros */
	}
	nobreak {
		if _rc {
			local rc = _rc
			if `rc'==902 & "$S_VYcomp"!="" {
				di in red "try using " _quote "available" _quote " option"
			}
			svy_pars cleanup
			exit `rc'
		}
	}

	svy_disp, `options'  /* display results */
end

program define Allratio /* markvar nh subvar varlist */
/*
   Main program for complete-case mean and ratio.
*/
	version 4.0
	local doit "`1'"
	local nh   "`2'"
	if "$S_VYnby"!="" { local subvar "`3'" }
	macro shift 3
	tempname r s

/* Preserve and drop. */

	preserve
	qui keep if `doit'
	keep `*' $S_VYstr $S_VYpsu $S_VYfpc $S_VYw `nh' `subvar'

/* Initialize stuff. */
		
	InitMat  /* initialize matrices */

	if "$S_VYcmd"=="mean" {
		local dy 1
		local x "one"
		local one 1
		local ro "*"  /* do for ratio only */
	}
	else {
		local dy 2
		local mo "*"  /* do for mean only */
	}

/* Compute covariance for meff. */

	AllMeff `subvar' `*'

/* Loop over variables and construct differences. */

	if "$S_VYnby"=="" {
		`mo' scalar `s' = $S_VYnpop
		local i 1
		local y 1
		while "``y''"!="" {
			GenSum `r' $S_VYw*``y''

			`ro' local x = `y' + 1
			`ro' GenSum `s' $S_VYw*``x''

			tempvar d
			if `s'==0 { /* set error */
				if $S_VYnpop==0 {
					matrix $S_VYerr[1,`i'] = 3
				}
				else matrix $S_VYerr[1,`i'] = 4
				qui gen byte `d' = 0
			}
			else {
				scalar `r' = `r'/`s'
				matrix $S_VYb[1,`i'] = `r'
				qui gen $S_VYtype `d' = (``y''-`r'*``x'')/`s'
			}
			local vars "`vars' `d'"
			`mo' drop ``y''
			local i = `i' + 1
			local y = `y' + `dy'
		}
		AllVsrs `vars'      /* compute SRS covariance */
		AllVar `nh' `vars'  /* compute covariance     */
		exit
	}

/* If here, there are subpops. */

	local i 1
	local y 1
	while "``y''"!="" {
		`ro' local x = `y' + 1
		local isub 1
		while `isub' <= $S_VYnby {
			local sub "(`subvar'==`isub')"

			GenSum `r' `sub'*$S_VYw*``y''
			
			`mo' scalar `s' = $S_VYnsub[1,`i']
			`ro' GenSum `s' `sub'*$S_VYw*``x''

			tempvar d
			if `s'==0 { /* set error */
				if $S_VYnsub[1,`i']==0 {
					matrix $S_VYerr[1,`i'] = 3
				}
				else matrix $S_VYerr[1,`i'] = 4
				qui gen byte `d' = 0
			}
			else {
				scalar `r' = `r'/`s'
				matrix $S_VYb[1,`i'] = `r'
				qui gen $S_VYtype `d' = `sub' /*
				*/	*(``y''-`r'*``x'')/`s'
			}
			local vars "`vars' `d'"
			local i = `i' + 1
			local isub = `isub' + 1
		}
		`mo' drop ``y''
		local y = `y' + `dy'
	}

	AllVsrs `vars'      /* compute SRS covariance */
	AllVar `nh' `vars'  /* compute covariance     */
end	

program define Alltotal /* markvar nh subvar varlist */
/*
   Main program for complete-case total.
*/
	version 4.0
	local doit "`1'"
	local nh   "`2'"
	if "$S_VYnby"!="" { local subvar "`3'" }
	macro shift 3

/* Preserve and drop. */

	preserve
	qui keep if `doit'
	keep `*' $S_VYstr $S_VYpsu $S_VYfpc $S_VYw `nh' `subvar'

	InitMat  /* initialize matrices */

/* Compute covariance for meff. */

	AllMeff `subvar' `*'

/* Compute estimates and covariance. */

	if "$S_VYnby"=="" {
		AllVsrsT `*'     /* compute SRS covariance */
		AllVar `nh' `*'  /* compute covariance     */
		exit
	}

/* If here, there are subpops. */

	local y 1
	while "``y''"!="" {
		local type : type ``y''
		local isub 1
		while `isub' <= $S_VYnby {
			tempvar t
			qui gen `type' `t' = (`subvar'==`isub')*``y''
			local vars "`vars' `t'"
			local isub = `isub' + 1
		}
		drop ``y''
		local y = `y' + 1
	}

	AllVsrsT `subvar' `vars'  /* compute SRS covariance */
	AllVar `nh' `vars'        /* compute covariance     */
end

program define OneByOne /* markvar subvar varlist */
/*
   Main program for variable-by-variable mean, ratio, and total.
*/
	version 4.0
	local doit   "`1'"
	local subvar "`2'"
	macro shift 2

/* Initialize stuff. */

	local prg = upper(substr("$S_VYcmd",1,1))

	if "$S_VYcmd"!="ratio" {
		local dy 1
		local ro "*"  /* do for ratio only */
	}
	else    local dy 2

	InitMat  /* initialize matrices */

/* Step through variables. */

	if "$S_VYnby"=="" {
		local i 1
		local y 1
		while "``y''"!="" {
			`ro' local x = `y' + 1
			`ro' local x "``x''"
			OneEst`prg' `doit' 1 ``y'' `x' `i'
			local i = `i' + 1
			local y = `y' + `dy'
		}
		exit
	}

/* If here, there are subpops. */

	local i 1
	local y 1
	while "``y''"!="" {
		`ro' local x = `y' + 1
		`ro' local x "``x''"
		local isub 1
		while `isub' <= $S_VYnby {
			OneEst`prg' `doit' (`subvar'==`isub') ``y'' `x' `i'
			local i = `i' + 1
			local isub = `isub' + 1
		}
		local y = `y' + `dy'
	}
end

/******************************* Subprograms *******************************/

program define InitMat /* initializes matrices */
	version 4.0

	if "$S_VYby"=="" { local dim "$S_VYnvl" }
	else		   local dim = $S_VYnby*$S_VYnvl

	matrix $S_VYb = J(1,`dim',0)
	matrix $S_VYerr  = J(1,`dim',0)

	if "$S_VYcomp"!="" { exit }

	matrix $S_VYv = J(1,`dim',0)
	matrix $S_VYdeff = J(1,`dim',0)
	matrix $S_VYmeff = J(1,`dim',0)
	matrix $S_VYnsub = J(1,`dim',0)
	matrix $S_VYosub = J(1,`dim',0)
	matrix $S_VYmstr = J(1,`dim',0)
	matrix $S_VYmpsu = J(1,`dim',0)
	if "$S_VYfpc"!="" { matrix $S_VYdeft = J(1,`dim',0) }
end

program define OneEstT /* markvar subpop variable index */
/*
   For total, computes #obs, population size, total, variance, meff,
   and deff for a single variable. 
*/
	version 4.0
	local doit "`1'"
	local sub  "`2'"
	local y    "`3'"
	local i    "`4'"
	tempname s w

/* Get number of obs. */

	qui count if `doit' & `y'!=.
	local nobs = _result(1)
	if `nobs'!=$S_VYnobs { 
		global S_VYmiss "missing"
	}

	qui count if `doit'*`sub' & `y'!=.
	if _result(1)==0 {
		matrix $S_VYerr[1,`i'] = 1
		exit
	}
	local o = _result(1)
	matrix $S_VYosub[1,`i'] = `o'

/* Compute subpopulation size. */

	GenSum `w' `doit'*`sub'*(`y'!=.)*$S_VYw
	matrix $S_VYnsub[1,`i'] = `w'

/* Compute variance for meff. */

	GenSum `s' `doit'*`sub'*`y'
	scalar `s' = `s'/$S_VYosub[1,`i']
	GenSum `s' `doit'*`sub'*(`y'-`s')^2
	matrix $S_VYmeff[1,`i'] = cond($S_VYosub[1,`i'] > 1, /*
	*/ `s'*$S_VYnsub[1,`i']^2/($S_VYosub[1,`i']*($S_VYosub[1,`i']-1)), 0)

/* Compute point estimate. */

	GenSum `s' `doit'*`sub'*$S_VYw*`y'
	matrix $S_VYb[1,`i'] = `s'

/* Compute variance. */

	OneVar `doit' `sub'*`y' `i'
	
/* Compute variance for deff and deft. */

	if "`sub'"!="1" & "$S_VYsrss"=="" {
			local o "`nobs'"
			GenSum `w' `doit'*(`y'!=.)*$S_VYw
	}
	
	if `w'==0 | `o'<=1 { exit }
	
	if "$S_VYsrss"=="" {
		GenSum `s' `doit'*$S_VYw*(`sub'*`y'-`s'/`w')^2
	}
	else    GenSum `s' `doit'*`sub'*$S_VYw*(`y'-`s'/`w')^2
	
	scalar `s' = `s'*`w'/(`o' - 1)
	
	if "$S_VYfpc"=="" {
		matrix $S_VYdeff[1,`i'] = `s'
		exit
	}
	
	matrix $S_VYdeft[1,`i'] = `s'
	matrix $S_VYdeff[1,`i'] = (1 - `o'/`w')*`s'
end

program define OneEstM /* markvar subpop variable index */
/*
   For mean, computes #obs, population size, mean, variance, meff,
   and deff for a single variable. 
*/
	version 4.0
	local doit "`1'"
	local sub  "`2'"
	local y    "`3'"
	local i    "`4'"
	tempname m w

/* Get number of obs. */

	qui count if `doit' & `y'!=.
	local nobs = _result(1)
	if `nobs'!=$S_VYnobs { 
		global S_VYmiss "missing"
	}

	qui count if `doit'*`sub' & `y'!=.
	if _result(1)==0 {
		matrix $S_VYerr[1,`i'] = 1
		exit
	}
	local o = _result(1)
	matrix $S_VYosub[1,`i'] = `o'

/* Compute variance for meff. */

	GenSum `m' `doit'*`sub'*`y'
	scalar `m' = `m'/$S_VYosub[1,`i']
	GenSum `m' `doit'*`sub'*(`y'-`m')^2
	matrix $S_VYmeff[1,`i'] = cond($S_VYosub[1,`i'] > 1, /*
	*/ `m'/($S_VYosub[1,`i']*($S_VYosub[1,`i']-1)), 0)

/* Compute subpopulation size. */

	GenSum `w' `doit'*`sub'*(`y'!=.)*$S_VYw
	if `w'==0 {
		matrix $S_VYerr[1,`i'] = 3
		exit
	}
	matrix $S_VYnsub[1,`i'] = `w'

/* Compute point estimate. */

	GenSum `m' `doit'*`sub'*$S_VYw*`y'
	scalar `m' = `m'/`w'
	matrix $S_VYb[1,`i'] = `m'

/* Compute variance. */

	OneVar `doit' `sub'*(`y'-`m')/`w' `i'
	
/* Compute variance for deff and deft. */

	GenSum `m' `doit'*`sub'*$S_VYw*((`y'-`m')/`w')^2

	if "`sub'"!="1" & "$S_VYsrss"=="" {
			local o "`nobs'"
			GenSum `w' `doit'*(`y'!=.)*$S_VYw
	}
	
	if `w'==0 | `o'<=1 { exit }
	
	scalar `m' = `m'*`w'/(`o' - 1)
	
	if "$S_VYfpc"=="" {
		matrix $S_VYdeff[1,`i'] = `m'
		exit
	}
	
	matrix $S_VYdeft[1,`i'] = `m'
	matrix $S_VYdeff[1,`i'] = (1 - `o'/`w')*`m'
end

program define OneEstR /* markvar subpop variable1 variable2 index */
/*
   For ratio, computes #obs, population size, ratio, variance, meff,
   and deff for a ratio of two variables.
*/
	version 4.0
	local doit "`1'"
	local sub  "`2'"
	local y    "`3'"
	local x    "`4'"
	local i    "`5'"
	tempname r s w
	tempvar doyx
	mark `doyx' if `doit'
	markout `doyx' `y' `x'

/* Get number of obs. */

	qui count if `doyx'
	local nobs = _result(1)
	if `nobs'!=$S_VYnobs { 
		global S_VYmiss "missing"
	}

	qui count if `doyx'*`sub'
	if _result(1)==0 {
		matrix $S_VYerr[1,`i'] = 1
		exit
	}
	local o = _result(1)
	matrix $S_VYosub[1,`i'] = `o'

/* Compute variance for meff. */

	GenSum `r' `doyx'*`sub'*`y'
	GenSum `s' `doyx'*`sub'*`x'
	if `s'!=0 {
		scalar `r' = `r'/`s'
		GenSum `s' `doyx'*`sub'*((`y'-`r'*`x')/`s')^2
		matrix $S_VYmeff[1,`i'] = cond($S_VYosub[1,`i'] > 1, /*
		*/ `s'*$S_VYosub[1,`i']/($S_VYosub[1,`i']-1), 0)
	}

/* Compute subpopulation size. */

	GenSum `w' `doyx'*`sub'*$S_VYw
	if `w'==0 {
		matrix $S_VYerr[1,`i'] = 3
		exit
	}
	matrix $S_VYnsub[1,`i'] = `w'

/* Compute point estimate. */

	GenSum `r' `doyx'*`sub'*$S_VYw*`y'
	GenSum `s' `doyx'*`sub'*$S_VYw*`x'
	if `s'==0 {
		if `w'== 0 {
			matrix $S_VYerr[1,`i'] = 3
		}
		else matrix $S_VYerr[1,`i'] = 4
		exit
	}
	scalar `r' = `r'/`s'
	matrix $S_VYb[1,`i'] = `r'

/* Compute variance. */

	OneVar `doit' `sub'*(`y'-`r'*`x')/`s' `i'
	
/* Compute variance for deff and deft. */

	GenSum `s' `doyx'*`sub'*$S_VYw*((`y'-`r'*`x')/`s')^2

	if "`sub'"!="1" & "$S_VYsrss"=="" {
			local o "`nobs'"
			GenSum `w' `doyx'*$S_VYw
	}
	
	if `w'==0 | `o'<=1 { exit }
	
	scalar `s' = `s'*`w'/(`o' - 1)
	
	if "$S_VYfpc"=="" {
		matrix $S_VYdeff[1,`i'] = `s'
		exit
	}
	
	matrix $S_VYdeft[1,`i'] = `s'
	matrix $S_VYdeff[1,`i'] = (1 - `o'/`w')*`s'
end

program define OneVar /* markvar variable index_of_variable */
/* 
   Computes variance for a single variable.
*/
	version 4.0
	local doit "`1'"
	local y    "(`2')"
	local i    "`3'"
	tempvar nh yhi yh
	tempname s

	if "$S_VYfpc"!="" {
		capture assert $S_VYfpc <= 1
		if _rc==0 { local fpc "(1-$S_VYfpc)*" }
		else        local fpc "(1-`nh'[_N]/$S_VYfpc)*"
	}

	quietly {

	/* Compute total #strata. */
	
		by `doit' $S_VYstr: gen byte `nh' = /*
		*/	cond(_n==_N, sum(`y'!=.)>0, .) if `doit'
		count if `nh' == 1
		matrix $S_VYmstr[1,`i'] = _result(1)
	
	/* `yhi' = (sum over PSU i) weight*variable */
	
		if "$S_VYpsu"!="" {
			by `doit' $S_VYstr $S_VYpsu: replace `nh' = /*
			*/	cond(_n==_N, sum(`y'!=.)>0, .) if `doit'

			by `doit' $S_VYstr $S_VYpsu: gen double `yhi' = /*
			*/	cond(_n==_N & `nh', sum($S_VYw*`y'), .) /*
			*/	if `doit'
		}
		else {
			replace `nh' = (`y'!=.) if `doit'
			gen double `yhi' = $S_VYw*`y' if `doit'
		}
		
		count if `nh' == 1
		matrix $S_VYmpsu[1,`i'] = _result(1) /* total #PSUs */

	/* `yh' = (sum over stratum h) `yhi' */

		if "$S_VYstr"!="" {
			by `doit' $S_VYstr: replace `nh' = sum(`nh') if `doit'
			capture by `doit' $S_VYstr: assert `nh'!=1 if _n==_N
			if _rc {
				matrix $S_VYerr[1,`i'] = 2
				exit
			}

			by `doit' $S_VYstr: gen double `yh' = sum(`yhi') /*
			*/	if `doit'

			by `doit' $S_VYstr: replace `yhi' = /*
			*/	`fpc'(`nh'[_N]/(`nh'[_N]-1)) /*
			*/	*(`yhi'-`yh'[_N]/`nh'[_N])^2 if `doit'
		}
		else {
			replace `nh' = sum(`nh')
			if `nh'[_N]==1 {
				matrix $S_VYerr[1,`i'] = 2
				exit
			}
			gen double `yh' = sum(`yhi')
			replace `yhi' = `fpc'(`nh'[_N]/(`nh'[_N]-1)) /*
			*/	*(`yhi'-`yh'[_N]/`nh'[_N])^2 if `doit'
		}
	}

	GenSum `s' `yhi'  /* sum to compute variance */

	matrix $S_VYv[1,`i'] = `s'
end

program define AllVar /* nh varlist */
/* 
   Computes complete-case covariance.
*/
	version 4.0
	local nh "`1'"
	macro shift

	if "$S_VYcmd"!="total" {
		local to "*"  /* do for total only */
	}
	`to' tempname s

	sort $S_VYstr $S_VYpsu $S_VYw

	local y 1
	while "``y''"!="" {
		`to' GenSum `s' $S_VYw*``y''  /* compute total */
		`to' matrix $S_VYb[1,`y'] = `s'
				
		GenDiff `nh' ``y''  /* generate differences */
		local y = `y' + 1
	}

	quietly {
		if "$S_VYpsu"!="" {
			by $S_VYstr $S_VYpsu: keep if _n==_N
		}
		
		if "$S_VYfpc"!="" {
			capture assert $S_VYfpc <= 1
			if _rc==0 { local fpc "(1-$S_VYfpc)*" }
			else        local fpc "(1-`nh'/$S_VYfpc)*"
		}

		matrix accum $S_VYv = `*' [iw=`fpc'`nh'/(`nh'-1)], /*
		*/	nocons
	}
end	

program define GenDiff /* nh varname */
/*
   Computes differences for complete-case covariance computation.
*/
	version 4.0
	local nh "`1'"
	local y  "`2'"
	tempvar yhi yh
	quietly {

	/* `yhi' = (sum over PSU i) weight*variable */

		if "$S_VYpsu"!="" {
			by $S_VYstr $S_VYpsu: gen double `yhi' = /*
			*/	cond(_n==_N, sum($S_VYw*`y'), .)
		}
		else gen double `yhi' = $S_VYw*`y'

		drop `y'

	/* `yh' = (sum over stratum h) `yhi' */

		if "$S_VYstr"!="" {
			by $S_VYstr: gen double `yh' = sum(`yhi')

			by $S_VYstr: gen $S_VYtype `y' = `yhi'-`yh'[_N]/`nh'
		}
		else {
			gen double `yh' = sum(`yhi')
			gen $S_VYtype `y' = `yhi'-`yh'[_N]/`nh'
		}
	}
end

program define AllVsrs /* varlist */
/*
   Computes Vsrs for mean and ratio.
*/
	version 4.0
	if $S_VYnpop <= 0 {
		local dim : word count `*'
		matrix $S_VYdeff = J(`dim',`dim',0)
		if "$S_VYfpc"=="" { exit }
		matrix $S_VYdeft = J(`dim',`dim',0)
		exit
	}
	
	qui matrix accum $S_VYdeff = `*' [iw=$S_VYw*$S_VYnobs/$S_VYnpop], /*
	*/	nocons
	
	if "$S_VYsrss"=="" {
		tempname s
		scalar `s' = $S_VYnpop^2/($S_VYnobs*($S_VYnobs - 1))
		
		if "$S_VYfpc"=="" { 
			matrix $S_VYdeff = `s'*$S_VYdeff
			exit
		}
		
		matrix $S_VYdeft = `s'*$S_VYdeff
		scalar `s' = 1 - $S_VYnobs/$S_VYnpop
		matrix $S_VYdeff = `s'*$S_VYdeft
		exit
	}

/* If here, we have subpops with srssubpop. */

	tempname At A I
	matrix `At' = J(1,$S_VYnby,0)
	matrix `A'  = J(1,$S_VYnby,0)
	local dim = colsof($S_VYdeff)/$S_VYnby
	matrix `I' = I(`dim')

	local isub 1
	while `isub' <= $S_VYnby {
		matrix `At'[1,`isub'] = cond($S_VYosub[1,`isub'] > 1, /*
		*/	$S_VYnsub[1,`isub']*$S_VYnpop /*
		*/	/($S_VYnobs*($S_VYosub[1,`isub'] - 1)), 0)

		if "$S_VYfpc"!="" {
			matrix `A'[1,`isub'] = cond($S_VYnsub[1,`isub'] > 0, /*
			*/	(1 - $S_VYosub[1,`isub'] /*
			*/	/$S_VYnsub[1,`isub'])*`At'[1,`isub'], 0)
		}
		local isub = `isub' + 1
	}
	
	matrix `At' = diag(`At')
	matrix `At' = `I' # `At'
	if "$S_VYfpc"=="" {
		matrix $S_VYdeff = `At'*$S_VYdeff
		exit
	}

	matrix $S_VYdeft = `At'*$S_VYdeff
	matrix `A' = diag(`A')
	matrix `A' = `I' # `A'
	matrix $S_VYdeff = `A'*$S_VYdeff
end

program define AllVsrsT /* [subvar] varlist */
/*
   Computes Vsrs for total.
*/
	version 4.0
	tempname s c0
	if "$S_VYnby"!="" {
		local subvar "`1'"
		macro shift
	}
	if $S_VYnpop <= 0 {
		local dim : word count `*'
		matrix $S_VYdeff = J(`dim',`dim',0)
		if "$S_VYfpc"=="" { exit }
		matrix $S_VYdeft = J(`dim',`dim',0)
		exit
	}
	if "$S_VYsrss"=="" {
		qui matrix accum $S_VYdeff = `*' /*
		*/ [iw=$S_VYw*$S_VYnobs/$S_VYnpop], nocons dev

	/* If y = constant, make Vsrs a hard zero. */
	
		local y 1
		while "``y''"!="" {
			capture assert ``y''==``y''[1]
			if _rc==0 { /* zero Vsrs */
				local dim = colsof($S_VYdeff)
				matrix `c0' = J(`dim',1,0)
				matrix substit $S_VYdeff[1,`y'] = `c0'
				matrix `c0' = `c0''
				matrix substit $S_VYdeff[`y',1] = `c0'
			}
			local y = `y' + 1
		}
		
		scalar `s' = $S_VYnpop^2/($S_VYnobs*($S_VYnobs - 1))
		
		if "$S_VYfpc"=="" {
			matrix $S_VYdeff = `s'*$S_VYdeff
			exit
		}
		
		matrix $S_VYdeft = `s'*$S_VYdeff
		scalar `s' = 1 - $S_VYnobs/$S_VYnpop
		matrix $S_VYdeff = `s'*$S_VYdeft
		exit
	}
	
/* If here, we have subpops with srssubpop. */

	tempname scale V Vnf Vf
	scalar `scale' = $S_VYnobs/$S_VYnpop
	local isub 1
	while `isub' <= $S_VYnby {
		qui matrix accum `V' = `*' [iw=`scale'*$S_VYw] /*
		*/	if `subvar'==`isub', nocons dev
		
		scalar `s' = cond($S_VYosub[1,`isub'] > 1, /*
		*/	$S_VYnsub[1,`isub'] /*
		*/	/(`scale'*($S_VYosub[1,`isub'] - 1)), 0)
		matrix `V' = `s'*`V'
		
		if `isub'== 1 { matrix `Vnf' = `V' }
		else matrix `Vnf' = `Vnf' + `V'
		
		if "$S_VYfpc"!="" {
			scalar `s' = cond($S_VYnsub[1,`isub'] > 0, /*
			*/	(1 - $S_VYosub[1,`isub'] /*
			*/	/$S_VYnsub[1,`isub']), 0)
			matrix `V' = `s'*`V'
			if `isub'== 1 { matrix `Vf' = `V' }
			else matrix `Vf' = `Vf' + `V'
		}
		local isub = `isub' + 1
	}
	
/* If y = constant, make Vsrs a hard zero. */
 
	local y 1
	while "``y''"!="" {
		local isub = mod(`y'-1,$S_VYnby) + 1
		capture assert ``y''==``y''[_n-1] /*
		*/	if `subvar'==`isub' & `subvar'[_n-1]==`isub'
		if _rc==0 { /* zero Vsrs */
			local dim = colsof(`Vnf')
			matrix `c0' = J(`dim',1,0)
			matrix substit `Vnf'[1,`y'] = `c0'
			cap matrix substit `Vf'[1,`y']  = `c0'
			matrix `c0' = `c0''
			matrix substit `Vnf'[`y',1] = `c0'
			cap matrix substit `Vf'[`y',1]  = `c0'
		}
		local y = `y' + 1
	}
	
	if "$S_VYfpc"=="" {
		matrix $S_VYdeff = `Vnf'
		exit
	}
	
	matrix $S_VYdeft = `Vnf'
	matrix $S_VYdeff = `Vf'
end

program define GenSum  /* scalarname expression */
/*
   Computes accurate sum: scalarname = sum(expression).
*/
	version 4.0
	local s "`1'"
	macro shift
	tempvar y e
	quietly {
		gen double `y' = sum(`*')
		gen double `e' = sum((`*')-(`y'-`y'[_n-1])) in 2/l
		scalar `s' = `y'[_N] + `e'[_N]
	}
end

program define AllMeff /* [subvar] varlist */
/*
   Main program for computing meff covariance for mean, ratio, and total.
*/
	version 4.0
	if "$S_VYcmd"=="ratio" {
		local cmd "rat"
	}
	else local cmd "mean"
	if "$S_VYnby"=="" {
		Meff`cmd' $S_VYmeff 1 `*'
		if "$S_VYcmd"!="total" { exit }
		tempname s
		scalar `s' = $S_VYnpop*$S_VYnpop
		matrix $S_VYmeff = `s'*$S_VYmeff
		exit
	}
	local subvar "`1'"
	macro shift
	tempname K V s
	matrix `K' = J($S_VYnby,$S_VYnby,0)
	local isub 1
	while `isub' <= $S_VYnby {
		Meff`cmd' `V' `subvar'==`isub' `*'
		if "$S_VYcmd"=="total" {
			scalar `s' = $S_VYnsub[1,`isub']^2
			matrix `V' = `s'*`V'
		} 
		matrix `K'[`isub',`isub'] = 1
		matrix `V' = `V' # `K'
		matrix `K'[`isub',`isub'] = 0
		if `isub'==1 {
			matrix $S_VYmeff = `V'
		}
		else matrix $S_VYmeff = $S_VYmeff + `V'
		local isub = `isub' + 1
	}
end

program define Meffrat /* matname expression varlist */
/*
   Subprogram for computing meff covariance for ratio.
*/
	version 4.0
	local V   "`1'"
	local exp "`2'"
	macro shift 2
	tempname m r s
	qui matrix accum `V' = `*' if `exp', mean(`m') dev nocons
	local cdim = rowsof(`V')
	local rdim = `cdim'/2
	matrix `r' = J(`rdim',`cdim',0)
	local i 1
	while `i' <= `rdim' {
		local j2 = 2*`i'
		local j1 = `j2' - 1
		if `m'[1,`j2']!=0 {
			matrix `r'[`i',`j1'] = 1/`m'[1,`j2']
			matrix `r'[`i',`j2'] = -`m'[1,`j1']/(`m'[1,`j2']^2)
		}
		local i = `i' + 1
	}
	matrix `V' = `V'*`r''
	matrix `V' = `r'*`V'
	scalar `s' = cond(_result(1) > 1, 1/(_result(1)*(_result(1)-1)), 0)
	matrix `V' = `s'*`V'
end

program define Meffmean /* matname expression varlist */
/*
   Subprogram for computing meff covariance for mean.
*/
	version 4.0
	local V   "`1'"
	local exp "`2'"
	macro shift 2
	tempname s
	qui matrix accum `V' = `*' if `exp', dev nocons
	scalar `s' = cond(_result(1) > 1, 1/(_result(1)*(_result(1)-1)), 0)
	matrix `V' = `s'*`V'
end

program define Labelmat /* list of matrices */
	version 4.0
	GetNames  /* get equation and column names */
	local M 1
	while "``M''"!="" {
		matrix coleq   ``M'' = $S_VYceq
		matrix colname ``M'' = $S_VYcnam
		if rowsof(``M'')==colsof(``M'') {
			matrix roweq   ``M'' = $S_VYceq
			matrix rowname ``M'' = $S_VYcnam
		}
		local M = `M' + 1
	}
end

program define GetNames
	version 4.0
	if "$S_VYby"=="" {
		if "$S_VYcmd"!="ratio" {
			global S_VYceq  "_"
			global S_VYcnam "$S_VYvl"
			exit
		}

	/* Do ratio with no by() subpops. */

		parse "$S_VYvl", parse(" ")
		local y 1
		while "``y''"!="" {
			local x = `y' + 1
			global S_VYceq "$S_VYceq ``y''"
			global S_VYcnam "$S_VYcnam ``x''"
			local y = `y' + 2
		}
		exit
	}

/* If here, there are by() subpops. */

	local nvar : word count $S_VYby
	if `nvar'==1 { /* there is one subpop label */
		local twoby "*"
		local cnam "$S_VYcnam"
		global S_VYcnam
	}
	else { /* use numbers (isub) for subpop labels */
		local isub 1
		while `isub' <= $S_VYnby {
			local cnam "`cnam' `isub'"
			local isub = `isub' + 1
		}
	}

	if "$S_VYcmd"!="ratio" {
		parse "$S_VYvl", parse(" ")
		local ro "*"  /* do for ratio only */
	}
	else local mo "*"  /* do for mean and total only */
	local y 1
	while `y' <= $S_VYnvl {
		global S_VYcnam "$S_VYcnam `cnam'"
		local isub 1
		while `isub' <= $S_VYnby {
			`mo' global S_VYceq "$S_VYceq ``y''"
			`ro' global S_VYceq "$S_VYceq `y'"
			local isub = `isub' + 1
		}
		local y = `y' + 1
	}
end
exit

error codes:

	 1 = no obs (only possible with available)
	 2 = nh is 1 (only possible with available; complete aborts earlier)
	 3 = sum of weights is 0
	 4 = denominator is 0 for ratio 
exit