R/spp.R

In bjw34032/waveslim: Basic Wavelet Routines for One-, Two-, and Three-Dimensional Signal Processing

#' Wavelet-based Maximum Likelihood Estimation for Seasonal Persistent
#' Processes
#' 
#' Parameter estimation for a seasonal persistent (seasonal long-memory)
#' process is performed via maximum likelihood on the wavelet coefficients.
#' 
#' The variance-covariance matrix of the original time series is approximated
#' by its wavelet-based equivalent.  A Whittle-type likelihood is then
#' constructed where the sums of squared wavelet coefficients are compared to
#' bandpass filtered version of the true spectral density function.
#' Minimization occurs for the fractional difference parameter \eqn{d} and the
#' Gegenbauer frequency \eqn{f_G}, while the innovations variance is
#' subsequently estimated.
#' 
#' @usage spp.mle(y, wf, J = log(length(y), 2) - 1, p = 0.01, frac = 1)
#' @usage spp2.mle(y, wf, J = log(length(y), 2) - 1, p = 0.01, dyadic = TRUE, frac = 1)
#' @aliases spp.mle spp2.mle
#' @param y Not necessarily dyadic length time series.
#' @param wf Name of the wavelet filter to use in the decomposition.  See
#' \code{\link{wave.filter}} for those wavelet filters available.
#' @param J Depth of the discrete wavelet packet transform.
#' @param p Level of significance for the white noise testing procedure.
#' @param dyadic Logical parameter indicating whether or not the original time
#' series is dyadic in length.
#' @param frac Fraction of the time series that should be used in constructing
#' the likelihood function.
#' @return List containing the maximum likelihood estimates (MLEs) of
#' \eqn{\delta}, \eqn{f_G} and \eqn{\sigma^2}, along with the value of the
#' likelihood for those estimates.
#' @author B. Whitcher
#' @seealso \code{\link{fdp.mle}}
#' @references Whitcher, B. (2004) Wavelet-based estimation for seasonal
#' long-memory processes, \emph{Technometrics}, \bold{46}, No. 2, 225-238.
#' @keywords ts
#' @export spp.mle
spp.mle <- function(y, wf, J=log(length(y),2)-1, p=0.01, frac=1)
{
  sppLL <- function(x, y) {
    delta <- x[1]
    fG <- x[2]

    ## cat("Parameters are: d =", delta, ", and f =", fG, fill=TRUE)
    
    y.dwpt <- y[[1]]
    y.basis <- y[[2]]
    n <- y[[3]]
    J <- y[[4]]
    
    ## Establish the limits of integration for the band-pass variances
    a <- unlist(apply(matrix(2^(1:J)-1), 1, seq, from=0, by=1)) / 
      2^(rep(1:J, 2^(1:J))) / 2
    b <- unlist(apply(matrix(2^(1:J)), 1, seq, from=1, by=1)) / 
      2^(rep(1:J, 2^(1:J))) / 2
    
    ## Define some useful parameters for the wavelet packet tree
                                        # n <- length(y)
    length.jn <- n / rep(2^(1:J), 2^(1:J))
    scale.jn <- rep(2^(1:J+1), 2^(1:J))
    
    ## Initialize various parameters for the reduced LL
    Basis <- (1:length(y.basis))[y.basis]
    bp.var <- numeric(length(Basis))
    delta.n <- 100

    ## Compute the band-pass variances according to \delta and f_G
    omega.diag <- NULL
    for(i in 1:sum(y.basis)) {
      jn <- Basis[i]
      bp.var[i] <- bandpass.spp(a[jn], b[jn], delta, fG)
      omega.diag <- c(omega.diag,
                      scale.jn[jn] * rep(bp.var[i], length.jn[jn]))
    }
    
    ## Compute reduced log-likelihood 
    rLL <- n * log(1/n * sum(y.dwpt^2 / omega.diag, na.rm=TRUE)) +
      sum(length.jn[y.basis] * log(scale.jn[y.basis] * bp.var))
    rLL
  }
  
  n <- length(y)
  x0 <- numeric(2)

  ## Perform discrete wavelet packet transform (DWPT) on Y
  y.dwpt <- dwpt(y, wf, n.levels=J)
  n <- length(y)
  if(frac < 1) {
    for(i in 1:length(y.dwpt)) {
      vec <- y.dwpt[[i]]
      ni <- length(vec)
      j <- rep(1:J, 2^(1:J))[i]
      vec[trunc(frac * n/2^j):ni] <- NA
      y.dwpt[[i]] <- vec
    }
  }
  y.basis <- as.logical(ortho.basis(portmanteau.test(y.dwpt, p)))
  y.dwpt <- as.matrix(unlist(y.dwpt[y.basis]))

  ## Compute initial estimate of the Gegenbauer frequency
  y.per <- per(y - mean(y))
  x0[2] <- (0:(n/2)/n)[max(y.per) == y.per]

  ## Compute initial estimate of the fractional difference parameter
  muJ <- (unlist(apply(matrix(2^(1:J)-1), 1, seq, from=0, by=1)) / 
          2^(rep(1:J, 2^(1:J))) + 
          unlist(apply(matrix(2^(1:J)), 1, seq, from=1, by=1)) /  
          2^(rep(1:J, 2^(1:J)))) / 4
  y.modwpt <- modwpt(y, wf=wf, n.levels=J)
  y.varJ <- rep(2^(1:J), 2^(1:J)) *
    unlist(lapply(y.modwpt,
                  FUN=function(x)sum(x*x,na.rm=TRUE)/length(x[!is.na(x)])))
  x0[1] <- min(-0.5 * lsfit(log(abs(muJ[y.basis] - x0[2])), 
                            log(y.varJ[y.basis]))$coef[2], 0.49)

  cat(paste("Initial parameters are: delta =", round(x0[1],4), 
            "freqG =", round(x0[2],4), "\n"))
  result <- optim(par=x0, fn=sppLL, method="L-BFGS-B",
                  lower=c(0.001,0.001), upper=c(0.499,0.499),
                  control=list(trace=0, fnscale=2),
                  y=list(y.dwpt, y.basis, n, J))
  return(result)
}

