pcamat.ado

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

		matrix rownames `bV' = `vfnames'
		matrix colnames `bV' = `vfnames'

		// diagonal block associated with variance matrices 
		// of eigenvectors (components)

		local ic = 1
		forvalues i = 1 / `nf' {
		    matrix `Vi' = J(`nvar',`nvar',0)
		    forvalues j = 1 / `nvar' {
		        if (`i'==`j')  continue
		      
		        scalar `t'  = (`Ev'[1,`i']*`Ev'[1,`j']) / ///
		                      (`Ev'[1,`i']-`Ev'[1,`j'])^2 
		        matrix `e'  = `L'[1...,`j']
		      
		        matrix `Vi' = `Vi' + `t' * `e'*`e'' 
		    }
		    matrix `bV'[`ic',`ic'] = `Vi'
		    local ic = `ic' + `nvar'
		}
		matrix drop `Vi' `e'

		// off-diagonal blocks associated with covariances of 
		// eigenvectors (components)

		local ic = 1
		forvalues i = 1 / `nf' {
		    matrix `ei' = `L'[1...,`i']
		    local jc = 1
		    forvalues j = 1 / `=`i'-1' {
		        scalar `t'  = (`Ev'[1,`i']*`Ev'[1,`j']) / ///
		                    (`Ev'[1,`i']-`Ev'[1,`j'])^2
		        matrix `ej' = `L'[1...,`j']
		        matrix `ee' = - `t' * `ej' * `ei''

		        matrix `bV'[`ic',`jc'] = `ee'
		        matrix `bV'[`jc',`ic'] = `ee''
		        local jc = `jc' + `nvar'
		    }
		    local ic = `ic' + `nvar'
		}

		// scale and attch name

		matrix `bV' = (1/`nobs') * `bV'

		capture matrix drop `ei'
		capture matrix drop `ej'
		capture matrix drop `ee'

	// additional statistics

		VarExplainedSE `nobs' `Ev' `ELC'

	// create (b,V) by concatenate eigenvalues and eigenvectors and 
	// their V-parts

		matrix `b' = `Ev'[1,1..`nf'] , `b'
		local names : colfullnames `b'

		matrix `bV' = ( `varEv'[1..`nf',1..`nf'] , J(`nf',`nr',0) \ ///
		                J(`nr',`nf',0)           , `bV'           )
		matrix colnames `bV' = `names'
		matrix rownames `bV' = `names'
	}

// ----------------------------------------------------------------------------
// post results
// ----------------------------------------------------------------------------
	
	matrix `L' = `L'[1...,1..`nf']
	matrix `Psi' = vecdiag( `C' - `L'*diag(`Ev'[1,1..`nf'])*`L'' )
	forvalues j = 1/`=colsof(`Psi')' {
		matrix `Psi'[1,`j'] = chop(`Psi'[1,`j'],1e-6)
	}
	matrix rownames `Psi' = Unexplained

	ereturn clear
	if "`normal'" != "" {
		ereturn post `b' `bV' , obs(`=`nobs'') properties(b V eigen)
		ereturn local  vce    "multivariate normality"

		ereturn scalar v_rho  = `ELC'[`nf',5]  // asymptotic var(rho)

		ereturn scalar chi2_i = `chi2_i'       // test for independence
		ereturn scalar df_i   = `df_i'
		ereturn scalar p_i    = `p_i'

		ereturn scalar chi2_s = `chi2_s'       // test for sphericity
		ereturn scalar df_s   = `df_s'
		ereturn scalar p_s    = `p_s'

		ereturn matrix Ev_bias  = `biasEv'     // bias of eigenvalues
		ereturn matrix Ev_stats = `ELC'        // stats on expl var
	}
	else {
		ereturn post , obs(`=`nobs'') properties(nob noV eigen)
	}

	ereturn matrix L  = `L'                        // eigvecs = components
	ereturn matrix Ev = `Ev'                       // eigenvalues
	ereturn matrix Psi  = `Psi'                    // unexplained corr/cov

	ereturn matrix C = `C'                         // covar/corr matrix
	ereturn local Ctype `type'

	if "`Means'" != "" {
		ereturn matrix means = `Means'
	}
	if "`Sds'" != "" {
		ereturn matrix sds   = `Sds'
	}

	ereturn scalar f     = `nf'                    // #comps=factors ret.
	ereturn scalar rho   = `Rho'                   // %explained variance
	
	ereturn scalar trace = `traceW'                // trace of C
	ereturn scalar lndet = `lndetW'                // ln(det) of C
	ereturn scalar cond  = `condW'                 // condition number(C)

	ereturn local predict    pca_p
	ereturn local rotate_cmd pca_rotate
	ereturn local estat_cmd  pca_estat

	ereturn local title   "Principal components"
	ereturn local cmd     pca                      // sic !

// and now display

	if "`display'" == "" {
		pca_display , `display_options'
	}
	_returnclear
end


program CheckEigenvalues, rclass
	args Ev tol matname mattype ignore normal

	tempname cond trace lndet
	local n = colsof(`Ev')

	local nneg = 0
	forvalues i = 1 / `n' {
		if `Ev'[1,`i'] <= 0 {
			local ++nneg
		}
	}
	
	if `nneg' == `n' {
		dis as err "{p}all eigenvalues of `mattype' `matname'" ///
			"are 0 or negative; this cannot be ignored{p_end}"
		exit 498
	}	

	if "`matname'" != "" & `nneg' > 0 {
		local color = cond("`ignore'"!="","txt","err") 
		if "`normal'" != "" {
			dis as `color' "standard errors assuming " ///
			    "multivariate normality cannot be computed"
		}	
		dis as `color' ///
		    "`mattype' `matname' has `nneg' nonpositive " ///
		    "eigenvals; it is not positive definite"
		matlist `Ev', noheader rowtitle(EigVals) left(4) twidth(8)
		
		if "`ignore'" != "" { 
			local normal 
		}
		else if "`ignore'" == "" {
			dis as err "specify option ignore to replace " /// 
			           "negative eigenvalues by 0"
			exit 198
		}
	}

