function [b, e] = MA(x, q)
%用Durbin算法計算MA(q)模型參數系數，返回e(1+B)中的B，e是sigma^2
L = 32; %round(exp((log(q)+log(N))/2))  % q<<L<<N
N = length(x);
[a, e] = AR2(x, L);
a = [1; a];
for i=1:q
    for j=1:q
        Raa(i, j) = 0;
        for n=1:L-abs(i-j)+1
            Raa(i,j) = Raa(i,j)+a(n)*a(n+abs(i-j));
        end
    end
    raa(i) = 0;
    for n=1:L-i+1
        raa(i) = raa(i)+a(n)*a(n+i);
    end
end
Raa = Raa/(L+1);
raa = raa/(L+1);
b = -1*pinv(Raa)*raa';