spp2.mle <- function(y, wf, J=log(length(y),2)-1, p=0.01,
                     dyadic=TRUE, frac=1)
{

  spp2LL <- function(x, y) {
    d1 <- x[1]
    f1 <- x[2]
    d2 <- x[3]
    f2 <- x[4]
    ## cat("Parameters are: d1 =", round(d1,6), ", and f1 =", round(f1,6),
    ##     ", d2 =", round(d2,6), ", and f2 =", round(f2,6), fill=TRUE)

    y.dwpt <- y[[1]]
    y.basis <- y[[2]]
    n <- y[[3]]
    J <- y[[4]]

    ## Establish the limits of integration for the band-pass variances
    a <- unlist(apply(matrix(2^(1:J)-1), 1, seq, from=0, by=1)) / 
      2^(rep(1:J, 2^(1:J))) / 2
    b <- unlist(apply(matrix(2^(1:J)), 1, seq, from=1, by=1)) / 
      2^(rep(1:J, 2^(1:J))) / 2
    
    ## Define some useful parameters for the wavelet packet tree
    length.jn <- n / rep(2^(1:J), 2^(1:J))
    scale.jn <- rep(2^(1:J+1), 2^(1:J))
    
    ## Initialize various parameters for the reduced LL
    Basis <- (1:length(y.basis))[y.basis]
    bp.var <- numeric(length(Basis))
    delta.n <- 100

    ## Compute the band-pass variances according to \delta and f_G
    omega.diag <- NULL
    for(i in 1:sum(y.basis)) {
      jn <- Basis[i]
      bp.var[i] <- bandpass.spp2(a[jn], b[jn], d1, f1, d2, f2)
      omega.diag <- c(omega.diag,
                      scale.jn[jn] * rep(bp.var[i], length.jn[jn]))
    }
    
    ## Compute reduced log-likelihood 
    n * log(1/n * sum(y.dwpt^2 / omega.diag, na.rm=TRUE)) +
      sum(length.jn[y.basis] * log(scale.jn[y.basis] * bp.var), na.rm=TRUE)
  }

  n <- length(y)
  x0 <- numeric(4)

  ## Perform discrete wavelet packet transform (DWPT) on Y
  y.dwpt <- dwpt(y, wf, n.levels=J)
  if(!dyadic) {
    for(i in 1:length(y.dwpt)) {
      vec <- y.dwpt[[i]]
      ni <- length(vec)
      j <- rep(1:J, 2^(1:J))[i]
      vec[trunc(frac * n/2^j):ni] <- NA
      y.dwpt[[i]] <- vec
    }
  }
  y.basis <- as.logical(ortho.basis(portmanteau.test(y.dwpt, p, type="other")))
  y.dwpt <- as.vector(unlist(y.dwpt[y.basis]))

  ## Compute initial estimate of the Gegenbauer frequencies
  if(dyadic)
    y.per <- per(y - mean(y))
  else
    y.per <- per(y[1:(frac*n)] - mean(y[1:(frac*n)]))
  freq.y <- (0:(frac*n %/% 2))/(frac*n)
  x0[2] <- freq.y[max(y.per) == y.per]
  x0[4] <- freq.y[max(y.per[freq.y > x0[2] + freq.y[10] | 
    freq.y < x0[2] - freq.y[10]]) == y.per]
  if(x0[2] > x0[4]) {
    xx <- x0[2]
    x0[2] <- x0[4]
    x0[4] <- xx
    rm(xx)
  }

  ## Compute initial estimate of the fractional difference parameters
  muJ <- (unlist(apply(matrix(2^(1:J)-1), 1, seq, from=0, by=1)) / 
      2^(rep(1:J, 2^(1:J))) + 
      unlist(apply(matrix(2^(1:J)), 1, seq, from=1, by=1)) /  
      2^(rep(1:J, 2^(1:J)))) / 4
  y.modwpt <- modwpt(y, wf=wf, n.levels=J)
  y.varJ <- rep(2^(1:J), 2^(1:J)) *
    unlist(lapply(y.modwpt,
                  FUN = function(x) sum(x*x,na.rm=TRUE)/length(x[!is.na(x)])))
  x0.mid <- (x0[2] + x0[4]) / 2
  muJ <- muJ[y.basis]
  y.varJ <- y.varJ[y.basis]
  x0[1] <- min(-0.5 * lsfit(log(abs(muJ[muJ < x0.mid] - x0[2])), 
                            log(y.varJ[muJ < x0.mid]))$coef[2], 0.49)
  x0[3] <- min(-0.5 * lsfit(log(abs(muJ[muJ > x0.mid] - x0[4])), 
                            log(y.varJ[muJ > x0.mid]))$coef[2], 0.49)

  cat(paste("Initial parameters: d1 = ", round(x0[1],4), 
            ", f1 = ", round(x0[2],4), ", d2 = ", round(x0[3],4), 
            ", f2 = ", round(x0[4],4), sep=""), fill=TRUE)
  result <- optim(par=x0, fn=spp2LL, method="L-BFGS-B",
                  lower=rep(0.001,4), upper=rep(0.499,4),
                  control=list(trace=1, fnscale=2),
                  y=list(y.dwpt, y.basis, n, J))
  return(result)
}