An Evaluation of ARFIMA (Autoregressive Fractional Integral Moving Average) Programs ^†

Abstract

Strong coupling between values at different times that exhibit properties of long range dependence, non-stationary, spiky signals cannot be processed by the conventional time series analysis. The autoregressive fractional integral moving average (ARFIMA) model, a fractional order signal processing technique, is the generalization of the conventional integer order models—autoregressive integral moving average (ARIMA) and autoregressive moving average (ARMA) model. Therefore, it has much wider applications since it could capture both short-range dependence and long range dependence. For now, several software programs have been developed to deal with ARFIMA processes. However, it is unfortunate to see that using different numerical tools for time series analysis usually gives quite different and sometimes radically different results. Users are often puzzled about which tool is suitable for a specific application. We performed a comprehensive survey and evaluation of available ARFIMA tools in the literature in the hope of benefiting researchers with different academic backgrounds. In this paper, four aspects of ARFIMA programs concerning simulation, fractional order difference filter, estimation and forecast are compared and evaluated, respectively, in various software platforms. Our informative comments can serve as useful selection guidelines.

1. Introduction

Humans are obsessed about their future so much that they worry more about their future more than enjoying the present. Time series modelling and analysis are scientific ways to predict the future. When dealing with empirical time series data, it usually comes to the classic book of Box and Jekin’s methodology for time series models in the 1970s, in which it introduced the autoregressive integrated moving average (ARIMA) models to forecast and predict the future behavior [1,2]. However, the ARIMA model as well as Poisson processes, Markov processes, autoregressive (AR), moving average (MA), autoregressive moving average (ARMA) and ARIMA processes, can only capture short-range dependence (SRD). They belong to the conventional integer order models [3].

In time series analysis, another traditional assumption is that the coupling between values at different time instants decreases rapidly as the time difference or distance increases. Long-range dependence (LRD), also called long memory or long-range persistence, is a phenomenon that may arise in the analysis of spatial or time series data. LRD was first highlighted in the hydrological data by the British hydrologist H. E. Hurst, and then the other statistics in econometrics, network traffic, linguistics and the Earth sciences, etc. LRD, which is characterized by the Hurst parameter, means that there is a strong coupling effect between values at different time separations. Thus, LRD also indicates that the decay of the autocorrelation function (ACF) is algebraic and slower than exponential decay so that the area under the function curve is infinite. This behavior can be also called inverse power-law delay. Different from the analytical results of linear integer-order differential equations, which are represented by the combination of exponential functions, the analytical results of the linear fractional-order differential equations are represented by the Mittag–Leffler function, which intrinsically exhibits a power-law asymptotic behavior [4,5,6].

Due to the increasing demand on modeling and analysis of LRD and self-similarity in time series, such as financial data, communications networks traffic data and underwater noise, the fractional order signal processing (FOSP) technique is becoming a booming research area. Moreover, fractional Fourier transform (FrFT), which is the generalization of the fast Fourier transform (FFT), has become one of the most valuable and frequently used techniques in the frequency domain of the fractional order systems [3].

Compared to the conventional integer order models, the ARFIMA model gives a better fit and result when dealing with the data which possess the LRD property. Sun et al. applied the ARFIMA model to analyze the data and predict the future levels of the elevation of Great Salt Lake (GSL) [7]. The results showed that the prediction results have a better performance compared to the conventional ARMA models. Li et al. examined four models for the GSL water level forecasting: ARMA, ARFIMA, autoregressive conditional heteroskedasticity (GARCH) and fractional integral autoregressive conditional heteroskedasticity (FIGARCH). They found that FIGARCH offers the best performance, indicating that conditional heteroscedasticity should be included in time series with high volatility [8]. Sheng and Chen proposed a new ARFIMA model with stable innovations to analyze the GSL data, and predicted the future levels. They also compared accuracy with previously published results [9]. Contreras-Reyes and Palma developed the statistical tools afmtools package in R for analyzing ARFIMA models. In addition, the implemented methods are illustrated with applications to some numerical examples and tree ring data base [10]. Baillie and Chung considered the estimation of both univariate and multivariate trend-stationary ARFIMA models, which generated a long memory autocorrelated process around a deterministic time trend. The model was found to be remarkably successful at representing annual temperature and width of tree ring time series data [11]. OxMetrics is an econometric software including the Ox programming language for econometrics and statistics, developed by Doornik and Hendry. Several papers and manuals are available for the ARFIMA model with OxMetrics [12,13,14].

