procrustes.ado

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

*! version 1.0.0  25feb2005
program	procrustes, byable(onecall)
	version 9

	if replay() {
		if "`e(cmd)'" != "procrustes" {
			error 301
		}
		if _by() {
			error 190
		}
		Display `0'
		exit
	}

	if _by() {
		by `_byvars'`_byrc0': Estimate `0'
	}
	else {
		Estimate `0'
	}
end


program Estimate, eclass byable(recall)

	#del ;
	syntax  anything [if] [in] [aw fw]
	[,
		TRansform(str) 
		noCONstant
		noRHo
		FORCE
		noFIt
		SCale  // undocumented option transform(ortho)
		LOG    // undocumented option transform(oblique)
	] ;
	#del cr
	
	local display_options `fit' 
	
	ParseTransform  `transform' 
	local transform `s(transform)' 
	
	gettoken ylist anything : anything, match(parens)
	gettoken xlist anything : anything, match(parens)
	if `"`anything'"' != "" { 
		dis as err `"unexpected input `anything'"' 
		exit 198
	}	
	if (`"`xlist'"' == "") | (`"`ylist'"' == "") {
		dis as err "two varlists expected"
		exit 198
	}
	unab xlist : `xlist' , name(varlist_x) min(1) 
	unab ylist : `ylist' , name(varlist_y) min(1) 
	if ("`:list dups xlist'" != "") & ("`force'" == "") {
		dis as err "duplicate variables in varlist-x"
		dis as err "specify option force to proceed anyway"
		exit 198
	}
	if ("`:list dups ylist'" != "") & ("`force'" == "") {
		dis as err "duplicate variables in varlist-y"
		dis as err "specify option force to proceed anyway"
		exit 198
	}
	if ("`:list xlist & ylist'" != "") & ("`force'" == "") {
		dis as err "varlist-x and varlist-y are not disjoint"
		dis as err "specify option force to proceed anyway"
		exit 198
	}
	confirm numeric var `xlist' `ylist'

	if "`weight'" != "" { 
		local wght [`weight'`exp'] 
	}	
	
	marksample touse, novarlist
	quietly summarize `touse' if `touse' `wght', meanonly
	local n0 = r(N)

	markout `touse' `xlist' `ylist'
	quietly summarize `touse' if `touse' `wght', meanonly
	local nobs = r(N)

	if (`nobs' == 0) error 2000
	if (`nobs' == 1) error 2001

	local nx : list sizeof xlist
	local ny : list sizeof ylist

// optimal transformation (A,rho,c)

	if "`transform'" == "orthogonal" {  
		Ortho "`ylist'" "`xlist'" `"`wght'"' "`touse'" /// 
		   "`scale'" "`constant'" "`rho'"
	}
	else if "`transform'" == "oblique" { 
		Oblique "`ylist'" "`xlist'" `"`wght'"' "`touse'" /// 
		   "`scale'" "`constant'" "`rho'" "`log'"
	}
	else if "`transform'" == "unrestricted" {
		local rho norho
		Unrestricted "`ylist'" "`xlist'" `"`wght'"' "`touse'" /// 
		   "`scale'" "`constant'" "`rho'"
	}
	else {
		_stata_internalerror	
	}

// extract results 

	tempname A c Rho
	
	matrix `A'   = r(A)
	matrix `c'   = r(c)
	scalar `Rho' = r(rho)
	
	local uniqueA  `r(uniqueA)' 

// degrees of freedom of model and residual

	tempname df_m df_r
	
	scalar `df_m' = `ny'*`nx' + ("`rho'"=="") + `ny'*("`constant'"=="")
	if "`transform'" == "orthogonal" { 
		scalar `df_m' = `df_m' - `nx'*(`nx'+1)/2 
	}
	else if "`transform'" == "oblique" { 
		scalar `df_m' = `df_m' - `nx'
	}
	scalar `df_r' = `nobs'*`ny' - `df_m'

// fit statistics per y-variable

	tempname b P rmse rss ss urmse ystats
	tempvar  res yhat

	matrix `ystats' = J(5,`ny',.)
	matrix colnames `ystats' = `ylist' 
	matrix rownames `ystats' = SS RSS RMSE Procrustes Corr_y_yhat
	
	local iy = 0
	foreach y of local ylist { 
		capture drop `yhat'
		capture drop `res' 
		local ++iy
	// ss
		SS `y' `touse' "`wght'" "`constant'" 
		matrix `ystats'[1,`iy'] = r(ss)
		
	// rss	
		matrix `b' = `A'[1...,`iy']'
		matrix score double `yhat' = `b' if `touse'
		quietly replace `yhat' = `Rho'*`yhat' + `c'[1,`iy'] if `touse'
		
		quietly gen double `res' = (`y'-`yhat')^2  if `touse' 
		quietly summ `res' if `touse' `wght' , meanonly 
		matrix `ystats'[2,`iy'] = r(N)*r(mean)
		
	// rmse	
		matrix `ystats'[3,`iy'] = sqrt(`ystats'[2,`iy']/(`df_r'/`ny'))
		
	// P = rss/ss	
		matrix `ystats'[4,`iy']	= `ystats'[2,`iy']/`ystats'[1,`iy']
		
	// pcorr(y,yhat)	
		quietly corr `y' `yhat' if `touse' `wght' 
		matrix `ystats'[5,`iy'] = r(rho)
	}
	
