#' 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)))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.