Skip to content

Instantly share code, notes, and snippets.

Created July 3, 2010 16:58
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save anonymous/462697 to your computer and use it in GitHub Desktop.
Save anonymous/462697 to your computer and use it in GitHub Desktop.
New England climate analysis, further addressing autocorrelation
# Attempting to replicate analysis in:
# http://wattsupwiththat.com/2010/06/29/waxman-malarkey-impact-zone-us-northeast/
# 7/3/2010: Nychka approach
# 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)
# Fit linear trends for full time series
fm <- function(ts) {
lm(ts ~ time(ts))
}
fm.annual <- fm(ts.annual)
fm.annual.1970 <- fm(ts.annual.1970)
fm.winter <- fm(ts.winter)
fm.winter.1970 <- fm(ts.winter.1970)
# Plot residuals
par(mfrow=c(2,2))
plot(residuals(fm.annual) ~ time(ts.annual),xlim=c(1895,2009))
plot(residuals(fm.annual.1970) ~ time(ts.annual.1970))
plot(residuals(fm.winter) ~ time(ts.winter),xlim=c(1895,2009))
plot(residuals(fm.winter.1970) ~ time(ts.winter.1970))
### Nychka analysis
nychka <- function(fm) {
output <- list()
# Calculate sample lag-1 autocorrelation
resids <- residuals(fm)
n <- length(resids)
rhohat <- cor(resids[2:n], resids[1:(n-1)])
output$n <- n
output$rhohat <- rhohat
# See Nychka Eq. 5
neff <- n * (1 - rhohat - 0.68 / sqrt(n)) / (1 + rhohat + 0.68 / sqrt(n))
output$neff <- neff
# Calculate scaling for SE, see Nychka Eq. 3, 6
sescale <- sqrt((n-2)/(neff-2))
output$sescale <- sescale
# Calculate adjusted SE for warming trend
seadj <- sescale * sqrt(vcov(fm)[2,2])
output$seadj <- seadj
# Calculate adjusted p-value
pval <- (1 - pt(coefficients(fm)[2] / seadj, neff)) * 2
output$pval <- pval
output
}
print(nychka(fm.annual))
print(nychka(fm.annual.1970))
print(nychka(fm.winter))
print(nychka(fm.winter.1970))
### MLE analysis
require(nlme)
mle.ar1 <- function(ts) {
fm.trend <- gls(ts ~ time(ts), corr=corAR1(), method="ML")
fm.notrend <- gls(ts ~ 1, corr=corAR1(), method="ML")
print(anova(fm.trend, fm.notrend))
}
mle.ar1(fm.annual)
mle.ar1(fm.annual.1970)
mle.ar1(fm.winter)
mle.ar1(fm.winter.1970)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment