function result = SolveMA (series, nq, feps, imax)
% PURPOSE: Solve for MA values using the Tunnicliffe-Wilson method
% -------------------------------------------------------
% Usage:  result = SolveMA (series, nq, feps, imax)
% Where:    series = a series matrix
%           nq = the order of q, feps = Termination tolerance on the function value
%           imax = Maximum number of iterations allowed
% -------------------------------------------------------
% RETURNS:
%           result = the parameters of the MA part.
% -------------------------------------------------------
% NOTES:    This iterative method is slow to converge.
% -------------------------------------------------------
% ACKNOWLEDGEMENT: The original ARFIMA code was written by Jurgen A Doornik
% -----------------------------------------------------

% The matlab port written and change by:
% Ying Ma
% yingma@email.com
    
    if nargin < 2
        error('Error:The nq  can must be empty.');
    elseif nargin < 3 
        feps = 1e-3;
        imax = 600;
    elseif nargin <4
        imax = 600;
    end

    if nq <= 0
        error('Error:The nq  can must be greater than zero.');
    end
    % solved by Tunnicliffe-Wilson method
    
    varv = std(series);
    varv = varv*varv;
    macf = acf(series, nq) * varv;
    
    ma  = zeros(nq+1,nq+1);
    vb  = zeros(nq+1,1);
    vb(1)= varv;
    
    for i=1:imax
        
        vg = -macf;
        
        for j=1:nq+1
            x1= vb(j:nq+1);
            x2= vb(1:nq-j+2);
            
            ma(j,1:nq+2-j)= x1';
            ma(j,j:nq+1) =  ma(j,j:nq+1) + x2';
            
            vg(j) = vg(j) + x1'*x2;
            
        end
        
        for k=1:nq+1
            vb(k) = vb(k) - vg'/(ma(:,k))';
        end
        
        var     = vb(1)^2;
        
        if i>=3 & max(abs(vg)/ varv) < feps
            break;
        end
        
    end
    
    vq  = vb(2:end)./vb(1);
    
    if nq ==1 & abs(vq(1))>1
        vq(1)  = 1/vq(1);
    end;
    result = FlipMA(vq);   
