asl / rssa

R package for Singular Spectrum Analysis
55 stars 27 forks source link

Example doesn't work? #226

Closed vdebuen closed 7 years ago

vdebuen commented 7 years ago

Hello

I've installed Rssa package on Windows and I've run the example code

# Construct the Toeplitz trajectory matrix for 'co2' series
h <- new.tmat(co2, L = 10)
# Print the number of columns and rows
print(trows(h)); print(tcols(h))

Results are both 10 columns and 10 rows while the length of co2 series is 468, but I expected a L × K toeplitz matrix with L = 10 and K = 468 - L +1 = 459 as defined on Definition 2 of section 4.2 Toeplitz and Hankel matrices of Computation- and space-efficient implementation of SSA, Anton Korobeynikov

Thanks in advance

neg99 commented 7 years ago

Thank you for this comment. There is a `copy-paste' misprint in the description of this function: new.tmat creates the Toeplitz "autocovariance" matrix, not a Toeplitz trajectory matrix.

Specific of SSA is that it uses Hankel trajectory matrices L x K and does not use Toeplitz matrices L x K, since they can be easily transformed one to the other. Thus, the trajectory L x K matrix X of a 1D series is always Hankel. There are two versions of SSA: Basic SSA and Toeplitz SSA. They differ in decomposition; in particular, in way of eigenvectors' calculation. In Basic SSA, eigenvectors are of X X^T, when in Toeplitz SSA eigenvectors are of a special Toeplitz matrix. This L x L Toeplitz matrix is constructed by new.tmat, see Basic Singular Spectrum Analysis and Forecasting with R, page 5, for the arxiv version of the CSDA paper.

Please note that the functions like new.tmat are internal in Rssa; it is recommended to use the interface functions such as ssa, reconstruct,....

vdebuen commented 7 years ago

Thanks a lot, Nina

Really, what I'm looking for is a method for near super-fast solving of medium-large symmetric Toeplitz least squares systems

y = X b + a; a ~ N(0,s^2)

I cannot find a direct super-fast method in CRAN but I will try Conjugate Gradient method using your function tmatmul to calculate X'X b in linear system

X'y = X'X b

asl commented 7 years ago

@vdebuen In addition, I would suggest you to look into Yule-Walker equations - effectively they allow one to solve a system of linear equations with Toeplitz matrix via several FFTs.

vdebuen commented 7 years ago

Thanks @asl and @neg99

Finally I've solved the problem using new.hmat and hmatmul over reversed vectors in a conjugate gradient cycle.

Best regards