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
For p=2, expression for y_tminus1 should be (see expressions after Eq.1 in the paper)
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.