######################################################################################
# ZW Ch 6 Example 2 (Dynamic regressions) Commented and Expanded.
######################################################################################

module(finmetrics)

#First, we recreate the "excess return" data set that was used in Example 1.

excessRet.ts <- seriesMerge(ret.ts,log(1+rf.30day))
excessRet.ts[,"MSFT"] <- excessRet.ts[,"MSFT"] - excessRet.ts[,"RF"]
excessRet.ts[,"SP500"] <- excessRet.ts[,"SP500"] - excessRet.ts[,"RF"]
excessRet.ts <- excessRet.ts[,1:2]


#
# Now, estimate an ADL version of CAPM for Microsoft using the 
#

adl.fit = OLS(MSFT~SP500+ar(1)+tslag(SP500),data=excessRet.ts)
summary(adl.fit)

#Extract the coefficients into a variable...

bhat = coef(adl.fit)
bhat

# Now look at the "long-range influence of a shock"

lr.beta = (bhat[2]+bhat[3])/(1-bhat[4])
lr.beta

#You might not have expected this out-put ... I did'n. Let's check the class

class(lr.beta)

#This is suggestive... but why this class?
#Let's look at bhat

class(bhat)

#We deduce that the lr.beta inherits the "name" from the coef vector. Note that the 
#name comes from the 2nd coordinate of the coef, each of which has a name.

###############################################
# Now for something a little more exciting...
# DYNAMIC REGRESSIONS....
# ....a name that you can sell, and even so it makes sense.
############################################################

# Part One: distributed lag regressions

dl.fit = OLS(diff(log(GDP))~FFR+tslag(FFR,1:12),data=policy.dat,
start="Jan 1990",in.format="%m %Y",na.rm=T)

#Just for humor lets apply the Information Criterion function to IC

IC(dl.fit,type="AIC")
IC(dl.fit,type="BIC")
IC(dl.fit, type="loglike")

#And how about two colinearity tests?
# First, a little background

help(collinearTest)

#Now the tests....

collinearTest(dl.fit,method="cn")
collinearTest(dl.fit,method="vif")

#OUCH... what we see here is a disgusting amount of colinearity.
#Could we have guessed this? 
#After the fact, it's understandable, but ex ante, people stumble into this all the time.

#Is there a cure? Even a herbalistic folk treatment?

########################################
# Part Two: pdl model for real gdp growth
###############################################

pdl.fit = OLS(diff(log(GDP))~pdl(FFR,d=2,q=12),
data=policy.dat,start="Jan 1990",
in.format="%m %Y",na.rm=T)
pdl.fit
summary(pdl.fit)

###############################################
# Part Three: Dyanmic CAPM predictions
###############################################

adl.fit = OLS(MSFT~SP500+ar(1)+tslag(SP500),data=excessRet.ts)

#Now, a leprechaun tells us the next three months' returns are -20%, 0% and, +20%,
#and we put this information into a data frame.

sp500.new = data.frame(c(-0.2,0, 0.2))

#and we tell S that this information is associated with SP500

colIds(sp500.new) = "SP500"

#The big puzzle: how does this leperchaun hint inform us about MSFT?

predict(adl.fit, n.predict=3, olddata=as.data.frame(excessRet.ts@data),
newdata=sp500.new)

#Originally I had "2" in stead of "0.2" for the 20% month and I got the
#silly answer that MSFT was forcast to go up 300 percent. Ironically, 
#depsite the coaching I give you, I did not check for a typo and almost
#slid off on a philosophical rocket. Now that the typo is fixed, the result
#look quite reasonable, though the forecast intervals are prehaps wider thatn
#on might wish.


############################################################################################
#The rest of this file contains "Advanced" material that will not be covered in 434.
#
# computing HC and HSCE values
#
ols.fit = OLS(MSFT~SP500, data=excessRet.ts)
avar.HC = vcov(ols.fit,correction="white")
summary(ols.fit,correction="white")

#
# Illustrate long-horizon regressions with standard errors computed using the
# NW correction. Use annual data in shiller.annual
#

colIds(shiller.annual)
# compute log of real data
ln.p = log(shiller.annual[,"real.price"])
colIds(ln.p) = "ln.p"
ln.d = log(shiller.annual[,"real.dividend"])
colIds(ln.d) = "ln.d"
ln.dpratio = ln.d - ln.p
colIds(ln.dpratio) = "ln.dpratio"
# compute cc real total returns - see CLM pg. 261
ln.r = diff(ln.p) + log(1+exp(ln.dpratio[-1,]))
colIds(ln.r) = "ln.r"
# create multi-year cc total returns
ln.r.2 = aggregateSeries(ln.r,moving=2,FUN=sum)
colIds(ln.r.2) = "ln.r.2"
ln.r.5 = aggregateSeries(ln.r,moving=5,FUN=sum)
colIds(ln.r.5) = "ln.r.5"
ln.r.10 = aggregateSeries(ln.r,moving=10,FUN=sum)
colIds(ln.r.10) = "ln.r.10"
stockdiv.ts = seriesMerge(ln.p,ln.d,ln.dpratio,
ln.r,ln.r.2,ln.r.5,ln.r.10,pos="union")

stockdiv.ts@title = "Log stock price and dividend data"
stockdiv.ts@documentation = 
c("Natural log of real stock prices and real dividends",
"Stock prices are on the Standard & Poor's composite index",
"source: Robert Shiller's homepage,",
"http://aida.econ.yale.edu/~shiller/data.htm")


smpl = (positions(stockdiv.ts) >= timeDate("1/1/1947")) &
(positions(stockdiv.ts) <= timeDate("1/1/1997"))

# regression over 1 year horizon
ols.fit = OLS(ln.r~tslag(ln.dpratio),data=stockdiv.ts,
start="1947",end="1995",in.format="%Y", na.rm=T)
ols.fit2 = OLS(ln.r.2~tslag(ln.dpratio),data=stockdiv.ts,
start="1947",end="1995",in.format="%Y", na.rm=T)
ols.fit5 = OLS(ln.r.5~tslag(ln.dpratio),data=stockdiv.ts,
start="1947",end="1995",in.format="%Y", na.rm=T)
ols.fit10 = OLS(ln.r.10~tslag(ln.dpratio),data=stockdiv.ts,
start="1947",end="1995",in.format="%Y", na.rm=T)

summary(ols.fit,correction="nw")
summary(ols.fit2,correction="nw")
summary(ols.fit5,correction="nw")
summary(ols.fit10,correction="nw")

summary(ols.fit10,correction="nw",
bandwidth=9,window="rectangular")

#
# testing for heteroskedasticity
#

heteroTest(ols.fit,method="white")

#
# Recursive least squares
#
args(RLS)
rls.fit = RLS(MSFT~SP500,data=excessRet.ts)
class(rls.fit)
names(rls.fit)
rls.fit
plot(rls.fit, ask=F)

cusumTest(rls.fit)