PCA4_TS: Principal component analysis for time serise
In HDTSA: High Dimensional Time Series Analysis Tools
PCA4_TS R Documentation
Principal component analysis for time serise
Description
PCA4_TS() seeks for a contemporaneous linear
transformation for a multivariate time series such that the transformed
series is segmented into several lower-dimensional subseries:
{\bf
y}_t={\bf Ax}_t,
where {\bf x}_t is an unobservable p \times 1
weakly stationary time series consisting of q\ (≥q 1) both
contemporaneously and serially uncorrelated subseries. See Chang, Guo and
Yao (2018).
Usage
PCA4_TS(
Y,
lag.k = 5,
thresh = FALSE,
tuning.vec = NULL,
K = 5,
prewhiten = TRUE,
permutation = c("max", "fdr"),
m = NULL,
beta,
just4pre = FALSE,
verbose = FALSE
)
Arguments
Y
{\bf Y} = \{{\bf y}_1, … , {\bf y}_n \}', a data matrix
with n rows and p columns, where n is the sample size and
p is the dimension of {\bf y}_t. The procedure will first
normalize {\bf y}_t as \widehat{{\bf V}}^{-1/2}{\bf y}_t, where
\widehat{{\bf V}} is an estimator for covariance of {\bf y}_t.
See details below for the selection of \widehat{{\bf V}}^{-1}.
lag.k
Time lag k_0 used to calculate the nonnegative definte
matrix \widehat{{\bf W}}_y:
\widehat{\mathbf{W}}_y\ =\
∑_{k=0}^{k_0}\widehat{\mathbf{Σ}}_y(k)\widehat{\mathbf{Σ}}_y(k)'=\mathbf{I}_p+∑_{k=1}^{k_0}\widehat{\mathbf{Σ}}_y(k)\widehat{\mathbf{Σ}}_y(k)',
where \widehat{\bf Σ}_y(k) is the sample autocovariance of
\widehat{{\bf V}}^{-1/2}{\bf y}_t at lag k. See (2.5) in
Chang, Guo and Yao (2018).
thresh
Logical. If FALSE (the default), no thresholding will
be applied to estimate \widehat{{\bf W}}_y. If TRUE, a
thresholding method will be applied first to estimate \widehat{{\bf
W}}_y, see (3.5) in Chang, Guo and Yao (2018).
tuning.vec
The value of the tuning parameter λ in the
thresholding level u = λ √{n^{-1}\log p}, where default
value is 2. If tuning.vec is a vector, then a cross validation
method proposed in Cai and Liu (2011) will be used to choose the best
tuning parameter λ.
K
The number of folders used in the cross validation for the
selection of λ, the default is 5. It is required when
thresh = TRUE.
prewhiten
Logical. If TRUE (the default), we prewhiten each
transformed component series of \hat{\bf z}_t [See Section 2.2.1 in
Chang, Guo and Yao (2018)] by fitting a univariate AR model with the order
between 0 and 5 determined by AIC. If FALSE, then prewhiten
procedure will not be performed to \hat{\bf z}_t.
permutation
The method of permutation procedure to assign the
components of \hat{\bf z}_t to different groups [See Section 2.2.1 in
Chang, Guo and Yao (2018)]. Option is 'max' (Maximum cross
correlation method) or 'fdr' (False discovery rate procedure based
on multiple tests), default is permutation = 'max'. See Sections
2.2.2 and 2.2.3 in Chang, Guo and Yao (2018) for more information.
m
A positive constant used in the permutation procedure [See (2.10) in
Chang, Guo and Yao (2018)]. If m is not specified, then default
option is m = 10.
beta
The error rate used in the permutation procedure when
permutation = 'fdr'.
just4pre
Logical. If TRUE, the procedure outputs \hat{\bf
z}_t, otherwise outputs \hat{\bf x}_t (the permutated version of
\hat{\bf z}_t).
verbose
Logical. If TRUE, the main results of the permutation
procedure will be output on the console. Otherwise, the result will not be
output.
Details
When p>n^{1/2}, the procedure use package clime to
estimate the precision matrix \widehat{{\bf V}}^{-1}, otherwise uses
function cov() to estimate \widehat{{\bf V}} and calculate its
inverse. When p>n^{1/2}, we recommend to use the thresholding method
to calculate \widehat{{\bf W}}_y, see more information in Chang, Guo
and Yao (2018).
Value
The output of the segment procedure is a list containing the
following components:
B
The p\times p transformation matrix
such that \hat{\bf z}_t = \widehat{\bf B}{\bf y}_t, where
\widehat{\bf B}=\widehat{\bf Γ}_y\widehat{{\bf V}}^{-1/2}.
Z
\hat{\bf Z}=\{\hat{\bf z}_1,…,\hat{\bf z}_n\}', the
transformed series with n rows and p columns.
The output of the permutation procedure is a list containing the
following components:
NoGroups
number of groups with at least two components series.
No_of_Members
The cardinalities of different groups.
Groups
The indices of the components in \hat{\bf z}_t that
belongs to a group.
References
Chang, J., Guo, B. & Yao, Q. (2018). Principal component
analysis for second-order stationary vector time series, The Annals of
Statistics, Vol. 46, pp. 2094–2124.
Cai, T. & Liu, W. (2011). Adaptive thresholding for sparse covariance
matrix estimation, Journal of the American Statistical Association, Vol.
106, pp. 672–684.
Cai, T., Liu, W., & Luo, X. (2011). A constrained l1 minimization
approach for sparse precision matrix estimation, Journal of the American
Statistical Association, Vol. 106, pp. 594–607.
Examples
## Example 1 (Example 5 of Chang Guo and Yao (2018)).
## p=6, x_t consists of 3 independent subseries with 3, 2 and 1 components.
p <- 6;n <- 1500
# Generate x_t
X <- mat.or.vec(p,n)
x <- arima.sim(model=list(ar=c(0.5, 0.3), ma=c(-0.9, 0.3, 1.2,1.3)),
n=n+2,sd=1)
for(i in 1:3) X[i,] <- x[i:(n+i-1)]
x <- arima.sim(model=list(ar=c(0.8,-0.5),ma=c(1,0.8,1.8) ),n=n+1,sd=1)
for(i in 4:5) X[i,] <- x[(i-3):(n+i-4)]
x <- arima.sim(model=list(ar=c(-0.7, -0.5), ma=c(-1, -0.8)),n=n,sd=1)
X[6,] <- x
# Generate y_t
A <- matrix(runif(p*p, -3, 3), ncol=p)
Y <- A%*%X
Y <- t(Y)
res <- PCA4_TS(Y, lag.k=5,permutation = "max")
res1=PCA4_TS(Y, lag.k=5,permutation = "fdr", beta=10^(-10))
# The transformed series z_t
Z <- res$Z
# Plot the cross correlogram of z_t and y_t
Y <- data.frame(Y);Z=data.frame(Z)
names(Y) <- c("Y1","Y2","Y3","Y4","Y5","Y6")
names(Z) <- c("Z1","Z2","Z3","Z4","Z5","Z6")
# The cross correlogram of y_t shows no block pattern
acfY <- acf(Y)
# The cross correlogram of z_t shows 3-2-1 block pattern
acfZ <- acf(Z)
## Example 2 (Example 6 of Chang Guo and Yao (2018)).
## p=20, x_t consists of 5 independent subseries with 6, 5, 4, 3 and 2 components.
p <- 20;n <- 3000
# Generate x_t
X <- mat.or.vec(p,n)
x <- arima.sim(model=list(ar=c(0.5, 0.3), ma=c(-0.9, 0.3, 1.2,1.3)),n.start=500,
n=n+5,sd=1)
for(i in 1:6) X[i,] <- x[i:(n+i-1)]
x <- arima.sim(model=list(ar=c(-0.4,0.5),ma=c(1,0.8,1.5,1.8)),n.start=500,n=n+4,sd=1)
for(i in 7:11) X[i,] <- x[(i-6):(n+i-7)]
x <- arima.sim(model=list(ar=c(0.85,-0.3),ma=c(1,0.5,1.2)), n.start=500,n=n+3,sd=1)
for(i in 12:15) X[i,] <- x[(i-11):(n+i-12)]
x <- arima.sim(model=list(ar=c(0.8,-0.5),ma=c(1,0.8,1.8)),n.start=500,n=n+2,sd=1)
for(i in 16:18) X[i,] <- x[(i-15):(n+i-16)]
x <- arima.sim(model=list(ar=c(-0.7, -0.5), ma=c(-1, -0.8)),n.start=500,n=n+1,sd=1)
for(i in 19:20) X[i,] <- x[(i-18):(n+i-19)]
# Generate y_t
A <- matrix(runif(p*p, -3, 3), ncol=p)
Y <- A%*%X
Y <- t(Y)
res <- PCA4_TS(Y, lag.k=5,permutation = "max")
res1 <- PCA4_TS(Y, lag.k=5,permutation = "fdr",beta=10^(-200))
# The transformed series z_t
Z <- res$Z
# Plot the cross correlogram of x_t and y_t
Y <- data.frame(Y);Z <- data.frame(Z)
namesY=NULL;namesZ=NULL
for(i in 1:p)
{
namesY <- c(namesY,paste0("Y",i))
namesZ <- c(namesZ,paste0("Z",i))
}
names(Y) <- namesY;names(Z) <- namesZ
# The cross correlogram of y_t shows no block pattern
acfY <- acf(Y, plot=FALSE)
plot(acfY, max.mfrow=6, xlab='', ylab='', mar=c(1.8,1.3,1.6,0.5),
oma=c(1,1.2,1.2,1), mgp=c(0.8,0.4,0),cex.main=1)
# The cross correlogram of z_t shows 6-5-4-3-2 block pattern
acfZ <- acf(Z, plot=FALSE)
plot(acfZ, max.mfrow=6, xlab='', ylab='', mar=c(1.8,1.3,1.6,0.5),
oma=c(1,1.2,1.2,1), mgp=c(0.8,0.4,0),cex.main=1)
# Identify the permutation mechanism
permutation <- res
permutation$Groups
HDTSA documentation built on Jan. 7, 2023, 5:26 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
| PCA4_TS | R Documentation |
Principal component analysis for time serise
Description
PCA4_TS() seeks for a contemporaneous linear
transformation for a multivariate time series such that the transformed
series is segmented into several lower-dimensional subseries:
{\bf y}_t={\bf Ax}_t,
where {\bf x}_t is an unobservable p \times 1 weakly stationary time series consisting of q\ (≥q 1) both contemporaneously and serially uncorrelated subseries. See Chang, Guo and Yao (2018).
Usage
PCA4_TS(
Y,
lag.k = 5,
thresh = FALSE,
tuning.vec = NULL,
K = 5,
prewhiten = TRUE,
permutation = c("max", "fdr"),
m = NULL,
beta,
just4pre = FALSE,
verbose = FALSE
)
Arguments
Y |
{\bf Y} = \{{\bf y}_1, … , {\bf y}_n \}', a data matrix with n rows and p columns, where n is the sample size and p is the dimension of {\bf y}_t. The procedure will first normalize {\bf y}_t as \widehat{{\bf V}}^{-1/2}{\bf y}_t, where \widehat{{\bf V}} is an estimator for covariance of {\bf y}_t. See details below for the selection of \widehat{{\bf V}}^{-1}. |
lag.k |
Time lag k_0 used to calculate the nonnegative definte matrix \widehat{{\bf W}}_y: \widehat{\mathbf{W}}_y\ =\ ∑_{k=0}^{k_0}\widehat{\mathbf{Σ}}_y(k)\widehat{\mathbf{Σ}}_y(k)'=\mathbf{I}_p+∑_{k=1}^{k_0}\widehat{\mathbf{Σ}}_y(k)\widehat{\mathbf{Σ}}_y(k)', where \widehat{\bf Σ}_y(k) is the sample autocovariance of \widehat{{\bf V}}^{-1/2}{\bf y}_t at lag k. See (2.5) in Chang, Guo and Yao (2018). |
thresh |
Logical. If |
tuning.vec |
The value of the tuning parameter λ in the
thresholding level u = λ √{n^{-1}\log p}, where default
value is 2. If |
K |
The number of folders used in the cross validation for the
selection of λ, the default is 5. It is required when
|
prewhiten |
Logical. If |
permutation |
The method of permutation procedure to assign the
components of \hat{\bf z}_t to different groups [See Section 2.2.1 in
Chang, Guo and Yao (2018)]. Option is |
m |
A positive constant used in the permutation procedure [See (2.10) in
Chang, Guo and Yao (2018)]. If m is not specified, then default
option is |
beta |
The error rate used in the permutation procedure when
|
just4pre |
Logical. If |
verbose |
Logical. If |
Details
When p>n^{1/2}, the procedure use package clime to
estimate the precision matrix \widehat{{\bf V}}^{-1}, otherwise uses
function cov() to estimate \widehat{{\bf V}} and calculate its
inverse. When p>n^{1/2}, we recommend to use the thresholding method
to calculate \widehat{{\bf W}}_y, see more information in Chang, Guo
and Yao (2018).
Value
The output of the segment procedure is a list containing the following components:
B |
The p\times p transformation matrix such that \hat{\bf z}_t = \widehat{\bf B}{\bf y}_t, where \widehat{\bf B}=\widehat{\bf Γ}_y\widehat{{\bf V}}^{-1/2}. |
Z |
\hat{\bf Z}=\{\hat{\bf z}_1,…,\hat{\bf z}_n\}', the transformed series with n rows and p columns. |
The output of the permutation procedure is a list containing the following components:
NoGroups |
number of groups with at least two components series. |
No_of_Members |
The cardinalities of different groups. |
Groups |
The indices of the components in \hat{\bf z}_t that belongs to a group. |
References
Chang, J., Guo, B. & Yao, Q. (2018). Principal component analysis for second-order stationary vector time series, The Annals of Statistics, Vol. 46, pp. 2094–2124.
Cai, T. & Liu, W. (2011). Adaptive thresholding for sparse covariance matrix estimation, Journal of the American Statistical Association, Vol. 106, pp. 672–684.
Cai, T., Liu, W., & Luo, X. (2011). A constrained l1 minimization approach for sparse precision matrix estimation, Journal of the American Statistical Association, Vol. 106, pp. 594–607.
Examples
## Example 1 (Example 5 of Chang Guo and Yao (2018)).
## p=6, x_t consists of 3 independent subseries with 3, 2 and 1 components.
p <- 6;n <- 1500
# Generate x_t
X <- mat.or.vec(p,n)
x <- arima.sim(model=list(ar=c(0.5, 0.3), ma=c(-0.9, 0.3, 1.2,1.3)),
n=n+2,sd=1)
for(i in 1:3) X[i,] <- x[i:(n+i-1)]
x <- arima.sim(model=list(ar=c(0.8,-0.5),ma=c(1,0.8,1.8) ),n=n+1,sd=1)
for(i in 4:5) X[i,] <- x[(i-3):(n+i-4)]
x <- arima.sim(model=list(ar=c(-0.7, -0.5), ma=c(-1, -0.8)),n=n,sd=1)
X[6,] <- x
# Generate y_t
A <- matrix(runif(p*p, -3, 3), ncol=p)
Y <- A%*%X
Y <- t(Y)
res <- PCA4_TS(Y, lag.k=5,permutation = "max")
res1=PCA4_TS(Y, lag.k=5,permutation = "fdr", beta=10^(-10))
# The transformed series z_t
Z <- res$Z
# Plot the cross correlogram of z_t and y_t
Y <- data.frame(Y);Z=data.frame(Z)
names(Y) <- c("Y1","Y2","Y3","Y4","Y5","Y6")
names(Z) <- c("Z1","Z2","Z3","Z4","Z5","Z6")
# The cross correlogram of y_t shows no block pattern
acfY <- acf(Y)
# The cross correlogram of z_t shows 3-2-1 block pattern
acfZ <- acf(Z)
## Example 2 (Example 6 of Chang Guo and Yao (2018)).
## p=20, x_t consists of 5 independent subseries with 6, 5, 4, 3 and 2 components.
p <- 20;n <- 3000
# Generate x_t
X <- mat.or.vec(p,n)
x <- arima.sim(model=list(ar=c(0.5, 0.3), ma=c(-0.9, 0.3, 1.2,1.3)),n.start=500,
n=n+5,sd=1)
for(i in 1:6) X[i,] <- x[i:(n+i-1)]
x <- arima.sim(model=list(ar=c(-0.4,0.5),ma=c(1,0.8,1.5,1.8)),n.start=500,n=n+4,sd=1)
for(i in 7:11) X[i,] <- x[(i-6):(n+i-7)]
x <- arima.sim(model=list(ar=c(0.85,-0.3),ma=c(1,0.5,1.2)), n.start=500,n=n+3,sd=1)
for(i in 12:15) X[i,] <- x[(i-11):(n+i-12)]
x <- arima.sim(model=list(ar=c(0.8,-0.5),ma=c(1,0.8,1.8)),n.start=500,n=n+2,sd=1)
for(i in 16:18) X[i,] <- x[(i-15):(n+i-16)]
x <- arima.sim(model=list(ar=c(-0.7, -0.5), ma=c(-1, -0.8)),n.start=500,n=n+1,sd=1)
for(i in 19:20) X[i,] <- x[(i-18):(n+i-19)]
# Generate y_t
A <- matrix(runif(p*p, -3, 3), ncol=p)
Y <- A%*%X
Y <- t(Y)
res <- PCA4_TS(Y, lag.k=5,permutation = "max")
res1 <- PCA4_TS(Y, lag.k=5,permutation = "fdr",beta=10^(-200))
# The transformed series z_t
Z <- res$Z
# Plot the cross correlogram of x_t and y_t
Y <- data.frame(Y);Z <- data.frame(Z)
namesY=NULL;namesZ=NULL
for(i in 1:p)
{
namesY <- c(namesY,paste0("Y",i))
namesZ <- c(namesZ,paste0("Z",i))
}
names(Y) <- namesY;names(Z) <- namesZ
# The cross correlogram of y_t shows no block pattern
acfY <- acf(Y, plot=FALSE)
plot(acfY, max.mfrow=6, xlab='', ylab='', mar=c(1.8,1.3,1.6,0.5),
oma=c(1,1.2,1.2,1), mgp=c(0.8,0.4,0),cex.main=1)
# The cross correlogram of z_t shows 6-5-4-3-2 block pattern
acfZ <- acf(Z, plot=FALSE)
plot(acfZ, max.mfrow=6, xlab='', ylab='', mar=c(1.8,1.3,1.6,0.5),
oma=c(1,1.2,1.2,1), mgp=c(0.8,0.4,0),cex.main=1)
# Identify the permutation mechanism
permutation <- res
permutation$Groups
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.