function ma=durbin_ma(x,p,q)
% Syntax is: ma=durbin_ma(x,p,q)
%   This program finds moving average with Durbin's method
%   INPUTS:
%   x : Given time series
%   q :  order of the moving average model
%	p :  order of the all-pole model used to approximation 1/B(z).
%	     This parameter should be at least 4*q as mentioned in the "Moving
%	     average parameter estimation paper" by Niclas Sandgren
b=durbin(x,p,q);
for i=1:1:q
    X(:,i)=x(q-i+1:end-i+1,1);
end
    ma=((X*b(2:end,1))/sum(b(2:end,1)));