WASP: An R package for complex system modelling and prediction

15 minute read

Published: June 24, 2023

The wavelet-based variance transformation method is used for system modelling and prediction. It refines predictor spectral representation using Wavelet Theory, which leads to improved model specifications and prediction accuracy. A supporting open-source software, Wavelet System Prediction (WASP), can be found under page of Software.

Required packages
DWT, MODWT and AT basic propertites
- Summary of various properties for the three DWT methods
- Illustration of three types of DWT methods
  - Daubechies 16 wavelet
  - Haar wavelet filter
Optimal variance transformation
- Preditive accuracy (RMSE)
- Transformed predictor variables
Stepwise variance transformation

Required packages

if(!require(SPEI)) devtools::install_github('sbegueria/SPEI@v1.7.1') # use 1.7.1
#> Loading required package: SPEI
#> Loading required package: lmomco
#> Loading required package: parallel
#> Loading required package: ggplot2
#> # Package SPEI (1.7) loaded [try SPEINews()].

library(WASP)
library(synthesis)
#> 
#> Attaching package: 'synthesis'
#> The following objects are masked from 'package:WASP':
#> 
#>     data.gen.ar1, data.gen.ar4, data.gen.ar9, data.gen.HL,
#>     data.gen.Rossler, data.gen.SW, data.gen.tar1, data.gen.tar2

library(ggplot2)
require(SPEI)
library(FNN)
#> 
#> Attaching package: 'FNN'
#> The following object is masked from 'package:WASP':
#> 
#>     knn
library(waveslim)
#> 
#> waveslim: Wavelet Method for 1/2/3D Signals (version = 1.8.4)
library(cowplot)
library(gridGraphics)
#> Loading required package: grid

DWT, MODWT and AT basic propertites

# data generation
x <- arima.sim(list(order = c(1, 0, 0), ar = 0.6), n = 512)
# x <- as.numeric(scale(data.gen.Rossler(time = seq(0, 50, length.out = 512))$x, scale=F))

# Daubechies wavelets
for (wf in c("haar", "d4", "d8", "d16")) {
  print(paste0("Wavelet filter: ", wf))
  #----------------------------------------------------------------------------
  # wavelet family, extension mode and package
  # wf <- "haar" # wavelet family D8 or db4
  boundary <- "periodic"
  if (wf != "haar") v <- as.integer(readr::parse_number(wf) / 2) else v <- 1

  # Maximum decomposition level J
  n <- length(x)
  J <- ceiling(log(n / (2 * v - 1)) / log(2)) - 1 # (Kaiser, 1994)

  cov <- rnorm(J + 1, sd = 2)
  Vr <- as.numeric(cov / norm(cov, type = "2") * sd(x))
  #----------------------------------------------------------------------------
  # DWT-MRA
  #print("-----------DWT-MRA-----------")
  x.mra <- waveslim::mra(x, wf = wf, J = J, method = "dwt", boundary = boundary)
  x.mra.m <- matrix(unlist(x.mra), ncol = J + 1)

  x.n <- scale(x.mra.m) %*% Vr
  var(x.n) - var(x)

  message(paste0("Additive decompostion: ", isTRUE(all.equal(as.numeric(x), rowSums(x.mra.m)))))
  message(paste0("Variance decompostion: ", isTRUE(all.equal(var(x), sum(apply(x.mra.m, 2, var))))))

  #----------------------------------------------------------------------------
  # MODWT
  #print("-----------MODWT-----------")
  x.modwt <- waveslim::modwt(x, wf = wf, n.levels = J, boundary = boundary)
  x.modwt.m <- matrix(unlist(x.modwt), ncol = J + 1)

  x.n <- scale(x.modwt.m) %*% Vr
  var(x.n) - var(x)

  message(paste0("Additive decompostion: ", isTRUE(all.equal(as.numeric(x), rowSums(x.modwt.m)))))
  message(paste0("Variance decompostion: ", isTRUE(all.equal(var(x), sum(apply(x.modwt.m, 2, var))))))

  #----------------------------------------------------------------------------
  # a trous
  #print("-----------AT-----------")
  x.at <- at.wd(x, wf = wf, J = J, boundary = boundary)
  x.at.m <- matrix(unlist(x.at), ncol = J + 1)

  # x.mra.modwt <- waveslim::mra(x,wf=wf, J=J, method="modwt", boundary=boundary)
  # x.mra.modwt <- matrix(unlist(x.mra.modwt), ncol=J+1)
  #
  # print(sum(abs(x.at.m-x.mra.modwt)))

  message(paste0("Additive decompostion: ", isTRUE(all.equal(as.numeric(x), rowSums(x.at.m)))))
  message(paste0("Variance decompostion: ", isTRUE(all.equal(var(x), sum(apply(x.at.m, 2, var))))))

  if (isTRUE(all.equal(x.at.m, x.modwt.m))) {
    message(paste0("AT and MODWT is equivalent using the", wf, "!"))
  }
}
#> [1] "Wavelet filter: haar"
#> Additive decompostion: TRUE
#> Variance decompostion: TRUE
#> Additive decompostion: TRUE
#> Variance decompostion: TRUE
#> Additive decompostion: TRUE
#> Variance decompostion: TRUE
#> AT and MODWT is equivalent using thehaar!
#> [1] "Wavelet filter: d4"
#> Additive decompostion: TRUE
#> Variance decompostion: TRUE
#> Additive decompostion: FALSE
#> Variance decompostion: TRUE
#> Additive decompostion: TRUE
#> Variance decompostion: FALSE
#> [1] "Wavelet filter: d8"
#> Additive decompostion: TRUE
#> Variance decompostion: TRUE
#> Additive decompostion: FALSE
#> Variance decompostion: TRUE
#> Additive decompostion: TRUE
#> Variance decompostion: FALSE
#> [1] "Wavelet filter: d16"
#> Additive decompostion: TRUE
#> Variance decompostion: TRUE
#> Additive decompostion: FALSE
#> Variance decompostion: TRUE
#> Additive decompostion: TRUE
#> Variance decompostion: FALSE