// overall statistics
	
	scalar `ss'  = 0 
	scalar `rss' = 0
	forvalues iy = 1/`ny' {
		scalar `ss'  = `ss'  + `ystats'[1,`iy']
		scalar `rss' = `rss' + `ystats'[2,`iy']
	}

	// procrustes statistic
	scalar `P'     = `rss' / `ss'
	scalar `urmse' = sqrt( `rss' / (`nobs'*`ny'))
	scalar `rmse'  = sqrt( `rss' / `df_r' )

// save results 

	ereturn post, obs(`nobs') esample(`touse') properties(nob noV)

	ereturn local ylist       `ylist'
	ereturn local xlist       `xlist'
	ereturn local wtype       `weight'
	ereturn local wexp        `"`exp'"' 

	ereturn local transform   `transform'
	ereturn local uniqueA     `uniqueA' 
	
	ereturn matrix A          = `A'
	ereturn matrix c          = `c'
	ereturn scalar rho        = `Rho'
	
	ereturn matrix ystats     = `ystats'  

	ereturn scalar N          = `nobs'
	ereturn scalar ss         = `ss' 
	ereturn scalar rss        = `rss'
	ereturn scalar P          = `P'
	ereturn scalar rmse       = `rmse'
	ereturn scalar urmse      = `urmse'
	ereturn scalar df_m       = `df_m'
	ereturn scalar df_r       = `df_r'

	ereturn local predict     procrustes_p
	ereturn local estat_cmd   procrustes_estat
	
	ereturn local cmd         procrustes

// display	

	Display , `display_options' 
end

// ----------------------------------------------------------------------------

program Display
	syntax [, noFIt ]

// header

	dis _n as txt "Procrustes analysis (`e(transform)')" ///
	    _col(46) "Number of observations" ///
	    _col(69) "= " as res %9.0f e(N)
	    
	dis _col(46) as txt "Model df (df_m)"    ///
	    _col(69) "= " as res %9.0f e(df_m)
	dis _col(46) as txt "Residual df (df_r)" ///
	    _col(69) "= " as res %9.0f e(df_r)

	if "`e(scale)'" != "" {
		dis as txt "X and Y are scaled to trace(XX') = trace(YY') = 1"
	}
	
	dis _col(46) as txt "SS(target)" ///
	    _col(69) "= " as res %9.0g e(ss)
	dis _col(46) as txt "RSS(target)" ///
	    _col(69) "= " as res %9.0g e(rss)
	dis _col(46) as txt "RMSE = root(RSS/df_r)" ///
	    _col(69) "= " as res %9.0g e(rmse)
	dis _col(46) as txt "Procrustes = RSS/SS" ///
	    _col(69) "= " as res %9.4f e(P)

// shift c

	dis as txt "Translation c"
	matlist e(c), left(4) border(top bot) 

// rotation A

	dis _n as txt "Rotation & reflection matrix A (`e(transform)')"
	matlist e(A), left(4) border(top bot) nohalf
	
	if "`e(uniqueA)'" == "0" {
		dis as txt _col(5) "(Warning: rotation A not unique)"
	}

// dilation factor

	if "`e(transform)'" != "unrestricted" { 
		dis _n as txt "Dilation factor " _n(2) ///
		    _col(14) "rho = " as res %9.4f e(rho)
	}	
	
// statistics

	if "`fit'" == "" { 
		dis _n as txt "Fit statistics by target variable"
		matlist e(ystats), row(Statistics) /// 
		   left(4) border(top bot) 
	}
end


// ============================================================================
// procrustes analysis:  y(i) = c + rho A' x(i) + e(i) , i = 1..nobs
//  in matrix notation:    Y  = 1c' + rho X A + E
//
//          orthogonal:  A is orthonormal A'A = AA' = I 
//             oblique:  A is normal      diag(A'A) = 1
//        unrestricted:  no restrictions on A, and rho = 1
//
// if norho  is specified,  rho is fixed at 1
// if nocons is specified,  cons = 0-vector
// ============================================================================
	
program Ortho, rclass
	args ylist xlist wght touse scale cons rho
	
	tempname A b c P Rho 
	tempname m traceX traceY traceW U V W x XtX YtX YtY Z 

// append zero x-vars so that # yvars = # xvars 

	local ny  : list sizeof ylist  
	local nx0 : list sizeof xlist 
	
	if `ny' < 2 {
		dis as err "at least two y-variables expected"
		exit 102
	}
	if `nx0' < 2 {
		dis as err "at least two x-variables expected"
		exit 102
	}
	if `nx0' > `ny' {
		dis as err "#x-variables should not exceed #y-variables"
		exit 198
	}
	local xxlist `xlist'
	if `nx0' < `ny' {
		tempvar zero
		gen byte `zero' = 0
		forvalues j = `=`nx0'+1' / `ny' {
			local xxnames `xxnames' zero
			local xxlist  `xxlist'  `zero'
		}
	}
	
