### R code for CSIRO half-day workshop `Wavelet Methods for Time Series Analysis'
###
### Part II: Wavelet Variance, Wavelet-Based Signal Extraction & Decorrelation of Time Series

library(wmtsa)

### load auxillary functions and data from downloaded file or from Internet

### print(load(file="/an/appropriate/folder/workshop.Rdata"))

con <- url("http://faculty.washington.edu/dbp/R-CODE/workshop.Rdata")
print(load(con))
close(con)

#############################################################################################
### atomic clock example
#############################################################################################

### time series in overhead II-35

lplot(clock/1000)
lplot(diff(clock))
lplot(diff(diff(clock)))

clock.wv.haar <- wavVar(atomclock,wave="haar")
clock.wv.haar
summary(clock.wv.haar)
plot(clock.wv.haar)
x.haar <- log10(clock.wv.haar$scales)  # for use below ...
y.haar <- log10(sqrt(clock.wv.haar$block$unbiased))

### create plot essentially same as right-hand side of overhead II-37 (latter shows sd
### rather than var) 

clock.wv.d4 <- wavVar(atomclock,wave="d4")
clock.wv.d4
summary(clock.wv.d4)
plot(clock.wv.d4)
x.d4 <- log10(clock.wv.d4$scales)  # for use below ...
y.d4 <- log10(sqrt(clock.wv.d4$block$unbiased))
lines(x.d4[1:4], lm(I(2*y.d4[1:4]) ~ I(x.d4[1:4]))$fitted.values, col="red")
lines(x.d4[5:8], lm(I(2*y.d4[5:8]) ~ I(x.d4[5:8]))$fitted.values, col="red")

clock.wv.d6 <- wavVar(atomclock,wave="d6")
clock.wv.d6
summary(clock.wv.d6)
plot(clock.wv.d6)
x.d6 <- log10(clock.wv.d6$scales)  # for use below ...
y.d6 <- log10(sqrt(clock.wv.d6$block$unbiased))

### create composite plot like left-hand side of II-37

plot(x.haar,y.haar,ylim=c(-3,2),pch="x")
lines(x.haar, lm(y.haar ~ x.haar)$fitted.values, col="red")
points(x.d4,y.d4,pch="o")
points(x.d6,y.d6,pch="+",col="blue")

#############################################################################################
### Nile River example of overhead II-39 ...
### first 94 observations cover years 622 to 715, while next 569 cover years 716 to 1284
#############################################################################################

lplot(nile.river)
nile.river.modwt <- wavMODWT(nile.river,wave="haar",n.levels=4)
plot(wavMRD(nile.river.modwt))
abline(v=94.5,lty="dotted",col="red")

nile.river.1.wv <- wavVar(nile.river[1:94],wave="haar")
nile.river.2.wv <- wavVar(nile.river[-(1:94)],wave="haar")

plot(nile.river.1.wv)
plot(nile.river.2.wv)

### create composite plot like right-hand side of II-39

x.2 <- log10(nile.river.2.wv$scales)
y.2 <- log10(nile.river.2.wv$block$unbiased)
plot(x.2,y.2,pch=19,col="blue", cex=1.5, xlim=c(0,1),ylim=c(-2,0),xlab="log10(scale)",ylab="log10(wavelet variances)")
conf <- nile.river.2.wv$confidence$n3
low <- log10(conf$low)
high <- log10(conf$high)
arrows(x.2, low, x.2, high, code = 3, col = "black", angle = 90, lwd = 2, length = 0.1)

x.1 <- log10(nile.river.1.wv$scales) + 0.05
y.1 <- log10(nile.river.1.wv$block$unbiased)
points(x.1,y.1,pch=19,col="red", cex=1.5)
conf <- nile.river.1.wv$confidence$n3
low <- log10(conf$low)
high <- log10(conf$high)
arrows(x.1, low, x.1, high, code = 3, col = "black", angle = 90, lwd = 2, length = 0.1)

#############################################################################################
### subtidal sea level example
#############################################################################################

