Algorithms for Linear Time Series Analysis: With R Package

Our ltsa package implements the Durbin-Levinson and Trench algorithms and provides a general approach to the problems of ﬁtting, forecasting and simulating linear time series models as well as ﬁtting regression models with linear time series errors. For computational eﬃciency both algorithms are implemented in C and interfaced to R . Examples are given which illustrate the eﬃciency and accuracy of the algorithms. We provide a second package FGN which illustrates the use of the ltsa package with fractional Gaussian noise (FGN). It is hoped that the ltsa will provide a base for further time series software.


Introduction
Let z t , t = 1, . . ., n, denote n successive observations from an ergodic covariance stationary Gaussian time series with autocovariance function (acvf) γ k = Cov (z t , z t−k ), k = 0, 1, . . ., n− 1 and mean µ.The general linear process (GLP) may be specified by its autocovariance sequence, γ k , or equivalently in terms of its autocorrelation sequence (acf), ρ k = γ k /γ 0 or coefficients ψ k in the infinite moving-average (MA), where a t ∼ NID (0, σ 2 a ) and ψ 2 1 + ψ 2 2 + . . .< ∞.Notice that both acf and infinite MA model specifications also require the innovation variance σ 2 a .The condition ψ 2 1 + ψ 2 2 + . . .< ∞ ensures that the acvf exists and that the GLP is stationary.For sufficiently large Q we may approximate using a MA of order Q, Most parametric time series models may be specified so that either the autocovariances, γ k , or the MA coefficients, ψ k , are functions of a small number of parameters, β.Both theory and experience suggests that exact maximum likelihood are preferable to other methods.Please see the Appendix for further discussion about the GLP.
where φ (t) = (φ t−1,1 , . . ., φ t−1,t−1 ) is determined by the linear equations Γ t φ (t) = (γ(1), . . ., γ(t)) , and the variance of the predictor is given by The covariance determinant is, The Durbin-Levinson algorithm (Golub and Loan 1996, Algorithm 5.7-1) provides an efficient algorithm for solving Equations ( 4) and (5) in O(n 2 ) flops.The Trench algorithm (Golub and Loan 1996;Trench 1964) may be derived by determining the parameters in an AR (n−1) with autocovariances γ 0 , . . ., γ n−1 using the Durbin-Levinson algorithm.After the AR coefficients have been found, Γ −1 n is obtained using the method given by Siddiqui (1958).The Trench algorithm evaluates the matrix inverse in O(n 2 ) flops.Our R function TrenchInverse uses the C interface to provide an efficient implementation of the Trench algorithm.On a typical current PC it takes less than a fraction of a second to evaluate this inverse for matrices with n = 1000.In the special case of ARMA time series, Zinde-Walsh (1988) obtained an explicit expression for computing the elements in Γ −1 n .This method is suitable for symbolic computation (McLeod 2005).
The Durbin-Levinson and Trench algorithms provide a convenient method for computing the exact likelihood function of a general linear Gaussian time series model (Li 1981;Sowell 1992;Brockwell and Davis 1991).In a suitable quantitative programming environment such as R the likelihood function may be explored graphically, optimized to obtain exact MLE or integrated for Bayesian estimates (McLeod and Quenneville 2001).For the well-known ARMA family of time series models there are many algorithms for the computation of the exact likelihood function which require only O(n) flops per function evaluation Box and Luceño (1997, Section 12B) as opposed to the O(n 2 ) flops required in the Durbin-Levinson or Trench algorithm.For long ARMA time series, a superfast likelihood algorithm requiring only O(1) flops per log-likelihood evaluation is available (McLeod and Zhang 2007).The advantage of the Durbin-Levinson or Trench algorithm is that they are more general and with current computing technology, MLE using this algorithm is sufficiently fast provided that n is not too large.When the ARMA model is extended to include long-memory alternatives such as the ARFIMA model Brockwell and Davis (1991, Section 13.2) and other ARMA long-memory extensions (Baillie 1996), the Durbin-Levinson likelihood algorithm is as computationally efficient as other commonly used exact likelihood methods such as the innovation algorithm or Kalman filter.A brief discussion and comparison with the superfast algorithm (Chen, Hurvich, and Lu 2006) is given in Section 2.3.
For linear parametric time series models, it is assumed that the acvf or MA coefficients, ψ k , k = 0, 1, . . ., are uniquely determined by p parameters, β = (β 1 , . . ., β p ).An example we will discuss in Section 3 is the fractional Gaussian noise time series model which is characterized by its autocorrelation function Hipel and McLeod (1994, page 340), In addition to fitting and forecasting linear time series models, there are many other applications for our algorithms in time series analysis including regression with autocorrelated error, an example of which is discussed in Section 3.4.Another example where the inverse matrix is needed is for power computations in intervention analysis (McLeod and Vingilis 2005).

Package overview
The R functions in the ltsa package are shown in Table 1.The TrenchInverse function is especially useful for efficient exact MLE for regression with autocorrelated error and for the mean, µ, as well for computing and updating forecasts.The log-likelihood may be computed using either the Durbin-Levinson or Trench algorithm.Residuals are useful for checking model adequacy and these may also be computed using the Durbin-Levinson algorithm.A linear time series may be simulated using the Durbin-Levinson recursion, the Davies-Harte algorithm (Davies and Harte 1987) or another method using the fast Fourier transform that is given in Section 2.6.In the next section we discuss in more detail TrenchInverse and in the following sections the remaining functions are discussed.The FGN package discussed in Section 3 describes how the ltsa functions may be used to develop a package for fractional Gaussian noise (FGN) time series modeling.

TrenchInverse
The TrenchInverse function in R is interfaced to a C function for maximum speed.R memory management C functions are used so there is complete compatibility with all R platforms.Our package has been tested in the Windows, Mac OS X and Debian Linux environments.
The function TrenchInverse in this package inverts a positive-definite Toeplitz matrix utilizing the interface provided in R to call a C function.If the matrix is not a positive definite, a suitable error message is returned.
The TrenchInverse function takes a single matrix argument Γ n .The built-in R function, toeplitz, may be used to create this matrix Γ n from the vector input (γ 0 , . . ., γ n−1 ).For maximum computational efficiency one could work with just this vector, and this is how in fact the underlying C algorithm is implemented.As a brief illustrative example of this package, in the code below, we subtract the product of a Toeplitz matrix and its inverse from the identity matrix and compute the largest absolute error in the result: R> phi <-0.8 R> n <-1000 R> r <-(1 / (1 -phi^2)) * phi^(0:(n-1)) R> G <-toeplitz(r) R> Gi <-TrenchInverse(G) R> id <-matrix(0, nrow=n, ncol=n) R> diag(id) < -1 R> max(abs(id -G%*%Gi)) [1] 6.661338e-16 We investigated the timings for the function TrenchInverse and compared them with the general purpose matrix inverse function solve in R. The timings reported in this section were done on a 3.6 GHz Pentium 4 PC running Windows XP.For these timings we first generated 25 uniform random values on the interval (−1, 1), denoted by φ k , k = 1, . . ., 25.These values were used to generate 25 Toeplitz matrices of the form, Γ n = (φ , for each n = 400, 800, 1200, 1600, 2000.The code for generating this table is included in the package documentation for TrenchInverse.We conclude from Table 2 that TrenchInverse is significantly faster than solve. In forecasting applications that we will discuss below in Section 2.5, the successive inverses of Γ n+k , k = 1, 2, . .., are needed.Then Γ −1 n+k may be more efficiently computed using the result for the inverse of partitioned matrices Graybill (1983, Section 8.3).When k = 1, we obtain, where a = ehh , e = 1/(γ 0 − h Γ −1 n h), h = (γ 1 , . . ., γ n ) , f = −eΓ −1 n h and I n is the n × n identity matrix.Then Equation ( 8) may be applied repeatedly to obtain Γ −1 n+k , k = 1, 2, . ... In Table 3, we compare the timing and accuracy of the updating approach described in Equation (8), implemented in ToeplitzInverseUpdate, with direct inversion in TrenchInverse.The function TrenchInverse is implemented using R interface to C, whereas ToeplitzInverseUpdate is implemented entirely in R. For the comparison we used the hyperbolic decay autocovariance function γ k = 1/(k + 1), k ≥ 0, and computed Γ −1 n+k , k = 1, 2, . . ., 100 and n = 100, 500, 1000, 2000.As a check, the absolute error for the difference between the two matrices was computed and it was found to be negligible.The CPU times are shown in Table 3.The updating algorithm is about 3 to 4 times as fast and this factor does not seem to change much with n.This might be expected since both algorithms require O(n 2 ) flops.

Log-likelihoood
In this section we discuss the functions DLLoglikelihood and TrenchLoglikelihood.
Assuming the mean is known and that z t has been mean-corrected, the log-likelihood function for parameters (β, σ 2 a ) given data z = (z 1 , . . ., z n ) may be written, after dropping the constant term, Letting M n = Γ n /σ 2 a and maximizing L over σ 2 a , the concentrated log-likelihood may be CPU time in seconds for direct computation using TrenchInverse and ToeplitzInverseUpdate.First the inverse of a matrix of order n is computed and then the inverses for matrices of orders n + 1, . . ., n + 100 are computed using TrenchInverse and ToeplitzInverseUpdate. written, where S(β) = z M −1 n z and g n = det(M n ).Note that L c is unchanged if we simply replace M n in S(β) and g n by the autocorrelation matrix R n = (ρ |i−j| ) n×n .Using the Durbin-Levinson algorithm, where ẑt and σ 2 k are given in Equation ( 3) and ( 5).The C functions we developed can evaluate S(β) as well as log(g n ) using the Durbin-Levinson method as well as the more direct Trench algorithm.These C functions are interfaced to DLLoglikelihood and TrenchLoglikelihood.Either can then be used with the optimization functions provided with R to obtain the MLE.In practice, the Durbin-Levinson is somewhat faster as can be seen in Table 4. Chen et al. (2006) have implemented a superfast method of solving linear Toeplitz systems of the form Γ n x = b, where Γ n is an order n symmetric positive-definite Toeplitz matrix, b is a known vector of length n and x is the vector of unknowns.With their method, x can be found in O(n log 5/2 n) flops.This provides an alternative and, in principle, computationally more efficient method of evaluating the loglikelhood.However, unless n is very large, the gain may not be significant.The time will also depend on the computer and the specific implementation.Timings they gave for their superfast algorithm, ML-PCG (S-PLUS), are reported in Table 4 and Table 5 along with timings for their S-PLUS version of the Durbin-Levinson method, ML-Levinson (S-PLUS).Apart from coding, their ML-Levinson (S-PLUS) is equivalent to our DLLoglikelihood, useC=TRUE.Our timings for the Durbin-Levison algorithm are much faster and this is no doubt due to language/machine differences.In Table 5 we used a Windows PC with a 3.6 Ghz Pentium processor and in Table 4 we used a newer PC running Debian Linnux since our Windows PC could not handle such large matrices.(Chen et al. 2006, Table 4) used a Sun Workstation running the Solaris OS.Tables 4 and 5 suggest that, for many purposes, the current implementation of our algorithms has a satisfactory performance with respect to computer time.
The superfast algorithm of Chen et al. (2006)  inverse.A further advantage of the Trench algorithm is that as a byproduct the exact value of the determinant is also obtained.

Estimation of the mean
Given the other parameters in the model, so that γ k , k = 0, . . ., n − 1 is specified, the best linear unbiased estiamte estimate (BLUE) for µ is given by (Beran 1994, Section 8.2), μ = where z = (z 1 , . . ., z n ) and 1 n is an n dimensional column vector whose entries are all equal to one.Notice that the autocovariances in Equation ( 12) may be replaced by autocorrelations.The function TrenchMean implements the computation in Equation ( 12).
An iterative algorithm may be used for the simultaneous joint MLE of µ and the other parameters β.
Step 3 Terminate when i+1 has converged or i > M .Otherwise set i ← i + 1 and return to Step 1 to perform the next iteration.
Convergence usually occurs in two or three iterations.
Another topic of interest is the question of the efficiency of the sample mean.The variance of the sample mean, z = (z 1 + . . .+ z n )/n, may be written, It may be shown that the exact finite sample efficiency of the sample mean is (Beran 1994, Section 8.2) Although the sample mean often provides an efficient estimate, situations exist it is not very efficient.We will discuss an example of this in Section 3.2.

Forecasting
The Trench algorithm is useful for the computation of exact finite-sample forecasts and their variances.Let z n (k) denote the minimum-mean-square-error linear predictor of z n+k , given the data z = (z 1 , . . ., z n ), the mean µ and autocovariances γ , = 0, . . ., n−1.Then Hamilton (1994, Section 4.3) or Hipel and McLeod (1994, Section 10.4.5), where g k = (γ n+k−1 , . . ., γ k ) and the variance for the forecast, For large n, V 1 .= σ 2 a .In Equation ( 15), the autocovariances may be replaced by autocorrelations, but in Equation ( 16) autocovariances must be used.Equations ( 15) and ( 16) may be vectorized to compute the k-step predictor for k = 1, 2, . . ., L, and this is implemented in TrenchForecast.The computation of the forecasts and their variances using Equations ( 15) and ( 16) requires the autocovariances γ 0 , . . ., γ n+k−1 .The prediction variance may also be computed in another way as described in Section 2.9, and this is implemented in the function PredictionVariance.
In practice we may also be interested in updating the forecasts given a new data values z n+k , k = 1, 2, . . ., L. This entails computing Γ −1 n+k for k = 1, 2, . . ., L given Γ −1 n and may be computed using the updating method described in Section 2.2.The function TrenchForecast implements the updating method as well as an option to compute the forecasts using direct matrix inversion.This is done in order to provide an optional check on the computations.
An example application of TrenchForecast to actual time series data is given in Section 3.3.
In the example below, we simulate a series of length n = 202 from an AR(1), then using the known parameter, we forecast starting at origin n = 200 for lead times k = 1, 2, 3, and updating the forecast origin to n = 201, 202.The result is a list with two 3 × 3 matrices.It is easy to check the forecasts and standard deviations in this case, R> n <-200 R> m <-2 R> maxLead <-3 R> N <-n +m R> phi <-0.9 R> r <-phi^seq(0, N + maxLead -w1) R> set.seed(19890824)R> z <-DLSimulate(N, r) R> out <-TrenchForecast(z, r, 0.0, n, maxLead) R> out$Forecasts Finally, it should be noted that often the Gaussian assumption may not be valid for the forecast errors.In the non-Gaussian case the prediction variances should not be used to set probability limits for the forecasts.Instead probability limits for the forecasts may be obtained by simulation or bootstrapping (McCullough 1994;Aronsson, Holst, Lindoff, and Svensson 2006).

Simulation
Simulation of time series is widely used in bootstrapping for statistical inference as well in the exploration of statistical properties of time series methods.Time series simulation is also important in engineering design and operational research (Dembo 1991;Hipel and McLeod 1994;Maceira and Damázio 2006).The Durbin-Levinson algorithm provides a convenient method for simulating a Gaussian time series, z 1 , . . ., z n , with autocovariances, γ 0 , . . ., γ n−1 (Hosking 1981;Hipel and McLeod 1994).Using the Durbin-Levinson algorithm to solve Equations ( 4) and ( 5), the series z 1 , . . ., z n may be generated for t = 2, . . ., n from where e t ∼ NID (0, σ 2 t−1 ) and z 1 ∼ NID (0, σ 2 0 ).The function DLSimulate implements this simulation method using an interface to C. As a check the algorithm was also implemented in R and may be invoked using the optional argument useC = FALSE.Davies and Harte (1987) gave an algorithm which only requires O(n log(n)) flops as compared with Durbin-Levinson's O(n 2 ) flops.This algorithm is implemented in R, in our function DHSimulate.However, the Davies-Harte algorithm requires a complicated non-negativity condition, and this condition may not always hold.For example, in Table 6 we generate a time series of length n = 5000 with fractional difference d = 0.45, and we found the Davies-Harte non-negativity condition failed, and so DLSimulate was needed.Another method is useful for simulating time series with innovations from a specified non-Gaussian distribution that also uses the Fast Fourier Transform (FFT).In this case we may approximate the linear time series model as a high-order MA The order Q may be quite large in some cases, and it may be chosen so that mean-square error difference, is made sufficiently small.It may be shown that by making E small we can make the Kullback-Leibler discrepancy between the exact model and the MA(Q) approximation negligible McLeod and Zhang (2007, Equation 13) An example R script for determining the approximation is given in the online documentation for DLAcfToAR.The sum involved in Equation ( 18) is efficiently evaluated using the R function convolve which uses the FFT method.Hence the simulation requires O(n log(n)) flops when n is a power of 2 and assuming n > Q.The built-in R function arima.simmay also be used, and it uses direct evaluation and drops the initial transients values.The number of initial transient values to drop is determined by the optional argument n.start.Only O(n) flops are required by arima.sim.In Table 6, we took Q = 1000 and n.start=1000 for SimGLP and arima.simrespectively.These two approximate methods were compared with the exact simulation methods, DLSimulate and DHSimulate, for the case of an hyperbolic decay time series with γ k = 1/ √ k + 1, k ≥ 0, for time series of lengths 100, 200, 500, 1,000, 5,000 and 10,000.For each n, the total time for 100 simulations was found.
Joint MLE for α, β and σ 2 a may be obtained using the following iterative algorithm.
Step 0 Initialization.Set i ← 0. Set R n to the identity matrix.Set 0 to a negative number with large absolute value.
Step 2 Taking ξ as the input time series, maximize L m using a nonlinear optimization function to obtain β(i) and i = L m ( β(i) ).
Step 3 If i − i−1 < 10 −3 , perform Step 4. Otherwise set i ← i + 1, and return to Step 1 to perform the next iteration.
Step 4 Compute the MLE for σ2 a .
The error tolerance is set to 10 −3 , since in terms of the log-likelihood the change in parameters is of negligible importance.With this error tolerance, convergence usually occurs in three or four iterations.An implementation of this method for multiple linear regression with FGN error is given in our FGN package.

Prediction residuals
The prediction residuals are defined by the difference between the observed value and the one-step ahead minimum-mean-square error forecast.These residuals may be standardized by dividing by their standard deviations.If the model is correct, these residuals should be approximately uncorrelated.It should be noted that asymptotically the prediction residuals are equivalent to the usual residuals, the estimated innovations, ât .Hence, the widely used Ljung-Box portmanteau test (Ljung and Box 1978) and other diagnostic checks (Li 2004) may be used to check the adequacy of the fitted model.

Acf to AR parameters and ARMA Acvf
Given the autocorrelations ρ 1 , . . ., ρ m , the function DLAcfToAR uses the Durbin-Levinson recursions to obtain the parameters φ 1 , . . ., φ m of the best linear predictor of order m as well the partial autocorrelations φ 1,1 , . . ., φ m,m and the minimum mean-square-errors σ 2 1 , . . ., σ 2 m corresponding to the k-th order predictor, k = 1, . . ., m.For our purposes DLAcfToAR provides more general and useful output than the built-in R function Acf2AR.
As an illustrative application, consider the computation R 2 (Nelson 1976), For a fitted model we may use estimates of the parameters.For sufficiently large m, where σ2 m is computed with DLAcfToAR using the fitted values, ρ1 , . . ., ρm ,and setting, without loss of generality, γ 0 = 1 in Equation ( 21).R 2 indicates the proportion of variability accounted for by the model and is sometimes called the coefficient of forecastibility.
In the brief example below, we show that for a FGN model with H = 0.84, R 2 .= 41%.
R> library("FGN") R> m <-10^4 R> r < -FGNAcf(1:m, 0.84) R> 1 -DLAcfToAR(r)[,3][m] 10000 0.4075724 (Box, Jenkins, and Reinsel 1994, Chapter 7) gave algorithms for forecasting ARMA models which work well when n is not too small and the model is invertible whereas the method given in Section 2.5 is always exact although not as computationally efficient or elegant their methods.(Box et al. 1994, Chapter 7) showed, Given autocorrelations ρ 1 , . . ., ρ n , we may fit the AR(n) linear predictor using DLAcfToAR, and then expand as a MA using the built-in R function ARMAtoAR.Our function PredictionVariance uses Equation ( 23) to compute V k , k = 1, . . ., L, where L is the maximum lead time.As an option, PredictionVariance also implements the exact method of Section 2.5.
Another useful function is tacvfARMA for the tacvf of the ARMA(p, q) which is similar to the built-in ARMAacf but provides autocovariances (McLeod 1975) instead of autocorrelations.Again this function generalizes the built-in R function in a useful way.The function tacvfARMA is needed to demonstrate that the output from our TrenchForecast normally agrees quite closely with the built-in predict.Arima for ARMA models -see Example 1 in the TrenchForecast documentation.This function is also useful in computing the variance of the sample mean using Equation (13).

Validation and checking of the algorithms
Many checks were done to ensure the accuracy of our results.Most of the checks described below are given in the help file documentation associated with the particular R functions.
The function TrenchInverse is easily checked by simply multiplying the output by the input to obtain the identity matrix.The exact concentrated log-likelihood function defined in Equation ( 10) and implemented in TrenchLoglikelihood was checked by implementing in R the Durbin-Levinson recursion to compute (10).This function, DLLoglikelihood, with the option useC=FALSE provides a slower alternative method for the evaluation of ( 10).The computation is also easily verified in the case of the Gaussian AR(1) model, . Given data z 1 , . . ., z n , Equation ( 10) for the concentrated log-likelihood reduces to where In the documentation for DLLoglikelihood we show numerically that these two methods of evaluating the concentrated log-likelihood for the AR(1) are equivalent.Also, for the AR(1), the exact forecast function may be written The output from TrenchForecast agrees with the results produced by this formula for the AR(1) case.Details are given in the package documentation for TrenchForecast.The function TrenchMean can be checked by computing the exact MLE for an AR(1) fitted to some test data and then using the function TrenchMean to compute the exact MLE for µ given the estimate φ obtained from arima.This estimate of µ closely agrees with that given by arima.An illustration of this check is provided in the documentation to the function TrenchMean.

FGN package
The FGN package illustrates the use of the ltsa package.This package provides modelling capabilities for the FGN time series defined by the autocorrelation function given in Equation (7).The principal functions available in this package are shown in Table 7 The functions FGNLL, FitFGN, FitRegressionFGN, GetFitFGN, and SimulateFGN are specific implementations of the general methods discussed in Section 2.5 for the case of the FGN model.
The function HurstK provides a nonparametric estimator of H which is described in detail in Hipel and McLeod (1994, page 232).
The simulation function SimulateFGN utilizes DHSimulate or DLSimulate.It was determined empirically that the Davies-Harte non-negativity condition holds for n ≥ 50 and 0 < H < 0.84.So in this case DHSimulate is used, and DLSimulate is used otherwise.Table 8 shows average time taken to simulate and fit an FGN model for various n on our PC Windows XP 3.6 GHz Pentium IV computer.
The output from FitFGN is an S3-class object "FitFGN" with methods implemented for the standard R generic functions: coef, plot, predict, print and summary.

Efficiency of the sample mean
It is shown in Beran (1994, Section 8.2) that the asymptotic efficiency of the sample mean is always greater than 98% in the persistent case, that is, when 1 2 < H < 1. Table 9, obtained by evaluating Equation ( 14), shows the exact small sample efficiency for various lengths n.The finite sample efficiency is in good agreement with the asymptotic limit in the persistent case.Note that when H = 1 2 , the series is simply Gaussian white noise, and when 0 < H < 1 2 the series is said to be antipersistent.Such antipersistent time series can arise when differencing is used to transform a nonstationary time series to a stationary one.In the strongly antipersistent case with H = 0.1, the efficiency of the small mean can be quite low as seen in Table 9.Since most annual geophysical time series exhibit persistence, the sample mean may be used for such data.term dependence.(Beran 1994, page 118) obtained Ĥ = 0.84 using the Whittle approximate MLE.The MLE using FitFGN is Ĥ = 0.831 which agrees with the Whittle estimate as well as with the Hurst K nonparametric estimate which was K = 0.825.The diagnostic plots produced by plot are shown in Figure 1. Figure 2 shows the normal probability of the residuals.These diagnostic plots confirm that the FGN is a reasonable model.An ARMA(2, 1) was fit using the R function arima and was found to fit the observed series very well.The exact likelihood for the fitted ARMA model was determined using the DLLoglikleihood function.Table 10 summarizes the fits in terms of L m .We see that the ARMA(2, 1) is slightly better in terms of the log-likelihood but requires more parameters than FGN so that FGN is better in terms of both the AIC and BIC criteria.
In Table 11, we compare the forecasts at origin time n = 663 for lead times k = 1, . . ., 5 for the FGN and ARMA models.The standard deviations of the forecasts were also computed.The differences between the models seem minor.In the case of the ARMA model we compared the built-in R function predict.Arima with our TrenchForecast for the fitted ARMA model.As expected there is almost no difference in this case.
A further forecasting experiment compares the quality of the forecasts using TrenchForecast with the built-in R function predict for the ARMA(2, 1) model.In addition, the forecast for the fitted FGN model was included to compare with the ARMA(2, 1).In each case, the model was fit to all the time series up to time t = n 1 + k for each k, k = 1, . . ., K, where we took n 1 = n − K, K = 100, and n = 663 is the length of the series NileMin.For each t, t = n 1 +1, . . ., n, we compute the forecasts z t ( ), = 1, 2, 3, and their errors e t ( ) = z t+ −z t ( ).The empirical root-mean-square errors (RMSE), were computed for = 1, 2, 3 and are shown in Table 12.It is interesting that the FGN outperformed ARMA and that at = 2 the TrenchForecast was more accurate than predict.Arima.The R script to produce of the time series is used for fitting the model, and the second portion is used to evaluate the out-of-sample forecast performance (Noakes, McLeod, and Hipel 1985;Jaditz and Sayers 1998).In the case of ARIMA models the built-in function predict can not do this, since it is necessary to update the fitted model for the next forecast.A script was written to fit ARMA(p, q) and FGN time series models to one part of the series and to compute the RMSE for the second part.This script is available in the online documentation for the dataset NileMin.Selecting p = 2, q = 1, we fit the ARMA(p, q) as well as the FGN to all but the last m = 100 observations, and then compared the forecasts with the actual values for the remaining 100 values.The RMSE for the forecasts at lead times = 1, 2, 3 was computed.
The forecasts for both models were computed using TrenchForecast.As shown in Table 13, the FGN model slightly outperforms the ARMA(2, 1) model.Notice that because the model is not updated, the forecast errors are larger in Table 13 than in Table 12.
Lead FGN ARMA(2, 1) 1 0.568 0.579 2 0.659 0.678 3 0.686 0.706 Table 13: RMSEs for forecasts for the last 100 values in the NileMin dataset.In this case, the model was fit only once to the first 563 values and then its forecasting performance was evaluated on the next 100 values.

Simulation and bootstrapping experiments
The computation results reported in this section were carried out using the Rmpi package (Yu 2002) on a Beowulf cluster with 48 CPUS running Debian Linux.This reduced the time needed for our computations by a factor of about 30.The R scripts used are available and provide a template for bootstrapping and simulating with Rmpi.

Bootstrapped forecasting experiment
The FGN model was also fit to the entire NileMin series, and then three independent bootstrap versions of the series were generated using Boot.The forecasting experiment discussed in the last paragraph of Section 3.3 was then repeated on each of the bootstrap series but in addition to fitting the FGN and ARMA(2, 1) models, we also used FGN model with the true parameter value H = 0.831 in order to allow us to investigate the effect of parameter estimation errors on the FGN forecast.
Then the bootstrapping was iterated B = 10 4 times for all three models.The average RMSE is shown in Table 15.As expected when the parameter H is known, the forecast is then, without question, the optimal RMSE forecast and as expected it did better than the other two methods using estimated parameters.The FGN only slightly outperforms the ARMA(2, 1) model at lead-times 2 and 3.The whole experiment was repeated six times, and each time results comparable to Table 15 were produced.If desired, the jackknife could be used to estimate the standard deviation or the experiment could be re-run as a paired comparison, and a suitable confidence interval for the difference given but we feel confident enough about the overall general conclusion.The R script used for Table 15

Nile river intervention analysis
In 1903 the Aswan Dam went into operation and it is of hydrological interest to quantify the effect of this dam on downstream annual riverflow (Hipel, Lennox, Unny, and McLeod 1975 so the FGN model may be preferred on these grounds.Figure 3 shows the confidence regions for both models.The larger region corresponds to FGN errors.This is due to the stronger form of dependence present in the FGN case.In conclusion, we see that the uncertainty in the effect of the intervention is much greater with FGN errors.

Concluding remarks
The ltsa package is implemented in R (R Development Core Team 2007) and C and is available from the Comprehensive R Archive Network at http://CRAN.R-project.org/.The FGN, available as well from CRAN, illustrated how this package can be used to provide a building block for other time series software.It is intended to develop a package for another paper which implements the ARFIMA model and some of its generalizations.Another possible extension is to the parametric model suggested by Bloomfield (1973).In this case a fast method is needed to evaluate the required autocovariances.It is hoped to report on these developments in another paper.
S-PLUS, Mathematica and MATLAB are some other computing environments that are very popular for teaching and research and it is hoped that others may find it useful to port the ltsa package to these and other computer environments.
It should be noted that the predict.Arima function in R can not be used to obtain the above results, since it only predicts for a single forecast origin time which is fixed to be at the end of the series.This hampers the use of R in forecasting experiments in time series where the data is divided into training and test samples.Examples of this technique are given in Sections 3.3 and 3.4.

3. 3 .
Fitting FGN model with ARMA comparisonThe Nile minima time series consists of n = 663 observations of the annual minimum flow from 622 AD to 1284 AD.It is an often cited example of a geophysical time series exhibiting long-

Figure 1 :
Figure 1: Diagnostic plots produced by plot for the fitted FGN model to the NileMin series.

Figure 2 :R>
Figure 2: Normal probability plot of the standardized prediction residuals of the fitted FGN model to the NileMin series.
Table10: Comparison of models fitted to the Nile minima time series.L m is the concentrated log-likelihood, BIC = −2L m + 2k, BIC = −2L m + k log(n), n = 663 is the series length and and k is the number of parameters.

Table 1 :
The principal functions in ltsa.

Table 2 :
CPU time in seconds for TrenchInverse and the R function solve.

Table 4 :
is an iterative method which is more complicated to program and has several other practical limitations.If only the inverse matrix is required then the Trench algorithm is always computationally more efficient, since the superfast algorithm only solves a set of linear equations and does not directly compute the matrix CPU time in seconds for d = 0.45, n = 10, 000 and n = 15, 000.

Table 5 :
Comparison of CPU time in seconds for n = 5000.

Table 6
compares the average time needed for 100 simulations for various series lengths, n. >From this table we see that DHSimulate is overall faster even though it is entirely implemented in R.

Table 6 :
CPU time in seconds for 100 simulations of a time series of length n.In general, if it is planned to do many simulations, DHSimulate may be preferred provided the Davies-Harte condition is met.As described in the documentation for DHSimulate, one can use the function DHSimulate itself to test, if the Davies-Harte condition is satisfied.Both DHSimulate and DLSimulate can also be used to generate non-Gaussian time series by setting the optional argument rand.gen to some other distribution.However, only SimGLP and arima.simallow one to specify the exact innovation distribution.See online documentation for examples.

Table 7 :
. The principal functions in FGN.

Table 8 :
Average CPU time to simulate and fit an FGN series of length n.

Table 9 :
Efficiency of the sample mean in FGN

Table 11 :
Table 12 is available with our paper.It is common practice to divide the data in two parts: a training sample and a test sample.The model is fit to the training sample and then its performance is evaluated on the test sample.In time series analysis, experiments are often reported in which an initial portion Forecasts and their standard deviations for NileMin.

Table 12 :
RMSEs for forecasts for the last 100 values in the NileMin dataset.In this case, the model was refit to each series for each new data value and the forecast for the updated model was computed.For comparison, the sample standard deviation of the series over the forecasting period was 0.733.

Table 14 :
Properties of the MLE for H.

Table 15 :
is provided for the interested reader.Average RMSE in B = 10 4 bootstrap iterations.

Table 16 :
). Intervention models fit the annual Nile riverflow series.The parameter β = φ or H correspond to the AR(1) and FGN models.