#Data is from Greene, Table F2.2
library(arm)
gasdata = read.csv("http://home.cc.umanitoba.ca/~godwinrt/7010/gas.csv")
attach(gasdata)

#View the break-point:
plot(YEAR,GAS)
lines(YEAR,GAS)
text(1973,1.4,"1973")
arrows(1973,1.35,1973,1.2,length = 0.1)

#Estimate the pooled model:
eq1 = lm(log(GASEXP/GASP/POP) ~ YEAR + log(INCOME/POP) + log(GASP) + log(PNC) + log(PUC))
#View the estimated coefficients:
coefplot(eq1,vertical=FALSE,var.las = 1,cex.var=1.2)
#Get the sum of squared residuals from the pooled (restricted) model:
sser = sum(eq1$residuals^2)

#Use only the first 21 observations (up to 1973):
preshock = gasdata[1:21,]
attach(preshock)
eq2 = lm(log(GASEXP/GASP/POP) ~ YEAR + log(INCOME/POP) + log(GASP) + log(PNC) + log(PUC))
coefplot(eq2,vertical=FALSE,var.las = 1,cex.var=1.2)
sseu1 = sum(eq2$residuals^2)

#Use only the last 31 observations (after 1973):
postshock = gasdata[22:52,]
attach(postshock)
eq3 = lm(log(GASEXP/GASP/POP) ~ YEAR + log(INCOME/POP) + log(GASP) + log(PNC) + log(PUC))
coefplot(eq3,vertical=FALSE,var.las = 1,cex.var=1.2)
sseu2 = sum(eq3$residuals^2)

#Calculate Chow test statistic:
chow = ((sser - sseu1 - sseu2)/6)/((sseu1 + sseu2)/(52 - 12))
#p-value:
1 - pf(chow,6,40)

#Estimate the model with dummy variables:
DUM = c(rep(0,21),rep(1,31))
attach(gasdata)
eq4 = lm(log(GASEXP/GASP/POP) ~ YEAR + log(INCOME/POP) + log(GASP) + log(PNC) + log(PUC) + DUM + DUM*YEAR + DUM*log(INCOME/POP) + DUM*log(GASP) + DUM*log(PNC) + DUM*log(PUC))
summary(eq4)
ssedum = sum(eq4$residuals^2)