R/func_whiten.R

In psd: Adaptive, Sine-Multitaper Power Spectral Density and Cross Spectrum Estimation

#' Prepare a series for spectral estimation
#' 
#' @description
#' Remove (optionally) mean, trend, and Auto Regressive (AR) model
#' from the original series.
#'
#' @details
#' The R-S multitapers do not exhibit the remarkable spectral-leakage 
#' suppression properties of the Thomson prolate tapers, 
#' so that in spectra with large dynamic range, 
#' power bleeds from the strong peaks into neighboring frequency bands of 
#' low amplitude -- spectral leakage. 
#' Prewhitening can ameliorate the problem, at least for red spectra 
#' [see Chapter 9, Percival and Walden (1993)]. 
#'
#' The value of the \code{AR.max} argument is made absolute, after which
#' this function has essentially two modes of operation (detailed below):
#' \describe{
#' \item{\code{AR.max} == 0}{Remove (optionally) a mean and/or linear trend.}
#' \item{\code{AR.max} > 0}{Remove an autoregressive model}
#' }
#' In the second case,
#' the time series is 
#' filtered in the time domain with a finite-impulse-response 
#' filter of \code{AR.max} terms. The filter is found by solving the Yule-Walker 
#' equations for 
#' which it is assumed the series was generated by an autoregressive process, up to
#' order \code{AR.max}.
# Next the multitaper spectral estimates are made on the filtered series. 
# Finally, the effect of the filter is removed by dividing out the spectrum of the 
# autoregressive process.
#
#' \subsection{Mean and trend (\code{AR.max == 0})}{
#'
#' Power spectral density estimates can become badly biased
#' (especially at lower frequencies) if a signal of the form
#' \eqn{f(x) = A x + B} is not removed from the series.  
#' If \code{detrend=TRUE} a model of this form is removed over the entire series using a
#' linear least-squares estimator; in this case a mean value is removed
#' regardless of the logical state of \code{demean}.
#' To remove \emph{only} a mean value, set \code{detrend=FALSE} and (obviously) \code{demean=TRUE}.
#' 
#' }
#'
#' \subsection{Auto Regressive (AR) innovations (\code{AR.max > 0})}{
#'
#' When an autoregressive model is removed from a non-stationary series, the residuals
#' are known as 'innovations', and may be stationary (or very-nearly stationary).  
#' This function fits an AR model [order at least 1, but up to and including AR(\code{AR.max})] to the series 
#' by solving the Yule-Walker equations; however, AIC is used to estimate the highest significant
#' order, which means that higher-order components may not necessarily be fit.
#' The resulting innovations can be used to better estimate the stationary component
#' of the original signal, and possibly in an interactive editing method.
#'
#' Note that the method used here--solving the Yule-Walker equations--is 
#' not a true maximum likelihood estimator; hence the AIC is calculated
#' based on the variance estimate (no determinant). From \code{?ar}:
#' \emph{In \code{ar.yw} the variance matrix of the innovations is 
#' computed from the fitted coefficients and the autocovariance of \code{x}.}
#'
#' A quick way to determine whether this may be needed for the series is to run
#' \code{acf} on the series, and see if significant non-zero lag correlations
#' are found.  A warning is produced if the fit returns an AR(0) fit, indicating
#' that AR prewhitening most likely inappropriate for the series, which
#' is apparently stationary (or very nearly so).  (The innovations could end up
#' having \emph{higher} variance than the input series in such a case.)
#'
#' \emph{Note that \code{AR.max} is restricted to the range \eqn{[1,N-1]} where
#' \eqn{N} is the series length.}
#'
#' }
#'
#' @section NA values:
#'
#' \code{NA} values are allowed.  If present, and \code{impute=TRUE}, 
#' the \code{na.locf} function in the package
#' \code{zoo} is used twice (with and without \code{fromLast} so that lead and
#' trailing \code{NA} values are also imputed).  The function name is an
#' acronym for "Last Observation Carried Forward", a very crude method
#' of imputation. 
#'
#' @name prewhiten
#' @export
#' @author A.J. Barbour and Robert L. Parker
#' @seealso \code{\link{psdcore}}, \code{\link{pspectrum}}
#' 
#' @param tser  vector; An object to prewhiten.
#' @param AR.max numeric; the maximum AR order to fit.
#' @param detrend  logical; Should a trend (and mean) be removed?
#' @param demean  logical; Should a mean value be removed?
#' @param impute  logical; Should NA values be imputed?
#' @param plot  logical; Should the results be plotted?
#' @param verbose  logical; Should messages be printed?
#' @param x.fsamp sampling frequency (for non \code{ts} objects)
#' @param x.start start time of observations (for non \code{ts} objects)
#' @param ... variables passed to \code{prewhiten.ts} (for non \code{ts} objects)
#' 
#' @return A list with the model fits (\code{lm} and \code{ar} objects),
#' the linear and AR prewhitened series (\code{ts} objects), and a logical
#' flag indicating whether the I/O has been imputed. This list includes:
#' \code{"lmdfit"}, \code{"ardfit"}, \code{"prew_lm"}, \code{"prew_ar"}, and \code{"imputed"}
#' @return \emph{Note that if \code{AR.max=0} the AR information will exist as \code{NULL}.}
#'
#' @example inst/Examples/rdex_prewhiten.R
prewhiten <- function(tser, ...) UseMethod("prewhiten")