lplot(subtidal)

### overhead II-42

st.modwt.01 <- wavMODWT(subtidal,n.levels=7)
plot(wavMRD(st.modwt.01))

### blow-up of portion of overhead II-42

st.modwt.02 <- wavMODWT(subtidal,n.levels=4)
plot(wavMRD(st.modwt.02))

### global look at wavelet variance ...

st.wv <- wavVar(subtidal)
plot(st.wv)

### overhead II-43 ... break up global wavelet variance across time ...

j <- 2
L <- 8
zap.me <- (2^j - 1)*(L-1)
sq.coeffs <- (st.modwt.01$data[[j]][-(1:zap.me)])^2
lplot(sq.coeffs)
lplot(log10(sq.coeffs))
dft.sq.coeffs <- dft(sq.coeffs)
N <- length(sq.coeffs)
rec.filter <- rep(0,N)
rec.filter[1:31] <- 1/61
rec.filter[(N-29):N] <- 1/61
sum(rec.filter)
dft.rec.filter <- dft(rec.filter)
tv.wv <- Re(inverse.dft(dft.sq.coeffs*dft.rec.filter))
lplot(log10(tv.wv))
abline(v=seq(0,9000,2*365.25),lty="dotted",col="red")

### do the same with a smoother that has the same effective bandwidth of the
### rectangular filer, but has a shape close to Gaussian ...

1/sum(rec.filter^2)  # 61 (bandwidth of rectangular filter)
small.rec.filter <- rep(0,N)
small.rec.filter[1:5] <- 1/10
small.rec.filter[(N-3):N] <- 1/10
small.rec.filter[6] <- 1/20
small.rec.filter[N-4] <- 1/20

sum(small.rec.filter)

dft.small.rec.filter <- dft(small.rec.filter)

1/sum(Re(inverse.dft(dft.small.rec.filter^35))^2)
irs <- Re(inverse.dft(dft.small.rec.filter^35))
lplot(irs[1:100])

tv.wv.smoother <- Re(inverse.dft(dft.sq.coeffs*dft.small.rec.filter^35))
par(mfrow=c(2,1))
lplot(log10(tv.wv.smoother))
lplot(log10(tv.wv))
par(mfrow=c(1,1))

### overhead II-57: sinusoid, bump & linear combination of sinusoid and bump and
###                 their NPESs

ts.1 <- 0.5 * cos(3*pi*(0:127)/32 + 0.08)
ts.2 <- rep(0,128)
ts.2[60:70] <- ts.1[60:70]
ts.3 <- ts.1/sqrt(120) + ts.2
lplot(ts.1)
lplot(ts.2)
lplot(ts.3)

td.npec.ts.1 <- cumsum(rev(sort(ts.1^2)))/ss(ts.1)
fd.npec.ts.1 <- cumsum(rev(sort((abs(dft(ts.1))^2)/128)))/ss(ts.1)
wd.npec.ts.1 <- cumsum(rev(sort(abs(unlist(wavDWT(ts.1)$data)^2))))/ss(ts.1)

plot(0:127,td.npec.ts.1,typ="l",ylim=c(0.9,1),main="D_1 (frequency domain signal)")
lines(0:127,fd.npec.ts.1,lty="dotted",col="red")
lines(0:127,wd.npec.ts.1,lty="dashed",col="blue")
legend(30, 0.94, c("identity","DFT","DWT"), lty=c("solid","dotted","dashed"), col=c("black","red","blue"))

td.npec.ts.2 <- cumsum(rev(sort(ts.2^2)))/ss(ts.2)
fd.npec.ts.2 <- cumsum(rev(sort((abs(dft(ts.2))^2)/128)))/ss(ts.2)
wd.npec.ts.2 <- cumsum(rev(sort(abs(unlist(wavDWT(ts.2)$data)^2))))/ss(ts.2)

