Title: | Spectral Density Estimation and Comparison for Functional Time Series |
---|---|
Description: | Functions for estimating spectral density operator of functional time series (FTS) and comparing the spectral density operator of two functional time series, in a way that allows detection of differences of the spectral density operator in frequencies and along the curve length. |
Authors: | Shahin Tavakoli [aut, cre] |
Maintainer: | Shahin Tavakoli <[email protected]> |
License: | GPL-2 |
Version: | 1.0.0 |
Built: | 2024-12-21 06:28:03 UTC |
Source: | CRAN |
The Epanechnikov weight function, with support in
Epanechnikov_kernel(x)
Epanechnikov_kernel(x)
x |
argument at which the function is evaluated |
ftsspec: collection of functions for estimating spectral density operator of functional time series (FTS) and comparing the spectral density operator of two functional time series, in a way that allows detection of differences of the spectral density operator in frequencies and along the curve length.
Tavakoli, Shahin and Panaretos, Victor M. "Detecting and Localizing Differences in Functional Time Series Dynamics: A Case Study in Molecular Biophysics", 2014, under revision
Generate the Filter of a multivariate MA process
Generate_filterMA(d.ts, d.n, MA.len = 3, ma.scale = rep(1, MA.len), a.smooth.coef = 0, seed = 1)
Generate_filterMA(d.ts, d.n, MA.len = 3, ma.scale = rep(1, MA.len), a.smooth.coef = 0, seed = 1)
d.ts |
dimension of the (output) time series |
d.n |
dimension of the noise that is filtered |
MA.len |
Length of the filter. Set to 3 by default. |
ma.scale |
scaling factor of each lag matrix. See details. |
a.smooth.coef |
A coefficient to shrink coefficients of filter. Set to 0 by default. |
seed |
The random seed used to generate the filter. Set to 1 by default. |
A d.ts x d.n x MA.len
array
Generates a filter (i.e. a d.ts x d.n x MA.len
array) for a moving
average process. The entries of the filter are generate randomly, but can be
reproduced by specifying the random seed seed
.
The ma.scale
parameter should be a vector of length MA.len
,
and corresponds to a scaling factor applied to each lag of the filter of the
MA process that is generated.
ma.scale1=c(-1.4,2.3,-2) a1=Generate_filterMA(10, 10, MA.len=3, ma.scale=ma.scale1, seed=10) str(a1) rm(a1)
ma.scale1=c(-1.4,2.3,-2) a1=Generate_filterMA(10, 10, MA.len=3, ma.scale=ma.scale1, seed=10) str(a1) rm(a1)
Get the square root of the covariance matrix associated to a noise type
Get_noise_sd(noise.type, d.n)
Get_noise_sd(noise.type, d.n)
noise.type |
the type of noise that is driving the MA process. See Details section. |
d.n |
dimension of the noise that is filtered |
SampleSpecDiffFreq
classPlotting function for SampleSpecDiffFreq
class
## S3 method for class 'SampleSpecDiffFreq' lines(x, method = NA, Kmax = 4, pch = 20, ...)
## S3 method for class 'SampleSpecDiffFreq' lines(x, method = NA, Kmax = 4, pch = 20, ...)
x |
object of the class |
method |
method used to adjust p-values |
Kmax |
maximum number of levels K for which the pvalues are plotted (used only if autok==0) |
pch |
the plot character to be used |
... |
additional parameters to be passed to plot() |
Compute the marginal p-values at each basis coefficients of for testing the equality of two spectral density kernels
Marginal_basis_pval(spec1, spec2, m, kappa.square, is.pi.multiple)
Marginal_basis_pval(spec1, spec2, m, kappa.square, is.pi.multiple)
spec1 |
The two sample spectral densities (at the same frequency |
spec2 |
The two sample spectral densities (at the same frequency |
m |
The number of Fourier frequencies over which the periodogram operator was smoothed. |
kappa.square |
the L2-norm of the weight function used to estimate the spectral density operator |
is.pi.multiple |
A logical variable, to specify if |
SampleSpec
Plotting method for object inheriting from class SampleSpec
## S3 method for class 'SampleSpec' plot(x, ...)
## S3 method for class 'SampleSpec' plot(x, ...)
x |
An object of the class |
... |
additional parameters to be passed to plot() |
SampleSpecDiffFreq
classPlotting function for SampleSpecDiffFreq
class
## S3 method for class 'SampleSpecDiffFreq' plot(x, method = NA, Kmax = 4, pch = 20, ...)
## S3 method for class 'SampleSpecDiffFreq' plot(x, method = NA, Kmax = 4, pch = 20, ...)
x |
object of the class |
method |
method used to adjust p-values |
Kmax |
maximum number of levels K for which the pvalues are plotted (used only if autok==0) |
pch |
the plot character to be used |
... |
additional parameters to be passed to plot() |
SampleSpecDiffFreqCurvelength
Plotting method for class SampleSpecDiffFreqCurvelength
## S3 method for class 'SampleSpecDiffFreqCurvelength' plot(x, ncolumns = 3, ...)
## S3 method for class 'SampleSpecDiffFreqCurvelength' plot(x, ncolumns = 3, ...)
x |
Object of the class |
ncolumns |
number of columns for the plots |
... |
additional parameters to be passed to |
SpecMA
Plotting method for object inheriting from class SpecMA
## S3 method for class 'SpecMA' plot(x, ...)
## S3 method for class 'SpecMA' plot(x, ...)
x |
A object of the class |
... |
additional parameters to be passed to plot() |
SampleSpecDiffFreqCurvelength
Printing method for class SampleSpecDiffFreqCurvelength
## S3 method for class 'SampleSpecDiffFreqCurvelength' print(x, ...)
## S3 method for class 'SampleSpecDiffFreqCurvelength' print(x, ...)
x |
Object of the class |
... |
Additional arguments for print |
Generic function to adjust pvalues
function to adjust pvalues for class SampleSpecDiffFreq
PvalAdjust(sample.spec.diff, method) ## S3 method for class 'SampleSpecDiffFreq' PvalAdjust(sample.spec.diff, method)
PvalAdjust(sample.spec.diff, method) ## S3 method for class 'SampleSpecDiffFreq' PvalAdjust(sample.spec.diff, method)
sample.spec.diff |
Object of the class
|
method |
method used to adjust p-values |
Simulate a new Moving Average (MA) vector time series and return the time series
Simulate_new_MA(a, T.len, noise.type, DEBUG = FALSE)
Simulate_new_MA(a, T.len, noise.type, DEBUG = FALSE)
a |
Array, returned by |
T.len |
Numeric, the length of the time series to generate |
noise.type |
the type of noise that is driving the MA process. See Details section. |
DEBUG |
Logical, for outputting information on the progress of the function |
A T.len x dim(a)[1]
matrix, where each column corresponds to a
coordinate of the vector time series
The function simulates a moving average process of dimension
dim(a)[1]
, defined by
noise.type
specifies the nature and internal correlation of the noise
that is driving the MA process. It can take the values
white-noise
the noise is Gaussian with covariance matrix identity
white-noise
the noise is Gaussian with diagonal covariance
matrix, whose j-th diagonal entry is
studentk
the coordinates of the noise are independent and have a student t distribution with 'k' degrees of freedom, standardized to have variance 1
ma.scale1=c(-1.4,2.3,-2) a1=Generate_filterMA(6, 6, MA.len=3, ma.scale=ma.scale1) X=Simulate_new_MA(a1, T.len=512, noise.type='wiener') plot.ts(X)
ma.scale1=c(-1.4,2.3,-2) a1=Generate_filterMA(6, 6, MA.len=3, ma.scale=ma.scale1) X=Simulate_new_MA(a1, T.len=512, noise.type='wiener') plot.ts(X)
This function estimates the spectral density operator of a Functional Time Series (FTS)
Spec(X, W = Epanechnikov_kernel, B.T = (dim(X)[1])^(-1/5), only.diag = FALSE, trace = FALSE, demean = TRUE, subgrid = FALSE, subgrid.density = 10, verbose = 0, subgrid.density.relative.to.bandwidth = TRUE)
Spec(X, W = Epanechnikov_kernel, B.T = (dim(X)[1])^(-1/5), only.diag = FALSE, trace = FALSE, demean = TRUE, subgrid = FALSE, subgrid.density = 10, verbose = 0, subgrid.density.relative.to.bandwidth = TRUE)
X |
A |
W |
The weight function used to smooth the periodogram operator. Set by default to be the Epanechnikov kernel |
B.T |
The bandwidth of frequencies over which the periodogram operator
is smoothed. If |
only.diag |
A logical variable to choose if the function only computes
the marginal spectral density of each basis coordinate
( |
trace |
A logical variable to choose if only the trace of the spectral
density operator is computed. |
demean |
A logical variable to choose if the FTS is centered before computing its spectral density operator. |
subgrid |
A logical variable to choose if the spectral density operator
is only returned for a subgrid of the Fourier frequencies, which can be
useful in large datasets to reduce memory usage. |
subgrid.density |
Only used if |
verbose |
A variable to show the progress of the computations. By
default, |
subgrid.density.relative.to.bandwidth |
logical parameter to specify if
|
A list containing the following elements:
The estimated spectral density operator. The first dimension corresponds to the different frequencies over which the spectral density operators are estimated.
The frequencies over which the spectral density is estimated.
The number of Fourier frequencies over which the periodogram operator was smoothed.
The equivalent Bandwidth used in the weight function W(), as defined in Bloomfield (1976, p.201).
The weight function used to smooth the periodogram operator.
The L2 norm of the weight function W.
spec.pgram function of R.
Bloomfield, P. (1976) "Fourier Analysis of Time Series: An Introduction", Wiley.
Panaretos, V. M. and Tavakoli, S., "Fourier Analysis of Functional Time Series", Ann. Statist. Volume 41, Number 2 (2013), 568-603.
ma.scale1=c(-1.4,2.3,-2) a1=Generate_filterMA(10, 10, MA.len=3, ma.scale=ma.scale1) X=Simulate_new_MA(a1, T.len=512, noise.type='wiener') ans=Spec(X, trace=FALSE, only.diag=FALSE) plot(ans) plot(Spec(X, trace=FALSE, only.diag=FALSE, subgrid=TRUE, subgrid.density=10, subgrid.density.relative.to.bandwidth=FALSE)) rm(ans)
ma.scale1=c(-1.4,2.3,-2) a1=Generate_filterMA(10, 10, MA.len=3, ma.scale=ma.scale1) X=Simulate_new_MA(a1, T.len=512, noise.type='wiener') ans=Spec(X, trace=FALSE, only.diag=FALSE) plot(ans) plot(Spec(X, trace=FALSE, only.diag=FALSE, subgrid=TRUE, subgrid.density=10, subgrid.density.relative.to.bandwidth=FALSE)) rm(ans)
A test for the null hypothesis that two spectral density operators (at the
same frequency ) are equal, using a pseudo-AIC criterion for the
choice of the truncation parameter. (used in
Spec_compare_localize_freq
)
Spec_compare_fixed_freq(spec1, spec2, is.pi.multiple, m, kappa.square, autok = 2, K.fixed = NA)
Spec_compare_fixed_freq(spec1, spec2, is.pi.multiple, m, kappa.square, autok = 2, K.fixed = NA)
spec1 , spec2
|
The two sample spectral densities (at the same frequency |
is.pi.multiple |
A logical variable, to specify if |
m |
The number of Fourier frequencies over which the periodogram operator was smoothed. |
kappa.square |
the L2-norm of the weight function used to estimate the spectral density operator |
autok |
A variable used to specify if (and which) pseudo-AIC criterion
is used to select the truncation parameter |
K.fixed |
The value of K used if |
Tavakoli, Shahin and Panaretos, Victor M. "Detecting and Localizing Differences in Functional Time Series Dynamics: A Case Study in Molecular Biophysics", 2014, under revision
Panaretos, Victor M., David Kraus, and John H. Maddocks. "Second-order comparison of Gaussian random functions and the geometry of DNA minicircles." Journal of the American Statistical Association 105.490 (2010): 670-682.
ma.scale2=ma.scale1=c(-1.4,2.3,-2) ma.scale2[3] = ma.scale1[3]+.3 a1=Generate_filterMA(10, 10, MA.len=3, ma.scale=ma.scale1) a2=Generate_filterMA(10, 10, MA.len=3, ma.scale=ma.scale2) X=Simulate_new_MA(a1, T.len=512, noise.type='wiener') Y=Simulate_new_MA(a2, T.len=512, noise.type='wiener') spec.X = Spec(X) spec.Y = Spec(Y) Spec_compare_fixed_freq(spec.X$spec[1,,], spec.Y$spec[1,,], is.pi.multiple=TRUE, spec.X$m, spec.X$kappa.square)
ma.scale2=ma.scale1=c(-1.4,2.3,-2) ma.scale2[3] = ma.scale1[3]+.3 a1=Generate_filterMA(10, 10, MA.len=3, ma.scale=ma.scale1) a2=Generate_filterMA(10, 10, MA.len=3, ma.scale=ma.scale2) X=Simulate_new_MA(a1, T.len=512, noise.type='wiener') Y=Simulate_new_MA(a2, T.len=512, noise.type='wiener') spec.X = Spec(X) spec.Y = Spec(Y) Spec_compare_fixed_freq(spec.X$spec[1,,], spec.Y$spec[1,,], is.pi.multiple=TRUE, spec.X$m, spec.X$kappa.square)
Compare the spectral density operator of two Functional Time Series and localize frequencies at which they differ.
Spec_compare_localize_freq(X, Y, B.T = (dim(X)[1])^(-1/5), W, autok = 2, subgrid.density, verbose = 0, demean = FALSE, K.fixed = NA, subgrid.density.relative.to.bandwidth)
Spec_compare_localize_freq(X, Y, B.T = (dim(X)[1])^(-1/5), W, autok = 2, subgrid.density, verbose = 0, demean = FALSE, K.fixed = NA, subgrid.density.relative.to.bandwidth)
X , Y
|
The |
B.T |
The bandwidth of frequencies over which the periodogram operator
is smoothed. If |
W |
The weight function used to smooth the periodogram operator. Set by default to be the Epanechnikov kernel |
autok |
A variable used to specify if (and which) pseudo-AIC criterion
is used to select the truncation parameter |
subgrid.density |
Only used if |
verbose |
A variable to show the progress of the computations. By
default, |
demean |
A logical variable to choose if the FTS is centered before computing its spectral density operator. |
K.fixed |
The value of K used if |
subgrid.density.relative.to.bandwidth |
logical parameter to specify if
|
X,Y
must be of equal size , where T.len is the length of the time series, and
is the number of basis functions. Each row corresponds to a time point, and each column
corresponds to the coefficient of the corresponding basis function of the FTS.
autok=0
returns the p-values for .
autok=1
uses the AIC
criterion of Tavakoli \& Panaretos
(2015), which is a generalization of the pseudo-AIC introduced in Panaretos
et al (2010).
autok=2
uses the AIC*
criterion of Tavakoli \& Panaretos
(2015), which is an extension of the AIC
criterion that takes into
account the difficulty associated with the estimation of eigenvalues of a
compact operator.
Tavakoli, Shahin and Panaretos, Victor M. "Detecting and Localizing Differences in Functional Time Series Dynamics: A Case Study in Molecular Biophysics", 2014, under revision
Panaretos, Victor M., David Kraus, and John H. Maddocks. "Second-order comparison of Gaussian random functions and the geometry of DNA minicircles." Journal of the American Statistical Association 105.490 (2010): 670-682.
ma.scale2=ma.scale1=c(-1.4,2.3,-2) ma.scale2[3] = ma.scale1[3]+.0 a1=Generate_filterMA(10, 10, MA.len=3, ma.scale=ma.scale1) a2=Generate_filterMA(10, 10, MA.len=3, ma.scale=ma.scale2) X=Simulate_new_MA(a1, T.len=512, noise.type='wiener') Y=Simulate_new_MA(a2, T.len=512, noise.type='wiener') ans0=Spec_compare_localize_freq(X, Y, W=Epanechnikov_kernel, autok=2, subgrid.density=10, verbose=0, demean=FALSE, subgrid.density.relative.to.bandwidth=TRUE) plot(ans0) plot(ans0, method='fdr') PvalAdjust(ans0, method='fdr') ## print FDR adjusted p-values abline(h=.05, lty=3) ans0=Spec_compare_localize_freq(X, Y, W=Epanechnikov_kernel, autok=0, subgrid.density=10, verbose=0, demean=FALSE, subgrid.density.relative.to.bandwidth=TRUE, K.fixed=4) ## fixed values of K plot(ans0) plot(ans0, 'fdr') plot(ans0, 'holm') PvalAdjust(ans0, method='fdr') rm(ans0)
ma.scale2=ma.scale1=c(-1.4,2.3,-2) ma.scale2[3] = ma.scale1[3]+.0 a1=Generate_filterMA(10, 10, MA.len=3, ma.scale=ma.scale1) a2=Generate_filterMA(10, 10, MA.len=3, ma.scale=ma.scale2) X=Simulate_new_MA(a1, T.len=512, noise.type='wiener') Y=Simulate_new_MA(a2, T.len=512, noise.type='wiener') ans0=Spec_compare_localize_freq(X, Y, W=Epanechnikov_kernel, autok=2, subgrid.density=10, verbose=0, demean=FALSE, subgrid.density.relative.to.bandwidth=TRUE) plot(ans0) plot(ans0, method='fdr') PvalAdjust(ans0, method='fdr') ## print FDR adjusted p-values abline(h=.05, lty=3) ans0=Spec_compare_localize_freq(X, Y, W=Epanechnikov_kernel, autok=0, subgrid.density=10, verbose=0, demean=FALSE, subgrid.density.relative.to.bandwidth=TRUE, K.fixed=4) ## fixed values of K plot(ans0) plot(ans0, 'fdr') plot(ans0, 'holm') PvalAdjust(ans0, method='fdr') rm(ans0)
Compare the spectral density operator of two Functional Time Series and localize frequencies at which they differ, and (spatial) regions where they differ
Spec_compare_localize_freq_curvelength(X, Y, B.T = (dim(X)[1])^(-1/5), W, alpha = 0.05, accept = 0, reject = 1, verbose = 0, demean = FALSE)
Spec_compare_localize_freq_curvelength(X, Y, B.T = (dim(X)[1])^(-1/5), W, alpha = 0.05, accept = 0, reject = 1, verbose = 0, demean = FALSE)
X |
The |
Y |
The |
B.T |
The bandwidth of frequencies over which the periodogram operator
is smoothed. If |
W |
The weight function used to smooth the periodogram operator. Set by default to be the Epanechnikov kernel |
alpha |
level for the test |
accept , reject
|
values for accepted, rejected regions |
verbose |
A variable to show the progress of the computations. By
default, |
demean |
A logical variable to choose if the FTS is centered before computing its spectral density operator. |
ma.scale2=ma.scale1=c(-1.4,2.3,-2) ma.scale2[3] = ma.scale1[3]+.4 a1=Generate_filterMA(10, 10, MA.len=3, ma.scale=ma.scale1) a2=Generate_filterMA(10, 10, MA.len=3, ma.scale=ma.scale2) X=Simulate_new_MA(a1, T.len=2^9, noise.type='wiener') Y=Simulate_new_MA(a2, T.len=2^9, noise.type='wiener') ans0=Spec_compare_localize_freq_curvelength(X, Y, W=Epanechnikov_kernel, alpha=.01, demean=TRUE) print(ans0) plot(ans0) rm(ma.scale1, ma.scale2, a1, a2, X, Y, ans0)
ma.scale2=ma.scale1=c(-1.4,2.3,-2) ma.scale2[3] = ma.scale1[3]+.4 a1=Generate_filterMA(10, 10, MA.len=3, ma.scale=ma.scale1) a2=Generate_filterMA(10, 10, MA.len=3, ma.scale=ma.scale2) X=Simulate_new_MA(a1, T.len=2^9, noise.type='wiener') Y=Simulate_new_MA(a2, T.len=2^9, noise.type='wiener') ans0=Spec_compare_localize_freq_curvelength(X, Y, W=Epanechnikov_kernel, alpha=.01, demean=TRUE) print(ans0) plot(ans0) rm(ma.scale1, ma.scale2, a1, a2, X, Y, ans0)
'Spectral density operator of a MA vector process' Object
SpecMA(a, nfreq = 2^9, noise.type)
SpecMA(a, nfreq = 2^9, noise.type)
a |
the filter of the moving average |
nfreq |
the number of frequencies between 0 and pi at which the spectral density has to be computed |
noise.type |
the type of noise that is driving the MA process. See
|
ma.scale1=c(-1.4,2.3,-2) a1=Generate_filterMA(6, 6, MA.len=3, ma.scale=ma.scale1) a1.spec=SpecMA(a1, nfreq=512, noise.type='wiener') plot(a1.spec) rm(a1, a1.spec)
ma.scale1=c(-1.4,2.3,-2) a1=Generate_filterMA(6, 6, MA.len=3, ma.scale=ma.scale1) a1.spec=SpecMA(a1, nfreq=512, noise.type='wiener') plot(a1.spec) rm(a1, a1.spec)