*! version 1.0.2  24apr1996
program define svylc
	version 4.0
 	if "`*'"=="" { error 198 }
	if substr("$S_E_cmd",1,3)!="svy" {
		error 301
	}
	capture di matrix(S_E_b[1,1])
	if _rc==0 {
		di in red "must run svy command with " _quote "complete" /*
		*/ _quote " option before using this command"
		exit 301
	}
	parse "`*'", parse(",")
	if "`1'"=="," {
		local options "SHOW"
		capture parse "`*'"
		if _rc { error 198 }
		if "`show'"!="" {
			di
			matrix mlout
			exit
		}
	}
	local formula "`1'"
	local equal = index("`formula'","=")
	if `equal'!=0 {
		di in red _quote "=" _quote " not allowed in expression"
		error 101
	}
	if "`2'"=="," {
		macro shift
		if "$S_E_cmd"=="svylogit" { local or "OR" }
		local options "Level(int $S_level) DEFF DEFT MEFF MEFT `or'"
		parse "`*'"
	}
	else local level = $S_level

	tempname estname x b V Vsrs Vmsp

	qui test `formula' = 0
	scalar `x' = _result(6)

	matrix `b' = get(_b)
	matrix `V' = get(VCE)

/* Save deff, deft, meff matrices and other S_E_. */

	matrix `Vsrs' = S_E_Vsrs
	matrix `Vmsp' = S_E_Vmsp 
	local cmd  "$S_E_cmd"
	local depn "$S_E_depv"
	local fpc  "$S_E_fpc"
	local npsu "$S_E_npsu"
	local nstr "$S_E_nstr"
	if "`fpc'"!="" {
		tempname Vswr
		matrix `Vswr' = S_E_Vswr
	}

	nobreak {
		estimates hold `estname'
		capture noisily break {

			Get_b_V `x' `b' `V' `formula'

			local bb = `b'[1,1]
			local VV = `V'[1,1]

		/* Display formula. */

			noisily test `formula' = 0, notest

		/* Convert Vsrs, Vmsp, Vswr to C*Vsrs*C', etc. */

			Get_CVC `Vsrs' `formula'
			local Deff = `VV'/`Vsrs'[1,1]

			Get_CVC `Vmsp' `formula'
			local Meft = sqrt(`VV'/`Vmsp'[1,1])

			if "`fpc'"!="" {
				Get_CVC `Vswr' `formula'
				local Deft = sqrt(`VV'/`Vswr'[1,1])
			}
			else local Deft = sqrt(`Deff')

		/* Post results for svy_dreg. */

			matrix post `b' `V' /* for matrix b = get(_b), etc. */

			if `Deff'!=. { matrix S_E_deff = (`Deff') }
			else	       matrix S_E_deff = (0)
			if `Deft'!=. { matrix S_E_deft = (`Deft') }
			else	       matrix S_E_deft = (0)
			if `Meft'!=. { matrix S_E_meft = (`Meft') }
			else	       matrix S_E_meft = (0)

			global S_E_npsu "`npsu'"
			global S_E_nstr "`nstr'"
			global S_E_depv "`depn'"
			global S_E_cmd  "`cmd'"
			noisily {
				di
				svy_dreg, nohead level(`level') prob ci /*
				*/ `deff' `deft' `meff' `meft' `or'
			}
		}
		local rc = _rc
		estimates unhold `estname'
		if `rc' { error `rc' }
	}
	global S_1 "`bb'"
	global S_2 = sqrt(`VV')
	global S_3 "`nstr'"
	global S_4 "`npsu'"
	global S_5 "`Deff'"
	global S_6 "`Deft'"
	global S_7 "`Meft'"
end

program define Get_CVC /* matrix_V `formula' */
	version 4.0
	local V "`1'"
	macro shift
	tempname w

	local dim = colsof(`V')
	local eqnames : coleq(`V')
	local cnames  : colnames(`V')
	matrix `w' = J(1,`dim',0)
	matrix coleq    `w' = `eqnames'
	matrix colnames `w' = `cnames'

	matrix post `w' `V'
	qui test `*' = 1
	scalar `w' = _result(6)

	matrix `V' = (0)
	if 1/`w'!=. { matrix `V'[1,1] = 1/`w' }
end

program define Get_b_V  /* test_stat b V formula */
/*
   Test statistic x = (Cb - q)*(CVC')^(-1)*(Cb - q).
*/
	version 4.0
	local x "`1'"  /* test statistic with q = 0 */
	local b "`2'"
	local V "`3'"
	macro shift 3
	tempname U w r

	local dim = colsof(`b')
	local eqnames : coleq(`b')
	local cnames  : colnames(`b')
	matrix `U' = `V'

	matrix `w' = J(1, `dim', 0)
	matrix coleq    `w' = `eqnames'
	matrix colnames `w' = `cnames'

	quietly {
	/* b = 0, true V, q = 1 gives test statistic = (CVC')^(-1). */

		matrix post `w' `V'
		test `*' = 1
		scalar `w' = _result(6)

	/* Var(Cb) = CVC' = 1/w. */

		matrix `V' = (0)
		if 1/`w'!=. { matrix `V'[1,1] = 1/`w' }
		matrix colnames `V' = (1)
		matrix rownames `V' = (1)

	/* With q = -1, get same result unless there is a constant. */

		test `*' = -1
		if abs(`w' - _result(6))/`w' > 1e-4 & `V'[1,1]!=0 {
			di in red "constant terms not allowed"
			exit 101
		}

	/* We eariler computed x = (Cb)^2*(CVC')^(-1).
   	   Thus, r = |Cb| = sqrt(x/w).
	*/
		scalar `r' = sqrt(`x'/`w')

		if `r'==. {
			di in red "point estimate equals " /*
			*/ "missing value"
			exit 504
		}

	/* Determine sign of r = Cb. */

		matrix post `b' `U'  /* same as original post. */
		test `*' = `r'
		scalar `w' = _result(6)
		test `*' = -`r'
		if _result(6) < `w' { scalar `r' = -`r' }

		matrix `b' = (0)
		matrix `b'[1,1] = `r'
		matrix colnames `b' = (1)
	}
end