## ----setup, include=FALSE, cache=FALSE------------------------------ # set global chunk options ## this code is only required to set up knitr for reproducible documentation ## it has nothing to do with the exercise itself, please ignore! opts_knit$set(aliases=c(h = 'fig.height', w = 'fig.width'), out.format='latex', use.highlight=TRUE) opts_chunk$set(fig.path='kgraph/exTSA-', fig.align='center', fig.show='hold', prompt=FALSE, fig.width=6, fig.height=6, out.width='.7\\linewidth', size='scriptsize') ## comment=NA, ## this does not comment out the output # options to be read by formatR options(replace.assign=TRUE,width=70) ## ----echo=FALSE, results='hide'------------------------------------- rm(list=ls()) ## ----child-data, child='ts_1.Rnw'----------------------------------- ## ----ts1-set-parent, echo=FALSE, cache=FALSE------------------------ set_parent('exTSA.Rnw') opts_chunk$set(fig.path='kgraph/ts1-') ## ----eval=FALSE----------------------------------------------------- ## file.show("./ds_tsa/anatolia_hati.txt") ## ----scan-hati------------------------------------------------------ gw <- scan("./ds_tsa/anatolia_hati.txt") str(gw) ## ------------------------------------------------------------------- gw <- ts(gw, start=1975, frequency=12) str(gw) attributes(gw) start(gw) end(gw) ## ----eval=FALSE----------------------------------------------------- ## print(gw) ## ----eval=FALSE----------------------------------------------------- ## time(gw) ## ----eval=FALSE----------------------------------------------------- ## cycle(gw) ## ------------------------------------------------------------------- frequency(gw) deltat(gw) ## ----fig.width=15,fig.height=5,out.width='\\textwidth'-------------- par(mfrow=c(1,3)) # pdf("AnatoliaWell1.pdf", width=10, height=5) plot(gw, ylab="Depth to water table (m)", main="Anatolia Well 1") # dev.off() plot(gw, type="o", pch=20, cex=0.6, col="blue", ylab="Depth to water table (m)", main="Anatolia Well 1") plot(gw, type="h", col="blue", ylab="Depth to water table (m)", main="Anatolia Well 1") par(mfrow=c(1,1)) ## ------------------------------------------------------------------- window(gw, start=c(1990,4), end=c(1993,3)) plot(window(gw, start=c(1990,4), end=c(1993,3)), type="o", ylab="Depth to water table (m)", main="Anatolia Well 1") ## ------------------------------------------------------------------- window(gw, start=c(1990,1), end=c(1990,1)) ## ------------------------------------------------------------------- window(gw, 1990, c(1992,12)) diff(window(gw, 1990, c(1992,12))) diff(window(gw, 1990, c(1992,12)), lag=2) diff(window(gw, 1990, c(1992,12)), lag=12) ## ------------------------------------------------------------------- diff(window(gw, 1990, c(1993,12)), lag=12, differences=1) diff(window(gw, 1990, c(1993,12)), lag=12, differences=2) diff(window(gw, 1990, c(1993,12)), lag=12, differences=3) ## ------------------------------------------------------------------- plot(diff(gw), ylab="One-month differences in groundwater depth (m)", main="Anatolia well 1") ## ------------------------------------------------------------------- i <- which.min(diff(gw)) diff(gw)[i] time(gw)[i] cycle(gw)[i] i <- which.max(diff(gw)) diff(gw)[i] time(gw)[i] cycle(gw)[i] i <- which(abs(diff(gw)) > 2) diff(gw)[i] time(gw)[i] cycle(gw)[i] ## ------------------------------------------------------------------- plot(window(gw, start=1986, end=1990), ylab="Groundwater depth (m)", main="Anatolia well 1", type="h", col="blue") lines(window(gw, start=1986, end=1990), lty=2) ## ------------------------------------------------------------------- window(gw, start=c(1987,3), end=c(1988,6)) ## ----eval=FALSE----------------------------------------------------- ## file.show("./ds_tsa/Tana_Daily.csv") ## ------------------------------------------------------------------- tana <- read.csv("./ds_tsa/Tana_Daily.csv", skip=1, header=T, colClasses=c(rep("integer",2), rep("character",12)), blank.lines.skip=T,na.strings=c("N.A","NA"," ")) str(tana) tana[1:35,] ## ------------------------------------------------------------------- sum(is.na(tana[,3:14])) sum(tana[,3:14]=="",na.rm=TRUE) ## ----show-unique-all------------------------------------------------ head(sort(unique(stack(tana[,3:14])$values))) tail(sort(unique(stack(tana[,3:14])$values))) ## ------------------------------------------------------------------- sum(tana[,3:14]=="TR", na.rm=TRUE) sum(tana[,3:14]=="tr", na.rm=TRUE) sum(tana[,3:14]=="0.01", na.rm=TRUE) sum(tana[,3:14]=="0", na.rm=TRUE) ## ------------------------------------------------------------------- require(car) for (i in 3:14) { tana[,i] <- recode(tana[,i], "c('TR','tr','0.01')='0'") } ## ----show-sorted-unique--------------------------------------------- head(sort(unique(stack(tana[,3:14])$values)),12) tail(sort(unique(stack(tana[,3:14])$values)),12) ## ------------------------------------------------------------------- sum(tana[,3:14]=="TR", na.rm=TRUE) sum(tana[,3:14]=="tr", na.rm=TRUE) sum(tana[,3:14]=="0.01", na.rm=TRUE) sum(tana[,3:14]=="0", na.rm=TRUE) ## ------------------------------------------------------------------- tana[tana$DATE==29,"FEB"] ## ------------------------------------------------------------------- tana.ppt <- NULL; month.days <- c(0,0,31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31) for (yr.first.row in seq(from=1, by=32, length=(2006 - 1981 + 1))) { for (month.col in 3:14) { tana.ppt <- c(tana.ppt, tana[yr.first.row:(yr.first.row + month.days[month.col]-1), month.col]) } }; str(tana.ppt) rm(month.days, yr.first.row, month.col) ## ------------------------------------------------------------------- length(tana.ppt)/365 ## ------------------------------------------------------------------- tana.ppt <- ts(tana.ppt, start=1981, frequency=365) str(tana.ppt) ## ----fig.width=12,fig.height=4, out.width='\\textwidth'------------- plot(tana.ppt, main="Lake Tana rainfall", ylab="mm") abline(h=100, col="gray") points(xy.coords(x=time(tana.ppt), y=100, recycle=T), pch=ifelse(is.na(tana.ppt),"l",""), col="red") grid() ## ------------------------------------------------------------------- yrs <- c(1982, 1983, 1988, 1989,1991,1998,1999, 2006); ymax <- 0 for (i in yrs) { ymax <- as.numeric(max(ymax, window(tana.ppt, start=i, end=i+1, extend=TRUE), na.rm=T)) } (ymax <- ceiling(ymax)) ## ----fig.height=11, fig.width=9, out.width='\\linewidth'------------ par(mfrow=c(4,2)) for (i in yrs) { plot(window(tana.ppt, start=i, end=i+1, extend=TRUE), type="h", ylab="mm", ylim=c(0,ymax)); title(main=paste("Lake Tana rainfall", i), sub=paste("Annual total:", sum(as.numeric(window(tana.ppt, start=i, end=i+1)), na.rm=T))) abline(h=ymax, col="gray") points(xy.coords(x=time(window(tana.ppt, start=i, end=i+1, extend=TRUE)), y=ymax, recycle=T), pch=ifelse(is.na(window(tana.ppt, start=i, end=i+1, extend=TRUE)),"l",""), col="red") grid() } par(mfrow=c(1,1)) ## ----fig.width=12,fig.height=6-------------------------------------- str(tana.ppt) sum(is.na(tana.ppt)) tana.ppt.c <- na.contiguous(tana.ppt) str(tana.ppt.c) frequency(tana.ppt.c) head(time(tana.ppt.c)) head(cycle(tana.ppt.c)) tail(time(tana.ppt.c)) tail(cycle(tana.ppt.c)) sum(is.na(tana.ppt.c)) plot(tana.ppt.c, main="Lake Tana rainfall", ylab="mm", sub="continuous record") ## ----child-data, child='ts_2.Rnw'----------------------------------- ## ----ts2-set-parent, echo=FALSE, cache=FALSE------------------------ set_parent('exTSA.Rnw') opts_chunk$set(fig.path='kgraph/ts2-') ## ------------------------------------------------------------------- summary(gw) ## ------------------------------------------------------------------- gw.f <- data.frame(gw, year=as.numeric(floor(time(gw))), cycle=as.numeric(cycle(gw)), time=time(gw)) str(gw.f) ## ------------------------------------------------------------------- head(by(gw.f$gw, IND=gw.f$year, FUN=summary)) ## ------------------------------------------------------------------- by(gw.f$gw, IND=gw.f$year, FUN=summary)[["1978"]] by(gw.f$gw, IND=gw.f$year, FUN=max)[["1978"]] by(gw.f$gw, IND=gw.f$year, FUN=min)[["1978"]] by(gw.f$gw, IND=gw.f$year, FUN=quantile, probs=.9)[["1978"]] ## ----fig.width=12,fig.height=6, out.width='\\textwidth'------------- boxplot(gw.f$gw ~ gw.f$year, main="Anatolia well 1", ylab="groundwater level (m)") grid() ## ------------------------------------------------------------------- (ann.mean <- aggregate(gw, nfrequency=1, FUN=mean)) time(ann.mean) ## ------------------------------------------------------------------- gw.f$in.yr <- as.numeric(gw - ann.mean[match(gw.f$year, time(ann.mean))]) str(gw.f) ## ----fig.width=12, fig.height=8, out.width='0.8\\textwidth'--------- boxplot(gw.f$in.yr ~ gw.f$cycle, xlab="Month", ylab="Deviation from annual mean (m)", main="Anatolia groundwater well 1") ## ----fig.width=12,fig.height=6, out.width='\\textwidth'------------- gw.q <- aggregate(window(gw, 1980, 1986), nfrequency=4, FUN=mean) str(gw.q) par(mfrow=c(1,2)) plot(gw.q, ylab="Depth to water table (m)", main="Anatolia Well 1, quarterly", type="b") plot(window(gw, 1980, 1986), ylab="Depth to water table (m)", main="Anatolia Well 1, monthly", type="b") par(mfrow=c(1,1)) ## ----out.width='0.75\\textwidth'------------------------------------ k <- c(.25,.5,.75,1,.75,.5,.25) (k <- k/sum(k)) fgw <- filter(gw, sides=2, k) plot.ts(gw, main="Anatolia Well 1, 7-month filter", ylab="groundwater depth (m)") lines(fgw, col="red") ## ----out.width='0.75\\textwidth'------------------------------------ fgw.2 <- filter(gw, sides=2, rep(1,7)/7) fgw.3 <- filter(gw, sides=2, rep(1,12)/12) plot.ts(gw, xlim=c(1990,1995), ylim=c(37,48), ylab="groundwater depth (m)") title(main="Anatolia Well 1, three filters") lines(fgw, col="blue", lty=2) lines(fgw.2, col="red") lines(fgw.3, col="magenta") text(1995,40,"1/4,1/2,1 filter", col="blue", pos=2) text(1995,39,"1,1,1 filter", col="red",pos=2) text(1995,38,"annual filter", col="magenta",pos=2) ## ----out.width='0.75\\textwidth'------------------------------------ plot.ts(gw, main="Anatolia Well 1, annual filter", ylab="groundwater depth (m)") lines(fgw.3, col="magenta", lwd=1.5) ## ----out.width='0.8\\textwidth'------------------------------------- plot.ts(gw, main="Anatolia well 1 with smoothers", ylab="groundwater depth (m)") lines(lowess(gw), col="blue") lines(lowess(gw, f=1), col="green") lines(lowess(gw, f=1/3), col="red") lines(lowess(gw, f=1/10), col="purple") text(1990, 36, "Smoothing parameter: 2/3 (default)", col="blue", pos=4) text(1990, 34, "Smoothing parameter: 1", col="green", pos=4) text(1990, 32, "Smoothing parameter: 1/3", col="red", pos=4) text(1990, 30, "Smoothing parameter: 1/10", col="purple", pos=4) ## ------------------------------------------------------------------- frequency(gw) ## ------------------------------------------------------------------- tapply(gw, cycle(gw), mean) ## ------------------------------------------------------------------- head(gw, 2*frequency(gw)) head(gw-rep(tapply(gw, cycle(gw), mean), length(gw)/frequency(gw)), 2*frequency(gw)) ## ----fig.width=10, fig.height=10, out.width='\\textwidth'----------- par(mfrow=c(2,1)) plot(gw, ylab="depth to groundwater", main="Original series") plot(gw-rep(tapply(gw, cycle(gw), mean), length(gw)/frequency(gw)), ylab="difference from cycle mean", main="Seasonally-corrected series") abline(h=0, lty=2) par(mfrow=c(2,1)) ## ------------------------------------------------------------------- gw.stl <- stl(gw, s.window="periodic") str(gw.stl) tmp <- stl(gw, s.window="periodic", t.window=19) unique(tmp$time.series[,"trend"]-gw.stl$time.series[,"trend"]) rm(tmp) ## ----out.width='\\textwidth'---------------------------------------- plot(gw.stl) ## ----fig.width=10,fig.height=4, out.width='\\textwidth'------------- plot(gw.stl$time.series[,"trend"], main="Anatolia well 1, trend", ylab="Groundwater level (m)") ## ----fig.width=10,fig.height=4, out.width='\\textwidth'------------- ts.plot(gw.stl$time.series, col=c("black","blue","red"), main="Anatolia well 1, decomposition", ylab="Groundwater level (m)") tmp <- attributes(gw.stl$time.series)$dimnames[[2]] for (i in 1:3) { text(1995, 24-(i*4), tmp[i], col=c("black","blue","red")[i], pos=4) } grid() ## ----out.width='\\textwidth'---------------------------------------- gw.stl <- stl(gw, s.window=2*frequency(gw)+1) plot(gw.stl) ## ----out.width='\\textwidth'---------------------------------------- gw.stl <- stl(gw, s.window=25, t.window=85) plot(gw.stl) ## ----fig.width=12,fig.height=12, out.width='\\textwidth'------------ lag.plot(gw, lags=12, main="Anatolia well 1", cex=0.3, pch=20, lty=1) ## ------------------------------------------------------------------- window(gw, 2000, 2001) lag(window(gw, 2000, 2001), 1) lag(window(gw, 2000, 2001), 2) lag(window(gw, 2000, 2001), -1) lag(window(gw, 2000, 2001), -2) ## ----fig.width=12,fig.height=12, out.width='\\textwidth'------------ lag.plot(window(gw, start=1990, end=1993), lags=12, main="Anatolia well 1") ## ------------------------------------------------------------------- print(acf(gw, plot=F)) acf(gw, main="Autocorrelation, groundwater levels, Anatolia well 1") ## ------------------------------------------------------------------- acf(gw, lag.max=60, main="Autocorrelation, groundwater levels, Anatolia well 1") ## ----plot-gw-stl-remainders, fig.width=8, fig.height=5, out.width='.75\\linewidth'---- gw.stl <- stl(gw, s.window=2*frequency(gw), t.window=84) gw.r <- gw.stl$time.series[,"remainder"] plot(gw.r, ylab="Remainder (m)") title(main="Anatolia well 1 level, remainder from trend & cycle") abline(h=0, lty=2, col="red") ## ----plot-gw-stl-remainders-acf------------------------------------- print(acf(gw.r, plot=F)) acf(gw.r, main="Autocorrelation of remainders, Anatolia well 1") ## ----plot-gw-stl-pacf----------------------------------------------- print(pacf(gw,plot=F)) pacf(gw, main="Partial autocorrelation, Anatolia well 1") ## ----plot-gw-stl-remainders-pacf------------------------------------ print(pacf(gw.r,plot=F)) pacf(gw.r, main="Parial autocorrelation of remainders") ## ------------------------------------------------------------------- s <- spectrum(gw, spans=c(5,7)); grid() str(s) head(s$spec, n=12) ## ------------------------------------------------------------------- frequency(gw) length(gw)/frequency(gw) head(s$freq,n=length(gw)/frequency(gw)) ## ------------------------------------------------------------------- ss <- sort(s$spec, decreasing=T, index.return=T) str(ss) hi <- ss$x>.15 which(hi) ss$x[hi] ss$ix[hi] sort(s$freq[ss$ix[hi]]) ## ----fig.width=12, fig.height=4, out.width='\\textwidth'------------ par(mfrow=c(1,3)) spectrum(gw, spans=c(5,7), log="no", main="Anatolia groundwater level", sub="annual cycle, smoothed") grid() spectrum(gw, spans=c(5,7), log="yes", main="Anatolia groundwater level", sub="annual cycle, smoothed") grid() spectrum(gw, spans=c(5,7), log="dB", main="Anatolia groundwater level", sub="annual cycle, smoothed") grid() par(mfrow=c(1,1)) ## ------------------------------------------------------------------- sp.gw <- spectrum(gw.r, span=c(5,7), log="dB", main="Periodogram, residuals", sub="annual cycle, smoothed") grid() ## ----fig.width=12, fig.height=6, out.width='\\textwidth'------------ par(mfrow=c(1,2)) spectrum(gw.r, span=c(5,7), log="no", main="Periodogram, residuals", sub="annual cycle, smoothed", xlim=c(0,2), type="h") grid() s <- spectrum(gw.r, span=c(5,7), log="dB", main="Periodogram, residuals", sub="annual cycle, smoothed", xlim=c(0,2)) grid() sp.gw$freq[which.max(s$spec)] frequency(gw)/(s$freq[which.max(s$spec)]) which.max(s$spec[16:30]) sp.gw$freq[which.max(s$spec[16:30])+15] frequency(gw)/(s$freq[which.max(s$spec[16:30])+15]) par(mfrow=c(1,1)) ## ----child-data, child='ts_3.Rnw'----------------------------------- ## ----ts3-set-parent, echo=FALSE, cache=FALSE------------------------ set_parent('exTSA.Rnw') opts_chunk$set(fig.path='kgraph/ts3-') ## ----out.width='0.9\\textwidth'------------------------------------- plot(stl(gw, s.window="periodic")$time.series, main="Well 1, periodic decomposition") ## ----out.width='0.9\\textwidth'------------------------------------- plot(stl(gw, s.window=25, t.window=85)$time.series, main="Anatolia well 1, s=25, t=85 decomposition") ## ------------------------------------------------------------------- m.gw <- lm(gw ~ time, data=gw.f) summary(m.gw) ## ----plot-gw-ols---------------------------------------------------- plot(gw.f$time, gw.f$gw, type="l", col="darkgreen", ylab="Groundwater level (m)", xlab="Year") title("Anatolia well 1, linear trend") lines(x=as.numeric(gw.f$time), y=fitted(m.gw), col="red") ## ------------------------------------------------------------------- summary(m.gw)$adj.r.squared coefficients(m.gw)["time"] ## ------------------------------------------------------------------- plot(m.gw, which=1) ## ----fit-reml------------------------------------------------------- library(nlme) (cor.value <- acf(gw, plot=FALSE)$acf[2]) summary(m.gw.gls <- gls(model=gw ~ time, data=gw.f, correlation=corAR1(value=cor.value, form = ~1))) coef(m.gw) coef(m.gw.gls) ## ----plot-gw-gls---------------------------------------------------- plot(gw.f$time, gw.f$gw, type="l", col="darkgreen", ylab="Groundwater level (m)", xlab="Year") title("Anatolia well 1, linear trend") abline(m.gw.gls, col="darkgreen") lines(x=as.numeric(gw.f$time), y=fitted(m.gw), col="red") legend("topleft", c("GLS","OLS"), lty=1, col=c("darkgreen","red")) ## ------------------------------------------------------------------- m.gw.3 <- lm(gw ~ I(time^3) + time, data=gw.f) summary.lm(m.gw.3) anova(m.gw.3, m.gw) ## ------------------------------------------------------------------- plot(gw.f$time, gw.f$gw, type="l", col="darkgreen", ylab="Groundwater level (m)", xlab="Year") title("Anatolia well 1, cubic trend") lines(x=as.numeric(gw.f$time), y=fitted(m.gw.3), col="red") ## ------------------------------------------------------------------- gw.f.78 <- subset(gw.f, gw.f$year > 1977) summary(m.gw.78 <- lm(gw ~ time, data=gw.f.78)) (cor.value <- acf(gw.f.78, plot=FALSE)$acf[2]) summary(m.gw.78.gls <- gls(model=gw ~ time, data=gw.f.78, correlation=corAR1(value=cor.value, form = ~1))) ## ----plot-gw-f-78, fig.width=8, fig.height=5, out.width='.75\\linewidth'---- plot(gw.f.78$time, gw.f.78$gw, type="l", col="darkgreen", ylab="Groundwater level (m)", xlab="Year") title("Anatolia well 1, linear trend since 1978") lines(x=as.numeric(gw.f.78$time), y=fitted(m.gw.78.gls), col="red") abline(m.gw.78, col="red") abline(m.gw.78.gls, col="darkgreen") legend("topleft", c("GLS","OLS"), lty=1, col=c("darkgreen","red")) ## ------------------------------------------------------------------- coef(m.gw.78) coef(m.gw.78.gls) ## ------------------------------------------------------------------- coefficients(m.gw.gls)["time"] coefficients(m.gw.78.gls)["time"] ## ------------------------------------------------------------------- predict(m.gw.78, data.frame(time=2010), interval="prediction", level=0.9) ## ------------------------------------------------------------------- gw.2050 <- predict.lm(m.gw.78, data.frame(time=2005:2050), interval="prediction", level=0.95) str(gw.2050) plot(gw.f$time, gw.f$gw, type="l", col="darkgreen", ylab="Groundwater level (m)", xlab="Year", xlim=c(1978,2050), ylim=c(floor(min(gw.f$gw)),ceiling(max(gw.2050[,"upr"])) )) title("Anatolia well 1, linear trend since 1978") grid() abline(v=2005, lty=2) lines(x=as.numeric(gw.f$time[gw.f$year > 1977]), y=fitted(m.gw.78), col="blue") lines(2005:2050, gw.2050[,"fit"]) lines(2005:2050, gw.2050[,"upr"], col="red", lty=2) lines(2005:2050, gw.2050[,"lwr"], col="red", lty=2) text(1990,60,"fit", col="blue") text(2030,60,"extrapolation") ## ------------------------------------------------------------------- gw.stl <- stl(gw, s.window=2*frequency(gw)+1) gw.f$nonseas <- gw.stl$time.series[,"trend"] + gw.stl$time.series[,"remainder"] str(gw.f) m.gw.nonseas <- lm(nonseas ~ time, data=subset(gw.f, gw.f$year > 1977)) summary(m.gw.nonseas) ## ----fig.height=8, fig.width=8-------------------------------------- plot(gw.f$time, gw.f$gw, type="l", col="darkgreen", ylab="Groundwater level (m)", xlab="Year") title("Anatolia well 1, linear trend since 1978, non-seasonal component") lines(x=as.numeric(gw.f$time), y=gw.f$nonseas, col="blue") lines(x=as.numeric(gw.f$time[gw.f$year > 1977]), y=fitted(m.gw.nonseas), col="red") text(1995, 38, col="darkgreen", pos=4, "time series") text(1995, 36, col="blue", pos=4, "non-seasonal component") text(1995, 34, col="red", pos=4, "linear trend") ## ------------------------------------------------------------------- (tmp <- summary(m.gw.nonseas)$adj.r.squared - summary(m.gw.78)$adj.r.squared) tmp/summary(m.gw.78)$adj.r.squared coefficients(m.gw.nonseas)["time"] coefficients(m.gw.78)["time"] ## ----fig.width=10,fig.height=5, out.width='\\textwidth'------------- par(mfrow=c(1,2)) plot(m.gw.nonseas, which=1:2) par(mfrow=c(1,1)) ## ------------------------------------------------------------------- require(Kendall) gw.1978 <- window(gw, start=c(1978,1), end=c(2005,12), extend=TRUE) SeasonalMannKendall(gw.1978) gw.nonseas.1978 <- ts(subset(gw.f, gw.f$year > 1977)$nonseas) MannKendall(gw.nonseas.1978) ## ------------------------------------------------------------------- MKslope <- function(ts) { f <- frequency(ts) n <- length(ts)/f d <- NULL for (j in n:2) for (k in (j-1):1) for (i in 1:f) d <- c(d, (ts[i + (j-1)*f] - ts[i + (k-1)*f])/(j-k)); hist(d, main="individual slope estimates", xlab="slope") print(summary(d)) return(median(na.omit(d))) } ## ------------------------------------------------------------------- print(MKslope(gw.1978)) coefficients(m.gw.nonseas)["time"] ## ------------------------------------------------------------------- data(GuelphP) plot(GuelphP, type="b", ylab="P concentration, mg l-1", main="Speed River water quality") ## ------------------------------------------------------------------- SeasonalMannKendall(GuelphP) print(guelphP.b <- MKslope(GuelphP)) ## ------------------------------------------------------------------- mean(gw.r) ## ------------------------------------------------------------------- lag.plot(gw.r, lags=1, main="Anatolia well 1 remainders, lag 1 scatterplot") ## ------------------------------------------------------------------- gw.r.0 <- gw.r[2:(length(gw.r)-1)] - mean(gw.r) gw.r.1 <- lag(gw.r,1)[1:(length(gw.r)-2)] - mean(gw.r) ## ------------------------------------------------------------------- plot(gw.r.0 ~ gw.r.1, xlab="Lag 1",ylab="Original") title(main="Anatolia well 1, remainders, series vs. lag-1 series") m.lag1 <- lm(gw.r.0 ~ gw.r.1) summary(m.lag1) abline(m.lag1) (alpha.1 <- cor(gw.r.0,gw.r.1)) ## ------------------------------------------------------------------- cor(gw[2:(length(gw)-1)] - mean(gw), lag(gw,1)[1:(length(gw)-2)] - mean(gw)) ## ------------------------------------------------------------------- acf(gw.r, lag.max=1, plot=F)$acf[2,,] ## ------------------------------------------------------------------- var(gw.r) (var.ar.1 <- (length(gw.r)-1)/(length(gw.r)-2) * (1 - alpha.1^2) * var(gw.r)) ## ----fig.width=12, fig.height=6, out.width='\\textwidth'------------ par(mfrow=c(1,2)) pacf(gw) pacf(gw.r) par(mfrow=c(1,1)) pacf(gw, lag.max=5, plot=F)$acf pacf(gw.r, lag.max=5, plot=F)$acf ## ------------------------------------------------------------------- (ar.gw.r <- ar(gw.r, method="yule-walker")) ar.gw.r$var.pred ar.gw.r$ar ## ------------------------------------------------------------------- (ar.gw.r.1 <- ar(gw.r, order.max=1)) ar.gw.r.1$ar ## ------------------------------------------------------------------- (arima.gw.r <- arima(gw.r, order=c(2,0,0))) ## ------------------------------------------------------------------- ar.gw.r$ar arima.gw.r$coef ## ------------------------------------------------------------------- gw.f.78 <- subset(gw.f, gw.f$year > 1977) coef(m.gw.78 <- lm(gw ~ time, data=gw.f.78)) ## ------------------------------------------------------------------- gw.1978 <- window(gw, c(1978,1), c(2004,12)) str(gw.1978) str(fitted(m.gw.78)) str(gw.1978.0 <- gw.1978 - fitted(m.gw.78)) plot(gw.1978.0, ylab="Deviations from linear trend, m", main="Anatolia well 1", pch=20) grid() abline(h=0, lty=2) ## ------------------------------------------------------------------- (par.gw.1978 <- arima(gw.1978.0, order=c(2,0,0), seasonal=c(0,0,0))) (par.gw.1978 <- arima(gw.1978.0, order=c(2,0,0), seasonal=c(1,0,0))) (par.gw.1978 <- arima(gw.1978.0, order=c(2,0,0), seasonal=c(2,0,0))) ## ----fig.width=12, fig.height=6, out.width='\\textwidth'------------ par(mfrow=c(1,2)) plot(gw, main="Anatolia well 1", ylab="Groundwater level (m)") plot(window(gw, 1990, 1994), main="Anatolia well 1", ylab="Groundwater level (m)") par(mfrow=c(1,2)) ## ----fig.width=12, fig.height=6, out.width='\\textwidth'------------ par(mfrow=c(1,2)) plot(diff(gw), main="Anatolia well 1", ylab="Groundwater level (m), delta-1") plot(diff(window(gw, 1990, 1994)), main="Anatolia well 1", ylab="Groundwater level (m), delta-1") par(mfrow=c(1,2)) ## ----fig.width=12, fig.height=6, out.width='\\textwidth'------------ par(mfrow=c(1,2)) plot(diff(diff(gw)), main="Anatolia well 1", ylab="Groundwater level (m), delta-2") plot(diff(diff(window(gw, 1990, 1994))), main="Anatolia well 1", ylab="Groundwater level (m), delta-2") par(mfrow=c(1,2)) ## ----fig.width=12, fig.height=6, out.width='\\textwidth'------------ par(mfrow=c(1,2)) plot(diff(diff(diff(gw))), main="Anatolia well 1", ylab="Groundwater level (m), delta-3") plot(diff(diff(diff(window(gw, 1990, 1994)))), main="Anatolia well 1", ylab="Groundwater level (m), delta-3") par(mfrow=c(1,2)) ## ----fig.height=12,fig.width=12, out.width='\\textwidth'------------ par(mfrow=c(2,2)) acf(gw) acf(diff(gw)) acf(diff(diff(gw))) acf(diff(diff(diff(gw)))) par(mfrow=c(1,1)) ## ----out.width='0.55\\textwidth'------------------------------------ pacf(gw.1978) ## ----fig.width=12,fig.height=6, out.width='\\textwidth'------------- par(mfrow=c(1,2)) pacf(diff(gw.1978)) pacf(diff(diff(gw.1978))) par(mfrow=c(1,1)) ## ----out.width='0.55\\textwidth'------------------------------------ acf(diff(diff(gw.1978))) ## ------------------------------------------------------------------- (m.ar <- arima(gw.1978, order=c(13,2,0))) ## ----m-ar-tsdiag, fig.width=8, fig.height=12------------------------ m.ar <- arima(gw.1978, c(3,2,0)) tsdiag(m.ar) ## ----m-ar-tsdig-2, fig.width=8, fig.height=12----------------------- m.ar <- arima(gw.1978, c(13,2,0)) tsdiag(m.ar) ## ----m-ar-predict--------------------------------------------------- p.ar <-predict(m.ar, n.ahead=12*(2013-2005)) str(p.ar) plot(p.ar$pred, ylim=c(35, 75), ylab="Predicted groundwater level, m", main="Anatolia well 1") lines(p.ar$pred+p.ar$se, col="red", lty=2) lines(p.ar$pred-p.ar$se, col="red", lty=2) grid() ## ----compute-AR1-coef, echo=FALSE----------------------------------- tmp <- exp(coef(m.gw.gls$modelStruct$corStruct)) gls.ar1.coeff <- round((tmp-1)/(tmp+1),4) ## ----child-data, child='ts_7.Rnw'----------------------------------- ## ----ts7-set-parent, echo=FALSE, cache=FALSE------------------------ set_parent('exTSA.Rnw') opts_chunk$set(fig.path='kgraph/ts7-') ## ----guelph-p-intervention, fig.height=5, fig.width=8--------------- require(Kendall) data(GuelphP) str(GuelphP) plot(GuelphP, type="b", ylab="P concentration, mg l-1", main="Speed River water quality") abline(v=1974+2/12, col="red", lwd=2) text(x=1974+2/12, y=1,"known intervention", col="red", pos=4) ## ------------------------------------------------------------------- (guelph.1 <- na.omit(as.vector(window(GuelphP, start=NULL, end=1974+1/12)))) (guelph.2 <- na.omit(as.vector(window(GuelphP, start=1974+2/12, end=NULL)))) mean(guelph.1); mean(guelph.2) median(guelph.1); median(guelph.2) sd(guelph.1); sd(guelph.2) ## ----child-data, child='ts_4.Rnw'----------------------------------- ## ----ts4-set-parent, echo=FALSE, cache=FALSE------------------------ set_parent('exTSA.Rnw') opts_chunk$set(fig.path='kgraph/ts4-') ## ----scan-alibe----------------------------------------------------- gw.2 <- ts(scan("./ds_tsa/anatolia_alibe.txt"), start=1975, frequency=12) plot(gw.2, ylab="Groundwater level (m)", main="Anatolia well 2") ## ----make-gw2------------------------------------------------------- gw2 <- ts.intersect(gw,gw.2) class(gw2) str(gw2) ## ----plot-gw2------------------------------------------------------- plot(gw2, plot.type="single", main="Anatolia wells 1 and 2", ylab="Groundwater depth (m)") lines(lowess(gw, f=1/3), col="red") lines(lowess(gw.2, f=1/3), col="red") ## ----summary-gw2---------------------------------------------------- summary(gw2) ## ----plot-gw2-multiple---------------------------------------------- plot(gw2, plot.type="multiple", main="Anatolia, two wells") ## ------------------------------------------------------------------- (cc <- ccf(gw,gw.2)) plot(cc, main="Cross-correlation, Anatolia well 1 vs. well 2") ## ------------------------------------------------------------------- str(cc) max(abs(cc$acf)) # check highest correlation (i <- which.max(abs(cc$acf))) cc$acf[i] cc$acf[i]^2 cc$lag[i] # check maximum lag 12*cc$lag[i] # maximum lag, in months ## ----compare-2-spectra, fig.height=5, fig.width=10, out.width='0.9\\textwidth'---- par(mfrow=c(1,3)) spectrum(gw2, spans=c(5,7), lty=1, col=c("black","red"), plot.type="marginal") spectrum(gw2, spans=c(5,7), lty=1, col=c("black","red"), plot.type="coherency") spectrum(gw2, spans=c(5,7), lty=1, col=c("black","red"), plot.type="phase") par(mfrow=c(1,1)) ## ----child-data, child='ts_5.Rnw'----------------------------------- ## ----ts5-set-parent, echo=FALSE, cache=FALSE------------------------ set_parent('exTSA.Rnw') opts_chunk$set(fig.path='kgraph/ts5-') ## ------------------------------------------------------------------- str(gw) ## ------------------------------------------------------------------- gwg <- gw ## ------------------------------------------------------------------- set.seed(0044) (ix <- sample(length(gw), size=5)) sort(ix) ## ------------------------------------------------------------------- summary(gwg) gwg[ix] <- NA summary(gwg) ## ----fig.width=10,fig.height=5, out.width='\\textwidth'------------- plot(gwg, main="Groundwater depth, with gaps", ylab="Groundwater depth (m)", sub="Dates with missing values shown as red bars") abline(h = min(gw), col = "gray") abline(v=time(gw)[ix], , col = "red", lwd=1) grid() ## ------------------------------------------------------------------- try(gwg.stl <- stl(gwg, s.window=25, t.window=85)) geterrmessage() ## ----try-stl-------------------------------------------------------- try(gwg.stl <- stl(gwg, s.window=25, t.window=85, na.action="na.exclude"), silent=TRUE) geterrmessage() ## ------------------------------------------------------------------- print(ix) gw[ix] gwg[ix] ## ------------------------------------------------------------------- time(gw)[ix] ## ------------------------------------------------------------------- gw.fill.linear <- approx(x=time(gwg), y=gwg, xout=time(gw)[ix]) str(gw.fill.linear) print(gw.fill.linear) ## ------------------------------------------------------------------- time(gw)[ix] gw.fill.linear$y gw[ix] summary((gw[ix] - gw.fill.linear$y)) ## ----fig.height=5, fig.width=10, out.width='\\textwidth'------------ plot(gwg, main="Gap-filled time series", sub="reconstructed values: red; true values: green", ylab="Groundwater depth (m)") points(gw.fill.linear$x, gw.fill.linear$y, col="red", cex=2) points(gw.fill.linear$x, gw.fill.linear$y, col="red", pch=20) points(time(gw)[ix], gw[ix], col="darkgreen", pch=20) ## ------------------------------------------------------------------- require(akima) ## ------------------------------------------------------------------- gw[ix] (gw.fill.aspline <- aspline(x=time(gwg), y=gwg, xout=time(gw)[ix]))$y (gw.fill.spline <- spline(x=time(gwg), y=gwg, xout=time(gw)[ix]))$y ## ------------------------------------------------------------------- summary(gw.fill.aspline$y - gw.fill.spline$y) summary(gw.fill.aspline$y - gw.fill.linear$y) summary(gw.fill.spline$y - gw.fill.linear$y) ## ----fig.height=6, fig.width=10, out.width='\\textwidth'------------ plot(gwg, main="Gap-filled time series", type="l", ylab="Groundwater depth (m)") points(gw.fill.aspline$x, gw.fill.aspline$y, col="red", cex=2) points(gw.fill.aspline$x, gw.fill.aspline$y, col="red", pch=20) points(gw.fill.spline$x, gw.fill.spline$y, col="blue", cex=2) points(gw.fill.spline$x, gw.fill.spline$y, col="blue", pch=20) points(gw.fill.linear$x, gw.fill.linear$y, col="brown", cex=2) points(gw.fill.linear$x, gw.fill.linear$y, col="brown", pch=20) points(time(gw)[ix], gw[ix], col="darkgreen", cex=2) points(time(gw)[ix], gw[ix], col="darkgreen", pch=20) text(2000, 35.5, "linear", col="brown", pos=2) text(2000, 34, "default spline", col="blue", pos=2) text(2000, 32.5, "Akima spline", col="red", pos=2) text(2000, 31, "true value", col="dark green", pos=2) ## ----fig.height=4, fig.width=6-------------------------------------- plot(window(gwg, start=1997, end=1997.5), main="Gap-filled time series", ylim=c(42,45), type="b", ylab="Groundwater depth (m)") points(gw.fill.aspline$x, gw.fill.aspline$y, col="red", cex=2) points(gw.fill.aspline$x, gw.fill.aspline$y, col="red", pch=20) points(gw.fill.spline$x, gw.fill.spline$y, col="blue", cex=2) points(gw.fill.spline$x, gw.fill.spline$y, col="blue", pch=20) points(gw.fill.linear$x, gw.fill.linear$y, col="brown", cex=2) points(gw.fill.linear$x, gw.fill.linear$y, col="brown", pch=20) points(time(gw)[ix], gw[ix], col="darkgreen", cex=2) points(time(gw)[ix], gw[ix], col="darkgreen", pch=20) text(1997.5, 43, "linear", col="brown", pos=2) text(1997.5, 42.7, "default spline", col="blue", pos=2) text(1997.5, 42.5, "Akima spline", col="red", pos=2) text(1997.5, 42.2, "true value", col="dark green", pos=2) ## ------------------------------------------------------------------- gwg.r <- gwg sum(is.na(gwg.r)) gwg.r[ix] <- gw.fill.aspline$y str(gwg.r) sum(is.na(gwg.r)) ## ------------------------------------------------------------------- gwg <- gw set.seed(0044) (six <- sort(ix <- sample(length(gw), size=length(gw)/5))) time(gw)[six] gwg[ix] <- NA summary(gwg) ## ----fig.width=10, fig.height=5, out.width='\\textwidth'------------ plot(gwg, ylab="Groundwater depth (m)", main="72 missing values") abline(v=time(gw)[six], , col = "red", lwd=1) grid() ## ------------------------------------------------------------------- gw.fill.linear <- approx(x=time(gwg), y=gwg, xout=time(gw)[six], rule=2) gw.fill.aspline <- aspline(x=time(gwg), y=gwg, xout=time(gw)[six]) summary(gw.fill.linear$y - gw[six]) summary(gw.fill.aspline$y - gw[six]) summary(gw.fill.aspline$y - gw.fill.linear$y) ## ----fig.height=6, fig.width=10, out.width='\\textwidth'------------ plot(gwg, main="Gap-filled time series", type="l", ylab="Groundwater depth (m)") points(gw.fill.aspline$x, gw.fill.aspline$y, col="red", cex=2) points(gw.fill.aspline$x, gw.fill.aspline$y, col="red", pch=20) points(gw.fill.linear$x, gw.fill.linear$y, col="brown", cex=2) points(gw.fill.linear$x, gw.fill.linear$y, col="brown", pch=20) points(time(gw)[ix], gw[ix], col="darkgreen", cex=2) points(time(gw)[ix], gw[ix], col="darkgreen", pch=20) text(2000, 34, "linear", col="brown", pos=2) text(2000, 32.5, "Akima spline", col="red", pos=2) text(2000, 31, "true value", col="dark green", pos=2) ## ----fig.height=6, fig.width=10, out.width='\\textwidth'------------ plot(window(gwg, start=1987, end=1989), main="Gap-filled time series", sub="True values in dark green", type="l", ylab="Groundwater depth (m)", ylim=c(39,47)) points(x=time(gw), y=as.numeric(gw), col="darkgreen", lty=2, pch=20, type="b") points(time(gw)[ix], gw[ix], col="darkgreen", cex=2) points(gw.fill.aspline$x, gw.fill.aspline$y, col="red", cex=2, type="b", lty=1) points(gw.fill.aspline$x, gw.fill.aspline$y, col="red", pch=20) points(gw.fill.linear$x, gw.fill.linear$y, col="brown", cex=2, type="b", lty=1) points(gw.fill.linear$x, gw.fill.linear$y, col="brown", pch=20) text(1989, 46, "linear", col="brown", pos=2) text(1989, 44.5, "Akima spline", col="red", pos=2) text(1989, 43, "true value", col="dark green", pos=2) ## ------------------------------------------------------------------- unique(floor(time(gw))) ## ------------------------------------------------------------------- gww <- window(gw, start=1990) gwg <- gww (six <- sort(which(floor(time(gwg))==1997))) gwg[six] <- NA summary(gwg) plot(gwg, main="Missing year 1997", ylab="Groundwater depth (m)") ## ------------------------------------------------------------------- gw.fill.linear <- approx(x=time(gwg), y=gwg, xout=time(gww)[six], rule=2) gw.fill.aspline <- aspline(x=time(gwg), y=gwg, xout=time(gww)[six]) summary(gw.fill.linear$y - gww[six]) summary(gw.fill.aspline$y - gww[six]) summary(gw.fill.aspline$y - gw.fill.linear$y) ## ----fig.height=6, fig.width=10, out.width='\\textwidth'------------ plot(window(gwg, start=1996, end=1999), main="Gap-filled time series", type="l", ylab="Groundwater depth (m)", ylim=c(43,49)) points(gw.fill.aspline$x, gw.fill.aspline$y, col="red", cex=2, lty=1, type="b") points(gw.fill.aspline$x, gw.fill.aspline$y, col="red", pch=20) points(gw.fill.linear$x, gw.fill.linear$y, col="brown", cex=2, lty=1, type="b") points(gw.fill.linear$x, gw.fill.linear$y, col="brown", pch=20) points(time(gww)[six], gww[six], col="darkgreen", cex=2, lty=1, type="b") points(time(gww)[six], gww[six], col="darkgreen", pch=20) text(1999, 45, "linear", col="brown", pos=2) text(1999, 44, "Akima spline", col="red", pos=2) text(1999, 43, "true value", col="dark green", pos=2) ## ------------------------------------------------------------------- tana.2<- read.csv("./ds_tsa/Tana_Daily_WeteAbay.csv", skip=1, header=T, colClasses=c(rep("integer",2), rep("character",12)), blank.lines.skip=T,na.strings=c("N.A","NA"," ")) str(tana.2) sum(is.na(tana.2[,3:14])) tana.3<- read.csv("./ds_tsa/Tana_Daily_Dangila.csv", skip=1, header=T, colClasses=c(rep("integer",2), rep("character",12)), blank.lines.skip=T,na.strings=c("N.A","NA"," ")) str(tana.3) sum(is.na(tana.3[,3:14])) ## ------------------------------------------------------------------- require(car) for (i in 3:14) { tana.2[,i] <- recode(tana.2[,i], "c('TR','tr','0.01')='0'") } for (i in 3:14) { tana.3[,i] <- recode(tana.3[,i], "c('TR','tr','0.01')='0'") } sum(c(tana.2[,3:14],tana.3[,3:14])=="TR", na.rm=TRUE) sum(c(tana.2[,3:14],tana.3[,3:14])=="tr", na.rm=TRUE) sum(c(tana.2[,3:14],tana.3[,3:14])=="0.01", na.rm=TRUE) sum(c(tana.2[,3:14],tana.3[,3:14])=="0", na.rm=TRUE) ## ------------------------------------------------------------------- sort(unique(tana$YEAR)) sort(unique(tana.2$Year)) sort(unique(tana.3$Year)) ## ------------------------------------------------------------------- tana.2[tana.2$DATE==29,"FEB"] tana.3[tana.3$DATE==29,"FEB"] ## ------------------------------------------------------------------- month.days <- c(0,0,31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31) tana.2.ppt <- NULL; for (yr.first.row in seq(from=1, by=32, length=(2006 - 2000 + 1))) { for (month.col in 3:14) { tana.2.ppt <- c(tana.2.ppt, tana.2[yr.first.row:(yr.first.row + month.days[month.col]-1), month.col]) } }; str(tana.2.ppt) length(tana.2.ppt)/365 tana.3.ppt <- NULL; for (yr.first.row in seq(from=1, by=32, length=(2006 - 1987 + 1))) { for (month.col in 3:14) { tana.3.ppt <- c(tana.3.ppt, tana.3[yr.first.row:(yr.first.row + month.days[month.col]-1), month.col]) } }; str(tana.3.ppt) length(tana.3.ppt)/365 rm(month.days, yr.first.row, month.col) ## ------------------------------------------------------------------- ## ------------------------------------------------------------------- tana.2.ppt <- ts(tana.2.ppt, start=2000, frequency=365) str(tana.2.ppt) tana.3.ppt <- ts(tana.3.ppt, start=1987, frequency=365) str(tana.3.ppt) ## ----fig.width=12,fig.height=4, out.width='\\textwidth'------------- plot(tana.ppt, main="Lake Tana rainfall, Station 1", ylab="mm", sub="Missing dates with red bars") abline(h=60, col="gray") points(xy.coords(x=time(tana.ppt), y=60, recycle=T), pch=ifelse(is.na(tana.ppt),"l",""), col="red") grid() ## ----fig.width=12,fig.height=4, out.width='\\textwidth'------------- plot(tana.2.ppt, main="Lake Tana rainfall, Station 2", ylab="mm", sub="Missing dates with red bars") abline(h=60, col="gray") points(xy.coords(x=time(tana.2.ppt), y=60, recycle=T), pch=ifelse(is.na(tana.2.ppt),"l",""), col="red") grid() ## ----fig.width=12,fig.height=4, out.width='\\textwidth'------------- plot(tana.3.ppt, main="Lake Tana rainfall, Station 3", ylab="mm", sub="Missing dates with red bars") abline(h=60, col="gray") points(xy.coords(x=time(tana.3.ppt), y=60, recycle=T), pch=ifelse(is.na(tana.3.ppt),"l",""), col="red") grid() ## ------------------------------------------------------------------- (miss.2 <- which(is.na(tana.2.ppt))) (time.miss.2 <- time(tana.2.ppt)[miss.2]) miss.3 <- which(is.na(tana.3.ppt)) time.miss.3 <- time(tana.3.ppt)[miss.3] miss.1 <- which(is.na(tana.ppt)) time.miss.1 <- time(tana.ppt)[miss.1] ## ------------------------------------------------------------------- length(miss.12 <- intersect(time.miss.1, time.miss.2)) length(miss.13 <- intersect(time.miss.1, time.miss.3)) length(miss.23 <- intersect(time.miss.2, time.miss.3)) ## ------------------------------------------------------------------- rm(miss.1, miss.2, miss.3, time.miss.1, time.miss.2, time.miss.3) ## ----fig.width=12, fig.height=8, out.width='\\textwidth'------------ t3 <- ts.intersect(tana.ppt, tana.2.ppt, tana.3.ppt) str(t3) class(t3) plot(t3, main="Lake Tana rainfall, 3 stations") ## ----child-data, child='ts_6.Rnw'----------------------------------- ## ----ts6-set-parent, echo=FALSE, cache=FALSE------------------------ set_parent('exTSA.Rnw') opts_chunk$set(fig.path='kgraph/ts6-') ## ----fig.height=16,fig.width=10, out.width='\\textwidth'------------ par(mfrow=c(4,1)) for (i in 1:3) { plot(arima.sim(model=list(ar=ar.gw.r.1$ar), n=12*15, rand.gen=rnorm, sd=sqrt(ar.gw.r.1$var.pred)), main=paste("Simulated AR(1) process",i), ylab="modelled") abline(h=0, lty=2) } plot(window(gw.r, 1989, 1989+15), main="Remainders 1989 -- 2004", ylab="actual") abline(h=0, lty=2) par(mfrow=c(1,1)) ## ----fig.height=16,fig.width=10, out.width='\\textwidth'------------ par(mfrow=c(4,1)) for (i in 1:3) { plot(arima.sim(model=list(ar=ar.gw.r$ar), n=12*15, rand.gen=rnorm, sd=sqrt(ar.gw.r$var.pred)), main=paste("Simulated AR(2) process",i), ylab="modelled") abline(h=0, lty=2) } plot(window(gw.r, 1989, 1989+15), main="Remainders 1989 -- 2004", ylab="actual") abline(h=0, lty=2) par(mfrow=c(1,1))