*! version 1.0.2  16apr1996
program define svy_reg /* markvar nh residuals varlist */
/*
   Driver program for survey-design estimation commands.

   Returns in S_VY* global macros:

	S_VYnobs
	S_VYnstr
	S_VYnpsu
	S_VYnpop
	S_VYnsub  [if subpop()]
	S_VYosub  [if subpop()]

   Other S_VY macros are used for argument passing. 
*/
	version 4.0
	local D  "`1'"
	macro shift
	local varlist "req ex"
	local weight "pweight iweight noprefix"
	local if "opt"
	local in "opt"
	#delimit ;
	local options "Var(string) STRata(string) PSU(string) FPC(string)
		       DEFF(string) DEFT(string) SUBpop(string) SRSsubpop
		       FLOAT noCONStant noMINUS" ;
	#delimit cr
	parse "`*'"
	parse "`varlist'", parse(" ")
	local res "`1'"
	macro shift

/* Set global macros and arguments. */

	global S_VYv    "`var'"
	global S_VYstr  "`strata'"
	global S_VYpsu  "`psu'"
	global S_VYfpc  "`fpc'"
	global S_VYdeff "`deff'"
	global S_VYdeft "`deft'"

	if "`float'"=="" { global S_VYtype "double" }
	else               global S_VYtype "float"

	if "`minus'"!="" { local minus +0 }
	else local minus -1

	if "`subpop'"!="" {
		tempname nsub
		global S_VYnsub "`nsub'"
	}
	if "`srssubp'"!="" & "`subpop'"!="" {
		local srssubp 1
	}
	else local srssubp 0
	
	tempname npop

	global S_VYnpop "`npop'"

/* Preserve. */

	preserve
	if "`if'`in'"!="" { qui keep `if' `in' }

	keep `*' `res' `exp' `strata' `psu' `fpc' `subpop'

/* Sort. */

	if "`strata'`psu'`exp'"!="" {
		sort `strata' `psu' `exp'
	}

	if "`exp'"!="" { global S_VYw "`exp'" }
	else		 global S_VYw 1

/* Generate `nh' variable and numbers of strata, psu, population. */

	tempvar nh z

	Gen_nh `nh' `subpop' `nsub'

/* Construct new varlist. */

	local x 1
	while "``x''"!="" {
		qui gen $S_VYtype `z' = `res'*``x''
		drop ``x''
		rename `z' ``x''
		local x = `x' + 1
	}

	if "`constan'"!="" { /* no constant term */
		local res
	}

	AllVsrs `srssubp' `*' `res'    /* compute SRS covariance    */
	AllVar `minus' `nh' `*' `res'  /* compute design-covariance */
	MakeVar `D'	               /* compute D*V*D             */

	global S_VYnpop = $S_VYnpop    /* make scalar a macro */
	if "$S_VYnsub"!="" {
		global S_VYnsub = $S_VYnsub    /* make scalar a macro */
	}
end

program define Gen_nh
	version 4.0
	local nh     "`1'"
	local subpop "`2'"
	local nsub   "`3'"
	if "$S_VYstr"!="" { local bystr "by $S_VYstr:" }

	quietly {

	/* Set total #obs. */

		global S_VYnobs = _N

	/* Compute total #strata. */
	
		`bystr' gen byte `nh' = (_n==1)
		count if `nh' == 1
		global S_VYnstr = _result(1)
		
	/* `nh' = #PSU in stratum h */
	
		if "$S_VYpsu"!="" {
			by $S_VYstr $S_VYpsu: replace `nh' = (_n==1)
		}
		else { /* observations are PSUs */
			replace `nh' = 1
		}

	/* Compute total #PSU. */

		count if `nh' == 1
		global S_VYnpsu = _result(1)

	/* Finish `nh' computation.  Note: automatic promotion of type. */

		`bystr' replace `nh' = sum(`nh')
		`bystr' replace `nh' = `nh'[_N]

	/* Check if `nh' >= 2. */ 

		capture assert `nh' >= 2
		if _rc {
			di in red "stratum with only one PSU detected"
			exit 499
		}

	/* Check fpc variable for sensible ranges. */

		if "$S_VYfpc"!="" {
			CheckFPC $S_VYstr $S_VYfpc `nh'
		}

	/* Compute sum of weights. */

		GenSum $S_VYnpop $S_VYw
		
		if "`subpop'"=="" { exit }
		
	/* If here, we have a subpopulation. */
	
		count if `subpop'
		global S_VYosub = _result(1)
		GenSum $S_VYnsub `subpop'*$S_VYw
	}
end

program define CheckFPC
	version 4.0
	local strata "`1'"
	local fpc    "`2'"
	local nh     "`3'"
	
	capture assert `fpc' >= 0
	if _rc {
		di in red "fpc() must be >= 0"
		exit 499
	}
	capture by `strata': assert abs((`fpc' - `fpc'[1]) /*
	*/	/max(`fpc'[1],1)) < 1e-5
	if _rc {
		di in red "all observations within a stratum must have " /*
		*/ "the same fpc() value"
		exit 499
	} 
	capture assert `fpc' >= `nh'
	if _rc {
		capture assert `fpc' <= 1
		if _rc {
			di in red "fpc() must be <= 1 if a rate " /*
			*/ "or >= # sampled PSUs per stratum" 
			exit 499
		}
	}
end

program define AllVar /* minus nh varlist */
/* 
   Computes survey-design covariance.
*/
	version 4.0
	local minus "`1'"
	local nh    "`2'"
	macro shift 2

	local y 1
	while "``y''"!="" {
		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'`minus')], /*
		*/	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 /* srssubpop varlist */
/*
   Computes Vsrs.
*/
	version 4.0
	local srssubp "`1'"
	macro shift
	tempname s f
	
	qui matrix accum $S_VYdeff = `*' [iw=$S_VYw*$S_VYnobs/$S_VYnpop], /*
	*/	nocons
	
	if `srssubp'==0 {
		scalar `s' = $S_VYnpop^2/($S_VYnobs*($S_VYnobs - 1))
		scalar `f' = 1 - $S_VYnobs/$S_VYnpop
	}
	else {
		scalar `s' = cond($S_VYosub > 1, /*
		*/ $S_VYnsub*$S_VYnpop/($S_VYnobs*($S_VYosub - 1)), 0)
		
		scalar `f' = 1 - $S_VYosub/$S_VYnsub
	}
		
	if "$S_VYfpc"=="" { 
		matrix $S_VYdeff = `s'*$S_VYdeff
		exit
	}
		
	matrix $S_VYdeft = `s'*$S_VYdeff
	matrix $S_VYdeff = `f'*$S_VYdeft
end

program define MakeVar  /* `D' */
/*
   Create covariance = D*V*D.
*/
	version 4.0
	local D "`1'"  /* note that `D' is destroyed */
	tempname s

/* Generate V = D*V*D. */
	
	matrix $S_VYv  = $S_VYv*`D'
	matrix $S_VYv  = `D'*$S_VYv

	matrix $S_VYdeff = $S_VYdeff*`D'
	matrix $S_VYdeff = `D'*$S_VYdeff
	
	if "$S_VYfpc"!="" {
		matrix $S_VYdeft = $S_VYdeft*`D'
		matrix $S_VYdeft = `D'*$S_VYdeft
	}

/* Make sure V is symmetric. */

	matrix `D' = $S_VYv'
	matrix $S_VYv = $S_VYv + `D'
	scalar `s' = 0.5
	matrix $S_VYv = `s'*$S_VYv
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