Before I go into the specific issues, I'd like to note that there is a minor confusion with respect to M not representing the AR process order. More specifically an AR process of order M=4 would use the M previous values of y1 rather than M+1 values (in the loop with index j). So, with M=4 your implementation is really that of an AR process of order $p = 5$.
Coefficients and stability
The problem is that I need to multiply the previous output y1 and current input x with a=0.1, otherwise it gives an unstable filter which leads to infinity. Please clarify which mistake I'm making while filtering x using AR filter or I'm doing it the wrong way?
The stability of the AR process can be determined from the location of the poles of the transfer function. To find those poles, we start off with the recurrence equation of the AR(p) process:
$$
y[n] = a x[n] + a \sum_{j=1}^{p} y[n-j]
$$
From which the transfer function in the $z$-plane is:
$$
\begin{align}
H(z) &= \frac{a}{1-a\sum_{j=1}^{p} z^{-j}} \\
&= \frac{az^{p}}{z^{p}-a\sum_{j=0}^{p-1} z^j}
\end{align}
$$
Hence the poles are given as the roots of the $z^{p}-a\sum_{j=0}^{p-1} z^j$ polynomial. In your specific case of $p = 5$, the pole location are given by the following plot:
As you may notice, there is a pole near $z=2$ (more precisely at approximately 1.96595) when $a=1$ which indicates that the system is not stable. On the other hand, when $a=0.1$, all the poles are located inside the unit circle and the system is stable. Beware that even with $a=0.1$, starting with $p=11$, the system also has a pole outside the unit-circle and starts to be unstable.
In other words, there are no specific mistaking in your filtering code (ie. it correctly implements an unstable filter when a=1).
Yule-Walker equation
But its not giving correct coefficients (or filter taps). Please somebody help me if there's something I'm missing in understanding the question fully or making some mistake somewhere in the code?
There are indeed a few issues with that part of the code. Perhaps the most significant issue is that the solution
$$
\begin{align}
\vec{a} &= \mathbf R^{-1} \vec r \\
&= \begin{bmatrix}
r_0 & r_1 &\cdots & r_{p-1} \\
r_1 & r_0 &\cdots & r_{p-2} \\
\vdots & \vdots & \ddots & \vdots \\
r_{p-1} & r_{p-2} &\cdots & r_0 \\
\end{bmatrix}^{-1}
\begin{bmatrix}
r_1 \\
r_2 \\
\vdots \\
r_p
\end{bmatrix} \tag{1}
\end{align}
$$
where
$$
\begin{align}
r_j &= \frac{c_j}{c_0} \\
&= \frac{E\left\{y_n y_{n-j}\right\}}{E\left\{y^2_n\right\}}
\end{align}
$$
applies to an AR process of the form
$$
y[n] = x[n] + \sum_{j=1}^p a_j y[n-j]
$$
where $x[n]$ is a zero-mean random variable. Unfortunately, the line
x = rand(1,1000);
Generates a uniformly distributed random variable in the range [0,1] with a mean of $\frac{1}{2}$. An example of random variable with zero-mean could be a uniformly distributed random variable in the range [-1,1]:
x = 2*rand(1,1000)-1;
For an input which does not have a mean of zero, or equivalently for an AR process of the form:
$$
y[n] = K + x'[n] + \sum_{j=1}^p a_j y[n-j]
$$
where $x'[n] = x[n] - K$ is zero-mean, the equations can be updated to:
$$
\begin{align}
\begin{bmatrix}
\mu \\
c_1 \\
c_2 \\
\vdots \\
c_p
\end{bmatrix} &=
\begin{bmatrix}
1 & \mu & \mu &\cdots & \mu \\
\mu & c_0 & c_1 &\cdots & c_{p-1} \\
\mu & c_1 & c_0 &\cdots & c_{p-2} \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
\mu & c_{p-1} & c_{p-2} &\cdots & c_0 \\
\end{bmatrix}
\begin{bmatrix}
K \\
a_1 \\
a_2 \\
\vdots \\
a_p
\end{bmatrix} \tag{2}
\end{align}
$$
where $\mu = E\left\{y_n\right\}$.
Implementation for zero-mean AR process
Finally there are a few small issues in the implementation of the equation. For one thing, the vector returned by autocorr(y2,p) contains p+1 values: from 0 to p lags. You need the first p values for the $\mathbf R$ matrix (leaving out lag p), and the last p values for the $\vec r$ vector (leaving out the zero lag value since $\vec r$ starts at $r_1$, not $r_0$).
So you can construct the Yule Walker equation with:
R = toeplitz(r(1:end-1)); % r(0), ..., r(p-1)
w = inv(R)*r(2:end); % as Yule Walker eq: w=R^-1*r
So putting it all together, you should get something like this:
x = 2*rand(1,1000)-1; % uniformly distributed over [-1,1]
a=0.1
p=5;
M=p-1;
for i = 1:length(x)
sum=0;
y1(i)=0;
for j = 1:M+1
if (i-j>0)
sum = sum + a*y1(i-j);
end
end
y1(i)=sum +a*x(i);
end
y2 = awgn(y1,25);
r = autocorr(y2,p)';
% Or if you don't have autocorr, replace with the following 2 lines:
% r = xcorr(y2,p)';
% r = r(p+1:end); % keep only non-negative lags
R = toeplitz(r(1:end-1)); % r(0), ..., r(p-1)
w = inv(R)*r(2:end); % as Yule Walker eq: w=R^-1*r
Implementation for non-zero-mean AR process
The solution of equation (2) can be constructed in a similar fashion as for equation (1). The main thing to note is that in this case the $\mathbf R$ matrix elements $c_0, \dots, c_{p-1}$ are not normalized by $c_0$ (which would have yielded $r_0, \dots, r_{p-1}$). In this case you may replace autocorr with xcorr, which returns a sum of product without normalization (you would then need some adjustment since is xcorr is double sided, meaning it contains 2*p-1 values, including the values for positive, zero and negative lags). This would give you the following:
x = rand(1,1000);
...
mu = mean(y2);
r = xcorr(y2,p)';
r = r(p+1:end); % keep only non-negative lags
R = toeplitz(r(1:end-1)); % r(0), ..., r(p-1)
A = [1 mu*ones(1,p);mu*ones(p,1) R];
b = [mu;r(2:end)];
w = inv(A)*b;
% separate the bias term K from the AR process weights:
K = w(1);
w = w(2:end);