Summary of various properties for the three DWT methods

Summary of various properties for the three DWT methods
Wavelet Method	Additive decomposition	Variance decomposition	No dependence on future data	Dyadic sample size
DWT-MRA	✓	✓		✓
MODWT		✓	✓
AT	✓		✓
Note: When Haar wavelet filter is used, MODWT and AT are equivalent and both of them preserves additive and variance decomposition.

Illustration of three types of DWT methods

p.list <- NULL
wf.opts <- c("d16", "haar")
for (k in seq_along(wf.opts)) {
  # data generation
  x <- arima.sim(list(order = c(1, 0, 0), ar = 0.6), n = 128)

  #----------------------------------------------------------------------------
  # wavelet family, extension mode and package
  wf <- wf.opts[k] # wavelet family D8 or db4
  boundary <- "periodic"
  if (wf != "haar") v <- as.integer(readr::parse_number(wf) / 2) else v <- 1

  # Maximum decomposition level J
  n <- length(x)
  J <- ceiling(log(n / (2 * v - 1)) / log(2)) - 1 # (Kaiser, 1994)

  limits.x <- c(0, n)
  limits.y <- c(-3, 3)
  #----------------------------------------------------------------------------
  # DWT-MRA
  x.mra <- waveslim::mra(x, wf = wf, J = J, method = "dwt", boundary = boundary)
  x.mra.m <- matrix(unlist(x.mra), ncol = J + 1)

  p1 <- mra.plot(x, x.mra.m, limits.x, limits.y,
    ylab = "X", col = "red", type = "details",
    main = paste0("DWT-MRA", "(", wf, ")"), ps = 12
  )
  # p1 <- recordPlot()

  #----------------------------------------------------------------------------
  # MODWT
  x.modwt <- waveslim::modwt(x, wf = wf, n.levels = J, boundary = boundary)
  x.modwt.m <- matrix(unlist(x.modwt), ncol = J + 1)

  p2 <- mra.plot(x, x.modwt.m, limits.x, limits.y,
    ylab = "X", col = "red", type = "coefs",
    main = paste0("MODWT", "(", wf, ")"), ps = 12
  )

  #----------------------------------------------------------------------------
  # a trous
  x.at <- at.wd(x, wf = wf, J = J, boundary = boundary)
  x.at.m <- matrix(unlist(x.at), ncol = J + 1)

  p3 <- mra.plot(x, x.at.m, limits.x, limits.y,
    ylab = "X", col = "red", type = "coefs",
    main = paste0("AT", "(", wf, ")"), ps = 12
  )

  p.list[[k]] <- list(p1, p2, p3)
}