plot(0:127,td.npec.ts.2,typ="l",ylim=c(0.9,1),main="D_2 (time domain signal)")
lines(0:127,fd.npec.ts.2,lty="dotted",col="red")
lines(0:127,wd.npec.ts.2,lty="dashed",col="blue")
legend(30, 0.94, c("identity","DFT","DWT"), lty=c("solid","dotted","dashed"), col=c("black","red","blue"))

td.npec.ts.3 <- cumsum(rev(sort(ts.3^2)))/ss(ts.3)
fd.npec.ts.3 <- cumsum(rev(sort((abs(dft(ts.3))^2)/128)))/ss(ts.3)
wd.npec.ts.3 <- cumsum(rev(sort(abs(unlist(wavDWT(ts.3)$data)^2))))/ss(ts.3)

plot(0:127,td.npec.ts.3,typ="l",ylim=c(0.9,1),main="D_3 (mixture signal)")
lines(0:127,fd.npec.ts.3,lty="dotted",col="red")
lines(0:127,wd.npec.ts.3,lty="dashed",col="blue")
legend(30, 0.94, c("identity","DFT","DWT"), lty=c("solid","dotted","dashed"), col=c("black","red","blue"))

### II-58: NMR example

nmr.dft <- dft(nmr)
temp <- rev(sort(abs(nmr.dft)))
thres.50  <- mean(temp[50:51])
thres.100 <- mean(temp[100:101])
thres.200 <- mean(temp[200:201])
thres.400 <- mean(temp[400:401])
nmr.dft[abs(nmr.dft) < thres.400] <- rep(0,1024-400)
nmr.dft.400 <- Re(inverse.dft(nmr.dft))
nmr.dft[abs(nmr.dft) < thres.200] <- rep(0,1024-200)
nmr.dft.200 <- Re(inverse.dft(nmr.dft))
nmr.dft[abs(nmr.dft) < thres.100] <- rep(0,1024-100)
nmr.dft.100 <- Re(inverse.dft(nmr.dft))
nmr.dft[abs(nmr.dft) < thres.50] <- rep(0,1024-50)
nmr.dft.50 <- Re(inverse.dft(nmr.dft))

lplot(nmr,ylim=c(-15,60),main="original series")
lplot(nmr.dft.50,ylim=c(-15,60),main="DFT, 50 coefficients")
lplot(nmr.dft.100,ylim=c(-15,60),main="DFT, 100 coefficients")
lplot(nmr.dft.200,ylim=c(-15,60),main="DFT, 200 coefficients")
lplot(nmr.dft.400,ylim=c(-15,60),main="DFT, 400 coefficients")

nmr.dwt.50  <- recon.ts.from.largest.LA8.coeffs(nmr,50)
nmr.dwt.100 <- recon.ts.from.largest.LA8.coeffs(nmr,100)
nmr.dwt.200 <- recon.ts.from.largest.LA8.coeffs(nmr,200)
nmr.dwt.400 <- recon.ts.from.largest.LA8.coeffs(nmr,400)

lplot(nmr,ylim=c(-15,60),main="original series")
lplot(nmr.dwt.50,ylim=c(-15,60),main="LA(8), 50 coefficients")
lplot(nmr.dwt.100,ylim=c(-15,60),main="LA(8), 100 coefficients")
lplot(nmr.dwt.200,ylim=c(-15,60),main="LA(8), 200 coefficients")
lplot(nmr.dwt.400,ylim=c(-15,60),main="LA(8), 400 coefficients")

### II-60, II-64, II-66: hard/soft/mid thresholding functions

xs<- seq(-3,3,0.001)
plot(xs,sapply(xs,hard.thresholding),typ="l",xlim=c(-3,3),ylim=c(-3,3),main="thresholding functions",xlab="original coefficient",ylab="thresholded coefficient")
lines(xs,sapply(xs,soft.thresholding),col="blue")
lines(xs,sapply(xs,mid.thresholding),col="red")
legend(-2.75,2.5,c("hard","soft","mid"),lty=rep("solid",3),col=c("black","blue","red"))

### II-73 and II-74: hard thresholding of NMR data using D(4) and LA(8) wavelets

