*! version 1.0.7  04aug1998
program define canon
	version 4.0
	parse "`*'", parse(" ,")
	local options "Level(integer $S_level)"
	if ("`1'"!="" & substr("`1'",1,1)!=",") {
		local eq1 `1'
		local eq2 `2'
		mac shift 2

		local options "`options' noConstant LC(int 1)"
		local weight  "aweight fweight"
		local if "opt"
		local in "opt"

		parse "`*'"

		tempname ABC A B C Ai Ci T L LL SD V SD1 SD2
		tempname sse t r2 tt ccor
		tempvar xx yy
		if ("`constan'"=="") { local dev "dev" }
		else { local nocons "nocons" }

		eq ? `eq1'
		local vl1 "$S_1"
		local eq1 "$S_3"    /* expand abbreviations */
		eq ? `eq2'
		local vl2 "$S_1"
		local eq2 "$S_3"

		qui mat accum `ABC' = `vl2' `vl1' `if' `in'  /*
			*/ [`weight'`exp'] , nocons `dev'
		local nobs = _result(1)
		local nobs1 = `nobs' - 1
		scalar `sse' = 1/(`nobs1')
		mat def `ABC' = `ABC' * `sse' /* convert to covariance */
		local v2e = rowsof(`ABC')
		local v1e : word count `vl2'
		local v2 = `v1e' + 1
		mat `SD' = vecdiag(`ABC')
		mat `ABC' = corr(`ABC')
		mat `A' = `ABC'[1..`v1e',1..`v1e']
		mat `B' = `ABC'[1..`v1e',`v2'..`v2e']
		mat `C' = `ABC'[`v2'..`v2e',`v2'..`v2e']
		mat `Ai' = syminv(`A')
		mat `Ci' = syminv(`C')
		local n2e = `v2e' - `v1e'
		if (`n2e' > `v1e') {
			local r1list "vl1"
			local r2list "vl2"
			* looking for left symeigenvectors of B' Ai B Ci
			mat `SD1' = `SD'[1,`v2'..`v2e']
			mat `SD2' = `SD'[1,1..`v1e']
			mat `T' = `B'' * `Ai'
			mat `T' = `T' * `B'
			mat `L' = cholesky(`Ci')
			mat `T' = `L'' * `T'
			mat `T' = `T' * `L'
			mat symeigen `V' `LL' = `T'
			mat `LL' = `LL'[1,1..`v1e']
			local n2s = `v1e' + 1
			local XXi `Ai'
			local YYi `Ci'
			local XY `B''
		}
		else {
			local r1list "vl2"
			local r2list "vl1"
			local x `eq1'
			local eq1 `eq2'
			local eq2 `x'
			* looking for left symeigenvectors of B Ci B' Ai
			mat `SD2' = `SD'[1,`v2'..`v2e']
			mat `SD1' = `SD'[1,1..`v1e']
			mat `T' = `B' * `Ci'
			mat `T' = `T' * `B''
			mat `L' = cholesky(`Ai')
			mat `T' = `L'' * `T'
			mat `T' = `T' * `L'
			mat symeigen `V' `LL' = `T'
			mat `LL' = `LL'[1,1..`n2e']
			local n2s = `n2e' + 1
			local XXi `Ci'
			local YYi `Ai'
			local XY `B'
		}
		if colsof(`SD1')==1 {
			local name : colname(`SD1')
			mat rownames `SD1' = `eq1':`name'
			mat colnames `SD1' = `eq1':`name'
		}
		else {
			mat colnames `SD1' = `eq1':
			mat `SD1' = diag(`SD1')
		}
		mat `SD1' = cholesky(`SD1')
		mat `SD1' = syminv(`SD1')
		mat `V' = `L' * `V'       /* Undo the similarity transform */

		mat `L' = `V''*`XY'
		mat `L' = `L'*`XXi'       /* The other linear combination */
		if colsof(`SD2')==1 {
			local name : colname(`SD2')
			mat rownames `SD2' = `eq2':`name'
			mat colnames `SD2' = `eq2':`name'
		}
		else {
			mat colnames `SD2' = `eq2':
			mat `SD2' = diag(`SD2')
		}
 		mat `SD2' = cholesky(`SD2')
		mat `SD2' = syminv(`SD2')
		mat `V' = `V'' * `SD1'
		mat `V' = `V'[`lc',.]     /* The linear combination we want */
		mat `L' = `L' * `SD2'
		mat `L' = `L'[`lc',.]     /* The linear combination we want */
		scalar `r2' = `LL'[1,`lc']
		scalar `t' = 1 / sqrt(`r2')
		mat `L' = `L' * `t'
		mat `T' = `L',`V'
		local nv : colnames(`T')
		local ne : coleq(`T')

		/* conditional variance calculation */
		mat `XXi' = `XXi' * `SD2'
		mat `XXi' = `SD2' * `XXi'
		scalar `tt' = (1 - `r2')/`r2'/(1/`sse'-colsof(`XXi'))
		mat `XXi' = `XXi' * `tt'
		mat `YYi' = `YYi' * `SD1'
		mat `YYi' = `SD1' * `YYi'
		scalar `tt' = (1 - `r2')/`r2'/(1/`sse'-colsof(`YYi'))
		mat `YYi' = `YYi' * `tt'
		mat rownames `ABC' = `nv'
		mat roweq `ABC' = `ne'
		mat colnames `ABC' = `nv'
		mat coleq `ABC' = `ne'
		mat `ABC' = `ABC' * 0
		mat subst `ABC'[1,1] = `XXi'
		mat subst `ABC'[`n2s',`n2s'] = `YYi'

						/* post results */
		matrix `ccor' = diag(`LL')
		matrix `ccor' = cholesky(`ccor')
		matrix `ccor' = vecdiag(`ccor')
		mat post `T' `ABC', dof(`nobs1')
		mat S_E_ccor = `ccor'
		mat drop `ccor'
		global S_E_nobs `nobs'
		global S_E_tdf = `nobs1'
		global S_E_lc "`lc'"
		global S_E_cmd "canon"
	}
	else {
		if ("$S_E_cmd"!="canon") { error 301 }
		parse "`*'"
		di
	}
	if (`level'<10 | `level'>99) { local level 95 }
	di _n in gr "Linear combinations for canonical correlation " $S_E_lc /*
	*/ _col(56) "Number of obs =" in ye %8.0f $S_E_nobs
	mat mlout , level(`level')
	di in gr _col(42) "(Std. Errors estimated conditionally)"
	di in gr "Canonical correlations:"
	mat list S_E_ccor, nohead nonames noblank format(%9.4f)
end
exit