Nowadays, there are lots of numerical tools available for the analysis of the ARFIMA processes since these applications are developed by different groups based on different algorithms and definitions of accuracies and procedures. As a consequence, the estimation and prediction results may be different or even conflicting with others. For the scholars or engineers who are going to do the modeling work of the ARFIMA processes, they might get confused as to which tool is more suitable to choose. Thus, we have evaluated techniques concerning the ARFIMA process so as to provide some guidelines when choosing appropriate methods to do the analysis. With this motivation, this paper briefly introduces their usage and algorithms, evaluates the accuracy, compares the performance, and provides informative comments for selection. Through such efforts, it is hoped that informative guidelines are provided to the readers when they face the problem of selecting a numerical tool for a specific application.

For one thing, many publications about fractional systems dynamics use their novel fractional order calculus ideas to represent with encouraging results [15]. However, in reality, when it comes to engineers with zero background, they do not even know which tool to start to use.

The rest of the paper is organized as the follows: Section 2 introduces the basic mathematics of LRD and the ARFIMA model. Section 3 gives a brief review and description on the software commonly used for the analysis of the ARFIMA processes. In Section 3, the quantitative performances of the tools are evaluated and compared in four primary categories—simulation, processing, estimation and prediction in the ARFIMA process. Conclusions are given in Section 4 and Section 5.

2. LRD and ARFIMA Model

When the hydrologist H.E. Hurst spent many years analyzing the records of elevation of the Nile river in the 1950s, he found a strange phenomena: the long-range recording of the elevation of the Nile river has much stronger coupling, and the autocorrelation function (ACF) decays slower than exponentially [16]. In order to quantify the level of coupling, the rescaled range (R/S) analysis method was provided to estimate the coupling level, which is now called the Hurst parameter. Furthermore, many valuable Hurst parameter estimators were provided to more accurately characterize the LRD time series [17]. Since then, the LRD or long memory phenomenon has attracted numerous research studies. Based on Hurst’s analysis, more suitable models, such as ARFIMA and fractional integral generalized autoregressive conditional heteroscedasticity (FIGARCH) were built to accurately analyze LRD processes.

The rescaled range (R/S) method is one of the time-domain analysis of Hurst parameter defined as follows [16]:

$$E{[\frac{R\left(n\right)}{S\left(n\right)}]}_{n\to \infty }=C{n}^H,$$

(1)

where $E(·)$ denotes the expected value of the observations, $R\left(n\right)$ is the range of the first n values, $S\left(n\right)$ is their standard deviation, and C is a constant. Whittle’s Maximum Likelihood Estimator (MLE) and wavelet analysis using periodogram based analysis in the frequency domain [18].

Autocorrelation function (ACF) analysis is one of the useful techniques for identifying trends and periodicities in the data, in a manner that is often more precise than can be obtained with simple visual inspections. In addition, LRD or long memory property can be defined by ACF.

Let $\left\{X\right(t);t\in (-\infty ,+\infty \left)\right\}$ and the ACF $\rho \left(k\right)$ is defined as:

$$\rho \left(k\right)=\frac{Cov(X_t,X_{t-k})}{Var\left(X_t\right)},$$

(2)

where $Cov(·)$ is the covariance and $Var(·)$ is the variance.