max.nmr <- max(nmr)
ws.nmr.d4.hard <- wavShrink(nmr,n.level=6,wav="d4",reflect=FALSE,xform="dwt")
ws.nmr.la8.hard <- wavShrink(nmr,n.level=6,reflect=FALSE,xform="dwt")

lplot(nmr,ylim=c(-15,60),main="original series")
lines(c(460,480),rep(max.nmr,2),col="red",lwd=2)

lplot(ws.nmr.d4.hard,ylim=c(-15,60),main="D(4), hard thresholding")
lines(c(460,480),rep(max.nmr,2),col="red",lwd=2)

lplot(ws.nmr.la8.hard,ylim=c(-15,60),main="LA(8), hard thresholding")
lines(c(460,480),rep(max.nmr,2),col="red",lwd=2)

### II-75 and II-76: soft/mid LA(8) thresholding of NMR data

ws.nmr.la8.soft <- wavShrink(nmr,n.level=6,reflect=FALSE,xform="dwt",shrink.fun="soft")
ws.nmr.la8.mid  <- wavShrink(nmr,n.level=6,reflect=FALSE,xform="dwt",shrink.fun="mid")

lplot(ws.nmr.la8.hard,ylim=c(-15,60),main="LA(8), hard thresholding")
lines(c(460,480),rep(max.nmr,2),col="red",lwd=2)

lplot(ws.nmr.la8.soft,ylim=c(-15,60),main="LA(8), soft thresholding")
lines(c(460,480),rep(max.nmr,2),col="red",lwd=2)

lplot(ws.nmr.la8.mid,ylim=c(-15,60),main="LA(8), mid thresholding")
lines(c(460,480),rep(max.nmr,2),col="red",lwd=2)

### II-78, II-79 and II-80: hard/soft/mid LA(8) MODWT-based thresholding of NMR data

ws.nmr.la8.hard.mod <- wavShrink(nmr,n.level=6,reflect=FALSE,xform="modwt",shrink.fun="hard")
ws.nmr.la8.soft.mod <- wavShrink(nmr,n.level=6,reflect=FALSE,xform="modwt",shrink.fun="soft")
ws.nmr.la8.mid.mod  <- wavShrink(nmr,n.level=6,reflect=FALSE,xform="modwt",shrink.fun="mid")

lplot(ws.nmr.la8.hard.mod,ylim=c(-15,60),main="LA(8), hard thresholding, MODWT")
lines(c(460,480),rep(max.nmr,2),col="red",lwd=2)

lplot(ws.nmr.la8.soft.mod,ylim=c(-15,60),main="LA(8), soft thresholding, MODWT")
lines(c(460,480),rep(max.nmr,2),col="red",lwd=2)

lplot(ws.nmr.la8.mid.mod,ylim=c(-15,60),main="LA(8), mid thresholding, MODWT")
lines(c(460,480),rep(max.nmr,2),col="red",lwd=2)

### II-88: conditional mean shrinkage rule (sparse signal model)

xs<- seq(-8,8,0.001)
plot(xs,sapply(xs,cm.shrinkage),typ="l",col="blue",xlim=c(-6,6),ylim=c(-6,6),main="conditional mean shrinkage functions",xlab="original coefficient",ylab="shrunk coefficient")
abline(a=0,b=1,lty="dotted")
lines(xs,sapply(xs,cm.shrinkage,G=10))
lines(xs,sapply(xs,cm.shrinkage,G=25),col="red")
legend(-5.5,5.5,c("sig^2_G = 5","sig^2_G = 10","sig^2_G = 25"),lty=rep("solid",3),col=c("blue","black","red"))

### II-91, II-92 and II-93: conditional mean shrinkage applied to NMR data ...

