function [vParameters, vMean, logvar] = DoInit(bRstrD, mSeries, nAR, nMA, vMean, mX, mZ, fid )
% PURPOSE: The initialization of the parameters for arfima model
% -----------------------------------------------------
% USAGE:     [vParameters, logvar] = DoInit(bRstrD, mSeries, nAR, nMA, vMean, nX, nZ, mX, mZ, fid )
%  where:    bRstrD = a flag for restricting the value of d 
%            mSeries = a series matrix
%            nAR, nMA = the order number of AR and MA parts
%            vMean   =  a vector of mean values
%            mX, mZ = the exogenous variable matrixs 
%            fid    = the handle of a file.
% -----------------------------------------------------
% RETURNS:   vParameters= a matrix with all estimated coefficients, in the initial stage
%            vMean      = a vector of mean values
%            logvar     = a log-variance
% -----------------------------------------------------             
% NOTE:      the function will be used for the initialization or estimation 
% -----------------------------------------------------
% 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 < 4
        error('Error:The number of parameters provided can not be sufficient');
    end
    
    if any([nAR nMA]<0)
        error('Error:AR or MA cannot be negative');
    end
    
    if nargin < 5 | isempty(vMean)
        vMean = mean(mSeries);
    end
    
    nX  = 0; nZ = 0;
    if nargin > 5 & ~isempty(mX);
        nX  = size(mX, 2);
    end
    
    if nargin > 6 & ~isempty(mZ);
        nZ  = sizes(mZ, 2);
    end
    
  
    if nargin < 8
        fid = 1;
    end
    
    fprintf(fid,strjust('-------Begin initial stage of ARFIMA-------\n','center'));
    
    [nT nY] =   size(mSeries); 
% first the mSeries data are filtered by mean value
    mu  = mSeries - ones(nT,1)*vMean;
    
    vParameters = zeros(1+nAR+nMA+nX+nZ,1);
    
% estimate the start beta coefficient value, now regress residuals on X
    fprintf(fid,strjust('-------First estimate the start beta coefficient value-------\n','center'));
    if nX > 0
        [vBeta, mu] = olsv(mu, mX);
        vParameters(2+nAR+nMA:1++nAR+nMA+nX) = vBeta;
        fprintf(fid,'Beta(n):%4.3f\n', vBeta);
    end
   
    
% find the starting value for d 
    if nAR+nMA > 0
        Trunc   = sqrt(nT);
    else
        Trunc   = nT;
    end
    
    [d, seasy, seprd] = EstimateGPH(mu,Trunc);
    
    fprintf(fid,strjust('-------Second determine starting value for d-------\n','center'));
    fprintf(fid,'starting value for d = %6.4f its significant at %4.3f\n',d,seprd);

    if bRstrD
% not adjust the value for d if significant below or equal to the level of 5%     
        if seprd > 0.05
            d = 0;
        end
% restrict d to stationary values
        if d > 0.45
            d   =   0.4;
        end
    
        if d < -0.45
            d   =   -0.4;
        end
        
        fprintf(fid,'modified value for d = %6.4f (because of not the stationary value below the confidence level or not stationary values)\n',d);
    end
    
    fprintf(fid,strjust('-------Third compute u_t = (1-L)^d y_t -----\n','center'));
    if abs(d) >1e-5;
        mu  = powerdiff(d,mu);
        mu(1,:) = zeros(1,nY);
    end
    
    vParameters(1)  =   d;
    
    fprintf(fid,strjust('-------Fourth estimate starting value for AR parameters -----\n','center'));
    if nAR
        mAR = SolveAR(mu, nAR, nMA);
        vParameters(2:1+nAR)  =   mAR;
        if ~IsStationary(mAR)
            mAR = SolveAR(mu,nAR, 0);
        end
        phi =   [1; -mAR]';
        mu  =   filter(phi,1,mu);
        mu(1:nAR,:)  =   zeros(nAR,size(mu,2));
        fprintf(fid,'phi(n):%4.3f\n', mAR);
    end
    
    
    
    fprintf(fid,strjust('-------Fifth estimate starting value for MA parameters -----\n','center'));
    if nMA
        mMA = SolveMA(mu, nMA);
        vParameters(2+nAR:1+nAR+nMA)  =   mMA;
        theta = [1; mMA]';
        mu  =   filter(1,theta,mu);
        fprintf(fid,'theta(n):%4.3f\n', mMA);
    end
    

     
    %now regress residuals on Z
    if nZ > 0
        [mGamma, mu] = olsv(mu, mZ);
        vParameters(2+nAR+nMA+nX:1++nAR+nMA+nX+nZ) = mGamma;
        
    end
    
    fprintf(fid,strjust('-------Finally estimate starting value for sigma parameters -----\n','center'));
    varv = std(mu);
    varv = varv*varv;
    fprintf(fid,'sigma2:%4.3f\n', varv);
    logvar = log(varv);
