function [mforcasts, mrmses, mes, mforclevel, mrmselevel, melevel] = Doforecast(bnaiveonly, nTBegin, nTForc, mSeries, fLogVar, vParameters, nAR, nMA, vMean, mX, mZ, vLevel, fid)
% PURPOSE: The forecasting function in the arfima model
% -----------------------------------------------------
% USAGE:     [mforcasts, mrmses, mes, mforclevel, mrmselevel, melevel] = Doforecast(Method, nTBegin, nTForc, mSeries, vMean, fLogVar, vParameters, nAR, nMA, nX, nZ, mX, mZ, fid)
%  where:    Method = Naive forecasts or best linear prediction
%            nTBegin, nTForc = The estimation period from nTBegin to nTForc
%            mSeries = a series matrix
%            vMean   =  a vector of mean values
%            fLogVar     = estimated standard errors 
%            vParameters = a matrix with all coefficients, in order:d, AR, MA, beta, gamma
%            nAR, nMA = the order number of AR, MA
%            mX, mZ = the exogenous variable matrixs 
%            fid    = the handle of a file.
% -----------------------------------------------------
% RETURNS:   mforcasts = a forecasts series matrix produced by the two method
%            mrmses = a root mean square error matrix
%            mes    = a error matrix
%            mforclevel = a level forecasts series matrix produced by the two method
%            mrmselevel = a level forecasts root mean square error matrix
%            melevel    = a level forecasts error matrix
% -----------------------------------------------------             
% NOTE:     
% -----------------------------------------------------

