# Change-point analysis of HadCRUT3 annual time series data

require(strucchange)

# Data import
dat.raw <- read.table("~/Downloads/annual")
dat <- data.frame(temp=ts(dat.raw[,2], start=1850))
dat <- window(cbind(temp=dat$temp, temp_lag1=lag(dat$temp, k = -1)), start=1851, end=2010)

# Structural change model

# Fit model
# Minimal segment size set to 5 years; overfit will be handled later
sc.alan <- breakpoints(temp ~ 1, h = 10, data=dat)
sc <- breakpoints(temp ~ time(temp), h = 5, data=dat)
sc.lag <- breakpoints(temp ~ time(temp) + temp_lag1 * time(temp), h = 5, data=dat)
plot(summary(sc.alan))
plot(summary(sc))
plot(summary(sc.lag))

# Alan's model
# By inspecting the BIC, we select a model with 7 breaks
print(breakpoints(sc.alan, breaks=7))
sc.alan <- lm(temp ~ breakfactor(sc.alan, breaks = 7), data=dat)
plot(dat[,1], ylab="temperature", main="Structural change model: no lag, 7 breakpoints")
lines(ts(fitted(sc.alan), start=1850), col=4)
lines(confint(sc.alan))

# No red noise correction
# By inspecting the BIC, we select a model with 4 breaks
print(breakpoints(sc, breaks = 4))
sc.4 <- lm(temp ~ time(temp) + breakfactor(sc, breaks = 4) * time(temp), data=dat)
plot(dat[,1],ylab="temperature",main="Structural change model: no lag, 4 breakpoints")
lines(ts(fitted(sc.4), start=1850), col=4)
lines(confint(sc))

# Red noise correction
# By inspecting the BIC, we select a model with 1 break
print(breakpoints(sc.lag, breaks=1))
sc.lag.1 <- lm(temp ~ time(temp) + breakfactor(sc.lag, breaks = 1) * time(temp), data=dat)
plot(dat[,1],ylab="temperature",main="Structural change model: 1 year lag, 1 breakpoint")
lines(ts(fitted(sc.lag.1), start=1851), col=4)
lines(confint(sc.lag))

# Compare models
print(anova(sc.alan, sc.4, sc.lag.1))