library("waveslim")
### read in oxygen isotope data from Antarctic ice core
### using one of the following two R commands:
temp <- scan("/Users/dbp/R/Lund/oxygen.dat")  ### will need to adjust pathname
temp <- scan("http://faculty.washington.edu/dbp/DATA/oxygen.dat")
### extract actual data from even indexed elements of temp
ox <- temp[2*(1:(length(temp)/2))]

### wavelet variance estimation ...
ox.modwt.haar.4 <- modwt(ox, "haar", 4)
ox.modwt.la8.4  <- modwt(ox, "la8", 4)
ox.modwt.haar.4.bw <- brick.wall(ox.modwt.haar.4, "haar")
ox.modwt.la8.4.bw  <- brick.wall(ox.modwt.la8.4,  "la8")
ox.wvar.haar       <- wave.variance(ox.modwt.haar.4.bw)
ox.wvar.haar
ox.wvar.la8        <- wave.variance(ox.modwt.la8.4.bw)
ox.wvar.la8
### verify a couple of these results ..
tW1 <- ox.modwt.la8.4$d1
### there are L-1 = 7 boundary coefficients at level 1
sum(tW1[-(1:7)]^2)/(length(tW1)-7)
###
tW4 <- ox.modwt.la8.4$d4
### there are L_4 -1 boundary coefficients at level 4,
### so need to compute L_4 for LA(8) wavelet
Lj <- function(L,j) (2^j - 1) * (L - 1) + 1
Lj(8,4)
### L_4 = 106, so L_4 -1 = 105 
sum(tW4[-(1:105)]^2)/(length(tW1)-105)
###
par(mfrow=c(1,1), las=1, mar=c(5,4,4,2)+.1)
matplot(2^(0:3), ox.wvar.haar[-5,], type="b", log="xy",
        xaxt="n", ylim=c(.01, 10.0), pch="*LU", lty=1, col=c(1,4,4),
        xlab="Wavelet Scale", ylab="", main = "Haar wavelet variance")
axis(side=1, at=2^(0:3))
###
quartz()  ### works on a Mac - need to adjust for Windows etc
par(mfrow=c(1,1), las=1, mar=c(5,4,4,2)+.1)
matplot(2^(0:3), ox.wvar.la8[-5,], type="b", log="xy",
        xaxt="n", ylim=c(.01, 10.0), pch="*LU", lty=1, col=c(1,4,4),
        xlab="Wavelet Scale", ylab="", main = "LA(8) wavelet variance")
axis(side=1, at=2^(0:3))

### read in Peter Craigmile's code for simulating FD processes
### using one of the following R commands:
source("/Users/dbp/R/Lund/Davies_Harte.R")   ### will need to adjust pathname
source("http://faculty.washington.edu/dbp/RCODE/Davies_Harte.R")
### generate a set of 1024 random normal deviates
rn.1024 <- rnorm(1024)
### define function for simulating an FD series and
### computing its wavelet variance, along with plots
wvar.of.FD.series <- function(delta, zs=rn.1024, plot.wv=T) {
  csum <- 0
  delta.title <- delta
  if(delta >= 0.5)
    {
      delta <- delta - 1.0
      csum = 1
    }
  sim.ts <- Davies.Harte.sim(DH.obj=Davies.Harte.sim.setup(512, fd.acvs, delta),
                             zs=zs,csum=csum)
  dev.set(2)
  par(mfrow=c(1,1))
  plot(sim.ts,typ="l",xlab="t",main=paste("FD(",delta.title,")",sep=""))
  ### plot(0:99,my.acf(sim.ts)[1:100],typ="h")
  if(plot.wv) {
    dev.set(3)
    par(mfrow=c(1,1), las=1, mar=c(5,4,4,2)+.1)
    modwt.d4    <- modwt(sim.ts, "d4", 7)
    modwt.d4.bw <- brick.wall(modwt.d4, "d4")
    sim.wvar.d4 <- wave.variance(modwt.d4.bw)
    matplot(2^(0:6), sim.wvar.d4[-8,], type="b", log="xy",
            xaxt="n", ylim=c(.001, 10.0), pch="*LU", lty=1, col=c(1,4,4),
            xlab="Wavelet Scale", ylab="", main = "D(4) wavelet variance")
    axis(side=1, at=2^(0:6))
  }
 } 
