function [result] = DoCovar(Method, fLogVar, vParameters, mSeries, vMean, nAR, nMA,  nX, nZ, mX, mZ, fid)
% PURPOSE: estimates the Variance-covariance matrix of the arfima model
% -----------------------------------------------------
% USAGE:     result = DoCovar(Method, fLogVar, vParameters, mSeries, vMean, nAR, nMA,  nX, nZ, mX, mZ, fid)
%  where:    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, nX, nZ = the order number of AR, MA, BETA, GAMMA 
%            mX, mZ = the exogenous variable matrixs 
%            fid    = the handle of a file.
% -----------------------------------------------------
% RETURNS:   result = the variance-covariance matrix for the EML and MPL estimates.
% -----------------------------------------------------             
% NOTE:      
% ----------------------------------------------------
% ACKNOWLEDGEMENT: The original ARFIMA code was written by Jurgen A Doornik
% -----------------------------------------------------

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


    if isempty(vParameters)
        error('error: no Parameter');
    end
    
    [nT, nY] = size(mSeries);
    
    if nargin < 8 
        nX  =   0;
        mX  =  [];
    end
    
    if nargin < 9 
        nZ  =   0;
        mZ  =  [];
    end
    
    if nargin < 12;
       fid = 1;
    end
    
    fprintf(fid,strjust('-------Compute covariance of the second estimation stage------\n','center'));
    
    mserr = zeros(size(vParameters, 1), 1);
    
    if nX  & ~cZ  
        vParams = vParameters(1:1+nAR+nMA);  %remove beta from params
    else
        vParams = vParameters;
    end
    
    if Method == 'EML'
        mv  =   hessian('LogEMLik', vParams, mSeries, nAR, nMA, vMean, nX, nZ, mX, mZ);
    elseif Method =='MPL'
        mv  =   hessian('LogMPLik', vParams, mSeries, nAR, nMA, vMean, nX, nZ, mX, mZ);
    end
    
    covarma =   inv(mv) / nT;
    mserr(1:size(covarma,1))  =   sqrt(diag(covarma));
    
    if  nX & ~nZ
        vsig    = acv(vParameters, nAR, nMA, nT);
        mv      = pacf(vsig, mX);
        covbeta = inv(mv'*mv / exp(fLogVar));
        mserr(size(covarma,1)+1:end)  =   sqrt(diag(covbeta));
    end

    result = real(mserr);
