Finding the $S$ matrix for algebraic Riccati equation from data?

Question

Assume that you have two vectors. The gaussian white noise with zero mean $e \in \mathbb{R}^{p\times N}$ and the white disturbance zero mean vector $\omega \in \mathbb{R}^{n\times N}$

If I create this matrix

$$A = \frac{1}{N}\begin{bmatrix} \omega \omega^T & \omega e^T \\ e\omega^T & ee^T \end{bmatrix}$$

And then I take the $Q = A(1:n, 1:n)$ I get the same results as if I compute the MATLAB commando

Q = cov(w')

Same if I compute $R = A(n+1:n+p, n+1:n+p)$

R = cov(e')

But what if I want to find the $S$ matrix for the Algebratic Riccati Equation? The simple solution is $S = A(1:n, n+1:n+p)$.

Question:

Is there any way I could use MATLAB commando cov for computing $S$ if I know $e$ and $\omega$?

How to reproduce this issue:

1. Install MATLAB or GNU Octave
2. Download MataveID and Matavecontrol and install it, or just stand in the same folder as the .m files together. Installation instructions are available at the page.
3. open cca.m file

Then change these lines

% Computing the covariance matrix
covariance = [w*w' w*e'; e*w' e*e']/(N-2); % <<-- Change this to 2 instead of 1

% Compute Q, R, S for the riccati equation
Q = covariance(1:n, 1:n) % No semicolon
R = covariance(n+1:n+p, n+1:n+p) % No semicolon
S = covariance(1:n, n+1:n+p) % No semicolon

% Add these
Q1 = cov(w')
R1 = cov(e')
%S1 = cov ?

Run these in the command window.

G = tf(1, [1 1.5 1]);
[u, t] = gensig('square', 10, 10, 100);
u = [u*5 u*2 -u 10*u -2*u];
t = linspace(0, 50, length(u));
y = lsim(G, u, t);
close
yn = y + randn(1, length(y)); % Random noise.
sysd = cca(u, yn, 50, t(2), 0, 2);

You are going to see this output

Q =

   3.8349e-09   2.2685e-08
   2.2685e-08   1.8307e-07

R = 1.0207
S =

  -5.0560e-05
  -2.4577e-04

Q1 =

   3.8348e-09   2.2685e-08
   2.2685e-08   1.8307e-07

R1 = 1.0207
>>

As you can see, Q = Q1 and R = R1, but how about S = S1 ? How can I use cov to reproduce S1 that is going to exactly as S ?

Answer 1

The solution was not what I was expecting. I just wanted to use cov command from MATLAB to find Q, R, S matrices for the Algebraic Riccati Equation from measurement data. The solution given by KBS in their first comment is very simple and simply consists of doing

covariance = cov([w' e']); % Old way = [w*w' w*e'; e*w' e*e']/(N-1);

% Compute Q, R, S for the Riccati equation
Q = covariance(1:n, 1:n);
R = covariance(n+1:n+p, n+1:n+p);
S = covariance(1:n, n+1:n+p);