New England climate analysis, addressing autocorrelation

ne_climate_2.R

# Attempting to replicate analysis in:
# http://wattsupwiththat.com/2010/06/29/waxman-malarkey-impact-zone-us-northeast/
# 7/2/2010: addressing autocorrelation

# Data import and cleanup
# Note that 2010 is excluded due to missing values
dat <- read.table("~/Desktop/drd964x.tmpst.txt", colClasses=c("character",rep("numeric", 12)))
colnames(dat) <- c("id", "jan", "feb", "mar", "apr", "may", "jun", "jul", "aug", "sep", "oct", "nov", "dec")
dat$year <- as.numeric(substr(dat$id, 7, 10))
dat <- dat[dat$year < 2010,]
dat$region <- substr(dat$id, 1, 3)
dat$annual <- apply(dat[,2:13], 1, mean)
dat.ne <- dat[dat$region == "101",]

# Ugliest possible way to average previous year's December
# and current year's January and February
dat.ne$last.dec <- (c(NA,dat.ne$dec))[1:115]
dat.ne$winter <- (dat.ne$jan + dat.ne$feb + dat.ne$last.dec) / 3

# Create time series for autocorrelation analysis
ts.annual <- ts(dat.ne$annual, start=c(1895, 1), frequency=1)
ts.winter <- ts(dat.ne$winter[-1], start=c(1896, 1), frequency=1)
ts.annual.1970 <- window(ts.annual, start=1970)
ts.winter.1970 <- window(ts.winter, start=1970)
par(mfrow=c(2,2))
acf(ts.annual, ci.type="ma")
pacf(ts.annual)
acf(ts.winter, ci.type="ma")
pacf(ts.winter)

# Create 8 year lagged series suitable for regression
ts.annual.8 <- c(rep(NA,8),ts.annual)[1:115]
ts.annual.1970.8 <- c(rep(NA,8),ts.annual.1970)[1:40]

# Regression for Northeast region, annual mean
print(summary(lm(ts.annual ~ time(ts.annual))))
print(summary(lm(ts.annual ~ time(ts.annual) + ts.annual.8)))

# Regression for Northeast region, annual mean, 1970--2009
print(summary(lm(ts.annual.1970 ~ time(ts.annual.1970))))
print(summary(lm(ts.annual.1970 ~ time(ts.annual.1970) + ts.annual.1970.8)))

# Regression for Northeast region, Winter mean
print(summary(lm(ts.winter ~ time(ts.winter))))

# Regression for Northeast region, Winter mean, 1970--2009
print(summary(lm(ts.winter.1970 ~ time(ts.winter.1970))))