###
wvar.of.FD.series(-0.5, plot.wv=F)
wvar.of.FD.series(-0.4, plot.wv=F)
wvar.of.FD.series(-0.3, plot.wv=F)
wvar.of.FD.series(-0.2, plot.wv=F)
wvar.of.FD.series(-0.1, plot.wv=F)
wvar.of.FD.series(0.0, plot.wv=F)   ### white noise ...
wvar.of.FD.series(0.1, plot.wv=F)
wvar.of.FD.series(0.2, plot.wv=F)
wvar.of.FD.series(0.3, plot.wv=F)
wvar.of.FD.series(0.4, plot.wv=F)
wvar.of.FD.series(0.5, plot.wv=F)
wvar.of.FD.series(0.6, plot.wv=F)
wvar.of.FD.series(0.7, plot.wv=F)
wvar.of.FD.series(0.8, plot.wv=F)
wvar.of.FD.series(0.9, plot.wv=F)
wvar.of.FD.series(1.0, plot.wv=F)
###
wvar.of.FD.series(-0.5)
wvar.of.FD.series(-0.4)
wvar.of.FD.series(-0.3)
wvar.of.FD.series(-0.2)
wvar.of.FD.series(-0.1)
wvar.of.FD.series(0.0)   ### white noise ...
wvar.of.FD.series(0.1)
wvar.of.FD.series(0.2)
wvar.of.FD.series(0.3)
wvar.of.FD.series(0.4)
wvar.of.FD.series(0.5)
wvar.of.FD.series(0.6)
wvar.of.FD.series(0.7)
wvar.of.FD.series(0.8)
wvar.of.FD.series(0.9)
wvar.of.FD.series(1.0)

### look at Nile River data
data(nile)
nile <- nile/100.0
dev.set(2)
plot(nile)
### Haar MODWT MRA out to level 4 ...
dev.set(3)
nile.mra <- mra(nile, "haar", 4, "modwt")
tS4 <- nile.mra$S4
tD4 <- nile.mra$D4
tD3 <- nile.mra$D3
tD2 <- nile.mra$D2
tD1 <- nile.mra$D1
stack.plot(data.frame(cbind(tS4,tD4,tD3,tD2,tD1,nile)))
### wavelet variance ...
nile.up.to.715 <- nile[1:94]
nile.after.715 <- nile[-(1:94)]
### before ... 
dev.set(2)
par(mfrow=c(1,1), las=1, mar=c(5,4,4,2)+.1)
modwt.before    <- modwt(nile.up.to.715, "haar", 4)
modwt.before.bw <- brick.wall(modwt.before, "haar")
wvar.before     <- wave.variance(modwt.before.bw)
matplot(2^(0:3), wvar.before[-5,], type="b", log="xy",
          xaxt="n", ylim=c(.01, 1.0), pch="*LU", lty=1, col=c(1,4,4),
          xlab="Wavelet Scale", ylab="", main = "Haar wavelet variance up to 715")
axis(side=1, at=2^(0:3))
### after ... 
dev.set(3)
par(mfrow=c(1,1), las=1, mar=c(5,4,4,2)+.1)
modwt.after    <- modwt(nile.after.715, "haar", 4)
modwt.after.bw <- brick.wall(modwt.after, "haar")
wvar.after     <- wave.variance(modwt.after.bw)
matplot(2^(0:3), wvar.after[-5,], type="b", log="xy",
          xaxt="n", ylim=c(.01, 1.0), pch="*LU", lty=1, col=c(1,4,4),
          xlab="Wavelet Scale", ylab="", main = "Haar wavelet variance after 715")
axis(side=1, at=2^(0:3))