Daubechies 16 wavelet

#----------------------------------------------------------------------------
# plot and save
cowplot::plot_grid(
  plotlist = p.list[[1]], ncol = 3, labels = c("(a)", "(b)", "(c)"),
  label_size = 12
)

Illustration of three types of DWT methods

Haar wavelet filter

#----------------------------------------------------------------------------
# plot and save
cowplot::plot_grid(
  plotlist = p.list[[2]], ncol = 3, labels = c("(a)", "(b)", "(c)"),
  label_size = 12
)

Illustration of three types of DWT methods

Optimal variance transformation

Preditive accuracy (RMSE)

if (FALSE) {
  ### Synthetic example
  # data generation
  set.seed(2020)
  sample <- 512
  # frequency, sampled from a given range
  fd <- c(3, 5, 10, 15, 25, 30, 55, 70, 95)
  # data <- WASP::data.gen.SW(nobs=sample,fp=25,fd=fd)
  data <- WASP::data.gen.SW(nobs = sample, fp = c(15, 25, 30), fd = fd)

  # ts = data.gen.Rossler(time = seq(0, 50, length.out = sample))
  # data <- list(x=ts$z, dp=cbind(ts$x, ts$y))
} else {
  ### Real-world example
  data("obs.mon")
  data("rain.mon")

  if (TRUE) { # SPI12 as response
    SPI.12 <- SPEI::spi(rain.mon[, 5], scale = 12)$fitted
    x <- window(SPI.12, start = c(1950, 1), end = c(2009, 12))
    dp <- window(obs.mon, start = c(1950, 1), end = c(2009, 12))
  } else { # rainfall as response
    x <- window(rain.mon[, 5], start = c(1950, 1), end = c(2009, 12))
    dp <- window(obs.mon, start = c(1950, 1), end = c(2009, 12))
  }
  data <- list(x = x, dp = dp)
  sample <- length(x)
}

# plot.ts(cbind(data$x,data$dp))

tab.list <- list()
mode.opts <- c("MRA", "MODWT", "AT")
for (mode in mode.opts) {
  print(mode)

  # cov.opt <- switch(2,"auto","pos","neg")
  if (mode == "MRA") {
    method <- switch(1,"dwt", "modwt")
  }

  # wavelet family, extension mode and package
  # wf <- switch(mode, "MRA"="haar", "MODWT"="haar", "AT"="haar")
  wf <- "haar"
  pad <- "zero"
  boundary <- "periodic"
  if (wf != "haar") v <- as.integer(readr::parse_number(wf) / 2) else v <- 1

  # Maximum decomposition level J
  n <- sample
  J <- ceiling(log(n / (2 * v - 1)) / log(2)) - 1 # (Kaiser, 1994)

  tab <- NULL
  for (cov.opt in c("auto", "pos", "neg")[1]) {
    # variance transform - calibration
    if (mode == "MRA") {
      dwt <- dwt.vt(data, wf, J, method, pad, boundary, cov.opt)
    } else if (mode == "MODWT") {
      dwt <- modwt.vt(data, wf, J, boundary, cov.opt)
    } else {
      dwt <- at.vt(data, wf, J, boundary, cov.opt)
    }

    # optimal prediction accuracy
    opti.rmse <- NULL
    dp.RMSE <- NULL
    dp.n.RMSE <- NULL
    S <- dwt$S
    ndim <- ncol(S)
    for (i in 1:ndim) {
      x <- dwt$x
      dp <- dwt$dp[, i]
      dp.n <- dwt$dp.n[, i]

      # ts.plot(cbind(dp,dp.n), col=1:2)

      dp.RMSE <- c(dp.RMSE, sqrt(mean(lm(x ~ dp)$residuals^2)))
      dp.n.RMSE <- c(dp.n.RMSE, sqrt(mean(lm(x ~ dp.n)$residuals^2)))

      # small difference due to the reconstruction
      opti.rmse <- c(opti.rmse, sqrt((n - 1) / n * (var(x) - sum(S[, i]^2) * var(dp) / var(dp.n))))
      # opti.rmse <- c(opti.rmse, sqrt((n-1)/n*(var(x)-sum(S[,i]^2))))
    }

    tab <- rbind(tab, data.frame(cov.opt, var=1:ndim, dp.RMSE, dp.n.RMSE, opti.rmse))
  }

  colnames(tab) <- c("Sign of covariance", "Variable", "Std", "VT", "Optimal")
  tab.list[[length(tab.list) + 1]] <- tab
}

