%% Time-Varying Systems and Computations
% Klaus Diepold, Matthias Kissel, Sven Gronauer, Nehir Güzelkaya
%
%% Stochastic Variables
% We are considering a time-varying stochastic variable <x_k>, where k is the time index. The variable is sampled from x_k \tilde{} \mathcal{N}(0, 1) for all k. This means that the distribution is time-invariant (we are looking at a stationary process).
% We are going to sample x_k and stack the results into a X \in \mathcal{R}^{K \times M} matrix. Here K is the number of timesteps, and M the number of samples we draw for each time step (more details are given in the following). We treat K and M as hyperparameters:

K = 300;
M = 200;

%% 1. Multiple Realizations available
% We now create a matrix containing realisations of the stochastic variable over time. For that, we stack different timesteps k=1,...,K and different measurements of the stochastic variable m=1,...,M into a X \in \mathcal{R}^{K \times M} matrix.

X = normrnd(0, 1, K, M);

% Since our stochastic variable is stationary and ergodic, X^TX should be a symmetric, positive (semi-) definite Toeplitz matrix. Let's checkt that:

R = X'*X;

display("Is symmetric: ");
display(issymmetric(R));

display("Is positive definite: ");
display(all(eig(R) >= 0))

Toeplitz_A = diag(ones(1,size(R,1)-1),-1);
Toeplitz_A(1, end) = 1;
Toeplitz_B = diag(ones(1,size(R,2)-1),-1);
Toeplitz_B(1, end) = -1;

R_displacement = Toeplitz_A * R - R * Toeplitz_B;
[U, S, V] = svd(R_displacement);
plot(diag(S));

% We see that the rank of the matrix with the displacements drops at the first singular values. This effect gets bigger with increasing matrix sizes (i.e. increasing the sample size and the number of timesteps).
%% 2. Only one realization available
% Now we try to recover X from only a single realization <x_k> for different timesteps. Therefore we have

x = normrnd(0, 1, K, 1);

% For that, we build a matrix \hat{X} based on shifted versions of x. We investigate 2 methods of embedding the shifted versions of x into X: Add zeros at the top and the bottom, or cut out the bottom.
%% 2.1 Add zeros at the top and the bottom

padding_size = 10;
padded_dim_size = 2*K-1+2*padding_size;
X_hat = zeros(padded_dim_size,K);
for k=1:K
    X_hat(padding_size+k:padding_size+k+K-1,k) = x;
end

% Now let's check if X_{hat}^TX_{hat} is symmetric positive definite and has toeplitz shape.

R = X_hat'*X_hat;

display("Is symmetric: ");
display(issymmetric(R));

display("Is positive definite: ");
display(all(eig(R) >= 0));

Toeplitz_A = diag(ones(1,size(R,1)-1),-1);
Toeplitz_A(1, end) = 1;
Toeplitz_B = diag(ones(1,size(R,2)-1),-1);
Toeplitz_B(1, end) = -1;

R_displacement = Toeplitz_A * R - R * Toeplitz_B;
[U, S, V] = svd(R_displacement);
plot(diag(S));

%% 2.2 Cut out the bottom

X_hat = zeros(K,K);
for k=1:K
    X_hat(k:end,k) = x(1:K-k+1);
end

% Now let's check if X_{hat}^TX_{hat} is symmetric positive definite and has toeplitz shape.

R = X_hat'*X_hat;

display("Is symmetric: ");
display(issymmetric(R));

display("Is positive definite: ");
display(all(eig(R) >= 0));

Toeplitz_A = diag(ones(1,size(R,1)-1),-1);
Toeplitz_A(1, end) = 1;
Toeplitz_B = diag(ones(1,size(R,2)-1),-1);
Toeplitz_B(1, end) = -1;

R_displacement = Toeplitz_A * R - R * Toeplitz_B;
[U, S, V] = svd(R_displacement);
plot(diag(S));