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...
- Divide into training and testing, about 80% in training
- Create sub vectors
nxi, using the testing data, with a delay of
nxdbetween the end of a vector and the start of the next one.
- Compute an autocorrelation matrix
R = sum(xi xi')
- Compute the unit normal eigenvectors of
- Select the
Npminor eigenvectors as
B(nxp, Nm) = R(1:Nm, nxp)'
- 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
X = B' C C = (B B')^(-1) B X
Since B consists of orthonormal vectors,
B B' = I C = B X
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,
Bp, are not normal basis vectors, so products
Bf Bf' and
Bp Bp' do not equal
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