# print(tab.list)

kable(tab.list[[1]][,-1], caption = "Optimal RMSE using DWT-based VT",
      booktabs = T, align = "c", digits = 3) %>%
kable_styling("striped", position = "center", full_width = FALSE)  %>%
collapse_rows(columns = 1, valign = "middle")

Optimal RMSE using DWT-based VT
Variable	Std	VT	Optimal
1	0.976	0.899	0.899
2	0.937	0.819	0.819
3	0.957	0.901	0.901
4	0.988	0.933	0.933
5	0.984	0.902	0.902
6	0.983	0.892	0.892
7	0.987	0.917	0.917

kable(tab.list[[2]][,-1], caption = "Optimal RMSE using MODWT/AT-based VT",
      booktabs = T, align = "c", digits = 3) %>%
kable_styling("striped", position = "center", full_width = FALSE)  %>%
collapse_rows(columns = 1, valign = "middle")

Optimal RMSE using MODWT/AT-based VT
Variable	Std	VT	Optimal
1	0.976	0.909	0.909
2	0.937	0.785	0.785
3	0.957	0.893	0.893
4	0.988	0.954	0.954
5	0.984	0.941	0.941
6	0.983	0.824	0.824
7	0.987	0.957	0.957

Transformed predictor variables

#-------------------------------------------------------------------
if (FALSE) {
  set.seed(2020)
  ### synthetic example - Rossler
  sample <- 10000
  s <- 0.1
  ts.list <- list()
  for (i in seq_along(s)) {
    ts.r <- data.gen.Rossler(a = 0.2, b = 0.2, w = 5.7, start = c(-2, -10, 0.2), time = seq(0, 50, length.out = sample))

    # add noise
    ts.r$x <- ts(ts.r$x + rnorm(n = sample, mean = 0, sd = s[i]))
    ts.r$y <- ts(ts.r$y + rnorm(n = sample, mean = 0, sd = s[i]))
    ts.r$z <- ts(ts.r$z + rnorm(n = sample, mean = 0, sd = s[i]))

    ts.list[[i]] <- ts.r
  }

  data.list <- lapply(ts.list, function(ts) list(x = ts$z, dp = cbind(ts$x, ts$y)))

  lab.names <- c("x", "y")
  xlim<- c(0,n); ylim <- c(-55, 55)
} else {

  ### Real-world example
  data("obs.mon")
  data("rain.mon")

  SPI.12 <- SPEI::spi(rain.mon[, 5], scale = 12)$fitted
  x <- window(SPI.12, start = c(1950, 1), end = c(2009, 12))
  dp <- window(obs.mon, start = c(1950, 1), end = c(2009, 12))

  data.list <- list(list(x = x, dp = dp))
  sample <- length(x)

  lab.names <- colnames(obs.mon)
  xlim<- NULL; ylim <- NULL
}

#-------------------------------------------------------------------
p.list <- list()
dp.list <- list()
if (wf != "haar") mode.opts <- c("MRA", "MODWT", "AT")[1:3] else mode.opts <- c("MRA", "MODWT","AT")[1:2]

