% ou1.m simulates Ornstein-Uhlenbeck by inverting the covariance matrix

beta = 10;                                  % decay rate of OU covariance
a    = 5;                                  % length of spatial interval
n    = 100;                                % number of spatial points

rho  = exp(-beta*a/n);                     % covariance at distance 1

S = zeros(n+1);                            % covariance matrix
S(1,:) = rho .^ (0:n);                 % first row

for i=2:(n+1),
  S(i,i:(n+1)) = S(1,1:(n+2-i));           % fill in upper right part of S
end

S = S + S' - eye(n+1);                     % fill in lower triangle

[U,lam] = eig(S);                          % eigenvalues and vectors

[lam,k] = sort(-diag(lam));                % sort in decreasing order
lam = -lam;                                % fix signs

Z = randn(n+1,1);                          % random amplitudes
X = U(:,k) * diag(sqrt(lam)) * Z;          % the Ornstein-Uhlenbeck field
t = (0:n)*a/n;                             % time points

clf
subplot(3,1,1)
plot(t, X);
title('Ornstein-Uhlenbeck field')

subplot(3,1,2)
plot(log(lam),'*');                             % plot eigenvalues
title('Logarithms of Eigenvalues')

fprintf('Press any key to see eigenvectors and partial sums')

pause

Y = zeros(n+1,1);

for i=1:(n+1),
  Y = Y + U(:,k(i)) * sqrt(lam(i)) * Z(i);

  subplot(3,1,2)
  plot(t, U(:,k) * diag(sqrt(lam)) * [Z(1:i); zeros(n+1-i,1)]);
  title(['Partial sum with ' int2str(i) ' terms']);
   
  subplot(3,1,3)
  plot(t, U(:,k(i)));
  title(['Eigenvector number ' int2str(i)]);

  drawnow;

  pause
end