function result= LogMPLik(vParameters, mSeries, nAR, nMA, vMean, mX, mZ)   
% PURPOSE: The modified profile likelihood function for arfima
% -----------------------------------------------------
% USAGE:     result = LogMPLik(vParameters, mSeries, nAR, nMA, vMean, nX, nZ, mX, mZ) 
%  where:    mSeries = a time series matrix
%            vParameters = a matrix with all coefficients, in order:d, AR, MA excluding Beta and Gamma
%            nAR, nMA = the order of AR and MA parts
%            vMean   =  a vector of mean values
%            nX, nZ = the order of Beta, Gamma
%            mX, mZ = the exogenous variable matrixs 
%            fid    = the handle of a file.
% -----------------------------------------------------
% RETURNS:   result = The time series likelihoods 
% -----------------------------------------------------             
% NOTE:      The function is used in the maximum likelihood procedure.
% -----------------------------------------------------
% ACKNOWLEDGEMENT: The original ARFIMA code was written by Jurgen A Doornik
% -----------------------------------------------------

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


    if nargin < 5
        error('Error:no sufficient parameters');
    end
    
    nX = 0; nZ = 0;
    if nargin > 5 & ~isempty(mX);
        nX  =   size(mX, 2);
    end
    
    if nargin > 6 & ~isempty(nZ);
        nZ  =   size(mZ, 2);
    end

    
    [nT nY] = size(mSeries);

    d   =   vParameters(1);

    if  abs(d) <= 1e-9 | ~nAR % if d not 0: checks the coefficients of AR, otherwise, no check
        mAR = [];
    else
        mAR = vParameters(2:nAR+1);  
    end
    
    if ~IsStationary(mAR,d)
        result = d*d + mAR'*mAR; % The penalty function is defined here
        return;
    end
    
    mXBeta  = [];
    
    mZGamma = [];
    
    logdetx = 0;
    
    global LogVar mBeta

    
    if nX & ~nZ % filter mX 
        [ve, logdet, vsig]  = acfilter(mSeries, vMean, vParameters, nAR, nMA, [], []);
     
        me = pacf(vsig, mX);
        
        mp    = me'*me;
        logdetx  = log(det(mp));
        
        minv  = inv(mp);
        mBeta   =   (ve'*me)*minv;
        drss    =   ve'*ve - mBeta * (me'*ve);
            
    else
        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
        
        [ve, logdet, vsig]  = acfilter(mSeries, vMean, vParameters, nAR, nMA, mXBeta, mZGamma);
        drss    =   ve'*ve;
        
        if nX
            me = pacf(vsig, mX);        
            mp    = me'*me;
            logdetx  = log(det(minv));
        end
        
    end
   
    if isempty(ve)| all(ve == 0)
        result = d*d + mAR'*mAR;
        return;
    end
    
    invt    =   1 / nT; 
    
    LogVar  = log(drss / (nT-nX));
    
    result  = (( 0.5 - invt) * logdet + (0.5 * (nT - nX) - 1)* log(drss*invt) + 0.5 *logdetx)*invt;
    