princomp and varimax

[Kraus]

Hello, I am an unfrequent user of RATS and therefore may be asking a stupid question. Anyway. What I am trying to do is to extract the principal components of a set of time series and then rotate them to simplify the interpretation. As I understand it, the two procedures princomp and varimax are to tools to do that. However, I don't get it how to use the output of the first as an input for the second. Any help is appreciated. Thank you.

[TomDoan]

This is an example from Tsay's book. You need to save the loadings from @PRINCOMP and feed them into VARIMAX. (That's also a QUARTIMAX, which is used in the TSAYP495 exampe).

Code: Select all

*
* Tsay, Analysis of Financial Time Series, 3rd edition
* Example 9.2 from pp 493-494
*
open data m-5clog-9008.txt
calendar(m) 1980:1
data(format=prn,org=columns) 1980:01 1998:12 ibm hpq intc jpm bac
*
* Compute multivariate Q statistics
*
dofor lags = 1 4 8
   @mvqstat(lags=lags)
   # ibm hpq intc jpm bac
end dofor
*
* Compute the (means subtracted) correlation matrix of the data
*
vcv(center,matrix=r)
# ibm hpq intc jpm bac
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" "JPM" "BAC" "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" "JPM" "BAC" "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)