*! 1.0.1 (Jan 03, 1997) Jeroen Weesie/ICS STB-37 smv7
program define pca
    version 5.0      
    local options "Level(integer $S_level)"
    if ("`1'" != "" & substr("`1'",1,1) != ",") {
        local varlist "opt"
        local if "opt"
        local in "opt"
        local weight "aweight fweight" 
        local options "Factor(int 1) Std `options'"
        parse "`*'"
        local nvar : word count `varlist'

      * options #factors 

        if `factor' < 1 { 
            di in re "#factors should be specified/positive"
            exit 198
        }
        else if `factor' > `nvar' {
            di in bl "#factors set at `nvar'"
            local factor `nvar'
        }

        if "`weight'" == "aweight" {
            di in bl "I am not quite sure that aweights mean here..."
        }
                
      * scratch
 
        tempname W Ve L Vi b b0 bV nobs iw e ee ei ej LC rho vrho Rho Vrho 

      * compute spectral/pca-decomposition of covariance matrix

        qui mat acc `W' = `varlist' `if' `in' [`weight'`exp'], dev nocons
        sca `nobs' = _result(1)
        sca `iw' = 1/(`nobs'-1)
        mat `W' = `iw' * `W'
        if "`std'" != "" { mat `W' = corr(`W') }
        mat syme `Ve' `L' = `W' 
        mat drop `W'

      * inference on percentage variance explained (see Kshirsagar (1972:454))

        mat `LC' = J(`nvar',4,0)
        mat `LC'[1,1] = `L'[1,1]
        mat `LC'[1,2] = `L'[1,1]^2
        local i 2
        while `i' <= `nvar' {
            mat `LC'[`i',1] = `LC'[`i'-1,1] + `L'[1,`i']
            mat `LC'[`i',2] = `LC'[`i'-1,2] + `L'[1,`i']^2
            local i = `i'+1
        }

        local i 1
        while `i' <= `nvar' {
            * rho = percentage of variance explained
            mat `LC'[`i',3] = `LC'[`i',1] / `LC'[`nvar',1]
            * vrho = asymptotic variance of rho
            mat `LC'[`i',4]=(2/`nobs')*((1-`LC'[`i',3])^2*`LC'[`i',2]  /*
             */ +`LC'[`i',3]^2*(`LC'[`nvar',2]-`LC'[`i',2]))/`LC'[`nvar',1]^2
            local i = `i'+1
        }
                
        di in gr "#factors    EigenValue    %Explained     Std. Err"
        di in gr _dup(50) "-"
        local i 1
        while `i' <= `nvar' {
            * rho = percentage of variance explained
            sca `rho' = `LC'[`i',1] / `LC'[`nvar',1]
            * vrho = asymptotic variance of rho
            sca `vrho' = (2/`nobs') * ((1-`rho')^2*`LC'[`i',2] + /*
                    */ `rho'^2*(`LC'[`nvar',2]-`LC'[`i',2])) / `LC'[`nvar',1]^2
            di in ye %6.0f `i' /*
                  */ %14.4f `L'[1,`i'] %14.4f `rho' %14.4f sqrt(`vrho')    
            if `i' == `factor' {
                sca `Rho'  = `rho'
                sca `Vrho' = `vrho'
            }
            local i = `i'+1
        }

      * create b : vectorize eigenvectors 
  
        local i 1
        while `i' <= `factor' {
            mat `b0' = `Ve'[.,`i']
            mat roweq `b0' = fac`i'
            mat `b' = `b' \ `b0'
            local i = `i'+1
        }
        mat drop `b0'

      * create properly label bV (AVAR of eigenvectors based on normality), 
      * initialized to zero  

        local nr = rowsof(`b')
        local vnames: rownames(`b')
        local eqnames: roweq(`b')

        mat `bV' = J(`nr',`nr',0)
        mat rownames `bV' = `vnames'
        mat colnames `bV' = `vnames'
        mat roweq `bV' = `eqnames'
        mat coleq `bV' = `eqnames'

      * first, define diagonal block associated with variance matrixes of 
      * eigenvectors (components)

        local i 1
        local ic 1
        while `i' <= `factor' {
            mat `Vi' = J(`nvar',`nvar',0)
            local j 1    
            while `j' <= `nvar' {
                if `i' != `j' {
                    mat `e' = `Ve'[.,`j']
                    mat `ee' = `e' * `e''
                    sca `iw' = `L'[1,`i']*`L'[1,`j']/(`L'[1,`i']-`L'[1,`j'])^2
                    mat `ee' = `iw' * `ee'
                    mat `Vi' = `Vi' + `ee'
                }
                local j = `j'+1
            }    
            mat subst `bV'[`ic',`ic'] = `Vi'
            local i = `i'+1
            local ic = `ic'+`nvar'    
        }
        mat drop `Vi'

      * next define the off-diagonal blocks associated with covariances
      * of the eigenvectors (components)
        
        local i 1
        local ic 1
        while `i' <= `factor' {
            mat `ei' = `Ve'[.,`i']
            local j 1
            local jc 1    
            while `j' < `i' {
                mat `ej' = `Ve'[.,`j']
                mat `ee' = `ej' * `ei''
                sca `iw' = -`L'[1,`i']*`L'[1,`j']/(`L'[1,`i']-`L'[1,`j'])^2
                mat `ee' = `iw' * `ee'
                mat subs `bV'[`ic',`jc'] = `ee'
                mat `ee' = `ee''
                mat subs `bV'[`jc',`ic'] = `ee'
                local j = `j'+1
                local jc = `jc'+`nvar'
            }
            local i = `i'+1
            local ic = `ic'+`nvar'
        }

      * save results

        mat `b' = `b''     
        sca `iw' = 1 / `nobs'
        mat `bV' = `iw' * `bV'
        mat post `b' `bV'

        global S_E_cmd  "pca"        /* name of command */
        global S_E_nobs = `nobs'     /* number of observations (cases) */
        global S_E_rho  = `Rho'      /* %explained variance */
        global S_E_vrho = `Vrho'     /* asymptotic variance of rho */
        global S_E_fact `factor'     /* #factors/components */
    }
    else {
        if ("$S_E_cmd" != "pca") { error 301 }
        parse "`*'"
    }

  * display as equations

    if (`level'<10 | `level'>99) { local level 95 }
    if "`std'" == "" { local clab "covariance" }
    else               local clab "correlation"
            
    #del ;    
    di _n       in gr "Principal components of" in ye " `clab' " in gr "matrix" 
       _col(52) in gr "Number of obs     =" in ye %8.0f $S_E_nobs _n
                in gr "(Std. errors based on multivariate normality)"
       _col(52) in gr "Number of factors =" in ye %8.0f $S_E_fact _n
       _col(52) in gr "Rho (%var expl.)  =" in ye %8.4f $S_E_rho  _n 
       _col(52) in gr "Std Err Rho       =" in ye %8.4f sqrt($S_E_vrho) _n ;
    #del cr                    
    mat mlout, level(`level')
    di in gr ""
end
exit