sábado, 12 de dezembro de 2009

Ainda sobre o aquecimento

Quem tiver interesse nos dados relacionados à temperatura mundial [usados em Koenker, R. and Schorfheide F. (1994) Quantile Spline Models for Global Temperature Change; Climate Change 28, 395–404], pode encontrá-los na "library" COBS do software R.

Caso você tenha o R instalado e a library disponínel, siga estes passos:

data(globtemp)

plot(globtemp, main = "Annual Global Temperature Deviations")


Para reproduzir uma das estimativas do texto do Koenker e Schorfheide:

a50 <- cobs(year, temp, knots.add = TRUE, degree = 1, constraint = "increase")
summary(a50)


a50 <- cobs(year, temp, nknots = 9, knots.add = TRUE, degree = 1,
constraint = "increase")
summary(a50)


a10 <- cobs(year, temp, nknots = length(a50$knots), knots = a50$knot,
degree = 1, tau = 0.1, constraint = "increase")
summary(a10)
a90 <- cobs(year, temp, nknots = length(a50$knots), knots = a50$knot,
degree = 1, tau = 0.9, constraint = "increase")
summary(a90)

which(hot.idx <- temp >= a90$fit)
which(cold.idx <- temp <= a10$fit)
normal.idx <- !hot.idx & !cold.idx

plot(year, temp, type = "n", ylab = "Temperature (C)", ylim = c(-.7,.6))
lines(predict(a50, year, interval = "both"), col = 2)
lines(predict(a10, year, interval = "both"), col = 3)
lines(predict(a90, year, interval = "both"), col = 3)
points(year, temp, pch = c(1,3)[2 - normal.idx])



text(year[hot.idx], temp[hot.idx] + .03, labels = year[hot.idx])
text(year[cold.idx],temp[cold.idx]- .03, labels = year[cold.idx])

Nenhum comentário: