Karhunen Loeve Expansion in Matlab

[confused_engineer]

Hello everyone.

I am trying to generate the KL expansion of a stochastic process. I use a Monte Carlo sampling method to generate the process, which involves two random variables and I compare it with it's theoretical mean for 50 values of time and they look quite similar. Then, I calculate the autocovariance matrix, from where I extract the eigenvalues and eigenvectors which are used to calculate the new variable Y which is the stochastic process S in a base where the random variables are not correlated. However, if I calculate the mean of Y it doesn't look like the the mean of S at all.

Can someone please tell me what am I doing wrong?
Thank you very much for reading.

clear
close all
V = unifrnd(1,2,1,10000);
A = betarnd(2,2,1,10000);
t=50;
for i=1:t
S(i,:)=V*i+0.5*A*i^2;
theoreticalmeanS(i)=3/2*i+1/4*i^2;
meanS(i)=mean(S(i));
end
for i=1:t
for j=1:t
AutocorrelationS(i,j)=(S(i,:)*S(j,:)')/10000;
AutocovarianceS(i,j)=(S(i,:)*S(j,:)')/10000-meanS(i)*meanS(j);
end
end

[eigenvectors, eigenvalues]=eig(AutocovarianceS);
eigenvalues_column=eig(AutocovarianceS);
Y=eigenvalues*S;Note: I am following the book of Alberto Leon Garcia (probability statistics and random processes for electrical engineering), pages 324-325.

 

[JMz]

Not sure if I can help here (I am extremely familiar with the core of MATLAB, but not with its stats toolbox, which you appear to be using), but let me ask a question: For the data sizes you're dealing with, is there a reason not to use SVD to give you the PCs and their weights directly, instead of going through the covariance intermediate and finding eigenvalues?

My reason for asking is partly just that, in theory and in my experience, SVD is more stable and well-behaved than EVD (and there are deep reasons for that). For instance, with round-off, you might have a complex eigenvalue for the covariance matrix; SVD is simply not capable of that.

 

[Stephen Tashi]

for i=1:t
S(i,:)=V*i+0.5*A*i^2;
theoreticalmeanS(i)=3/2*i+1/4*i^2;
meanS(i)=mean(S(i));
end

What does $S(i,k)$ represent? The $k$ value of the $i$th realization of the process? - or the $i$th value of the $k$th realization of the process? (If the process takes place in time, which index represents time?)