% The matlab port written by:
% Ying Ma
% yingma@email.com
    
    if nargin < 8
        error ('The number of parameters provided is not enough');
    end
    
    if isempty(vParameters)
        error('The parameter  cannot be empty.')
    end
    
    
    if nTForc <= 0 | nTBegin <= 1
        return;
    end
    
    [nT nY] = size(mSeries);
    
    if nT < nTBegin-1
        error('The beginning time cannot be bigger than sample time.');
    end
   
    if  nTBegin ==-1
        nTBegin = nT + 1;
    else
        nT  = nTBegin -1;
    end
    
    if nargin < 9 | isempty(vMean)
        vMean = mean(mSeries);
    end
    
    nX = 0; nZ = 0;
    if nargin > 9 & ~isempty(mX);
        nX  =   size(mX, 2);
    end
    
    if nargin > 10 & ~isempty(nZ);
        nZ  =   size(mZ, 2);
    end
    

    if  nX
        if nTBegin+nTForc-1 > size(mX,1)
            error('No enough regressors over forecast period');
        end
        mX1 =  mX(nTBegin:nTBegin+nTForc-1);
        mX  =  mX(1:nTBegin-1);
    end
    
    if nargin < 12
        vLevel = [];
    end
    
    if nargin < 13
        fid = 1;
    end
    
    fprintf(fid,strjust('-------Begin the state of prediction over  forecast period ------\n','center'));
    
    sigma2 = exp(fLogVar);
    nTsample =   min(nTForc, size(mSeries,1)-nT);
    
    mXForc =  ones(nTForc,1)*vMean;
    mXBeta =  ones(nT,1)*vMean;
    
    %------ get XBeta over fitted and forecast period
    if nX
        mXBetaForc = mX1 * vParameters(2+nAR+nMA:1+nAR+nMA+nX);
        mXBeta  = mXBeta + mX * vParameters(2+nAR+nMA:1+nAR+nMA+nX);
        mXForc = mXForc + mXBetaForc;
    end
    
    %-------AR coefs computation 
    
    %-------AR coefs of fractional part
    arcoefs    = powerdiff(vParameters(1), [1; zeros(nT + nTForc-1, 1)]);
    
    %-------AR coefs of AR part
    if nAR
        arcoefstemp = conv(arcoefs, [1; -vParameters(2:1+nAR)]);
        arcoefs = TrimMatrix(arcoefstemp,nT+nTForc);
    end
    
    %-------AR coefs of MA part
    if nMA
        arcoefstemp = deconv([arcoefs;zeros(nMA,1)], [1; vParameters(2+nAR:1+nAR+nMA)]);
        arcoefs = TrimMatrix(arcoefstemp,nT+nTForc);
    end
    
    %-------MA coefs of MA part
    macoefstemp =   deconv([1;zeros(2*(nTForc-1),1)], arcoefs(1:nTForc));
    macoefs = TrimMatrix(macoefstemp,nTForc);
    
    %-------Store arcoefs in reverse order - speed up computation 
    arcoefs = -flipud(arcoefs(2:end));
    
    %-------Naive forecasts
    
    my   =  mSeries(1:nT) - mXBeta;
        
    my   =  [my; zeros(nTForc,1)];
    
    for i=1:nTForc
         my(nT+i) = arcoefs(nTForc-i+1:end)'* my(1:nT+i-1);
    end
    
    mnaive  =   my(nT+1:end) + mXForc;
        
    mrmsen  =   sqrt(sigma2* cumsum((macoefs.^2)));
    
    
    %-------The best linear prediction    
    if ~bnaiveonly
        
        nTT     =   nT+nTForc;
        vsig    =   acv( vParameters, nAR, nMA, nTT);
    
        mgamma  = zeros(nT, nTForc);
    
        for k = 1:nTForc
           mgamma(:,k)  = vsig(k:nT+k-1); 
        end
    %-------Flip each columns of mgamma in reverse order
        
        mgamma  = flipud(mgamma);
    
        vsig    = vsig(1:nT);

        mq  =   inv(toeplitz(vsig))*mgamma;
    
        mforc = mq'*(mSeries(1:nT) - mXBeta) + mXForc;
    
        mrmse = sqrt(sigma2)*sqrt(abs(vsig(1)-diag(mgamma'*mq)));

    end
    
    if bnaiveonly
        mforcasts = mnaive;
        mrmses = mrmsen;
    else
        mforcasts = mforc;
        mrmses = mrmse;
        
    end
    
    if nargout >2
        if  nTsample
            mes      =   mSeries(nTBegin:nTBegin+nTsample-1) - mforcasts(1:nTsample);
        end
    end
    
    mforclevel = [];
    mrmselevel = [];
    melevel    = [];
    
    %-------integrate the forecasts to levels
    
    if  ~isempty(vLevel)
        [nR nC] = size(vLevel);
        if nR < nC
            vLevel = vLevel';
            [nR nC] = size(vLevel);
        end
        
        vinitlevel = zeros(nR,1);
       
        for k = 1:nR
            vinitlevel(nR - k +1) = vLevel(nR);
            vLevel = [0; diff(vLevel)];
        end
        
        mforcvar = [];
        melevel  = [];
        
        %------------- compute the forecast errors covariance matrix
        if ~bnaiveonly
            mforcvar = toeplitz(vsig(1:nTForc)) - mgamma'*mq;
        end
        
        mact = mSeries(nTBegin:nTBegin+nTsample-1);
        
        for j=1:nR
            macoefs = cumsum(macoefs);
            mnaive  = vinitlevel(j) + cumsum(mnaive);
            mrmsen  = sqrt(cumsum(sigma2*macoefs.^2));
            
            if nTsample
                mact = vinitlevel(j) + cumsum(mact);
            end
            
            if ~bnaiveonly
                mforc = vinitlevel(j) + cumsum(mforc);
                
                %-----------compute  variance of cumsums of forecast errors 
                [nrcum nccum]  = size(mforc);
                mcum = lagmatrix(ones(nrcum,1),0:nrcum - 1);
                
                mforcvar = mcum * mforcvar * mcum';
                mrmse   = sqrt(sigma2 *diag(mforcvar));
                
            end
        end
        
        if bnaiveonly
            mforclevel = mnaive;
            mrmselevel = mrmsen;
        else
            mforclevel = mforc;
            mrmselevel = mrmse;
        end
        
        if nargout > 2
            if  nTsample
                melevel      =   mforclevel - mact;
            end
        end
        
    end
    
    
    