A stationary time series defined over $t=0,1,2,3\cdots$ is said to be long memory if ${\sum }_{k=0}^{\infty }\left|{\rho }_{\left(k\right)}\right|$ diverges, where ${\rho }_{\left(k\right)}$ is the ACF of the process. Otherwise, the time series is said to be short memory or SRD. Another definition of long memory if for some frequency, $f\in [0,0.5]$, the power spectrum $P\left(f\right)$, becomes unbounded.

The power spectrum $P\left(f\right)$ is defined by:

$$P\left(f\right)={\int }_{-\infty }^{\infty }{e}^{-2\pi ifk}\rho \left(k\right)dk,$$

(3)

where $-\infty <f<\infty$, $i=\sqrt{-1}$ and $\rho \left(k\right)$ is the ACF.

The spectral density $S\left(f\right)$ is a normalized form of $P\left(f\right)$, defined by:

$$S\left(f\right)=\frac{P\left(f\right)}{{\sigma }^2}={\int }_{-\infty }^{\infty }{e}^{-2\pi ifk}\rho \left(k\right)dk.$$

(4)

If the spectrum becomes unbounded, then the ACF are not absolutely summable [19]. Therefore, ACF is defined as time domain analysis, while power spectrum density (PSD) is used for the frequency domain analysis.

The ACF of the stationary SRD stochastic models, such as the ARMA processes and Markov processes, is absolutely summable, while the correlations function ${\rho }_k$ is not absolutely summable for the processes with long-range dependence [19]. Signals with long-range correlations, which are characterized by inverse power-law decaying autocorrelation function, occur ubiquitously in nature and many man-made systems. Because of the strong coupling and the slow decaying autocorrelation, these processes are also said to be long memory processes. Typical examples of LRD signals include financial time series, underwater noise, electroencephalography (EEG) signal, etc. The level of the dependence or coupling of LRD processes can be indicated or measured by the estimated Hurst parameter, or the Hurst exponent [16]. The value of the Hurst Exponent varies between 0 and 1. If H = 0.5, the time series has no statistical dependence. If H < 0.5, the time series is a negatively correlated process or an anti-persistent process. If H > 0.5, the time series is a positively correlated process [20]. The LRD processes are also closely related to fractional calculus. In order to capture the property of coupling or hyperbolic decaying autocorrelation, fractional calculus based LRD models have been suggested, such as ARFIMA and FIGARCH models [21,22]. The ARFIMA model is a generalization of ARMA model, which is a typical fractional order system.

2.1. Autoregressive (AR) Model

The notation AR(p) refers to the autoregressive model of order p. The AR(p) model is written as [2]:

$$X_t=c+{\displaystyle \sum _{i=1}^p}{\varphi }_i{X}_{t-i}+{\epsilon }_t,$$

(5)

where ${\varphi }_1$, ⋯, ${\varphi }_p$ are autoregressive parameters, c is a constant, and the random variable ${\epsilon }_t$ is the white noise. Some constraints are necessary on the values of the parameters so that the model remains stationary. For example, processes in the AR(1) model with $|{\varphi }_1|\ge 1$ are not stationary. In statistics and signal processing, an autoregressive (AR) model is a representation of a type of random process; as such, it describes certain time-varying processes in nature, economics, etc.

2.2. Moving Average (MA) Model

The notation MA(q) refers to the moving average model of order q [2]:

$$X_t=\mu +{\displaystyle \sum _{i=1}^q}{\theta }_i{\epsilon }_{t-i}+{\epsilon }_t,$$

(6)

where the ${\theta }_1$, ⋯, ${\theta }_q$ are the moving average parameters of the model, $\mu$ is the expectation of $X_t$ (often assumed to equal 0), and the ${ϵ}_t$, ${ϵ}_{t-1}$,… are again, white noise error terms. The moving average (MA) smooths a time series, which can produce cyclic and a trend like plots even when the original data are themselves independent random events with fixed mean. This characteristic lessens its usefulness as a control mechanism.