#' @rdname prewhiten
#' @export
prewhiten.default <- function(tser, x.fsamp=1, x.start=c(1, 1), ...){
  Xts <- stats::ts(tser, frequency=x.fsamp, start=x.start)
  prewhiten.ts(Xts, ...)
}

#' @rdname prewhiten
#' @export
prewhiten.ts <- function(tser, AR.max=0L, detrend=TRUE, demean=TRUE, impute=TRUE, plot=TRUE, verbose=TRUE, ...){
  
  # prelims
  stopifnot(is.ts(tser))
  sps <- stats::frequency(tser)
  tstart <- stats::start(tser)
  n.o <- NROW(tser)
  ttime <- sps*n.o
  
  # NA action on input series
  NAFUN <- function(tso){
    stopifnot(is.ts(tso))
    # twice, to catch leading or trailing NA
    if (any(is.na(tso))){
      # forward l-o-c-f, then reverse
      #as.ts(zoo::na.locf(zoo::na.locf(as.zoo(tso), na.rm=FALSE), fromLast=TRUE, na.rm=FALSE))
      tso[] <- na_locf(tso)
    } else {
      tso
    }
  }
  
  tser <- NAFUN(tser)
  
  lmdfit <- ardfit <- NULL
  tser_prew_lm <- tser_prew_ar <- NULL
  
  ## Mean and Trend
  AR.max <- abs(AR.max)
  if (AR.max >= 0){
    

    X <- if (detrend){
      if (verbose) message("\tdetrending (and demeaning)")
      as.matrix(stats::residuals(stats::lm(tser~base::seq_len(n.o))))
      
    } else if (demean) {
      if (verbose) message("\tdemeaning")
      as.matrix(stats::residuals(stats::lm(tser~rep.int(1, n.o))))
      
    } else {
      if (verbose) message("\tnothing was done to the timeseries object")
      tser
    }
    
    # TS object of residuals or equivalent
    tser_prew_lm <- stats::ts(as.matrix(X), frequency=sps, start=tstart)
    
  }
  
  ## AR innovations
  if (AR.max >= 1) {
    
    # AR.max cannot be greater than N-1 in length
    AR.max <- as.integer(max(1, min(n.o-1, AR.max)))
    if (verbose){ message("\tautoregressive model fit (returning innovations)") }
    
    # solve the Yule-Walker equations
    ardfit <- stats::ar.yw(tser_prew_lm, aic=TRUE, order.max=AR.max, demean=TRUE)

    # returns a TS object
    tser_prew_ar <- ts(as.matrix(ardfit[['resid']]), frequency=sps, start=tstart)
    if (impute) tser_prew_ar <- NAFUN(tser_prew_ar)
    
  }
  
  if (plot){
    
    opar <- graphics::par(no.readonly = TRUE)
    on.exit(graphics::par(opar))
    graphics::par(las=0, xpd=FALSE)
    
    pltprew <- if (!is.null(ardfit)){
      # AR + linear
      nAR <- ardfit[['order']]
      ftyp <- sprintf("AR(%i) model fit", nAR)
      if (nAR==0) warning("AR(0) was the highest model found!")
      stats::ts.union(x=tser, x_lm=tser_prew_lm, x_ar=tser_prew_ar)
    } else if (is.null(lmdfit)) {
      # nothing
      ftyp <- ""
      tser
    } else {
      # Linear only
      ftyp <- "linear fit"
      stats::ts.union(x=tser, x_lm=tser_prew_lm)
    }
    
    PANELFUN <- function(x, col = col, bg = bg, pch = pch, type = type, ...){
      graphics::lines(x, col = col, bg = bg, pch = pch, type = type, ...)
      graphics::abline(h=c(0, mean(x)), lty=c(1,3), lwd=2, col=c("dark grey","red"))
    }
    
    plot(pltprew, 
         main="Raw and prewhitened series",
         cex.lab=0.7, cex.axis=0.7, xy.labels=FALSE,
         xaxs="i", yax.flip=TRUE, panel=PANELFUN)
    graphics::mtext(sprintf("%s  (demean %s | detrend %s )", ftyp, demean, detrend), line=0.7, cex=0.7)
    
  }
  
  return(invisible(list(lmdfit=lmdfit, 
                        ardfit=ardfit, 
                        prew_lm=tser_prew_lm, 
                        prew_ar=tser_prew_ar, 
                        imputed=impute)))
}