Estimate the process (cross) spectral density function via nonparametric models.
SDF(x, method="direct", taper.=NULL, window=NULL,
n.taper=5, overlap=0.5, blocksize=NULL,
single.sided=TRUE, sampling.interval=NULL,
center=TRUE, recenter=FALSE, npad=2*numRows(x))a vector or matrix containing uniformly-sampled real-valued time series.
If a matrix, each column should contain a different time series.
an integer representing the number
of points (width) of each block in the WOSA estimator scheme.
Default: floor(N/4) where N is the number of
samples in each series.
a logical value. If TRUE, the mean of each
time series is recentered prior to estimating
the SDF. Default: TRUE.
a character string denoting the method to use in estimating the SDF.
Choices are "direct", "lag window", "wosa" (Welch's Overlapped Segment Averaging),
"multitaper". See DETAILS for more information. Default: "direct".
an integer defining the number of tapers
to use in a multitaper scheme. This value is overwritten if
the taper input is of class taper. Default: 5.
an integer representing the total length of each
time series to analyze after padding with zeros. This argument
allows the user to control the spectral resolution of the SDF
estimates: the normalized frequency interval is
\(\Delta f = 1 / \hbox{npad}\).
This argument must be set such that
\(\hbox{npad} > 2\).
Default: 2*numRows(x).
a numeric value on \([0,1]\) denoting the fraction
of window overlap for the WOSA estimator. Default: 0.5.
a logical value. If TRUE, the mean of each
time series is recentered after (posssibly) tapering the series prior to estimating
the SDF. Default: FALSE.
a numeric value representing the interval
between samples in the input time series x. Default: NULL, which
serves as a flag to obtain the sampling interval via the deltat
function. If x is a list, the default sampling interval is deltat(x[[1]]).
If x is an atomic vector (ala isVectorAtomic), then the default
samplign interval is established ala deltat(x). Finally, if the
input series is a matrix, the sampling interval of the first series (assumed to
be in the first column) is obtained ala deltat(x[,1]).
a logical value. If TRUE, a single-sided
SDF estimate is returned corresponding to the normalized frequency
range of \([0,1/2]\). Otherwise, a double-sided SDF estimate
corresponding to the normalized frequency interval \([-1/2,1/2]\)
is returned. Default: TRUE.
an object of class taper or a character string denoting the primary taper.
If an object of class taper, the length of the taper is checked to ensure
compatitbility with the input x.
See DETAILS for more information. The default values are a function
of the method as follows:
normalized rectangular taper
normalized Parzen window with a cutoff at \(N/2\) where \(N\) is the length of the time series.
normalized Hanning taper
normalized Hanning taper
an object of class taper or a character string denoting the (secondary) window
for the lag window estimator.
If an object of class taper, the length of the taper is checked to ensure
compatitbility with the input x. See DETAILS for more information.
Default: Normalized Hanning window.
an object of class SDF.
converts the (cross-)SDF estimate(s) as a matrix. Optional arguments
are passed directly to the matrix function during the conversion.
plots the (cross-)SDF estimate(s). Optional arguments are:
a character string defining the scaling to perform on the (common) frequency vector
of the SDF estimates. See the scaleData function for supported choices. Default: "linear".
a character string defining the scaling to perform on the SDF estimates.
See the scaleData function for supported choices. Default: "linear".
a single character defining the plot type (ala the par function) of the
SDF plots. Default: ifelse(numRows(x) > 100, "l", "h").
a character string representing the x-axis label. Default: "FREQUENCY (Hz)".
a (vector of) character string(s), one per (cross-)SDF estimate,
representing the y-axis label(s). Default: in the multivariate case, the strings
"Sij" are used for the y-axis labels, where i and j are the indices of the
different variables. For example, if the user supplies a 2-column matrix for x,
the labels "S11", "S12", and "S22" are used to label the y-axes of the corresponding
(cross-)SDF plots. In the univariate case, the default string "SDF" prepended with a string
describing the type of SDF performed (such as "Multitaper") is used to label the y-axis.
a logical value. If TRUE, the SDF value at normalized frequency \(f=0\)
is plotted for each SDF. This frequency is associated with the sample mean
of the corresponding time series. A relatively large mean value dominates
the spectral patterns in a plot and thus the corresponding frequency is typically not plotted.
Default: !attr(x,"center").
an integer defining the maximum number of SDF plots to place onto a single graph.
Default: 3.
a post processing function to apply to the SDF values prior to plotting. Supported
functions are Mod, Im, Re and Arg. See each of these functions for details.
If the SDF is purely real (no cross-SDF is calculated), this argument is coerced to the Mod function.
Default: Mod.
A logical value. If TRUE, the plot is added using the current
par() layout. Otherwise a new plot is produced. Default: FALSE.
additional plot parameters passed directly to the genPlot function used
to plot the SDF estimates.
prints the object. Available options are:
text justification ala prettPrintList. Default: "left".
header separator ala prettyPrintList. Default: ":".
Additional print arguments sent directly to the prettyPrintList function.
Let \(X_t\) be a uniformly sampled real-valued time series of length \(N\), Let an estimate of the process spectral density function be denoted as \(\hat{S}_X(f)\) where \(f\) are frequencies on the interval \([-1/(2\Delta t),1/(2\Delta t)]\) where \(\Delta t\) is the sampling interval. The supported SDF estimators are:
The direct SDF estimator is defined as
\(\hat{S}_X^{(d)}(f) = | \sum_{t=0}^{N-1} h_t X_t e^{-i2\pi f t}|^2\),
where \(\{h_t\}\) is a data taper normalized such that
\(\sum_{t=0}^{N-1} h_t^2 = 1\). If
\(h_t=1/\sqrt{N}\) then we obtain the definition
of the periodogram
\(\hat{S}_X^{(p)}(f) = \frac{1}{N} | \sum_{t=0}^{N-1} X_t e^{-i2\pi f t}|^2\).
See the taper function for more details on supported window types.
The lag window SDF estimator is defined as
\(\hat{S}_X^{(lw)}(f) = \sum_{\tau=-(N-1)}^{N-1} w_\tau
\hat{s}_{X,\tau}^{(d)} e^{-i2\pi f \tau}\),
where \(\hat{s}_{X,\tau}^{(d)}\) is the autocovariance
sequence estimator corresponding to some
direct spectral estimator (often the periodogram)
and \(w_\tau\) is a lag window (popular choices
are the Parzen, Papoulis, and Daniell windows). See the
taper function for more details.
Welch's Overlapped Segment Averaging SDF estimator is defined as
$$ \hat S^{(wosa)} = {1\over N_B} \sum_{j=0}^{N_B-1} \hat S^{(d)}_{jN_O} (f) $$ where $$ \hat S^{(d)}_{l}(f) \equiv \left| \sum_{t=0}^{N_S-1} h_t X_{t+l} e^{-i2\pi ft} \right|^2, \enskip 0 \le l \le N - N_S; $$ Here, \(N_O\) is a positive integer that controls how much overlap there is between segments and that must satisfy both \(N_O \le N_S\) and \(N_O(N_B-1) = N-N_S\), while \(\{ h_t \}\) is a data taper appropriate for a series of length \(N_S\) (i.e., \(\sum_{t=0}^{N_S-1} h_t^2 = 1\)).
A multitaper spectral estimator is given by
$$\hat S^{(mt)}_X(f)=
{1\over K}
\sum_{k=0}^{K-1}
\left| \sum_{t=0}^{N-1} h_{k,t} X_t e^{-i2\pi ft} \right|^2,
$$
where
\(S(k,f) = {|\sum_{t=0}^{N-1} h_{k,t} X_t \exp(-i 2 \pi f t)|}^2\)
and \(\{ h_{k,t} \}$, $k=0,\ldots,K-1\),
is a set of \(K\) orthonormal data tapers.
$$\sum_{t=0}^{N-1} h_{k,t} h_{k',t}
=
\left\{
\begin{array}{ll}
1,& \mbox{if $k=k'$;}\\
0,& \mbox{otherwise}
\end{array}
\right.
$$
Popular choices for multitapers include sinusoidal tapers and
discrete prolate spheroidal sequences (DPSS). See the
taper function for more details.
Cross spectral density function estimation:
If the input x is a matrix, where each column contains
a different time series, then the results are returned in
a matrix whose columns correspond to all possible unique combinations
of cross-SDF estimates. For example, if x has three columns,
then the output will be a matrix whose columns are
\(\{S_{11},S_{12},S_{13},S_{22},S_{23},S_{33}\}\)
where \(S_{ij}\) is the cross-SDF estimate of the ith
and jth column of x. All cross-spectral density
function estimates are returned as complex-valued series to maintain
the phase relationships between components. For all \(S_{ij}\)
where \(i=j\), however, the imaginary portions will be zero (up to a
numerical noise limit).
Percival, Donald B. and Constantine, William L. B. (2005) ``Exact Simulation of Gaussian Time Series from Nonparametric Spectral Estimates with Application to Bootstrapping", Journal of Computational and Graphical Statistics, accepted for publication.
D.B. Percival and A. Walden (1993), Spectral Analysis for Physical Applications: Multitaper and Conventional Univariate Techniques, Cambridge University Press, Cambridge, UK.
# NOT RUN {
## calculate various SDF estimates for the
## sunspots series. remove mean component for a
## better comparison.
require(ifultools)
data <- as.numeric(sunspots)
methods <- c("direct","wosa","multitaper",
"lag window")
S <- lapply(methods, function(x, data) SDF(data, method=x), data)
x <- attr(S[[1]], "frequency")[-1]
y <- lapply(S,function(x) decibel(as.vector(x)[-1]))
names(y) <- methods
## create a stack plot of the data
stackPlot(x, y, col=1:4)
## calculate the cross-spectrum of the same
## series: all spectra should be the same in
## this case
SDF(cbind(data,data), method="lag")
## calculate the SDF using npad=31
SDF(data, npad=31, method="multitaper")
# }
Run the code above in your browser using DataLab