canon.ado

是一个经济学管理应用软件 很难找的 但是经济学学生又必须用到

*! version 1.1.11  06apr2005
program define canon, eclass byable(recall)
	version 9
	local options `"Level(cilevel) STDCoef COEFMatrix Format(string) NOTESTs TESTs(numlist ascending integer >0)"'
	
	if _caller() < 9 {
		canon_8 `*'
		exit
	}

	if !replay() {
		if substr("`1'",1,1)=="(" {
			GetEq `0'
			local eq1 `s(name)'
			local vl1 `s(varlist)'
			local nl1 : word count `vl1'
			
			GetEq `s(rest)'
			local eq2 `s(name)'
			local vl2 `s(varlist)'
			local nl2 : word count `vl2'
			
			local 0 `"`s(rest)'"'
			sret clear
			if "`eq1'" == "" {
				local eq1 "u"
			}
			if "`eq2'"== "" { 
				local eq2 "v"
			}
		}
		else {
		/* this seems to be antiquated syntax for specifying equations
		   via the eq command, from Stata 5 */
			local iseq True
			gettoken eq1 0 : 0, parse(" ,")
			gettoken eq2 0 : 0, parse(" ,")
		}

		syntax [aw fw] [if] [in] [, `options' noConstant	///
			LC(numlist int min=1 max=1) 			///
			FIrst(numlist int min=1	max=1) ]
	
		local constan `constant'

		if "`constan'"=="" { 
			local dev "dev"
		}
		else 	local nocons "nocons" 

		if `"`stdcoef'"'!="" & `"`coefmatrix'"'!="" {
			di as err "cannot specify stdcoef and coefmatrix " ///
				"options together"
			error(198)
		}
		if `"`stdcoef'"'=="" & `"`coefmatrix'"'=="" & `"`format'"' != "" {
			di as err "format() may only be specified with stdcoef or coefmatrix"
			exit(198)
		}
		if "`format'" =="" {
			local format format(%8.4f)
		} 
		else {
			quietly di `format' 0
			local format format(`format')
		}
	
		marksample touse

		if "`iseq'"=="True" {
		/* this seems to be antiquated syntax for specifying equations
		   via the eq command, from Stata 5 */
			eq ? `eq1'
			local vl1 `"`r(eq)'"'
			local eq1 `"`r(eqname)'"'
			local nl1 : word count `vl1'
			eq ? `eq2'
			local vl2 `"`r(eq)'"'
			local eq2 `"`r(eqname)'"'
			local nl2 : word count `vl2'
		}

		markout `touse' `vl1' `vl2'

		if "`lc'"!="" & "`first'"!="" {
			di as err "lc and first may not be specified together"
			exit 198
		}
		if "`notests'"!="" & "`tests'"!="" {
	di as err "test and notests may not be specified together"
			exit 198
		}
		local ncc = min(`nl1', `nl2')
		if "`lc'"!="" {
			if `lc' > `ncc' | `lc' < 1 {
				di as err "lc() must be one "	///
				"of the `ncc' "	///
				"canonical correlations"
				exit 198
			}
		}
		if "`first'" !="" {
			if `first' > `ncc' {
				di as err "first() cannot be greater "	///
				"than `ncc', the number of "	///
				"canonical correlations"
				exit 198
			}
			if `first' < 1 {
				di as err "first() must be greater " ///
				"than 0"
				exit 198
			}
		}
		if "`tests'"!="" {
			local ltest : word count `tests'
			local ltest : word `ltest' of `tests'
			if (`ltest' > `ncc') {
				if `ncc' == 1 {
		di as err "test() may only test the one canonical correlation"
					exit 125
				}
				else {
	di as err "test() may only test the `ncc' canonical correlations"
					exit 125
				}
			}
		}


		tempname bigABC ABC A B C Ai Ci T TR TC TCR L LL SD V SD1 SD2
		tempname sse t r2 tt ccor foo D1 D2
		tempname corr_ind corr_dep corr_mixed
		tempvar xx yy

		/* `ABC' is almost the covariance matrix for the two
		    variable sets, all it needs is to divide by nobs-1 */
		qui capture mat accum `ABC' = `vl2' `vl1' if `touse' /*
			*/ [`weight'`exp'] , nocons `dev'
		if _rc == 2000 {
			di as error "insufficient observations"
			exit 2000
		}
		else{
			if _rc {
			 	return _rc
			}
		}
		capture matrix `foo' = cholesky(`ABC')
		if _rc == 506 {
			di as error "matrix not positive definite"
			di as error "failure for overall model"
			exit 506
		} 
		local nobs = r(N)
		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 p `v1e'
		local q = `v2e' - `v1e' 
		local v2 = `v1e' + 1
		mat `SD' = vecdiag(`ABC')
		mat `ABC' = corr(`ABC')  /* correlation transform */
		mat `A' = `ABC'[1..`v1e',1..`v1e'] /* first quadrant 
				correlations of first variable set */
		matrix `corr_dep' = `A'
		mat `B' = `ABC'[1..`v1e',`v2'..`v2e'] /* off diagonal 
				correlations between first and second 
				variable set */
		matrix `corr_mixed' = `B'
		mat `C' = `ABC'[`v2'..`v2e',`v2'..`v2e'] /* last quadarant 
				correlations of last variable set */
		matrix `corr_ind' = `C'
		mat `Ai' = syminv(`A')
		mat `Ci' = syminv(`C')
		local n2e = `v2e' - `v1e'
		if (`n2e' > `v1e') { /* first varlist is shorter */
			local r1list "vl1"
			local r2list "vl2"
		 	* looking for left symeigenvectors of B' Ai B Ci
			mat `SD1' = `SD'[1,`v2'..`v2e'] /* variances vars1*/
			mat `SD2' = `SD'[1,1..`v1e']  /* variances vars2 */
			mat `D1' = `SD1'
			mat `D2' = `SD2'
			mat `T' = (`B'' * `Ai') * `B'
			mat `L' = cholesky(`Ci')
			mat `T' = (`L'' * `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 { /* second varlist is shorter */
			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 `D1' = `SD1'
			mat `D2' = `SD2'
			mat `T' = (`B' * `Ci') * `B''
			mat `L' = cholesky(`Ai')
			mat `T' = (`L'' * `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'
			local flip True
		}
		mat colnames `SD1' = `eq1':
		mat `SD1' = syminv(cholesky(diag(`SD1')))
		mat `V' = `L' * `V'       /* Undo the similarity transform */

		mat `L' = (`V''*`XY')*`XXi' 	/* the other lin. comb. */
		mat colnames `SD2' = `eq2':
		mat `SD2' = syminv(cholesky(diag(`SD2')))
		mat `V' = `V'' * `SD1'
		local lcs = `n2s' - 1
		if "`lc'"!="" {
			local beglc `lc'
			local endlc `lc'
		}
		else {
			local beglc 1
			local endlc `lcs'
			if "`first'"!="" {
				local endlc `first'
			}
		}
		mat `L' = `L' * `SD2'
		tempname zerom zeromT
		forvalues kc = `beglc'/`endlc' {
			tempname V`kc' L`kc' T`kc' TC`kc' XXi2 YYi2
			mat `V`kc'' = `V'[`kc',1...]  
			mat `L`kc'' = `L'[`kc',1...]  
		/* The linear combination we want above */	
			scalar `r2' = `LL'[1,`kc']
			scalar `t' = 1 / sqrt(`r2')
			mat `L`kc'' = `L`kc'' * `t'
			mat `T`kc'' = `L`kc'', `V`kc''
			local nv : colnames(`T`kc'')
			local ne : coleq(`T`kc'')
			local ne2
			foreach equa of local ne {
				local ne2  `ne2' `equa'`kc'
			}
		/* conditional variance calculation */
			mat `XXi2' = `SD2' * (`XXi' * `SD2')
			scalar `tt' = (1 - `r2')/`r2'/(1/`sse'-colsof(`XXi'))
			mat `XXi2' = `XXi2' * `tt'
			mat `YYi2' = `SD1' * (`YYi' * `SD1')
			scalar `tt' = (1 - `r2')/`r2'/(1/`sse'-colsof(`YYi'))
			mat `YYi2' = `YYi2' * `tt'
			mat rownames `ABC' = `nv'
			mat roweq `ABC' = `ne2'
			mat colnames `ABC' = `nv'
			mat coleq `ABC' = `ne2'
			mat coleq `T`kc'' = `ne2'
			mat `ABC' = `ABC' * 0
			mat subst `ABC'[1,1] = `XXi2'
			mat subst `ABC'[`n2s',`n2s'] = `YYi2'

						/* post results */
			matrix `ccor' = vecdiag(cholesky(diag(`LL')))
			tempname tp1 tp2
			forvalues eqn = 1/2 {
				local mqn = `eqn'
				if "`flip'"!="" {
					local mqn = mod(`eqn',2) + 1
				}
				mat `tp2' = `T`kc''[1, "`eq`mqn''`kc':"]
				matrix `tp1' = cholesky(diag(`D`mqn''))*`tp2''
				local tp1r = rowsof(`tp1')
				local tp1c = colsof(`tp1')
				local pos
				forvalues ii = 1/`tp1r' {
					forvalues jj = 1/`tp1c' {
						if `tp1'[`ii', `jj'] > 0 {
							local pos pos
						}
					}
				}
				if "`pos'" =="" {
					matrix `tp1' = -1*`tp1'
				}
				local mynms : colfullnames `tp2'
				mat roweq `tp1' = `mynms'
				if `eqn' == 1 {
					mat `TC`kc'' = `tp1''
				}
				else {
					mat `TC`kc'' = (`TC`kc'', `tp1'')
				}
			}
			if "`flip'"=="" {		/* sic */
				mat `T`kc'' = (`T`kc''[1,"`eq1'`kc':"], `T`kc''[1,"`eq2'`kc':"])
				tempname Z
				mat `Z' = `ABC'["`eq1'`kc':", "`eq2'`kc':"] * 0
				mat `ABC' = ( /*
				*/ (`ABC'["`eq1'`kc':","`eq1'`kc':"], `Z'  ) \ /*
				*/ (`Z'', `ABC'["`eq2'`kc':","`eq2'`kc':"] ))
			}
			if `kc'==`beglc' {
				mat `bigABC' = `ABC'
				mat `T' = `T`kc''
				mat `TR' = (`T`kc'')'
				mat `TC' = `TC`kc''
				mat `TCR' = (`TC`kc'')'
				mat `zerom' = `ABC'*0
				mat `zeromT' = `zerom'
			}
			else {
				local nenew : coleq(`ABC')
				mat coleq `zerom' = `nenew'
				mat roweq `zeromT' = `nenew'
			mat `bigABC' = `bigABC', `zerom' \ `zeromT', `ABC'
				mat `T' = `T' ,  `T`kc''
				mat `TR' = `TR', (`T`kc'')'
				mat `TC' = `TC' ,  `TC`kc''
				mat `TCR' = `TCR', (`TC`kc'')'
				mat `zerom' = `zerom' \ `ABC'*0
				mat `zeromT' = `zeromT', `ABC'*0
			}
		}
		/* get lables right, no equations here, just variables */
		mat roweq `TCR' = "_:"
		local myrownms : rownames `TCR'
		matrix roweq `TR' = "_:"
		matrix rownames `TR' = `myrownms'
		
		/* Now lets compute the canonical loadings 
		   We have A = cov matrix for first set of vars
		   We have C = cov matrix for second set of vars
		   We have SD1 and SD2 which are diagonal matrices
		   of one over standard deviations of the two
		   equations sets.  SD1 goes with eq1 and SD2 with
		   eq2, but this may have been reversed due to 
		   the flip ... beware.
		   We have the coefficients in TR -- we'll
		   have to subdivide to get the individual coeffs.

		   I'm using _Applied Multivariate Statistical Analysis_
		   3rd Ed. by R. A. Johnson and D. W. Wichern, 
		   Prentice Hall, 1992 as a reference here.
		   Section 10.3, pp 466-467 in my edition
		*/
		local colA : colnames `A' /* a is correlation matrix (rho
		or R) for a variable set, C is for the other set,
		and B is the cross correlations.  Since there can be a flip
		depending on which variable set is shortest, all
		I'm doing here is figuring out what goes where.  */
		if "`colA'" == "``r1list''" {
			local SDA `SD1'
			local SDC `SD2'
			local rAlist ``r1list''
			local rClist ``r2list''
		}
		else{
			local SDA  `SD2'
			local SDC  `SD1' 
			local rAlist ``r2list''
			local rClist ``r1list''
		}
		/* now extract submatricies from TR.  TR is the matrix
		of coefficients from the original variables to make
		the canonical variables.  coefA, coefC are what
		Johnson and Wichern are calling A and B in U=AX(1)
		and V=BX(2).*/
		local wc : word count `rAlist'
		local w1 : word 1 of `rAlist'
		local last : word `wc' of `rAlist'
		tempname coefA coefC stdcoefA stdcoefC
		matrix `coefA' = `TR'["`w1'".."`last'", 1...]
		matrix `stdcoefA' = `TCR'["`w1'".."`last'", 1...]
		local wc : word count `rClist'
		local w1 : word 1 of `rClist'
		local last : word `wc' of `rClist'
		matrix `coefC' = `TR'["`w1'".."`last'", 1...]
		matrix `stdcoefC' = `TCR'["`w1'".."`last'",1...]
		tempname canloadA canloadC SDAi SDCi canloadM canloadM2 
		/* SDAi is what Johnson and Wichern call V^(1/2) 
		   SDA is V^(-1/2).  We need SDAi to convert
		   our correlation matrix to covariance