pcdpca: Compute periodically correlacted DPCA filter coefficients

Description Usage Arguments Value References See Also Examples

View source: R/pcdpca.R

Description

For a given periodically correlated multivariate process X eigendecompose it's spectral density and use an inverse fourier transform to get coefficients of the optimal filters.

Usage

1
pcdpca(X, period = NULL, q = 30, freq = (-1000:1000/1000) * pi)

Arguments

X

multivariate stationary time series

period

period of the periodic time series

q

window for spectral density estimation as in spectral.density

freq

frequency grid to estimate on as in spectral.density

Value

principal components series

References

Kidzinski, Kokoszka, Jouzdani Dynamic principal components of periodically correlated functional time series Research report, 2016

See Also

pcdpca.inverse, pcdpca.scores

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
## Prepare some process
library(fda)
library(freqdom)

MSE = function(X,Y=0){ sum((X-Y)**2) / nrow(X) }

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 = dpca(as.matrix(X[train,]),
   q=3,
   freq=pi*(-150:150/150),
   Ndpc=1)  # finds the optimal filter
Y.est = freqdom::filter.process(X, XI.est$filters )
Xdpca.est = freqdom::filter.process(Y.est, t(rev(XI.est$filters)) )    # deconvolution

## Periodically correlated PCA ##
XI.est.pc = pcdpca(as.matrix(X[train,]),
   q=3,
   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")

Example output

Loading required package: splines
Loading required package: Matrix

Attaching package: 'fda'

The following object is masked from 'package:graphics':

    matplot

Loading required package: mvtnorm
Loading required package: matrixcalc

Attaching package: 'freqdom'

The following object is masked from 'package:base':

    %*%

[1] "done"
NMSE PCA = 0.7630073
NMSE DPCA = 0.7866877
NMSE PCDPCA = 0.5405852

pcdpca documentation built on May 2, 2019, 3:38 p.m.