How to use minor component analysis to forecast time series?

Question

I would like to forecast a time series based on minor eigenvalues of the autocorrelation matrix. But I seem to be doing something wrong. I have the following steps implemented in MATLAB, but the forecast result is not comparable to the raw data. The paper that I'm trying to implement is "Ship-Motion Prediction: Algorithms and Simulation Results", X. Zhao, et. al, ICASSP 2004.

Here are the analysis steps for Minor Component Analysis. Can someone indicate a reference, or identify where the error is...

1. Given x(t)
2. Divide into training and testing, about 80% in training
3. Create sub vectors xi of length nxi, using the testing data, with a delay of nxd between the end of a vector and the start of the next one.
4. Compute an autocorrelation matrix R(nxp,nxp) as R = sum(xi xi')
5. Compute the unit normal eigenvectors of R as matrix U(nxp,nxp)
6. Select the Np minor eigenvectors as B(nxp, Nm) = R(1:Nm, nxp)'
7. Divide the minor eigenvectors into past and future B(nxp,Nm) = [Bp(Np,Nm), Bf(Nf,Nm)], where Np+Nf = nxp

Here's where I'm starting getting lost. I would typically calculate the coefficient vector C for basis vector matrix B:

X = B' C
C = (B B')^(-1) B X

Since B consists of orthonormal vectors,

B B' = I
C = B X

But, since Nf of the X's are in the future, they can't be used. So, we have (in more detail, since somewhere in these steps the ability to recreate the raw data is lost ...):

[ Xp ] = [Bp Bf]' C = [ Bp'] C
[ Xf ]                [ Bf']

then, as before

C = [Bp Bf] [Xp; Xf] = Bp Xp + Bf Xf

Note, the above equation is similar to equation 2 in the report, but in the report the sum is equal to zero instead of 'C'. I'm not sure why this is different.

Expanding the expression for past data,

Xp = Bp' C = Bp'[Bp Xp + Bf Xf]

Next solve this equation for the future values, Xf. Note that the partial basis vectors, Bf and Bp, are not normal basis vectors, so products Bf Bf' and Bp Bp' do not equal I.

Xf = (Bf' Bf)^(-1) Bf' ((Bp Bp')^(-1) Bp - Bp) Xp

In the paper, this result is equation 3, and is given by

Xf = (Bf' Bf)^(-1) Bf' (-Bp) Xp

So, I don't know where the error is, why the (Bp Bp')^(-1) Bp term is missing in the paper. Also, the additional term does not appear to be numerically stable. It's possible that adding sensor noise, or using SVD for basis vectors instead of eigenvectors, would improve the numerical stability.

Here's a sample MATLAB code to generate data similar to Figure 1 in the paper. I have implemented the MCA analysis function also, but didn't include it here since it's not correct, but the MCA function follows directly from the above derivation, with the option to leave the additional term out of Xf. But, the MCA forecast result doesn't look anything like Figure 1, or, in other words, the MCA forecast is not comparable with the data at the 20 second (40 data point) horizon, but the paper indicates that the MCA 20 second forecast has good agreement with the data.

function quest_1
npoints=5000+1000;
tspan=0:.5:npoints/2;
wn=28/500*2*pi;
zeta=0.02;
zeta_nl = 0.015;
ft=tspan;
f=(rand(size(ft))-.5)/4;
options = odeset('RelTol',1e-4,'AbsTol',[1e-4 1e-4]);
[T,Y] = ode45(@model,tspan,[-.4 -.2],options, wn, zeta, ft, f, zeta_nl);
figure;
plot(T,Y(:,1))
figure;
plot(T-2500,Y(:,1))
ax=axis;
axis([0,600,ax(3:4)])
end

function dy = model(t,y,wn, zeta, ft, f, zeta_nl)
dy = zeros(2,1);    % a column vector
dy(1) = y(2);
ff = interp1(ft,f,t);
dy(2) = -2*zeta*wn*y(2) - wn*wn*y(1) + ff; % ddy + 2 zeta wn dy + wn^2 y = f
dy(2) = dy(2) - zeta_nl*y(1)*y(1)*y(2);
end

Answer 1

The paper has a few errors and ambiguities.

First of all they confuse the integer constant n (lower case) which is the number of samples in their data window (the dimension of the vector Xi), and N (upper case) which is the number of data windows in their entire data set. In their example n = 800 and N = 420. The dimension of the eigenvectors is n, not N. They start mixing them up when describing the B matrix. B is P x n (not P x N), where P < n (you choose a subset of the eigenvectors), and B2 is P x (n-n1).

Secondly, the eigenvectors used to construct B need to be transposed and stacked as the rows of B (not the columns). B should be a "fat" matrix, i.e. more columns than rows.

I think your step 6 should be something like B(Nm,nxp) = U(nxp,1:Nm)'