*
* Example 9.3 from pp 431-432
*
open data m-bnd.txt
calendar(m) 1942
data(format=free,org=columns) 1942:1 1999:12 us30yr us20yr us10yr us5yr us1yr
*
vcv(center)
# us30yr us20yr us10yr us5yr us1yr
*
compute r=%cvtocorr(%sigma)
*
@prinfactors(ncomp=2,noprint,loadings=beta,commune=commune) r
compute betastar=beta
@varimax(standardize) betastar
*
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) "30 Years" "20 Years" "10 Years" "5 Years" "1 Year" "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)
*
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) "30 Years" "20 Years" "10 Years" "5 Years" "1 Year" "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)