function a = AR_estimate( x, M )
% do AR estimation
% a = AR_estimate( x, M )
% input
%   x : ( T * 1 ) column data vector
%   M : AR dimension
% output
%   a : AR factor
%

T = length(x);
t = M:(T-1);
A = zeros(M);
b = zeros(M,1);

for i=1:M
  for j=1:M
    A(i,j) = x(t-i+1)'*x(t-j+1);
  end
end

for j=1:M
  b(j) = x(t+1)'*x(t-j+1);
end

a = inv(A)*b;