// unscaled cross-products [X Y]'[X Y] 

	if "`cons'" == "" {
		local Dev dev
	}
	quietly matrix accum `Z' = `xxlist' `ylist' `wght' /// 
	   if `touse' , `Dev' means(`m') nocons

	local ny   : list sizeof ylist 
	local p    = `ny'   
	local pp1  = `p'+1

	matrix `XtX' = `Z'[1..`p', 1..`p']
	matrix `YtY' = `Z'[`pp1'..., `pp1'...]
	matrix `YtX' = `Z'[`pp1'..., 1..`p']
	matrix drop `Z'

	_m2scalar `traceX' = trace(`XtX')
	_m2scalar `traceY' = trace(`YtY')
	
	if `traceX' < 1e-8 {
		dis as err "no variation in x-variables" 
		exit 198
	}	
	if `traceY' < 1e-8 {
		dis as err "no variation in y-variables" 
		exit 198
	}	
	
	if "`scale'" != "" {
		matrix `XtX' = `XtX' / `traceX'
		matrix `YtY' = `YtY' / `traceY'
		matrix `YtX' = `YtX' / sqrt(`traceX'*`traceY')
		scalar `traceX' = 1
		scalar `traceY' = 1
	}

	matrix svd `U' `W' `V' = `YtX'

	matrix `A' = `V' * `U''
	_m2scalar `traceW' = `W' * J(colsof(`W'),1,1)
	
	UniqueA `W' 
	local uniqueA = `r(unique)'

	if "`rho'" == "" {
		scalar `Rho' = `traceW'/`traceX'
		// equivalent, 
		// matlist trace(`A'*`YtX') / trace(`XtX') 
	}
	else {
		scalar `Rho' = 1
	}

	if "`cons'" == "" {
		matrix `c' = `m'[1,`pp1'...] - `Rho'*`m'[1,1..`p']*`A'
	}
	else {
		matrix `c' = J(1,`ny',0)
		matrix colnames `c' = `ylist'
		matrix rownames `c' = _cons
	}
	
	if `nx0' < `ny' {
		matrix `A' = `A'[1..`nx0',1...]
	}
	
	return matrix A       = `A'
	return matrix c       = `c'
	return scalar rho     = `Rho'
	return local  uniqueA = `uniqueA' 
end


// ----------------------------------------------------------------------------
// Oblique Procrustes -- we may treat each y-variable separately if rho==1