## Read in dataset (from Brooks): the following works irrespective
## of a) or b)

## a) estimate from raw data
apt <- read.table(file="APT.txt", h=T, sep="\t")

## define model formula (an APT model)
fm <- I(diff(log(Microsoft))[-1]*100 - (USTB3M/12)[-(1:2)]) ~
       I(diff(log(SANDP))[-1]*100 - (USTB3M/12)[-(1:2)]) +
       diff(Industrial.production)[-1] + diff(CONSUMER.CREDIT)[-1] +
       diff(diff(log(CPI))*100) + diff(M1MONEY.SUPPLY)[-1] +
       diff(BAA.AAA.SPREAD)[-1] + diff(USTB10Y-USTB3M)[-1]

## define model formula for restricted model (CAPM)
fm0 <- I(diff(log(Microsoft))[-1]*100 - (USTB3M/12)[-(1:2)]) ~
        I(diff(log(SANDP))[-1]*100 - (USTB3M/12)[-(1:2)])
#####

## b) ready-made dataset in R format (easier)
load("apt.rda")

## define model formula (an APT model)
fm <- ermsoft ~ ersandp + dprod + dcredit + dinflation + 
        dmoney + dspread + rterm

## define model formula for restricted model (CAPM)
fm0 <- ermsoft ~ ersandp 
#####


## estimate model by OLS
aptmod <- lm(fm, apt)

summary(aptmod)

## estimate restricted model
capmod <- lm(fm0, apt)

## evaluate joint restriction
library(lmtest)
waldtest(aptmod, capmod)

## do an F-test for the joint restriction APT--> CAPM "by hand"
myFtest <- (sum(resid(capmod)^2)-sum(resid(aptmod)^2)) / sum(resid(aptmod)^2) *
           (length(resid(aptmod))-length(coef(aptmod))) /
           (length(coef(aptmod))-length(coef(capmod)))

mypval <- 1 - pf(myFtest, df1=length(coef(aptmod))-length(coef(capmod)),
                 df2=length(resid(aptmod))-length(coef(aptmod)))

myFtest
mypval

## test for parameter stability:
## see if there are breaks in Mar-00 (peak of the Dot.com Bubble)
t.bar <- which(apt$Date == "Mar-00")

## restricted model is now aptmod; unrestricted model is now:
aptmod1 <- lm(fm, apt[1:t.bar, ])
aptmod2 <- lm(fm, apt[(t.bar+1):dim(apt)[[1]], ])

RRSS <- sum(resid(aptmod)^2)
URSS <- sum(resid(aptmod1)^2) + sum(resid(aptmod2)^2)
T. <- length(resid(aptmod))
k <- length(coef(aptmod))

myChowTest <- (RRSS-URSS) / URSS * (T.-2*k)/k
myChowTest

## critical value:
qf(0.95, df1=k, df2=(T.-2*k))

## p.value:
1-pf(myChowTest, df1=k, df2=(T.-2*k))


## a different crash:
## (backwards) predictive failure test w.r.t. Oct-87

t0 <- which(apt$Date == "Oct-87")

aptmod0 <- lm(fm, apt[(t0+1):dim(apt)[[1]], ])
RSS1 <- sum(resid(aptmod0)^2)
T1 <- length(resid(aptmod0))
T2 <- T. - T1

myPFtest <- (RRSS-RSS1) / RSS1 * (T1-k)/T2
myPFtest

## critical value:
qf(0.95, df1=T2, df2=(T1-k))

## p.value:
1-pf(myPFtest, df1=T2, df2=(T1-k))

## draw your conclusions: which of the two crashes denotes a
## structural breakpoint?