for (mode in mode.opts) {
  cov.opt <- switch(1,"auto","pos","neg")
  flag <- switch(1,"biased", "unbiased")
  
  if (mode == "MRA") {
    method <- switch(1,"dwt","modwt")
  }

  # wavelet family, extension mode and package
  # wf <- switch(mode, "MRA"="haar", "MODWT"="haar", "AT"="haar")
  wf <- "d16"
  pad <- "zero"
  boundary <- "periodic"
  if (wf != "haar") v <- as.integer(readr::parse_number(wf) / 2) else v <- 1

  # Maximum decomposition level J
  n <- sample
  J <- ceiling(log(n / (2 * v - 1)) / log(2)) - 1 # (Kaiser, 1994)
  # J <- floor(log(n/(2*v-1))/log(2))

  # variance transform - calibration
  if (mode == "MRA") {
    dwt.list <- lapply(data.list, function(x) dwt.vt(x, wf, J, method, pad, boundary, cov.opt, flag))
  } else if (mode == "MODWT") {
    dwt.list <- lapply(data.list, function(x) modwt.vt(x, wf, J, boundary, cov.opt, flag))
  } else {
    dwt.list <- lapply(data.list, function(x) at.vt(x, wf, J, boundary, cov.opt, flag))
  }

  for (j in 1:length(dwt.list)) {
    dwt <- dwt.list[[j]]

    par(
      mfrow = c(ncol(dwt$dp), 1), mar = c(0, 2.5, 2, 1),
      oma = c(2, 1, 0, 0), # move plot to the right and up
      mgp = c(1.5, 0.5, 0), # move axis labels closer to axis
      pty = "m", bg = "transparent",
      ps = 12
    )

    # plot(dwt$x, type="l", xlab=NA, ylab="SPI12", col="red")
    # plot(dwt$x, type="l", xlab=NA, ylab="Rain", col="red")
    for (i in 1:ncol(dwt$dp)) {
      ts.plot(cbind(dwt$dp[, i], dwt$dp.n[, i]),
        xlab = NA, ylab = paste0(lab.names[i]),
        xlim = xlim, ylim = ylim,
        col = c("black", "blue"), lwd = c(1, 2)
      )
    }

    p.list[[length(p.list) + 1]] <- recordPlot()

    dp.list[[length(dp.list) + 1]] <- dwt$dp.n
  }
}

#----------------------------------------------------------------------------
# plot and save
fig <- cowplot::plot_grid(plotlist = p.list, nrow = 1, labels = c("(a)", "(b)", "(c)"))
fig

Orignal and VT predictors. (a): DWT-MRA (b): MODWT/AT

Stepwise variance transformation

#-------------------------------------------------------------------
### Real-world example
data("obs.mon")
data("rain.mon")
op <- par()
station.id <- 5
lab.names <- colnames(obs.mon)[c(1, 3, 4, 5, 7)]

if (TRUE) { # SPI12 as response
  SPI.12 <- SPEI::spi(rain.mon, scale = 12)$fitted
  x <- window(SPI.12, start = c(1950, 1), end = c(2009, 12))
  dp <- window(obs.mon[, lab.names], start = c(1950, 1), end = c(2009, 12))
} else { # rainfall as response
  x <- window(rain.mon, start = c(1950, 1), end = c(2009, 12))
  dp <- window(obs.mon[, lab.names], start = c(1950, 1), end = c(2009, 12))
}

data.list <- lapply(station.id, function(id) list(x = x[, id], dp = dp))


ylim <- data.frame(
  GPH = c(700, 900), TDP700 = c(5, 25), TDP500 = c(5, 25), EPT = c(300, 330),
  UWND = c(-5, 25), VWND = c(-5, 10), MSLP = c(-1, 1)
)[c(1, 3, 4, 5, 7)]