2.3. ARIMA and ARFIMA Model

The above AR and MA models can be generalized as follows [2]:

$$(1-{\displaystyle \sum _{i=1}^p}{\varphi }_i{B}^i){(1-B)}^d(X_t-\mu )=(1+{\displaystyle \sum _{i=1}^q}{\theta }_i{B}^i){ϵ}_t.$$

(7)

The above ${(1-B)}^d$ is called a difference operator ${\mathsf{\nabla }}^d$. The ARMA or ARIMA models can only capture the SRD property, since d is confined in the range of integer order. Therefore, in order to capture the LRD property of the fractional systems, the ARFIMA(p,d,q) model is thereby proposed accordingly. In fact, the operator can be defined in a natural way by using binomial expansion for any real number d with Gamma function:

$${(1-B)}^d={\displaystyle \sum _{k=0}^{\infty }}\left(\genfrac{}{}{0pt}{}{d}{k}\right){(-B)}^k={\displaystyle \sum _{k=0}^{\infty }}\frac{\mathsf{\Gamma }(d+1)}{\mathsf{\Gamma }(k+1)\mathsf{\Gamma }(d+1-k)}{(-B)}^k.$$

(8)

Many authors suggested the use of the fractionally ARIMA model by using a fractional difference operator rather than an integer one could better take into account this phenomenon of LRD [23]. Hosking et al. defined an extension of the ARIMA model, which allows for the possibility of stationary long-memory models [24]. Thus, the general form of ARIMA$(p,q,d)$ process $X_t$ in Equation (7)—the ARFIMA(p,d,q) process is defined as:

$$\mathsf{\Phi }\left(B\right){(1-B)}^d{X}_t=\mathsf{\Theta }\left(B\right){\epsilon }_t,$$

(9)

where $d\in (-0.5,0.5)$, and ${(1-B)}^d$ is defined as the fractional difference operator in Equation (8). ARFIMA$(p,d,q)$ processes are widely used in modeling LRD time series, where p is the autoregressive order, q is the moving average order and d is the level of differencing [25]. The larger the value of d, the more closely it approximates a simple integrated series, and it may approximate a general integrated series better than a mixed fractional difference and ARMA model.

Figure 1 presents the discrete ARFIMA process that can be described as the output of the fractional-order system driven by a discrete white Gaussian noise (wGn). The ARFIMA$(p,d,q)$ process is the natural generalization of the standard ARIMA or ARMA processes. In a fractionally differenced model, the difference coefficient d is a parameter to be estimated first [26]. The intensity of self-similar of ARFIMA is measured by a parameter d. For the finite variance process with fractional Gaussian noise, d has a closed relation with Hurst parameter H [3,26,27]:

$$d=H-1/2.$$

(10)

In addition, for the infinite variance process with fractional $\alpha$-stable noise, d is related with Hurst and characteristic exponent $\alpha$ [18,22]:

$$d=H-1/\alpha .$$

(11)

In this way, the parameter d may be chosen to model long-time effects, whereas p and q may be selected to model relatively short-time effects.

3. Review and Evaluation

ARFIMA$(p,d,q)$ processes are widely used in modeling LRD time series, especially for the high frequency trading data, network traffic and hydrology dataset, etc. In practice, several time series exhibit LRD in their observations, leading to the development of a number of estimation and prediction methodologies to account for the slowly decaying autocorrelations. The ARFIMA process is one of the best-known classes of long-memory models. As introduced in Section 1, most statistical analysis programs have been embedded with ARFIMA models. A summary of the current software dealing with ARFIMA model analysis is as follows:

• MATLAB applications
  MATLAB^® (Matrix Laboratory) is a multi-paradigm numerical computing environment and fourth-generation programming language developed by MathWorks (Natick, MA 01760-2098, USA). The MATLAB applications are interactive applications written to perform technical computing tasks with the MATLAB scripting language from MATLAB File Exchange, through additional MATLAB products, and by users.
