*
* @PRINFACTORS( options ) sigma
*
* @PRINFACTORS does a principal components based factor analysis of an
* input covariance or correlation matrix.
*
* Parameters:
*   sigma     input SYMMETRIC array (covariance or correlation matrix)
*
* Options:
*   NCOMPS=number of components [number of variables]
*   PRINT/[NOPRINT]  Controls whether the output will be displayed.
*   SIGMAHAT=output sigma matrix [not used]
*   LOADINGS=RECTANGULAR of loadings of components onto variables [not used]
*   VECTORS=RECTANGULAR of eigenvectors [not used]
*   VALUES=VECTOR of eigenvalues [not used]
*   COMMUNALITIES=VECTOR of communalities for variables [not used]
*
* Revision Schedule
*   05/2005 Written by Tom Doan, Estima
*   08/2006 Corrections to the comments describing the procedure syntax.
*   03/2007 Define %VARIANCE. Add COMMUNALITIES, VECTORS and VALUES option.
*           Sign flip loading vectors which are primarily negative.
*   09/2010 Use LOCAL REPORT
*
procedure prinfactors sigma
type symm sigma
*
option integer ncomps   0
option switch  print    0
option symm    *sigmahat
option rect    *loadings
option rect    *vectors
option vect    *values
option vect    *communalities
*
local rect 		 eigenvec
local vect 		 eigenval princoef scalings
local int  		 n nvars icomp
local real 		 totalvar running
local report 	preport
local string 	ltitle
*
compute nvars=%rows(sigma)
*
* Do an eigen decomposition of the input matrix.
*
eigen sigma eigenval eigenvec
if %defined(values)
   compute values=eigenval
if %defined(vectors)
   compute vectors=eigenvec
*
* Figure out how many components need to be computed
*
if .not.ncomps
   compute n = nvars
else
   compute n = ncomps

compute ltitle="Principal Components Analysis"
if print {
   report(use=preport,action=define,title=ltitle)
   report(use=preport,atrow=1,atcol=2,span) ltitle
   report(use=preport,atcol=1,atrow=2,fillby=cols) "Eigenvalue" "Proportion" $
       "Cumulative" "Eigenvector"
}
*
* Get the total variation for normalizing the contributions
*
compute totalvar=%sum(eigenval)
*
* Initialize the running sum of the covariances
*
if %defined(sigmahat)
   compute sigmahat=%zeros(nvars,nvars)
*
* And the running sum of the variance explained
*
if %defined(communalities)
   compute communalities=%zeros(nvars,1)

compute running=0.0
do icomp=1,n
   compute running=running+eigenval(icomp)
   compute princoef=%xcol(eigenvec,icomp)
   if %defined(communalities)
      compute communalities=communalities+eigenval(icomp)*princoef.^2
   if print
      report(use=preport,col=new,atrow=2,fillby=cols) eigenval(icomp) $
         eigenval(icomp)/totalvar running/totalvar princoef
   if %defined(sigmahat)
      compute sigmahat=sigmahat+%outerxx(princoef)*eigenval(icomp)
end do icomp
if print
  report(use=preport,action=show)
*
* Replace the diagonal elements of sigmahat with the input elements
*
if %defined(sigmahat) {
   do icomp=1,nvars
      compute sigmahat(icomp,icomp)=sigma(icomp,icomp)
   end do i
}
compute %variance   =totalvar
*
* Compute the loadings if requested. Sign flip column which is
* predominantly negative.
*
if %defined(loadings) {
   dim loadings(nvars,n)
   dim scalings(n)
   do j=1,n
      compute scalings(j)=sqrt(eigenval(j))*%sign(%sum(%xcol(eigenvec,j)))
   end do j
   ewise loadings(i,j)=eigenvec(i,j)*scalings(j)
}
*
end
