*
* Example 9.2.1, pp 408-412
*
open data m-fac9003.txt
calendar(m) 1990
data(format=free,org=columns) 1990:01 2003:12 aa age cat f fdx gm hpq kmb mel nyt pg trb txn sp
*
dec vect[strings] vl(13)
ewise vl(i)=%l(i)
*
* Regress all the returns on "constant" only to get the means, sample standard
* errors and covariance matrix of the excess returns. Because we'll need it later,
* pull the top 13x13 corner out. (A TABLE instruction would work fine to get the
* information in table 9.1 - we're using this because we need the covariance matrix
* as well).
*
sweep(coeffs=cconst)
# aa to sp
# constant
compute raw=%sigma
compute [vector] means  =tr(%xrow(cconst,1))
compute [vector] stderrs=%sqrt(%xdiag(raw)*%nobs/(%nobs-1))
compute covexcess=%xsubmat(raw,1,13,1,13)
*
report(action=define)
report(atcol=1,atrow=2,fillby=cols) vl "SP500"
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)
*
* Now regress all the returns on constant and the sp500. The betas will be in the
* second row of the cbeta matrix. The sigmas can be pull out of the diagonal of
* the %sigma array. Because the book shows these with degrees of freedom
* correction, we need to rescale. (%sigma just uses a T divisor, not T-2).
*
sweep(coeffs=cbeta)
# aa to txn
# constant sp
*
compute [vector] betas    =tr(%xrow(cbeta,2))
compute [vector] sigmas   =%sqrt(%xdiag(%sigma)*%nobs/(%nobs-2))
compute [vector] rsquareds=1.0-%xdiag(%sigma)./%xdiag(covexcess)
*
report(action=define)
report(atcol=1,atrow=2,fillby=cols) vl
report(atcol=2,atrow=1,fillby=cols) "beta.hat" betas
report(atcol=3,atrow=1,fillby=cols) "sigma"    sigmas
report(atcol=4,atrow=1,fillby=cols) "r.square" rsquareds
report(action=format,picture="*.###")
report(action=show)
*
set betacopy 1 13 = betas(t)
set rsqrcopy 1 13 = rsquareds(t)
spgraph(vfields=2)
graph(style=bar,xlabels=vl)
# betacopy
graph(style=bar,xlabels=vl)
# rsqrcopy
spgraph(done)
*
* Compute the model-based covariance matrix
*
compute covmarket=raw(14,14)*%outerxx(betas)+%diag(%xdiag(%sigma))
compute [rect] corrmarket=%cvtocorr(covmarket)
report(action=define)
report(atrow=1,atcol=2,fillby=rows) vl
report(atrow=2,atcol=1,fillby=cols) vl
report(atrow=2,atcol=2) corrmarket
report(action=format,picture="*.##")
report(action=show)
*
compute [rect] correxcess=%cvtocorr(%xsubmat(raw,1,13,1,13))
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 the GMVP's for the market model and for the excess returns
* Note that the direct calculation will not impose non-negativity. If
* you want to restrict yourself to non-negative weights, you can use
* the instruction LQPROG to solve the quadratic program
*
compute ones=%fill(13,1,1.0)
compute omegamarket=inv(covmarket)*ones/%qform(inv(covmarket),ones)
compute omegaexcess=inv(covexcess)*ones/%qform(inv(covexcess),ones)
*
report(action=define)
report(atrow=2,atcol=1,fillby=cols) vl
report(atrow=1,atcol=2,fillby=cols) "GMVP-Market" omegamarket
report(atrow=1,atcol=3,fillby=cols) "GMVP-Excess" omegaexcess
report(action=format,picture="*.####")
report(action=show)
*
* Show correlation matrix of residuals
*
compute [rect] corrresids=%cvtocorr(%sigma)
report(action=define)
report(atrow=1,atcol=2,fillby=rows) vl
report(atrow=2,atcol=1,fillby=cols) vl
report(atrow=2,atcol=2) corrresids
report(action=format,picture="*.##")
report(action=show)