• SAS software
  SAS (Statistical Analysis System) is a software suite developed by SAS Institute (Cary, NC 27513-2414, USA) for advanced analytics, multivariate analyses, business intelligence, data management, and predictive analytics.
• R packages
  R packages and projects are contributed by RStudio (Boston, MA 02210, USA) team on CRAN (Comprehensive R Archive Network). R users are doing some of the most innovative and important work in science, education, and industry. It is a daily inspiration and challenge to keep up with the community and all it is accomplishing.
• OxMetrics
  Ox^TM is an object-oriented matrix language with a comprehensive mathematical and statistical function library developed by Timberlake Consultants Limited (Richmond, Surrey TW9 3GA, UK). Many packages were written for Ox including software mainly for econometric modelling. The Ox packages for time series analysis and forecasting.

MATLAB codes are open-source applications where we could download, view and revise the codes if possible while other three are packaged and embedded in the software modules. In the following evaluation parts, we could clearly see the differences between them even with the same inputs. Four primary embedded functions concerning simulation, fractional difference filter, parameter estimation and forecast, are tested and evaluated for the ARFIMA processes in

Table 1. Numerical tools for the ARFIMA process.

Procedures	MATLAB	R	SAS	OxMetrics
Simulation	✔${}^{*}$	✔	✔	✗
Fractional Difference	✔	✔	✔	✔
Parameter Estimation	✔	✔	✔	✔
Forecast	✔	✔	✗	✔

✔${}^{*}$ means it can simulate ARFIMA processes, but cannot choose or define the initial seeds.

. It should be noted that the first two functions can be regarded as the forward problem solving systems, while the latter two are developed for the backward problem solving systems which are much more significant. In view of the above, this section can be divided into four parts.

3.1. Simulation

On the website of MATLAB Central, there are two files that can simulate ARFIMA processes. They are developed by Fatichi [28] and Caballero [29]. However, users cannot choose initial random seeds, that is, it can only simulate one certain series of ARFIMA process. The ARFIMA$(p,d,q)$ estimator is developed by Inzelt, which is used for a linear stationary ARFIMA$(p,d,q)$ process [30].

R is a freely available language and environment for statistical computing and graphics, which provides a wide variety of statistical and graphical techniques: linear and nonlinear modelling, statistical tests, time series analysis, classification, clustering, etc. Like the MATLAB Central, CRAN is a platform that stores identical, up-to-date, versions of code and documentation for R. There are several major packages concerning ARFIMA process according to the authors’ survey in

Table 2. Comparison of ARFIMA packages in R.

Package	Author	Release Date	Typical Functions	Requirements
fractal	William Constantine et al. [31]	2016-05-21	hurstSpec	R (≥ 3.0.2)
fracdiff	Martin Maechler et al. [32]	2012-12-02	fracdiff	longmemo, urca
afmtools	Javier E. Contreras-Reyes et al. [33]	2012-12-28	arfima.whittle	R (≥ 2.6.0), polynom fracdiff, hypergeo, sandwich, longmemo
ArfimaMLM	Patrick Kraft et al. [34]	2015-01-21	arfimaMLM	R (≥ 3.0.0), fractal
arfima	Justin Q. Veenstra et al. [35]	2015-12-31	arfima	R (≥ 2.14.0), ltsa

.

The first two packages are used for the processing of ARFIMA processes, including Hurst fitting, calculation and fractional order difference and so on, while the latter two are mainly used for the parameter estimation of ARFIMA. The last package arfima is the most comprehensive tool that could simulate, estimate and predict the results of ARFIMA processes. In the paper, we use the last one package to compare with the other software.

SAS and R could also generate the ARFIMA process by defining the order of AR(p) and MA(q), setting the parameters $\varphi$, $\theta$ and d, respectively. In addition, the number of the initial random seeds could/should be set for the stochastic process. Random seeds are defined by the internal algorithms, which make the initial stochastic process a difference. Therefore, it may be a big difference if picking arbitrary seeds. In order to illustrate the above problems, we have generated the ARFIMA$(1,0.4,1)$ process with $d=0.4$, $\varphi =0.5$, $\theta =-0.1$ and $\sigma =1$. Then, we set 100 different initial random seeds with 3000 observations and do the same estimation. It should be kept in mind that, even with the same simulation software that generates the processes, the estimation results could be a big difference in Figure 2. However, from the perspective of the sample-path analysis for the stochastic processes, this could be the advantage compared to the MATLAB ARFIMA applications, which can only generate one certain series (path). Furthermore, we have also found that the SAS software is somewhat better or “conservative”, while R software is more “aggressive” in Figure 2. We could check the comparisons below with dashed lines showing true values of parameters.

It should be also noted that, even with a certain series of the initial starting random seeds, the estimation results could also have quite a few variations. For example, we have set the fractional order d from 0 to 0.5, and do the simulation and estimation accordingly in MATLAB. It can be seen that the estimation $\widehat{d}$ is jumping up and down around the true values (red line) in Figure 3.

Here are some comments of this subsection:

• Estimation results also depend on the initial random seeds, even the series that are from their own simulations.
• The test results may be different if not enough points/observations are generated. More than 300 points are preferred.
• Estimation results may not be accurate if they only use one method. R should be more desirable to try first.

3.2. Fractional Order Difference Filter

Many time series signals contain trends, i.e., they are non-stationary. It is usually preferable to specify and remove the trends explicitly to get the smoothed or stationary data for the further analysis and modeling. According to the theory of Box–Jekins, an ARIMA model can be viewed as a “filter” that tries to separate the signal from the noise, and the signal is then extrapolated into the future to obtain forecasts [2]. Since the beginning of the 1980s, the long memory ARFIMA model has been introduced and investigated by many scholars especially for the parameter estimation problems. Shumway and Stoffer gave a brief overview of “long memory ARMA” models in [36]. This type of long memory model might be considered to use when the ACF of the series tapers slowly to 0 and spectral densities are unbounded at $f=0$. Jensen et al. derived an algorithm for the calculation of fractional differences based on circular convolutions method in [37]. In fact, there are a lot of estimation methods concerning fractional difference algorithms.

In some instances, however, we may see a persistent pattern of non-zero correlations that begins with a first lag correlation that is not close to 1. In these cases, models that incorporate “fractional differencing” may be useful. Therefore, differencing the time series data by using the approximated binomial expression of the long-memory filter is a prerequisite to estimates of the memory parameter in the ARFIMA$(p,d,q)$ model. The user should not only set numeric vector of p and q, but also specify the order of the fractional difference filter. By passing through fractional order difference filter, the ARFIMA series will yield residuals that are uncorrelated and normally distributed with constant variance in Figure 4.

The sample ACF and partial autocorrelation function (PACF) are useful qualitative tools to assess the presence of autocorrelation at individual lags. The Ljung-Box Q-test is a more quantitative way to test for autocorrelation at multiple lags jointly [38]. The Ljung-Box test statistic is given by:

$$Q\left(L\right)=N(N+2){\displaystyle \sum _{k=1}^L}\frac{{\widehat{\rho }}_k^2}{N-k},$$

(12)

where N is the sample size, L is the number of autocorrelation lags, and $\rho \left(k\right)$ is the sample autocorrelation at lag k. Under the null hypothesis, the asymptotic distribution of Q is chi-square with L degrees of freedom. If we use lbqtest function in the MATLAB Econometrics Toolbox, it returns the rejection decision and p-value for the hypothesis test. Similar functions Box.test of stats and ljung.wge of tswge are also available in the R package. p-values indicate the strength at which the test rejects the null hypothesis. If all of the p-values are larger than 0.01, there is strong evidence to accept the hypothesis that the residuals are not autocorrelated.

Thus, we have generated an ARFIMA(0,0.4,0) process in Figure 5 and use fractional order difference filter with the order $d=0.4$ to filter the LRD property in Figure 6. It is obvious that, by passing through the fractional order difference filter, the slowly decaying property of LRD has been eliminated.

In order to evaluate residuals, p-values are used to quantify the goodness of fitting in

Table 3. Fractional filters.

Software	MATLAB	R	SAS	Ox
Function	d_filter	diffseries	fdif	fracdiff
p-values with 1,5,10,15 lags	0.0710	0.09998	0.1062	0.0862
	0.2253	0.2395	0.2414	0.2114
	0.5850	0.5320	0.5198	0.5898
	0.5330	0.4571	0.4473	0.5473

.

Here are some comments of this subsection: all of the four programs above have fractional order operators to filter the LRD process successfully. In general, d is the parameter to be estimated first [26]. If we use the calculation defined by the Hurst method in Equations (10) and (11), d could probably be the fractional one. Therefore, the fractional order filter would be the primary tool to eliminate the LRD property or the heavy-tailedness in order to get the stationary series.

Meanwhile, however, the fractional order d is closely related to the Hurst parameter in Equations (10) and (11). There are more than ten methods to estimate Hurst parameters, R/S method, aggregated variance method, absolute value method, periodogram method, whittle method, Higuchi’s method, etc. These methods are mainly useful as simple diagnostic tools for LRD time series. These Hurst estimators have been introduced to analyze the LRD time series in [17,39,40]. Therefore, the results of Hurst estimators can be different if applying different methods. In addition, from Equation (8), it is interesting to note that there are infinite factorial series in the expansion of binomial expansion. In practice, we usually take the first three factorials for approximation. That is to say, the accuracy of differencing is also determined by how many factorials are used for approximation. Consequently, these different methods make the subsequent estimations differ from each other in the following sections.

3.3. Parameter Estimation

From Figure 2, we could see that, even though R and SAS can both simulate the ARFIMA processes, the properties of these processes are not the same mainly because of the distinctive random seeds defined by different software. Therefore, when reviewing and evaluating the accuracy of above software, the same ARFIMA series should be guaranteed first. Herein, we have proposed to use the following steps to compare the results in Figure 7. In addition, OxMetrics is the software that cannot generate ARFIMA simulation, but it can estimate and forecast ARFIMA-FIGARCH processes. MATLAB cannot generate multiple ARFIMA series for the same parameter combinations. We have thus used R and SAS to provide the ARFIMA series for the inputs of estimations.

Since we have received simulation results, the parameters of ARFIMA processes can be estimated and compared with true values (parameter setting values). First, we have used the simulation data from R software and have then used these three programs to do the estimation in Figure 8. Second, we have used SAS to do the same simulation and have then used the other three to do the estimation in Figure 9. Without loss of generality, we pick 10 groups of 3000 observations to see who could capture the accuracy.

Here are some comments from this subsection: from the above plots, it is very interesting to find that the estimation results of ARFIMA simulations are relatively accurate when they come from the same simulation data set. However, OxMetrics and MATLAB estimate the negative values of $\theta$.

In order to further test whether the Ox and MATLAB can only return the negative values of $\theta$ or if they just return the inverse values. The parameters are set to the inverse values, accordingly, with $\varphi =-0.5$ and $\theta =0.2$; while, in the previous test, they are $\varphi =0.5$ and $\theta =-0.1$, respectively. The result presented in Figure 10 validates the comments above.

3.4. Forecast

