*
* Example 9.3.1, pp 416-420
*
open data m-barra-9003.txt
calendar(m) 1990
data(format=free,org=columns) 1990:01 2003:12 age c mwd mer dell hpq ibm aa cat pg
*
dec vect[strings] vl(10)
ewise vl(i)=%l(i)
*
sweep(coeffs=cconst,series=demean)
# age to pg
# constant
compute [vector] means  =tr(%xrow(cconst,1))
compute [vector] stderrs=%sqrt(%xdiag(%sigma)*%nobs/(%nobs-1))
compute covexcess=%sigma
*
report(action=define)
report(atcol=1,atrow=2,fillby=cols) vl
report(atcol=2,atrow=1,fillby=cols) "mean"     means
report(atcol=3,atrow=1,fillby=cols) "std. err" stderrs
report(action=format,picture="*.##")
report(action=show)
*
dec vect group(10)
compute group=||1.0,1.0,1.0,1.0,2.0,2.0,2.0,3.0,3.0,3.0||
dec rect b(10,3)
ewise b(i,j)=(group(i)==j)
*
make x
# demean
*
*  OLS estimates
*
compute fprj=b*inv(tr(b)*b)
compute f=x*fprj
compute u=x-f*tr(b)
*
*  GLS estimates
*
compute uu=%innerxx(u)/(%nobs-1)
compute dd=inv(%diag(%xdiag(uu)))
compute fprj=dd*b*inv(tr(b)*dd*b)
compute f=x*fprj
compute u=x-f*tr(b)
compute uu=%innerxx(u)/(%nobs-1)
*
*  Residual covariance matrix
*
compute dg=%diag(%xdiag(uu))
compute ff=%innerxx(f)
compute ff=(ff-%nobs*%innerxx(%fill(1,%nobs,1.0/%nobs)*f))/(%nobs-1)
compute covmodel=b*ff*tr(b)+dg
*
*  Display the correlation matrices
*
compute [rect] correxcess=%cvtocorr(covexcess)
report(action=define)
report(atrow=1,atcol=2,fillby=rows) vl
report(atrow=2,atcol=1,fillby=cols) vl
report(atrow=2,atcol=2) correxcess
report(action=format,picture="*.#")
report(action=show)
*
compute [rect] corrmodel=%cvtocorr(covmodel)
report(action=define)
report(atrow=1,atcol=2,fillby=rows) vl
report(atrow=2,atcol=1,fillby=cols) vl
report(atrow=2,atcol=2) corrmodel
report(action=format,picture="*.#")
report(action=show)
*
* Extract the factors and graph
*
set financial = f(t,1)
set hightech  = f(t,2)
set other     = f(t,3)
*
spgraph(footer="Figure 9.3 Estimated factor realizations of a BARRA industrial factor model",vfields=3)
graph(header="(a) Financial sector")
# financial
graph(header="(b) High-tech sector")
# hightech
graph(header="(c) Other sector")
# other
spgraph(done)