pcamat.ado

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

*! version 2.0.0  29mar2005
program pcamat
	version 8

	if replay() {
		if "`e(cmd)'" != "pca" {
			error 301
		}
		pca_display `0'
	}
	else {
		Estimate `0'
	}
end


program Estimate, eclass

	#del ;
	syntax  anything(name=Cin id="covariance or correlation matrix"),
 	        n(numlist max=1 >=0 integer)
	[
	     // display options
	     	// MEans		not allowed with pcamat
	     	noDISPLAY
	        Level(passthru)
	        BLanks(passthru)
	        
	     // matrix data specification
 	        SDS(str)
 	        MEANS2(str)
 	        NAMes(passthru)
 	        SHape(passthru)
 	        FORCE	// undocumented (do not enforce same name stripes)
 	        
	     // model
 	        COMponents(numlist integer max=1 >0)
 	        FActors(numlist integer max=1 >0)  // synonym of components()
 	        MINEigen(numlist max=1 >=0)
 	        CORrelation
 	        COVariance
 	        
	     // SE/VCE
	     	noVCE
	     	VCE2(string)
 	        IGNORE
 	        tol(numlist max=1 >=0)
 	        
 	     // not to be documented -- used internally to speed up pca
 	        matrixtype(str)
	] ;
	#del cr
	local display_options  `blanks' `level' `means' `std' `vce'

	if `:list sizeof Cin' > 1 {
		dis as err "exactly one matrix name expected"
		exit 103
	}
	confirm matrix `Cin'

	ParseVCE `vce2' 
	if "`s(vce)'" == "normal" { 
		local normal normal
		local vce_arg vce(normal)
	}
	
	if "`correlation'" != "" & "`covariance'" != "" {
		dis as err "options correlation and covariance are exclusive"
		exit 198
	}	
	if "`covariance'" == "" { 
		local correlation correlation
	}	
	if "`normal'" != "" & "`correlation'" != "" {
		dis as txt /// 
		   "(with PCA/correlation, SEs and tests are approximate)"
	}

	if "`normal'" == "" {
		if ("`ignore'" != "") { 
			Invalid ignore
		}	
		if ("`level'"  != "") {
			Invalid level()
		}	
	}