#-------------------------------------------------------------------
p.list <- list()
RMSE <- NULL
mode.opts <- c("MRA", "MODWT", "AT")[1:2]
for (mode in mode.opts) {
  cov.opt <- switch(1,
    "auto",
    "pos",
    "neg"
  )
  if (mode == "MRA") {
    method <- switch(1,
      "dwt",
      "modwt"
    )
  }

  # wavelet family, extension mode and package
  wf <- switch(mode,"MRA" = "d4","MODWT" = "haar","AT" = "haar")
  pad <- "zero"
  boundary <- "periodic"
  if (wf != "haar") v <- as.integer(readr::parse_number(wf) / 2) else v <- 1

  # Maximum decomposition level J
  n <- nrow(x)
  J <- ceiling(log(n / (2 * v - 1)) / log(2)) - 1 # (Kaiser, 1994)

  # high order variance transformation
  dwt.list <- lapply(data.list, function(data) stepwise.VT(data, mode = mode, wf = wf, J=J))

  for (j in seq_len(length(dwt.list))) {
    dwt <- dwt.list[[j]]
    cpy <- dwt$cpy

    MSE <- NULL
    for (i in seq_len(length(cpy))) {
      m1 <- sqrt(FNN::knn.reg(train = dwt$dp[, 1:i], y = dwt$x)$PRESS / n)
      m2 <- sqrt(FNN::knn.reg(train = dwt$dp.n[, 1:i], y = dwt$x)$PRESS / n)

      MSE <- rbind(MSE, c(m1, m2))
    }

    RMSE <- cbind(RMSE, MSE)

    par(
      mfrow = c(length(cpy), 1), mar = c(0, 4, 2, 1),
      oma = c(2, 1, 0, 0), # move plot to the right and up
      mgp = c(1.5, 0.5, 0), # move axis labels closer to axis
      pty = "m", bg = "transparent",
      ps = 8
    )

    # plot(dwt$x, type="l", xlab=NA, ylab="SPI12", ylim=c(-3,3),col="red")
    # plot(dwt$x, type="l", xlab=NA, ylab="Rain", col="red")
    for (i in seq_len(length(cpy))) {
      ts.plot(dwt$dp[, i], dwt$dp.n[, i],
        xlab = NA, ylab = paste0(lab.names[cpy[i]]), # ylim=ylim[,i],
        col = c("black", "blue"), lwd = c(1, 2)
      )
    }

    p.list[[length(p.list) + 1]] <- recordPlot()
  }
}
#> [1] "Variance difference between transformed\n                        and original series by percentage: 17.2263230629905"
par(op)
#-------------------------------------------------------------------
# plot and save
cowplot::plot_grid(plotlist = p.list, nrow = 1, labels = c("(a)", "(b)", "(c)"))

Orignal and SVT predictors. (a): DWT-MRA (b): MODWT/AT

#-------------------------------------------------------------------
# RMSE when more predictors are included
tab1 <- round(RMSE, 3)
tab1 <- cbind(1:nrow(tab1), tab1)
colnames(tab1) <- c("No. of Predictors", rep(c("Original", "Transformed"), length(mode.opts)))

kable(tab1, caption = "Comparison of prediction accuracy using Std and SVT", booktabs = T) %>%
  kable_styling(latex_options = c("HOLD_position"), position = "center", full_width = FALSE)  %>%
  #  add_header_above(c(" " = 1, "DWT-MRA" = 2, "MODWT" = 2, "AT" = 2))
  add_header_above(c(" " = 1, "DWT-MRA" = 2, "MODWT/AT" = 2))

Comparison of prediction accuracy using Std and SVT
	DWT-MRA		MODWT/AT
No. of Predictors	Original	Transformed	Original	Transformed
1	1.114	1.007	1.112	0.995
2	1.093	0.876	1.100	0.934
3	1.089	0.743	1.050	0.821
4	1.086	0.744	1.053	0.824

Share on

Twitter Facebook LinkedIn 微博微信

Hailong Shi

Required packages

DWT, MODWT and AT basic propertites

Summary of various properties for the three DWT methods

Illustration of three types of DWT methods

Daubechies 16 wavelet

Haar wavelet filter

Optimal variance transformation

Preditive accuracy (RMSE)

Transformed predictor variables

Stepwise variance transformation

Share on