// varlist : negative eigenvalues must be due to rouding errors
//  matrix : option ignore specifies chopping negative eigenvalues

	forvalues i = `=`n'-`nneg'+1' / `n' {
		matrix `Ev'[1,`i'] = 0
	}

// compute some statistics

	if `nneg' == 0 {
		scalar `cond'  = sqrt(`Ev'[1,1]/`Ev'[1,`n'])
		scalar `lndet' = 0
		forvalues i = 1 / `n' {
			scalar `lndet' = `lndet' + ln(`Ev'[1,`i'])
		}
	}
	else {
		scalar `cond'  = .
		scalar `lndet' = .  
	}
	scalar `trace' = 0 
	forvalues i = 1 / `n' {
		scalar `trace' = `trace' + `Ev'[1,`i']
	}	
	scalar `trace' = chop(`trace',abs(`trace')*1e-6)

// check for (almost) zero or equal eigenvalues

	if "`normal'" != "" {
		// this is only reached if all eigvals>0
		local stop = 0
		if `Ev'[1,`n'] < `tol' * `Ev'[1,1] {
			local msg "matrix is nearly singular"
			local stop = 1
		}
		else {
			forvalues i = 2 / `n' {
				if `Ev'[1,`i'-1]-`Ev'[1,`i']<`tol'*`Ev'[1,1] {
					local stop = 1
				}
			}
			if `stop' {
				local msg /// 
				      "matrix has (almost) equal eigenvalues"
			}
		}
		
		if `stop' {
			dis _n as txt "`msg'" 
			dis as txt `"tolerance for "closeness": "' as res `tol'
			
			matlist `Ev', rowtitle(EigVals) left(4) twidth(8)
			if "`ignore'" == "" { 
				error 498
			}	
		}
	}

	return local  normal  `normal'
	return scalar cond  = `cond'    // condition number of C
	return scalar trace = `trace'   // trace(C)
	return scalar lndet = `lndet'   // ln(det(C))
end


// inference on percentage and cumulative percentage of variance
// explained by the "largest eigenvectors". See Kshirsagar (1972:454)
//
program VarExplainedSE
	args nobs   /// input:  number of obs
	     Ev	    /// input:  sorted eigenvalues, 1xnv ROWVECTOR
	     ELC     // output: nv x 5 matrix

	tempname EvC rho sEv sEv2

	local nvar = colsof(`Ev')
	matrix `EvC'  = J(`nvar',2,0)
	matrix `EvC'[1,1] = `Ev'[1,1]
	matrix `EvC'[1,2] = `Ev'[1,1]^2
	forvalues i = 2 / `nvar' {
		matrix `EvC'[`i',1] = `EvC'[`i'-1,1] + `Ev'[1,`i']
		matrix `EvC'[`i',2] = `EvC'[`i'-1,2] + `Ev'[1,`i']^2
	}
	scalar `sEv'  = `EvC'[`nvar',1]
	scalar `sEv2' = `EvC'[`nvar',2]

	matrix `ELC' = J(`nvar',5,.)
	forvalues i = 1 / `nvar' {
		matrix `ELC'[`i',1] = `Ev'[1,`i']
		matrix `ELC'[`i',2] = `Ev'[1,`i'] / `sEv'
		matrix `ELC'[`i',3] = ///
		   (2*`Ev'[1,`i']^2 / (`nobs'*`sEv'^2)) * ///
		   (1 - (2*`Ev'[1,`i']/`sEv') + (`sEv2'/`sEv'^2))

		// rho = cumulative percentage of variance explained
		scalar `rho' = `EvC'[`i',1] / `sEv'
		matrix `ELC'[`i',4] = `rho'

		// asymptotic variance of rho
		matrix `ELC'[`i',5] = /// 
		   (2/`nobs') * ( (1-`rho')^2 * `EvC'[`i',2]  ///
		   + `rho'^2 * (`sEv2' - `EvC'[`i',2])) / (`sEv'^2)
	}
	matrix rownames `ELC' = `:colnames `Ev''
	matrix colnames `ELC' = /// 
	   Eigenvalue Proportion var_Prop Cumulative var_Cum
end


program Invalid
	args optname
	
	dis as err "option `optname' only valid " /// 
	           "in combination with option vce(normal)"
	exit 198
end	


program ParseVCE, sclass
	local 0 ,`0' 
	syntax [, NORmal NONe ] 
	
	local arg `normal' `none' 
	if `:list sizeof arg' > 1 {
		exclusive_opts "`arg'" "vce()" 
	}
	
	sreturn clear
	sreturn local vce `arg' 
end	
exit