Matlab: Unable to perform inverse operation - How to create non-singular square matrix

Question

I am trying Matrix operations. The data is a time series model - Autoregressive model, AR(2) where the model order p =2 represented by the variable Y excited by white Gaussian noise, epsilon. I am stuck with conceptual questions for which I will be grateful for answers. The code computes the autocorrelation matrix after multiplying the 1 sampled time delayed vector Y with the transpose of the 1 samples time delayed vector. The general formula for correlation is E(Y*Y'). In this case I need to do E(Y(1:end)*Y(1:end)'). The answer is a scalar if done this way. But, based on the reply by Walter Roberts mathematical expectation I am able to get a 9 by 9 matrix. However, I am not getting a non-singular matrix.

How do I create a singular matrix for the autocorrelation matrix such that the inverse can be computed? Please help

THE CODE:

clc;
clear all;
 var_eps = 1;
    epsilon = sqrt(var_eps)*randn(5000,1); % Gaussian signal exciting the AR model

   Y(1) = 0.0;
    Y(2) = 0.0;
    for n= 3:5000
    Y(n)=  0.1950*Y(n-1) -0.9500*Y(n-2)+ epsilon(n); %AR(2) model
    end

  y_tminus1 = Y(1:end-1).';
    mult = y_tminus1*y_tminus1';  %This creates a square matrix
    autocorr = xcorr2(mult);  %To perform autocorrelation of 1 sampled lag time series with itself(1 sampled lag)
  inverse_autocorr = inv(autocorr);  **%PROBLEM**        


%Warning: Matrix is singular to working precision. 
  trace_inv=trace(inverse_autocorr);

UPDATE: I am facing this problem since I am trying to implement the terms of the matrix in Expression of Eq(20) and the RHS of Eq(24)

UPDATE : 01 Dec

In this code, the issues, the analytical form of variables

CRB_LHS = (b_transpose_b/N)*tr(inv(E[y(t-1)y(t-1)']))

CRB_RHS = sigma2_v*tr(inv(sum_t =1 to N E[y(t-1)y(t-1)']))

where y(t-1) =

[y(t-1), y(t-2),..,y(t-p)]'

The CRB_LHS is calculated from Eq(19) of the paper

`Y_ARMA(n)=  0.1950*Y_ARMA(n-1)- 0.1950*Y_ARMA(n-2) + b*x_gauss(n);` where

 x_gauss = N(0,1)

The CRB_RHS is calculated from Eq(1)

Y_AR(n) = 0.1950*Y_AR(n-1) -0.9500*Y_AR(n-2) + x_chaos(n); where

x_chaos(n)= f(x(n-1),x(n-2),...,x(n-d)); d = 2

I have assumed the coefficients b = 1 for the white noise in ARMA model of Eq(19).

For each SNR level, snrbd = -3:1:2 I am simulating M number of time series (M independent runs). There are 2 doubts regarding the implementation : (1) for the RHS of Eq(25), in the code the Expectation is taken over the independent runs (according to the suggestion in the answer, provided I understood correctly) - the expression shows that the sum is taken over all the data points & Expectation over M samples for (t-1). This is confusing to implement. (2) The code throws error due to the inverse. Unable to solve this.

The result should be CRB_LHS > CRB_RHS for each snr.

clc;
clear all;
N = 10;
M = 100;  % number of independent runs over which the 

  snrdb = -3:1:2;
for snr_levels = 1:length(snrdb)



for monte_carlo = 1: M
     x_gauss = randn(1,N);
     x_chaos(1) = rand();
% generate the chaotic code
    for i =1 : N
        x_chaos(i+1) = 4*x_chaos(i)*(1-x_chaos(i));
    end

    x_chaos = x_chaos-mean(x_chaos);
    Y_AR(1) = 0.0;
    Y_AR(2) = 0.0;

   Y_ARMA(1) = 0.0;
    Y_ARMA(2) = 0.0;
    for n= 3:N
    Y_ARMA(n)=  0.1950*Y_ARMA(n-1)- 0.1950*Y_ARMA(n-2) + x_gauss(n); %Eq(19) model
    Y_AR(n) =   0.1950*Y_AR(n-1) -0.9500*Y_AR(n-2)  + x_chaos(n); %Eq(1)
    end
    % signalPower_Y_AR = var(Y_AR);
 signalPower_Y_AR = 1;

    sigma2_v = signalPower_Y_AR .*(10^(-snrdb(snr_levels))) ;
    v = sqrt(sigma2_v)*randn(1,length(Y_AR));
    z = Y_AR + v;  %Eq(3)

Y_LHS = Y_ARMA(end-2:end-1).';
Y_RHS = z(end-2:end-1).';

A1(:,monte_carlo) = Y_LHS;
B1(monte_carlo,:) = Y_LHS.';

A2(:,monte_carlo) = Y_RHS;
B2(monte_carlo,:) = Y_RHS.';
end


dimension = length(Y_LHS);
sum_of_products_LHS = zeros(dimension,dimension);
sum_prod_RHS = zeros(dimension,dimension);

for runs = 1:M
A = A1(:,runs);
B = B1(runs,:);
mult_LHS = A*B;

C = A2(:,runs);
D = B2(runs,:);
mult_RHS = C*D;

sum_of_products_LHS = sum_of_products_LHS+ mult_LHS;
sum_of_products_RHS = sum_prod_RHS + mult_RHS;
end 

  b_transpose_b = 1;

Expectation_LHS = mean(sum_of_products_LHS);
% Inverse_LHS = inv(Expectation_LHS);
% trace_LHS = tr(Inverse_LHS);

Expectation_RHS = mean(sum_of_products_RHS);

% Inverse_RHS = inv(Expectation_RHS);
% trace_RHS = tr(Inverse_RHS);
%MANUALLY MAKING A SQUARE MATRIX
size_Inverse = 7;
InverseMatrix_LHS = eye(size_Inverse,size_Inverse);
InverseMatrix_RHS = eye(size_Inverse,size_Inverse);

for i = 1:size_Inverse
    InverseMatrix_LHS(i:size_Inverse,i) = Expectation_LHS(1:size_Inverse-i+1);
    InverseMatrix_LHS(i,i:size_Inverse) = Expectation_LHS(1:size_Inverse-i+1);

    InverseMatrix_RHS(i:size_Inverse,i) = Expectation_RHS(1:size_Inverse-i+1);
    InverseMatrix_RHS(i,i:size_Inverse) = Expectation_RHS(1:size_Inverse-i+1);

end
  trace_LHS = tr(InverseMatrix_LHS);
  CRLB_RHS(snr_levels)=  (b_transpose_b/N).*trace_RHS;

  trace_RHS = tr(InverseMatrix_RHS);
  CRLB_RHS(snr_levels)= sigma2_v*trace_RHS;

 end

Answer 1

For p=2, expression for y_tminus1 should be (see expressions after Eq.1 in the paper)

y_tminus1 = Y(end-2:end-1).' 

Your mult is OK. For Eq. 20 you need to take expectation E(mult). For this you need to generate multiple paths and take an average over them. For the RHS of Eq. 25, you need to use y_tminus1 for every step of your ARMA process, sum corresponding mult matrices, take expectation over ensemble and only after that take an inverse. Maybe, try to adjust your code along these lines and i'll correct it.