cap prog drop coldiag
program define coldiag
/*
This program calculates condition number, singular values,
and pi matrix of variance decomposition portions according
to Belsey, Kuh and Welsh 1980
Adapted from program by Jeroen Weesie/ICS  STB-39 dm49
code was modelled after -cond- in Matlab 4.0.
Updated for version 6.0 7/6/00
*/
     version 6.0
	syntax varlist [if] [in][, Full NOScale]
	tempname B v U w V  maxw minw cc  Phi phi Pi SV Y
	qui mat accum `B' = `varlist' `if' `in', noc
	local nrB=rowsof(`B')
*	matrix list `B'
	matrix `v'=diag(vecdiag(`B'))
	local i 1
	if "`noscale'"=="" { 
* 	scale A to unit length	
	while  `i'<=`nrB' {
		matrix `B'[`i',`i']=1
		local j=`i'+1
		while  `j'<=`nrB' {
			matrix `B'[`i',`j']=`B'[`i',`j']/(sqrt(`v'[`i',`i'])*sqrt(`v'[`j',`j']))
			matrix `B'[`j',`i']=`B'[`i',`j']
			local j=`j'+1
		}
		local i=`i'+1
	}
	}

	mat svd `U' `w' `V' = `B'
	local nw = colsof(`w')
	local i 1
	while `i'<=`nw' {
		matrix `w'[1,`i']=sqrt(`w'[1,`i'])
		local i=`i'+1
	}

     scalar `maxw' = 0
     scalar `minw' = 1E32
     local i 1
     while `i' <= `nw' {
     	scalar `maxw' = max(`maxw', `w'[1,`i']) 
          scalar `minw' = min(`minw', `w'[1,`i']) 
          local i = `i'+1
    }
 	matrix `SV'=J(`nw',1,0)   
    scalar `cc' = .
    if (`minw' > 0) { scalar `cc' = `maxw' / `minw' }
	disp _newline(1)
	di in gr "condition number = " in ye %12.2f scalar(`cc')

	if "`full'" == "" {
		exit
	}

	disp _s(2)
	di in gr "Singular values:"	
	local i 1
	while `i' <= `nw' {
		di in gr "   " %3.0f `i' ":" in ye /*
	*/ 		%12.2f `maxw'/`w'[1,`i']
		matrix `SV'[`i',1]= round(`maxw'/`w'[1,`i'],.01)
		local i = `i'+1
    }

	cap mat `Phi'=`V'
	mat `phi'=J(1,`nw',0)
	local i 1
	while `i'<=`nw' {
		local j 1
		while `j'<=`nw' {
			matrix `Phi'[`i',`j']= `Phi'[`i',`j']^2/`w'[1,`j']^2
			matrix `phi'[1,`i']= `phi'[1,`i']+ `Phi'[`i',`j']
			local j=`j'+1
		}
		local i=`i'+1
	}
	matrix `Pi'=J(`nw',`nw',0)
	local i 1
	while `i'<=`nw' {
		local j 1
		while `j'<=`nw' {
			matrix `Pi'[`j',`i']= round(`Phi'[`i',`j']/`phi'[1,`i'],.001)
			local j=`j'+1
		}
		local i=`i'+1
	}
	matrix `Y'=`SV',`Pi'
	matrix colnames `Y'=SV `varlist'
	matrix list `Y', noheader

end