// get the matrix for PCA

	if "`matrixtype'" == "" {
		tempname C
		if `"`means2'"' != "" { 
			local mopt means(`means2')
		}
		if "`sds'" != "" {
			local sopt sds(`sds')
		}	
		
		_getcovcorr `Cin', `covariance' `correlation' ///
		   `sopt' `mopt' `names' `shape' `force'
		   
		matrix `C' = r(C)
		local type `r(Ctype)'
		if "`r(means)'" == "matrix" {
			tempname Means
			matrix `Means' = r(means)
		}
		if "`r(sds)'" == "matrix" {
			tempname Sds
			matrix `Sds' = r(sds)
		}
	}
	else {
		local C `Cin'
		local type `matrixtype'
		local Means `means2' 
		local Sds   `sds' 
	}
	
	local nvar = rowsof(`C')
	local varlist : rownames `C'

//  number of components (factors) to be extracted

	if `"`components'"' != "" & `"`factors'"' != "" {
		dis as err "components() and factors() are synonyms"
		exit 198
	}
	if `"`factors'"' != "" {
		local optname factors()
		local nf = `factors'
	}
	if "`components'" != "" {
		local optname components()
		local nf = `components'
	}
	if "`nf'" == "" {
		local nf = `nvar'
	}
	else if !inrange(`nf',1,`nvar') {
		dis as err "`optname' should be between 1 and `nvar'"
		exit 125
	}
	if "`tol'" == "" {
		local tol = 1e-5
	}
	if "`mineigen'" == "" {
		local mineigen = `tol'
	}

// ---------------------------------------------------------------------------
// do actual PCA of C -- a spectral decomposition
// notes on symeigen
//   returns (eigenvectors,eigenvalues), sorted on eigvals (decreasing order)
//   returns eigenvectors e signed positively, i.e., with e'1 > 0
// ---------------------------------------------------------------------------

	tempname b bV Ev L nobs condW lndetW Rho traceW Psi

	scalar `nobs' = `n'
	quietly matrix symeigen `L' `Ev' = `C'

	CheckEigenvalues `Ev' `tol' "`matrix'" "`type'" "`ignore'" "`normal'"
	local  normal   `r(normal)'
	scalar `traceW' = r(trace)
	scalar `lndetW' = r(lndet)
	scalar `condW'  = r(cond)

	forvalues i = 1 / `nvar' {
		local facnames `facnames' Comp`i'
	}
	matrix colnames `Ev' = `facnames'   // eigenvalues
	matrix coleq    `Ev' = Eigenvalues
	matrix colnames `L'  = `facnames'   // eigenvectors

// retain factors

	forvalues i = `nf'(-1)1 {
		if `Ev'[1,`i'] < `mineigen' {
			local --nf
		}
	}
	
	if `nf' == 0 {
		dis as txt "(no components retained)"
		
		ereturn post , obs(`=`nobs'') properties(nob noV eigen)
		
		ereturn scalar f   = 0
		ereturn scalar rho = 0 
		ereturn matrix Ev  = `Ev'
	
		ereturn matrix C = `C'            // covar/corr matrix
		ereturn local Ctype `type'
		ereturn scalar trace = `traceW'   // trace of C
		if "`Means'" != "" {
			ereturn matrix means = `Means'
		}
		if "`Sds'" != "" {
			ereturn matrix sds   = `Sds'
		}
		
		ereturn local title      "Principal components/`type'"
		ereturn local predict    pca_p
		ereturn local rotate_cmd pca_rotate
		ereturn local estat_cmd  pca_estat
		
		ereturn local cmd  pca

		if "`display'" == "" {
			pca_display , `display_options'
		}
		exit
	}

	scalar `Rho' = 0
	forvalues i = 1 / `nf' {
		scalar `Rho' = `Rho' + `Ev'[1,`i']
	}
	scalar `Rho' = `Rho' / `traceW'

// ---------------------------------------------------------------------------
// VCE assuming multivariate normality
// ---------------------------------------------------------------------------

	if "`normal'" != "" {

		tempname b0 st1 st2 biasEv varEv ELC
		tempname chi2_i df_i p_i chi2_s df_s p_s 
		tempname e ee ei ej lsum llsum t Vi

		forvalues i = 1 / `nf' {
			matrix `b0' = `L'[1...,`i']'
			matrix coleq `b0' = Comp`i'
			matrix `b' = nullmat(`b'), `b0'
		}
		matrix drop `b0'

	// test for independence (Basilevsky: 187)

		if "`covariance'" != "" {
			tempname cC
			matrix `cC' = corr(`C')
		}
		else {
			local cC `C'
		}

		scalar `chi2_i' = -(`nobs'-(2*`nvar'+5)/6) * ln(det(`cC'))
		scalar `df_i'   = `nvar'*(`nvar'-1)/2
		scalar `p_i'    = chi2tail(`df_i', `chi2_i')

	// test for sphericity (Basilevsky: 192)

		scalar `chi2_s' = ///
		        -(`nobs' - (2*`nvar'^2+`nvar'+2)/(6*`nvar')) * ///
		         (`lndetW' - `nvar'*ln(`traceW'/`nvar'))
		scalar `df_s'   = (`nvar'+2)*(`nvar'-1)/2
		scalar `p_s'    = chi2tail(`df_i', `chi2_s')

	// compute bias(Ev) and var(Ev) of eigenvalues Ev
	// (asymptotic values based on normality)

		matrix `biasEv' = J(1,     `nvar',0)    // bias(Ev)
		matrix `varEv'  = J(`nvar',`nvar',0)    // var(Ev)
		matrix colnames `biasEv' = `facnames'
		matrix rownames `biasEv' = bias

		// estimate V(lambda) with terms of order 1/n and 1/n^2
		//
		// if this estimator that are not posdef, we switch to the
		// estimator that uses terms up to order 1/n only
		
		local EV_order2 = 1
		forvalues i = 1/`nvar' {
		    scalar `st1' = 0
		    scalar `st2' = 0
		    forvalues j = 1 / `nvar' {
		        if (`j' == `i')  continue

		        scalar `t' = `Ev'[1,`j'] / (`Ev'[1,`i']-`Ev'[1,`j'])
		        scalar `st1' = `st1' +  `t'
		        scalar `st2' = `st2' + (`t')^2

		        // cov(L_i,L_j)
		        matrix `varEv'[`i',`j'] = /// 
		           (2/`nobs'^2) * (`Ev'[1,`i']*`t')^2

		    }
		    matrix `biasEv'[1,`i']  = (`Ev'[1,`i']/`nobs') * `st1'
		    matrix `varEv'[`i',`i'] = ///
		       (2/`nobs')*`Ev'[1,`i']^2 * (1-`st2'/`nobs')
		      
		    if (`varEv'[`i',`i'] <= 0) local EV_order2 = 0
		}
		
		if (`EV_order2' == 0) { 	
		    // estimate var(Ev) with order 1/n term only
		    matrix `varEv'  = J(`nvar',`nvar',0)
		    forvalues i = 1/`nvar' {
		        matrix `varEv'[`i',`i'] =  (2/`nobs')*`Ev'[1,`i']^2
	            }
		}

	// Compute variance bV of eigvecs, asymptotic based on normality

		local nr = colsof(`b')
		local vfnames : colfullnames(`b')
		matrix `bV' = J(`nr',`nr',0)