## ----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))