*
* Example 9.2 from pp 430-431
*
open data m-5cln.dat
calendar(m) 1990
data(format=free,org=columns) 1990:1 1999:12 ibm hwp intc mer mwd
*
* Compute multivariate Q statistics
*
dofor lags = 1 4 8
   @mvqstat(lags=lags)
   # ibm hwp intc mer mwd
end dofor
*
* Compute the (means subtracted) correlation matrix of the data
*
vcv(center,matrix=r)
# ibm hwp intc mer mwd
compute r=%cvtocorr(r)
*
* Compute the principal components estimates of the factors
*
@prinfactors(ncomp=2,noprint,loadings=beta,communalities=commune) r
*
compute betastar=beta
@varimax(standardize) betastar
*
report(action=define)
report(atrow=1,atcol=1) "Variable" "f1" "f2" "f1*" "f2*" "Communalities"
report(atrow=2,atcol=2,tocol=5,span,align=centered) "Principal Component Method"
report(atrow=3,atcol=1,fillby=cols) "IBM" "HPQ" "INTC" "MER" "MWD" "Variance"
report(atrow=3,atcol=2) beta
report(atrow=3,atcol=4) betastar
report(atrow=3,atcol=6,fillby=cols) commune
report(atrow=8,atcol=2) %normsqr(%xcol(beta,1)) %normsqr(%xcol(beta,2)) %normsqr(%xcol(betastar,1)) %normsqr(%xcol(betastar,2))
report(action=format,picture="*.###")
report(action=show)
*
* Maximum likelihood factor analysis. Note that this uses a different
* normalization for the factors than the one used in the text. A model with m
* factors needs m(m-1)/2 restrictions in order to be identified. There are many
* ways to accomplish this. That's why choosing a rotation is important, since the
* normalization chosen by the software may have little meaning in a particular
* situation.
*
compute nfactor=2,n=5
dec symm   llead(nfactor,nfactor)
dec rect   lrem(n-nfactor,nfactor)
dec vect   d(n)
nonlin llead lrem d
compute llead=%identity(nfactor)
compute lrem =%const(0.0)
compute d=%sqrt(%xdiag(r))
*
declare real test
find(method=bfgs,noprint) maximum test
   compute lambda=llead~~lrem
   compute sigma=%outerxx(lambda)+%diag(d.*d)
   compute test=%logdensitycv(sigma,r,%nobs)
end find
*
dim commune(n)
ewise commune(i)=%normsqr(%xrow(lambda,i))
compute lambdastar=lambda
@varimax(standardize) lambdastar
report(action=define)
report(atrow=1,atcol=1) "Variable" "f1" "f2" "f1*" "f2*" "Communalities"
report(atrow=2,atcol=2,tocol=5,span,align=centered) "Maximum Likelihood Method"
report(atrow=3,atcol=1,fillby=cols) "IBM" "HPQ" "INTC" "MER" "MWD" "Variance"
report(atrow=3,atcol=2) lambda
report(atrow=3,atcol=4) lambdastar
report(atrow=3,atcol=6,fillby=cols) commune
report(atrow=8,atcol=2) %normsqr(%xcol(lambda,1)) %normsqr(%xcol(lambda,2)) %normsqr(%xcol(lambdastar,1)) %normsqr(%xcol(lambdastar,2))
report(action=format,picture="*.###")
report(action=show)