function [vParameters, fLogvar, vMean, vResidual] = DoEstimation(Method, mSeries, vParameters, nAR, nMA, vMean, mX, mZ, fid) 
% PURPOSE: estimates coefficients, with estimated residual variance for EML and MPL, using BFGS
% -----------------------------------------------------
% USAGE:     [vParameters, fLogvar, vResidual] = DoEstimation(Method, mSeries, vParameters, nAR, nMA, vMean, nX, nZ, mX, mZ, fid)
%  where:    Method  = a string stands for the methods:EML or MPL 
%            mSeries = a series matrix
%            vParameters = a matrix with all coefficients, in order:d, AR, MA, beta, gamma
%            vMean   =  a vector of mean values
%            nAR, nMA = the order number of AR, MA
%            mX, mZ = the exogenous variable matrixs 
%            fid    = the handle of a file.
% -----------------------------------------------------
% RETURNS:   vParameters = a matrix with all coefficients, in order:d, AR, MA, beta, gamma after estimation stage using BFGS convergence
%            fLogVar     = estimated standard errors 
%            vResidual   = the residual variance.
% -----------------------------------------------------             
% NOTE:      
% -----------------------------------------------------

% The matlab port written by:
% Ying Ma
% yingma@email.com

    if nargin < 5
        error('no sufficient parameter for the function');
    end
    
    if Method ~='EML' | Method ~='MPL'
        error('Now only two methods can be chosen from');
    end
    
    if nargin < 6 | isempty(vMean)
        vMean = mean(mSeries);
    end
    
    if any([nAR nMA]) < 0
        error('AR or MA cannot be negative');
    end
    
    nX  =  0; nZ = 0;
    if nargin > 6 & ~isempty(mX);
        nX  =   size(mX, 2);
    end
    
    if nargin > 7 & ~isempty(mZ);
        nZ  =   size(mZ, 2);
    end
    
    if ~nX
        mX = [];
    end
    
    if ~nZ
        mZ = [];
    end
    
    if nargin < 9
        fid = 1;
    end
    
    [nT nY] = size(mSeries);
    
    loglik  = 0;
 
    global LogVar mBeta
    
    LogVar  = 0;
    mBeta   = [];
    
    fprintf(fid,strjust('-------Begin the second estimation stage of ARFIMA-------\n','center'));
    
    if nX  & ~nZ  
        vParams = vParameters(1:1+nAR+nMA);  %remove beta from params
    else
        vParams = vParameters;
    end

    
    if nAR+nMA == 0  % if no parameter need to be estimated
        if Method   ==   'EML'
            fprintf(fid,strjust('-------Exact maximum likelihood function-------\n','center'));
            loglik = LogEMLik(vParLik, mSeries, 0, 0, vMean, mX, mZ);
        elseif Method == 'MPL'
            fprintf(fid,strjust('-------Modified profile likelihood function-------\n','center'));
            loglik = LogMPLik(vParLik, mSeries, 0, 0, vMean, mX, mZ);
        end
    else
        if Method == 'EML' | Method == 'MPL'
            options  =  optimset('fminunc');
            options  =  optimset(options , 'Diagnostics' , 'on');
            options  =  optimset(options , 'Display'     , 'iter');
            options  =  optimset(options , 'LargeScale'  , 'off');
            options  =  optimset(options , 'MaxFunEvals' , 400*(1+nAR+nMA));
            options  =  optimset(options , 'MaxIter'     , 600);
            options  =  optimset(options , 'Hessian'     , 'off');  
            options  =  optimset(options , 'DiffMinChange',  1e-5);
            %info.maxit = 600; 
            
            %info.prt   = 1;
            
            if Method   ==   'EML'
                [vParams, loglik, exitflag, output] = fminunc('LogEMLik', vParams, options, mSeries,  nAR, nMA, vMean, mX, mZ);
                %result = MaxBFGS('LogEMLik',vParams,info, mSeries,  nAR, nMA, vMean, nX, nZ, mX, mZ);
            elseif Method == 'MPL'
                [vParams, loglik, exitflag, output] = fminunc('LogMPLik', vParams, options, mSeries,  nAR, nMA, vMean, mX, mZ);
                %result = MaxBFGS('LogMPLik',vParams,info, mSeries,  nAR, nMA, vMean, nX, nZ, mX, mZ);
            end
            
        end
    end
    
    if nX  & ~nZ  
        vParameters = [vParams; mBeta ];  %add beta from params
    else
        vParameters = vParams ;
    end
    
    fLogvar =   LogVar;
    
    clear global LogVar mBeta 
      
    loglik = -loglik - (1 + log(2*pi)) / 2;	 
    loglik = loglik * nT;
        
    if nMA ~= 0
        vParameters(2+nAR:1+nAR+nMA) = FlipMA(vParameters(2+nAR:1+nAR+nMA));
    end;
   
    if nargout > 2
        mXBeta  = [];
        mZGamma = [];
        if nX 
            mbeta   = vParameters(2+nAR+nMA:1+nAR+nMA+nX);
            mXBeta  = mX*mbeta;
        end
    
        if nZ 
            mgamma  = vParameters(2+nAR+nMA+nX:1+nAR+nMA+nX+nZ);
            mZGamma = mZ*mgamma;
        end
        vResidual = acfilter(mSeries, vMean, vParameters, nAR, nMA, mXBeta, mZGamma);
    end