Simulation data could only be the auxiliary part of these software, since it can never be a powerful and useful tool for the ARFIMA process analysis if it cannot retrieve the estimation parameters from the real data with LRD. Moreover, the last and the most significant part of the ARFIMA process is to forecast and thereby predict the future behavior. Therefore, mean absolute percentage error (MAPE) values are used for the evaluation of the forecast results for the data from real life. The error square of the prediction results from different methods with the increasing number of predictions are illustrated in Figure 11:

$${\sigma }_t^2={\displaystyle \sum _{t=1}^n}{(\widehat{y_t}-y_t)}^2.$$

(13)

Data description: Centered annual pinus longaeva tree ring width measurements at Mammoth Creek, Utah, USA from 0 A.D (anno Domini) to 1989 A.D with 1990 sampling points in time series [41,42]. The data can be divided into two parts: the first part with 1900 observations are used to estimate ARFIMA parameters, and the second part with 90 points are used for comparison with the prediction results from the fitted ARFIMA models. Finally, the results with the implemented methods that are applied to real-life time series are summarized in

Table 4. Parameter estimations and forecast comparisons.

Number	Parameters	R	SAS	OxMetrices	MATLAB
1	mu	0.9833	N/A (Not Applicable)	0.98799	0.9878
2	d	0.1670	0.1479624	0.282087	0.2313
3	ar	0.9070119	0.8939677	−0.254265	0.6473
4	ma	0.8603811	0.8318787	0.18698	0.6393
5	sigma	0.1078173	0.1073417	0.1246	0.1163
6	p value Lag1	0.9195	N/A	0.7709458	0.9101
7	p value Lag5	0.6369	N/A	0.341324	0.6959
8	p value Lag10	0.8659	N/A	0.4367925	0.9037
9	p value Lag15	0.6491	N/A	0.6229542	0.6776
10	LogLikelihood	2117.224	1851.5512	−570.599	1162.527
11	MAPE	28.95	N/A	29.36	29.02

.

Here are some comments of this subsection:

• d is the parameter to be estimated first when doing ARFIMA model fitting. Therefore, if the estimation of d is different for a certain time series, the following estimations for AR($\mathsf{\Phi }$) and MA($\mathsf{\Theta }$) will be different.
• The ideal length (horizon) of predictions is within 30 steps. With the increasing steps of forecast, prediction errors are adding up. If a long range prediction series is required, R and MATLAB should be priorities for their smaller prediction errors.
• Compared with other forecast results with true values in Table 4, R produces the minimum prediction errors and MAPE.

4. Summary of Selection Guidelines

Qualitative analysis as well as quantitative evaluations of the selected ARFIMA tools have been conducted in the previous sections. In order to make it easier for researchers from different backgrounds, we summarize the selection guidelines for the ARFIMA process modeling and analysis.

• R and SAS software are priorities for the simulation of ARFIMA process, since they could define the initial seeds. R is one of the desirable tools for the estimation of ARFIMA process, since it has more than five packages including Hurst estimators, ACF plot, Quantile-Quantile (QQ) plot, white noise test and some LRD examples.
• Estimation results of the ARFIMA process may be different if the number of observations is not large enough. Therefore, more than one estimation method should be used in order to guarantee the accuracy.
• d is the parameter to be estimated first. All of this software could use fractional difference functions to filter the trend and thereafter stationarize time series data.
• The ideal length (horizon) of predictions is within 30 steps. If a long range prediction series is required, R and MATLAB are the priorities for their smaller prediction errors.

5. Conclusions

Compared to the conventional integer order models that can only capture SRD, the ARFIMA model gives a better fitting result, especially for the data with the LRD property. Nowadays, some programs have been integrated with ARFIMA solutions. However, the final results of estimation and prediction could be different or even conflicting if choosing different methods. Therefore, a comprehensive review and evaluation of the numerical tools for the ARFIMA process is presented in the paper so as to provide some guidelines when choosing appropriate methods to do the time series analysis of LRD data. Through such efforts, it is hoped that an informative guidance is provided to the readers when they face the problem of selecting a numerical tool for a specific application.