nmr.dwt <- wavDWT(nmr,n.level=6)
sigma.eps.2 <- mad(nmr.dwt$data[[1]],center=0)^2
wav.coeffs <- head(unlist(nmr.dwt$data),length(nmr)-length(nmr.dwt$data[[7]]))
sigma.W.2 <- mean(wav.coeffs^2)
p.1 <- 0.9
sigma.G.2.1 <- (sigma.W.2-sigma.eps.2)/(1-p.1)
p.2 <- 0.95
sigma.G.2.2 <- (sigma.W.2-sigma.eps.2)/(1-p.2)
p.3 <- 0.99
sigma.G.2.3 <- (sigma.W.2-sigma.eps.2)/(1-p.3)

nmr.signal.1 <- my.wavShrink(nmr.dwt,G=sigma.G.2.1,n=sigma.eps.2,p=p.1)
nmr.signal.2 <- my.wavShrink(nmr.dwt,G=sigma.G.2.2,n=sigma.eps.2,p=p.2)
nmr.signal.3 <- my.wavShrink(nmr.dwt,G=sigma.G.2.3,n=sigma.eps.2,p=p.3)
max.nmr <- max(nmr)

lplot(nmr.signal.1,ylim=c(-15,60),main="conditional mean shrinkage, p = 0.9")
lines(c(460,480),rep(max.nmr,2),col="red",lwd=2)

lplot(nmr.signal.2,ylim=c(-15,60),main="conditional mean shrinkage, p = 0.95")
lines(c(460,480),rep(max.nmr,2),col="red",lwd=2)

lplot(nmr.signal.3,ylim=c(-15,60),main="conditional mean shrinkage, p = 0.99")
lines(c(460,480),rep(max.nmr,2),col="red",lwd=2)

### exact simulation of FD process ...

delta <- 0.4
CE.stuff.1024 <- sim.circ.embed.setup(1024,fd.acvs,delta)
set.seed(42)
fd.ts <- sim.circ.embed.from.setup(CE.stuff.1024)

augmented.acf.plot <- function(dwt,j)
  {
    w.coeffs <- dwt$data[[j]]
    N.j <- length(w.coeffs)
    max.lag <- min(32,length(w.coeffs)-1)
    acf(w.coeffs,lag.max=max.lag,demean=FALSE,main=paste("level",j))
    taus <- 1:max.lag
    lines(taus,2*sqrt(N.j-taus)/N.j,lty="dashed",col="red")
    lines(taus,-2*sqrt(N.j-taus)/N.j,lty="dashed",col="red")
  }

### create plots similar to left- and right-hand sides of overhead II-98

lplot(fd.ts)
acf(fd.ts,lag.max=32,demean=FALSE)

### create plots similar to those comprising overhead II-99

dwt.fd.ts <- wavDWT(fd.ts,n.level=7)
plot(dwt.fd.ts)
augmented.acf.plot(dwt.fd.ts,1)
augmented.acf.plot(dwt.fd.ts,2)
augmented.acf.plot(dwt.fd.ts,3)
augmented.acf.plot(dwt.fd.ts,4)
augmented.acf.plot(dwt.fd.ts,5)
augmented.acf.plot(dwt.fd.ts,6)
augmented.acf.plot(dwt.fd.ts,7)

### overhead II-100

modwt.fd.ts <- wavMODWT(fd.ts,n.level=7)
plot(modwt.fd.ts)
augmented.acf.plot(modwt.fd.ts,1)
augmented.acf.plot(modwt.fd.ts,2)
augmented.acf.plot(modwt.fd.ts,3)
augmented.acf.plot(modwt.fd.ts,4)
augmented.acf.plot(modwt.fd.ts,5)
augmented.acf.plot(modwt.fd.ts,6)
augmented.acf.plot(modwt.fd.ts,7)

### plots similar to what is shown in overhead II-112

ts.sim.01 <- wavSimulateSeries()  # default is FD(0.4)
lplot(ts.sim.01)
as.vector(acf(ts.sim.01,plot=FALSE,demean=FALSE)$acf[1:10])

ts.sim.02 <- wavSimulateSeries(delta=0)  # white noise
lplot(ts.sim.02)
as.vector(acf(ts.sim.02,plot=FALSE,demean=FALSE)$acf[1:10])

