*
* Example 11.2, pp 533-535
*
open data m-fac9003.txt
calendar(m) 1990
data(format=free,org=columns) 1990:1 2003:12 aa age cat f fdx gm hpq kmb mel nyt pg trb txn sp
set gm = gm*.01
set sp = sp*.01
*
* Fixed coefficients model done with linreg
*
linreg gm
# constant sp
*
* Done with DLM. sv=1.0 with VAR=CONCENTRATE concentrates out the variance of the
* measurement equation. FREE=2 has no effect on the estimates, but corrects the
* likelihood function so it matches LINREG's by computing the likelihood
* conditional on the (freely estimated) coefficients.
*
dlm(y=gm,c=||1.0,sp||,sv=1.0,var=concentrate,exact,free=2,type=smoothed) / xstates vstates
*
* With the smoothed estimates, the state means and variances are the same for all
* time periods. Note that the %VARIANCE (the concentrated estimate of the
* measurement equation variance) isn't corrected for degrees of freedom, while the
* calculation used in the text is.
*
disp xstates(2003:12)
disp tr(%sqrt(%xdiag(vstates(2003:12))*%variance))
*
* Now estimate the component variances. Note that we're pegging sv=1.0 and using
* var=concentrate to  concentrate out the unknown measurement variance. The other
* variance estimates will have to be multiplied by %VARIANCE to get their true
* estimated values. In component models like this, it can be very difficult to
* estimate all the variances directly because they are so highly correlated. In
* this case, the estimates given in the book aren't really the correct likelihood
* maximizers. The likelihood is a bit higher with a somewhat higher variance on
* beta.
*
nonlin leps leta
compute leps=leta=0.0
dlm(y=gm,c=||1.0,sp||,sv=1.0,var=concentrate,sw=%diag(||exp(leta),exp(leps)||),$
   exact,method=bfgs,type=smoothed) / xstates vstates
disp "Component Std Dev"
disp "Measurement" @20 sqrt(%variance)
disp "Alpha" @20 sqrt(%variance*exp(leta))
disp "Beta" @20 sqrt(%variance*exp(leps))
*
set alpha = xstates(t)(1)
set beta  = xstates(t)(2)
set expret = %dot(||1.0,sp||,xstates)
*
spgraph(hfields=2,vfields=2,footer="Figure 11.5 Time plots from a time-varying CAPM")
graph(hlabel="Excess return")
# gm
graph(hlabel="Expected return")
# expret
graph(hlabel="Alpha")
# alpha
graph(hlabel="Beta")
# beta
spgraph(done)