https://github.com/cran/pcdpca
Raw File
Tip revision: ce9d67e15a5402f4d12e62fe9882ada23b12971e authored by Lukasz Kidzinski on 26 November 2016, 23:06:38 UTC
version 0.2.1
Tip revision: ce9d67e
pcdpca.Rd
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/pcdpca.R
\name{pcdpca}
\alias{pcdpca}
\title{Compute periodically correlacted DPCA filter coefficients}
\usage{
pcdpca(X, V = NULL, lags = -10:10, q = NULL, weights = NULL,
  freq = NULL, period = NULL)
}
\arguments{
\item{X}{multivariate stationary time series}

\item{V}{correlation structure between coefficients of vectors (default diagonal)}

\item{lags}{requested filter coefficients}

\item{q}{window for spectral density estimation as in \code{\link{spectral.density}}}

\item{weights}{as in \code{\link{spectral.density}}}

\item{freq}{frequency grid to estimate on as in \code{\link{spectral.density}}}

\item{period}{period of the periodic time series}
}
\value{
principal components series
}
\description{
For a given periodically correlated multivariate process \code{X} eigendecompose it's spectral density
and use an inverse fourier transform to get coefficients of the optimal filters.
}
\examples{
## Prepare some process
library(fda)
library(freqdom)
d = 7
n = 100
A = t(t(matrix(rnorm(d*n),ncol=d,nrow=n))*7:1)
B = t(t(matrix(rnorm(d*n),ncol=d,nrow=n))*7:1)
C = t(t(matrix(rnorm(d*n),ncol=d,nrow=n))*7:1)

X = matrix(0,ncol=d,nrow=3*n)
X[3*(1:n) - 1,] = A
X[3*(1:n) - 2,] = A + B
X[3*(1:n) ,] = 2*A - B + C

basis = create.fourier.basis(nbasis=7)
X.fd = fd(t(Re(X)),basis=basis)
plot(X.fd)

## Hold out some datapoints
train = 1:(50*3)
test = (50*3) : (3*n)

## Static PCA ##
PR = prcomp(as.matrix(X[train,]))
Y1 = as.matrix(X) \%*\% PR$rotation
Y1[,-1] = 0
Xpca.est = Y1 \%*\% t(PR$rotation)

## Dynamic PCA ##
XI.est = dprcomp(as.matrix(X[train,]),
   q=3,lags=-2:2,weights="Bartlett",
   freq=pi*(-150:150/150))  # finds the optimal filter
Y.est = XI.est \%c\% X  # applies the filter
Y.est[,-1] = 0 # forces the use of only one component
Xdpca.est = t(rev(XI.est)) \%c\% Y.est    # deconvolution

## Periodically correlated PCA ##
XI.est.pc = pcdpca(as.matrix(X[train,]),
   q=3,lags=-2:2,weights="Bartlett",
   freq=pi*(-150:150/150),period=3)  # finds the optimal filter
Y.est.pc = pcdpca.scores(X, XI.est.pc)  # applies the filter
Y.est.pc[,-1] = 0 # forces the use of only one component
Xpcdpca.est = pcdpca.inverse(Y.est.pc, XI.est.pc)  # deconvolution

## Results
cat("NMSE PCA = ")
r0 = MSE(X[test,],Xpca.est[test,]) / MSE(X[test,],0)
cat(r0)
cat("\\nNMSE DPCA = ")
r1 = MSE(X[test,],Xdpca.est[test,]) / MSE(X[test,],0)
cat(r1)
cat("\\nNMSE PCDPCA = ")
r2 = MSE(X[test,],Xpcdpca.est[test,]) / MSE(X[test,],0)
cat(r2)
cat("\\n")
}
\references{
Kidzinski, Kokoszka, Jouzdani
Dynamic principal components of periodically correlated functional time series
Research report, 2016
}
\seealso{
\code{\link{pcdpca.inverse}}, \code{\link{pcdpca.scores}}
}

back to top