*! version 1.0.1  15jun1995                     (sg37.1: STB-26)
* Computes orthogonal polynomials using the Christoffel-Darboux
* recurrence formula
program define orthpoly /* variable */
	version 4.0
	local varlist "req ex max(1)"
	local weight "aweight fweight"
	local if "opt"
	local in "opt"
	local options "Degree(integer 1) Generate(string) Poly(string)"
	parse "`*'"
	local x "`varlist'"
	local wght "`weight'"
	local wexp "`exp'"

	if "`generat'"=="" & "`poly'"=="" {
		di in red "generate() or poly() required"
		exit 100
	}
	if `degree' < 1 {
		di in red "degree() must be positive"
		exit 499
	}

	tempvar doit t
	tempname a b c
	mark `doit' [`wght'`wexp'] `if' `in'
	markout `doit' `x' 

	if "`generat'" ~= "" {
		parse "`generat'", parse("*")
		if "`2'"=="*" & "`3'"=="" {
			local generat "`1'1-`1'`degree'"
		}
		local varlist "req new"
		local weight
		local exp
		local if
		local in
		parse "`generat'"
		drop `varlist'
		local nvar : word count `varlist'
		if `nvar' ~= `degree' {
			di in red /*
		*/ "generate() must specify degree = `degree' new variables"
			if `nvar' < `degree' { error 102 }
			else                 { error 103 }
		}
	}

	quietly {

		if "`poly'" ~= "" {
			tempname p p1
			local dd = `degree' + 1
			matrix `p' = J(`dd', `dd', 0)
			matrix `p'[1, 1] = 1
		}

		gen double `t' = .
		local q0 1
		local i 1
		while `i' <= `degree' {
			local j = `i' - 1
			replace `t' = `x'*`q`j''^2 if `doit'
			summarize `t' [`wght'`wexp'] if `doit'
			scalar  `b' = _result(3)
			replace `t' = `x'*`t' if `doit'
			summarize `t' [`wght'`wexp'] if `doit'
			scalar  `a' = _result(3)
			if `i' > 1 {
				local k = `i' - 2
				replace `t' = `x'*`q`j''*`q`k'' if `doit'
				summarize `t' [`wght'`wexp'] if `doit'
				scalar  `c' = _result(3)
			}
			else {
				local k 0
				scalar `c' = 0
			}

			scalar `a' = 1/sqrt(`a' - `b'*`b' - `c'*`c')
			scalar `b' = -`a'*`b'
			scalar `c' =  `a'*`c'

			tempvar q`i'
			gen double `q`i'' = (`a'*`x' + `b')*`q`j'' /*
			*/		    - `c'*`q`k'' if `doit'

			local i = `i' + 1

			if "`poly'" ~= "" {
				local j = `j' + 1
				local k = `k' + 1
				local pnames "`pnames' deg`j'"
				matrix `p'[`i',1] =   `b'*`p'[`j',1] /*
				*/		    - `c'*`p'[`k',1]
				local h 2
				while `h' <= `i' {
					matrix `p'[`i',`h'] =        /*
					*/	  `a'*`p'[`j',`h'-1] /*
					*/	+ `b'*`p'[`j',`h']   /*
					*/	- `c'*`p'[`k',`h']
					local h = `h' + 1
				}
			}
		} /* end of i loop */
	} /* end of quietly */

	if "`generat'" ~= "" {
		parse "`varlist'", parse(" ")
		capture {
			local i 1
			while `i' <= `degree' {
				rename `q`i'' ``i''
				label var ``i'' "deg=`i' orth. poly. for `x'"
				local i = `i' + 1
			}
		}
		if _rc {
			local rc = _rc
			local i 1
			while `i' <= `degree' {
				capture drop ``i''
				local i = `i' + 1
			}
			error `rc'
		}
	}

	if "`poly'" ~= "" {
		matrix `p1' = `p'[1,.]  /* interchange top row and bottom */
		matrix `p' = `p'[2...,.]
		matrix `p' = `p' \ `p1'
		matrix `p1' = `p'[.,1]  /* interchange first column and last */
		matrix `p' = `p'[.,2...]
		matrix `poly' = `p' , `p1'
		matrix colnames `poly' = `pnames' _cons
		matrix rownames `poly' = `pnames' _cons
	}
end