ts.sim.03 <- wavSimulateSeries(delta=-0.4)  # antipersistent
lplot(ts.sim.03)
as.vector(acf(ts.sim.03,plot=FALSE,demean=FALSE)$acf[1:10])

### simulate a pure power law process rather than an FD process

ppl.sdf <- function(f,Cs=1,alpha=-0.25) Cs*f^alpha
ppl.var <- function(Cs=1,alpha=-0.25) Cs/((2^alpha) * (alpha + 1))

ts.sim.04 <- wavSimulateSeries(sdf=ppl.sdf,var=ppl.var)  # default is PPL(-0.25)
lplot(ts.sim.04)
as.vector(acf(ts.sim.04,plot=FALSE,demean=FALSE)$acf[1:10])

ts.sim.05 <- wavSimulateSeries(sdf=ppl.sdf,var=ppl.var,alpha=-0.999)
lplot(ts.sim.05)
as.vector(acf(ts.sim.05,plot=FALSE,demean=FALSE)$acf[1:10])

ts.sim.06 <- wavSimulateSeries(sdf=ppl.sdf,var=ppl.var,alpha=0)  # white noise
lplot(ts.sim.06)
as.vector(acf(ts.sim.06,plot=FALSE,demean=FALSE)$acf[1:10])

ts.sim.07 <- wavSimulateSeries(sdf=ppl.sdf,var=ppl.var,alpha=1)  # antipersistent
lplot(ts.sim.07)
as.vector(acf(ts.sim.07,plot=FALSE,demean=FALSE)$acf[1:10])

#############################################################################################
### bootstrapping Nile River ...
### first 94 observations cover years 622 to 715
### next 569 cover years 716 to 1284
#############################################################################################

nile.512 <- tail(nile.river,512)
lplot(nile.512,ylim=c(8,15))

log2(512)  # 9
dwt.nile.512 <- wavDWT(nile.512,n.level=6)

lplot(wavDWTbootstrap(dwt.nile.512)[[1]],ylim=c(8,15))

bootstrapped.nile.512 <- wavDWTbootstrap(dwt.nile.512,n.realization=100)

for(n in 1:length(bootstrapped.nile.512)) lplot(bootstrapped.nile.512[[n]],ylim=c(8,15))

acvs.lag.1 <- function(ts)
  {
    ts <- ts - mean(ts)
    sum(ts[-1]*ts[-length(ts)])/ss(ts)
  }

bootstrapped.samples <- sapply(bootstrapped.nile.512,acvs.lag.1)
sd(bootstrapped.samples)
plot(density(bootstrapped.samples))
abline(v=acvs.lag.1(nile.512),lty="dotted",col="red")

### NOTE: the wmtsa library has a function call wavBootstrap.  This is
### based upon an algorithm described in Percival, Sardy and Davison (2001)
### and differs from what is described in overhead II-22.
###   bootstrapped.nile.512 <- wavBootstrap(nile.512,n.realization=100)

#############################################################################################
### blue curve on overhead II-129
### MLE of FD parameters ... wavFDPBlock in wmtsa library should compute this, but is
### giving strange results for Nile River example, so will use alternative code instead
### (only works for N that is a power of 2)
#############################################################################################

dwt.nile.512 <- wavDWT(nile.512-mean(nile.512))

deltas <- seq(0.01,0.49,0.01)
reduced.ll.nile <- sapply(deltas,reduced.ll.approx.MLE.DWT,dwt.nile.512)
lplot(deltas,reduced.ll.nile)

optimize(reduced.ll.approx.MLE.DWT,c(0,0.49),dwt.nile.512)  # 0.4531676

#############################################################################################
### homogeneity of variance ...
#############################################################################################

dwt.nile <- wavDWT(nile.river,wavelet="haar",n.level=4)
dwt.nile.512 <- wavDWT(nile.512,wavelet="haar",n.level=4)

wavVarTest(dwt.nile)  # close to - but not exactly the same as -- values in overhead II-130
wavVarTest(dwt.nile.512)