How to plot the residual analysis plot manually

I would like to compute the residual analys plot manually, for instance the output of the following command:

Z = iddata(u, y, dt); % u are user input data, y are user output data
sys = oe(Z, [nb nf nk]) % nb nf nk is the user model structure to be identified
resid(Z, sys); % plot the built-in plot

What I tried:

[residuals, residuals_corr] = resid(Z, sys);

m = 25; % default number of lags
xx = -m : m;
Ree = residuals_corr(:,1,1); % autocorellation of the residuals
Reu = residuals_corr(:,2,1); % crosscorrelation between residuals and inputs
Ree0 = Ree(1); % variance of the residuals (autocorrellation at lag 0)
Ruu0 = residuals_corr(1,2,2); % variance of the inputs (autocorrellation at lag 0)
Ree_norm = Ree./Ree0; % normalization
Reu_norm = Reu./((Ree0*Ruu0)^(0.5)); % normalization
% alternative
Ree_norm2 = xcorr(residuals.OutputData,m,'normalized');
Reu_norm2 = xcorr(residuals.OutputData,residuals.InputData,m,'normalized');

% white test
p = 0.99;
CIy = norminv(p+(1-p)/2)*1/sqrt(N); % N is the number of obsevations in the input and output data, e.g. the length

% my plot
figure; t = tiledlayout(1, 2, "TileSpacing", "compact", "Padding", "compact");
title(t,'Residual Analysis', 'FontSize', 18);

nexttile(); hold on;
fill([-m m m -m], [CIy CIy -CIy -CIy], 'b', 'FaceAlpha', 0.1)
%stem(xx, [flip(Ree_norm); Ree_norm(2:end)]) does not work
stem(xx, Ree_norm2, "filled", "LineWidth", 2, "Color", "b") % works the same
ylim([-0.2 1]); xlabel("Lag $\tau$", "Interpreter", "latex", 'FontSize', 18);
ylabel("$\hat{R}_{\epsilon}(\tau)$", "Interpreter", "latex", 'FontSize', 18)

nexttile(); hold on;
fill([-m m m -m], [CIy CIy -CIy -CIy], 'b', 'FaceAlpha', 0.1)
%stem(xx, [flip(Reu_norm); Reu_norm(2:end)]) % does not work
stem(xx, Reu_norm2, "filled", "LineWidth", 2, "Color", "b") % works the same
ylim([-0.2 1]); xlabel("Lag $\tau$", "Interpreter", "latex", 'FontSize', 18);
ylabel("$\hat{R}_{\epsilon u}(\tau)$", "Interpreter", "latex", 'FontSize', 18)

The approach using the outputs of the resid function works the same while computing the crosscorrelation myself from input and output data fails at obtaing the same plot of the resid function. What am I missing?

[Jan Kappen]

I can't reproduce the data from residuals_corr, but via xcorr:

This reproduces the first cross correlation plot.

[E,R]=  resid(zData, sys);
x = xcorr(E.y(:,1),E.u(:,1), 25, "normalized"), stem(-25:25, x)

[ArchLucs]

That is the same result I got, isn't it. I don't understand how to reproduce the same data from the output of the built-in function.

Another question: do you know how to compute the confident band for the cross-correlation between noise and input?

Answer 1

[Milan Bansal]

Hi ArchLucs

I understand that you are not able to plot the same plot for residual auto-correlation and cross-correlation using the input and output data manually, as plotted by the "resid" function.

In the given code, it looks like the calculation used for normalization is not correct. You may calculate the auto-correlation and cross-correlation and perform normalization manually instead of using the correlation outputs of the "resid" function.

Please refer to the following code snippet to calculate and plot the residual auto-correlation and cross-correlation manually without using the built-in function "xcorr".

% Get the residuals
[E, R] = resid(Z, sys);
x = E.InputData; 
y = E.OutputData;

m = 25; % default number of lags
xx = -m : m;

%% 1.) Calculate Auto-Correlation
n = length(y); % Length of the data series
yZeroMean = y - mean(y); % Zero-mean series

% Preallocate array for autocorrelation values
acf = zeros(2*n-1, 1);
lags = -n+1:n-1;

% Calculate autocorrelation for each lag
for k = 1:length(lags)
    lag = lags(k);
    if lag < 0
        productSum = sum(yZeroMean(1-lag:end) .* yZeroMean(1:end+lag));
    else
        productSum = sum(yZeroMean(1:end-lag) .* yZeroMean(1+lag:end));
    end
    acf(k) = productSum;
end
% Normalize the autocorrelation values
acf = acf / (n * var(y, 1)); % Normalize by n*variance for consistency with xcorr 'coeff'

acflen = lenght(acf);
stem(xx, acf((acfLen+1)/2 -m:(acfLen+1)/2 + m ));
xlabel('Lag');
ylabel('Autocorrelation');


%% 2.) Calculate Cross-Correlation
nx = length(x);
ny = length(y);
n = max(nx, ny); % Use the longer length if they differ

% Ensure both series are zero-mean
xZeroMean = x - mean(x);
yZeroMean = y - mean(y);

% Preallocate array for cross-correlation values
% The total number of lags is nx+ny-1 (from -ny+1 to nx-1 for equal lengths)
ccf = zeros(nx+ny-1, 1);
lags = -ny+1:nx-1;

% Calculate cross-correlation for each lag
for k = 1:length(lags)
    lag = lags(k);
    if lag < 0
        productSum = sum(xZeroMean(1:end+lag) .* yZeroMean(1-lag:end));
    else
        productSum = sum(xZeroMean(1+lag:end) .* yZeroMean(1:end-lag));
    end
    ccf(k) = productSum;
end

% Normalize the cross-correlation values
ccf = ccf / sqrt(sum(xZeroMean.^2) * sum(yZeroMean.^2)); % Normalize for consistency

ccfLen = length(ccf);
stem(xx, ccf((ccfLen+1)/2 -m:(ccfLen+1)/2 + m ));
xlabel('Lag');
ylabel('Cross-Correlation');

Please refer to the following documentation link to learn more about "resid" function.

https://www.mathworks.com/help/ident/ref/resid.html

Please refer to the following documentation link to learn more about "xcorr" function.

https://www.mathworks.com/help/matlab/ref/xcorr.html

Hope this helps