Improved Subset Autoregression: With R Package

The FitAR R ( R Development Core Team 2008) package that is available on the Comprehensive R Archive Network is described. This package provides a comprehensive approach to ﬁtting autoregressive and subset autoregressive time series. For long time series with complicated autocorrelation behavior, such as the monthly sunspot numbers, subset autoregression may prove more feasible and/or parsimonious than using AR or ARMA models. The two principal functions in this package are SelectModel and FitAR for automatic model selection and model ﬁtting respectively. In addition to the regular autoregressive model and the usual subset autoregressive models (Tong 1977), these functions implement a new family of models. This new family of subset autoregressive models is obtained by using the partial autocorrelations as parameters and then selecting a subset of these parameters. Further properties and results for these models are discussed in McLeod and Zhang (2006). The advantages of this approach are that not only is an eﬃcient algorithm for exact maximum likelihood implemented but that eﬃcient methods are derived for selecting high-order subset models that may occur in massive datasets containing long time series. A new improved extended BIC criterion, UBIC, developed by Chen and Chen (2008) is implemented for subset model selection. A complete suite of model building functions for each of the three types of autoregressive models described above are included in the package. The package includes functions for time series plots, diagnostic testing and plotting, bootstrapping, simulation, forecasting, Box-Cox analysis, spectral density estimation and other useful time series procedures. As well as methods for standard generic functions including print , plot , predict and others, some new generic functions and methods are supplied that make it easier to work with the output from FitAR for bootstrapping, simulation, spectral density estimation and Box-Cox analysis.


Introduction
The family of AR(p) models for a time series z t , t = 1, 2, . . ., may be written in operator notation, φ(B)(z t − µ) = a t , where φ(B) = 1 − φ 1 B − . . .− φ p B p , µ is the mean of the series and a t ∼ NID (0, σ 2 a ).It is assumed that the parameters are in the admissible region so that all roots of the equation φ(B) = 0 lie outside the unit circle.The usual family of subset AR(p) models is obtained by taking a subset of the parameters φ 1 , . . ., φ p , and may be written as φ(B)(z t − µ) = a t , where φ(B) = 1 − φ i 1 B − . . .− φ im B im .This family of subset AR models may be denoted by ARp(i 1 , . . ., i m ), where 0 ≤ i 1 ≤ i 2 . . .≤ i m ≤ p1 .

Time series plots
Although time series plots can easily be produced with the built-in R function plot or plot method for 'ts' objects, there are several advantages to using our function TimeSeriesPlot.As pointed out in Cleveland (1993), many stationary time series are best visualized with a relatively low aspect-ratio and it is awkward to do this with standard R graphics.The function TimeSeriesPlot adjusts the aspect-ratio using par settings.In the lattice package (Sarkar 2008), the function xyplot also allows one to control the aspect-ratio2 .In general, Cleveland (1993) recommends a procedure, banking to 45 o .This algorithm determines the aspect-ratio so as to make the average slope 45 o .For time series this is only a rough guide since, for stationary time series, the resulting aspect-ratio will often be too close to zero, and for time series containing strong trends, it makes more sense to bank to the underlying trend.When strong trends are present, it is usually reasonable to use aspect = 1 but the user should be prepared to experiment with what aspect-ratio is most effective.
Figure 1 and Figure 2 compare the plot method for 'ts' objects and TimeSeriesPlot using the built-in lynx time series.Figure 2 with the lower aspect-ratio shows the asymmetry in the population cycles more clearly.
When plotting long time series the resolution may be improved by using a multipanel display which can be easily created using xyplot.Our function TimeSeriesPlot implements a lattice display for long time series.Figure 3 illustrates the multipanel time series display for the Ninemile treering time series.

Partial autocorrelation plot
The R function, pacf, plots the usual partial autocorrelations and their 95% confidence limits.Our function PacfPlot3 is more useful for identifying ARz subset models.Given a maximum order P for the model, the estimated partial autocorrelations ζ1 , . . ., ζP are computed along with their estimated standard errors.The estimated standard errors are computed from the large-sample covariance matrix of ζ1 , . . ., ζP given in McLeod and Zhang (2006).The plot produced by PacfPlot shows the 95% confidence interval for each individual of ζi , i = 1, . . ., P .From this plot, the nonzero subset of partial autocorrelations ζ 1 , . . ., ζ P for the ARz model may be selected.The use of PacfPlot is illustrated in Figure 4 for the natural logarithms of the built-in lynx time series.This plot shows that 95% confidence intervals for ζ1 , ζ2 , ζ4 , ζ7 , ζ10 , and ζ11 do not include zero and hence the subset partial autocorrelations at lags 1, 2, 4, 7, 10 and 11 are significantly different from zero at the 5% level.Thus the ARz (1,2,4,7,10,11) model is suggested.Since lags 4 and 10 are barely significant at 5%, an ARz(1, 2, 7, 11) model may also be considered.

Automatic model selection
It is often preferable to use an automatic model selection method.We consider automatic methods for the non-subset and the subset AR models.In all model selection algorithms, we will assume that the mean parameter is included in the model and the mean is estimated by the sample mean.

Non-subset models
The approximate concentrated log-likelihood for an AR(p) model may be written, where n is the series length.When p = 0, L = 0.The AIC model selection criterion,  is evaluated for p = 0, . . ., P , where P is the maximum order4 .Similarly for the BIC, BIC = −2L + p log n. (5) Initially κ 0 candidate models are selected using the AIC/BIC evaluated using the approximate likelihood in Equation 3.These models are then refit using exact likelihood method described in Section 3 and new AIC/BIC values are computed.From these models the best κ 1 < κ 0 are output as the final selection.As noted by Duong (1984) and Choi (1992), it is often useful to consider not just the best model but also close runner-ups.So in this case we would choose κ 1 > 1.In the function SelectModel, κ 0 and κ 1 correspond to the arguments Candidates and Best which have default settings, Candidates = 5 and Best = 3.
The example below illustrates why this two-step procedure is needed.We consider the highorder AR models fit by Percival and Walden (1993, Chapter 10) to the Willamette river data.Using the Akaike FPE, which is asymptotically equivalent to the AIC (Nishii 1984), Percival and Walden (1993, Section 10.15) find that the best model is an AR(38) which agrees with our function SelectModel: R> SelectModel(log(Willamette), lag.max = 150, Criterion = "AIC") p AIC AICx 1 38 -2474.058 -485.9111 2 27 -2473.823 -484.9443 3 39 -2472.308 -485.0843The AIC-column is based on the exact log-likelihood given in Equation 8 and the AICxcolumn is based on the approximate likelihood in Equation 3. In this case both methods agreed that the best model was with p = 38 but the exact method selected p = 27 as second best whereas the approximate method would have selected p = 39.
The default order selection criterion in SelectModel is BIC since in forecasting applications this often provides the most suitable model (Koehler and Murphree 1988;Granger and Jeon 2004;Hipel and McLeod 1994).For spectral density estimation using non-subset models, many authors prefer the AIC (Choi 1992;Percival and Walden 1993;Hipel and McLeod 1994).

Subset models
In the case of subset model selection of high-order AR, the model space can become very large.When the model space is large like this, the usual AIC/BIC criteria may select overparameterized models (Chen and Chen 2008;Broman and Speed 2002).Chen and Chen (2008) have developed the UBIC criterion for this large model space problem.Therefore, for subset model selection, we have implemented UBIC as the default in SelectModel5 .For comparison, the usual AIC/BIC criteria are also implemented.The computation of these criteria is outlined next.
The approximate concentrated log-likelihood function for an ARz(i 1 , . . ., i m ) evaluated at ζ1 , . . ., ζm may be written (McLeod and Zhang 2006, Equation 15), where n is the series length.The UBIC criterion (Chen and Chen 2008) may be written, where P is the maximum possible lag.Any value of P such that P ≥ i m produces an equivalent ordering of the UBIC.
The usual AIC/BIC criteria defined in Equation 4and Equation 5 may also be used but these criteria are not recommended for the subset problem.For the model order selection problem in the non-subset case, the UBIC reduces to the BIC.
The UBIC or AIC/BIC subset can be found simply by ordering the parameters in ascending order of magnitude, | ζi 1 | ≤ | ζi 2 | ≤ . . .≤ | ζi P |, and then evaluating successively L starting with L = 0 for the null model containing only the mean parameter 6 .As in the non-subset model selection algorithm, a second pass is made after an initial number of candidate models are found, in which each model is refit using exact maximum likelihood.For the log(lynx) data, we find using the UBIC, R> SelectModel(log(lynx), lag.max = 15, ARModel = "ARz", Best = 1) [1] 1 2 7 10 11 and this agrees well with the model suggested by the PacfPlot.To summarize the procedure used in this example, the best 5 models were found using the approximate likelihood to estimate the UBIC.Then for each of these 5 models, the exact likelihood was used to determine the UBIC for each model.Then the best model, the ARz(1, 2, 7, 10, 11) was found.For this data, the model space is not very large and consequently, it may be shown that the BIC selects the ARz(1, 2, 7, 10, 11).
Model identification for subset ARp models is supported using the R subset regression package leaps (Lumley and Miller 2004).Using this package the best ARp subset model is identified using a two-step procedure.In the first step, the leaps function is used to find the ARp(i 1 , . . ., i m ) subsets which when fit using least squares have minimum residual variance for each i m , m = 0, 1, . . ., P parameters.For the next step, all P models are refit using our least squares AR fitting algorithm, GetFitARpLS, and the exact log-likelihood is determined using LoglikelihoodAR.Based on these log-likelihoods, the UBIC/BIC/ AIC is re-evaluated and the best model is selected.
As shown below, the model selected by Tong (1977) for the lynx series, ARp(1, 2, 4, 10, 11), is selected using SelectModel, R> SelectModel(log(lynx), lag.max = 15, ARModel = "ARp", Criterion = "BIC", + Best = 1)) [1] 1 2 4 10 11 however, the UBIC selects the more parsimonious ARp(1, 2, 9, 12).When Best > 1, the output from SelectModel is a list and is defined as a 'SelectModel' S3 class object.This output may be graphically viewed using the plot function.For the lynx dataset, Figure 5 compares the 3 best ARp and ARz subset models selected using the UBIC and BIC.In Figure 5, the left vertical axis shows the scale in terms of BIC while the right one uses relative plausibility, defined by exp{(BIC 1 − BIC 0 )/2}, where BIC 0 and BIC 1 denote the BIC values for the best model and the alternative respectively.Likelihood and plausibility for time series models were introduced by Akaike (1978) and Akaike (1979).As with the relative likelihood (Fisher 1959;Royall 1997;Sprott 2000), alternative models with relative plausibility less than 1% may be said to be implausible.As shown in Figure 5, the BIC and UBIC select the same best three ARz models although for the ARp a more parsimonious model is selected using the UBIC than with the BIC.
The best subset for ARz models is easily found even when P is large since only sorting is needed.On the other hand, for the ARp best subset selection, the leaps package (Lumley and Miller 2004), available on the Comprehensive R Archive Network, is used and this package implements an exponential-time branch-and-bound algorithm.The function used, leaps, is interfaced to Fortran for maximum efficiency.But still it is limited to a maximum order of around P = 50 or so due to the rapidly increasing computing time needed as P increases.Thus it is possible to fit much higher order subset models using the ARz model.This is illustrated with the Zurich monthly sunspots 1749 − 1983 available in R as a built-in dataset, sunspots.Taking a maximum model order of P = 200 and using a square-root transformation, the best non-subset AR and subset ARz models were determined.In the case of AR model selection the AR(27) and AR(21) models were selected using the AIC and BIC criteria.For the subset selection, AIC , BIC and UBIC were used and the number of parameters not including the mean in the selected models were 55, 18 and 8 respectively for AIC , BIC and UBIC. Figure 6 compares the spectral density functions for the estimated models.The AIC subset model spectral density was very jagged and is not shown in Figure 6 7 .

Estimation
The log-likelihood function given n time series observations z 1 , . . ., z n from a covariance sta-7 See Example 3 in the documentation for sdfplot for an R script to estimate the spectral density function for the monthly sunspot series by using the AIC criterion to fit an ARz model.tionary Gaussian time series with mean µ and autocovariance function γ k , k = 0, 1, . .., may be determined from the multivariate normal distribution.Letting M n = Γ n /σ 2 a where Γ n is the covariance matrix of n successive observations and maximizing the log-likelihood function over σ 2 a , the concentrated log-likelihood may be written, after dropping constants, where S(φ) = z M −1 n z, g n = det(M n ) and σ2 a = S(φ)/n where φ = (φ 1 , . . ., φ p ).Note that L given in Equation 6 is an affine large-sample approximation to the exact log-likelihood given in Equation 8.
Assuming that the mean is known, the exact log-likelihood function may be computed in O(1) flops independent of n in repeated function evaluations after an initial setup computation is performed (McLeod and Zhang 2006, Section 2.1).This is done by conditioning on the sufficient statistics, D i,j = z i z j + • • • + z n−j z n−i , where we have assumed a mean-zero time series.The R function optim is then used to maximize the log-likelihood function.Initial starting values are obtained using Burg estimates.This exact maximum likelihood estimation algorithm is implemented in our functions FitAR and GetARFit.The function GetARFit may be preferred for bootstrapping applications while FitAR is intended as the main function for model fitting when diagnostic checking and forecasting are important.
FitAR also computes the covariance matrix of the estimates, residuals, fitted values and some diagnostic checks.The function FitAR returns a S3 class object 'FitAR' and methods functions print.FitAR, summary.FitAR, coef.FitAR, residual.FitAR, and plot.FitAR are supplied.In addition some new generic functions Boot and sdfplot are defined and methods are given for class 'FitAR', 'ts' and other objects of interest.More methods are supplied for these functions in McLeod, Yu, and Krougly (2007).
FitAR and getFitAR use several other functions which may be of interest in some applications and are described briefly in Appendix A. In Appendix B, we discuss the perils of fitting the ARp model using exact MLE.This provides another reason for preferring the ARz model.

Exact MLE for mean
Assuming that (φ 1 , . . ., φ p ) are known, the exact MLE for mean µ is given by μ where 1 n denotes the n dimensional column vector with all entries equal to 1, z = (z 1 , . . ., z n ) and Γ −1 n denotes the inverse of the covariance matrix of n successive observations.Since μ does not depend on σ 2 a , we may assume without loss of generality that σ 2 a = 1.Direct evaluation of Equation 9 using the exact inverse matrix derived by Siddiqui (1958) would require O(n 2 ) flops.A more efficient algorithm is given in McLeod and Zhang (2008) for evaluating Equation 9 in O(n) flops and this is implemented in our R function GetARMeanMLE.As a check for the correctness of GetARMeanMLE, in the example section of the help documentation for GetARMeanMLE, we provide a R script which compares the results obtained by computing the MLE of the mean µ directly and using GetARMeanMLE.
When (φ 1 , . . ., φ p ) is unknown, the iterative algorithm given in McLeod and Zhang (2008, Section 2.2) may be used to obtain the simultaneous joint MLEs of µ and the other model parameters.This iterative algorithm is implemented in FitAR with the option MeanMLEQ = TRUE.In the documentation for FitAR, an example is given to demonstrate the agreement with results obtained by the built-in function arima.Some simulations were carried out to compare the timings for GetFitAR, FitAR and the builtin function ar for fitting nonsubset AR(p), p = 1, 2 and p = 20, 40.The main purpose of these timings is just to indicate the relative computer times for the three different algorithms.As far as accuracy is concerned, there is no meaningful difference in accuracy between FitAR and ar when the latter function converges.As noted below, ar does not always work and sometimes it fails completely.It seems reasonable to recommend that users avoid the exact MLE option in ar.Between GetFitAR and FitAR there are some slight differences in the estimates but overall there does not seem to be any notable difference in accuracy for the models in this experiment.
For each model, the admissible parameters were randomly selected by randomly selecting partial autocorrelations uniformly from −1 to 1 and then reparameterizing through Equation 2. For each n and each p, 100 series were simulated and the models estimated using the functions GetFitAR, FitAR and ar.The average time over all 100 simulations is reported for each parameter combination.These timings are illustrative of the relative performance of the functions but are highly system dependent.A 3.6 GHz Pentium PC with R Version 2.70 running Windows XP was used.
Various lengths of series were used from n = 50, 100, 200, 500, 1000.The timings do not vary much once n and p are fixed but the timing for a single fitting is too short to get an accurate value for low order models, so 100 simulations were done for each of these parameter settings.The average CPU times required for simulations and fitting are given below in Table 3.During fitting, ar reported one error in optim which is shown below: Error in optim(init[mask], arma0f, method = "BFGS", hessian = TRUE, control = optim.control):non-finite finite-difference value [3] This error was trapped using the built-in R function try.Many warning messages about convergence were also generated by ar.On the other hand GetFitAR or FitAR worked correctly in all cases without generating any messages.The experiment was repeated, this time with p = 20, 40 and n = 1000, 2000, 5000.Once again ar generated error and warning messages.With GetFitAR, the sample mean was used but for FitAR the iterative algorithm was used to compute the exact MLE for the mean as well as the other parameters.With ar the option method = "mle" was used.All timings reported in this paper are for a Pentium 4 processor at 3 GHz PC running Windows XP with 2 GB RAM.The CPU times in Table 3 show that FitAR is much faster than ar for p = 40.If the sample mean is used then GetFitAR is always very fast and for most AR models that seem to occur in practice, there is not much difference between using the sample mean and the exact MLE. for data z t , t = 1, . . ., n, from a normal linear regression and ANOVA.This method was adapted to ARIMA time series in Hipel and McLeod (1977, Equation 10, p.571).For autoregressive models, the concentrated loglikelihood function can be written,

Box-Cox analysis
For fixed λ we can maximize over φ to obtain L(λ).The function L(λ) may be maximized numerically with the optimize function to obtain the MLE, λ.It is convenient to plot the relative likelihood function R(λ) = L(λ)/L( λ).The relative likelihood, R(λ), itself provides a quantification of the relative plausibility of various values of λ (Sprott 2000, Section 2.4 and Section 4.5).It may be shown that a 95% confidence interval based on the likelihood-ratio test corresponds to R(λ) ≥ 0.1465.A horizonal line is drawn to indicate a 95% confidence interval for λ.
The function BoxCox is implemented as a generic function with methods for classes 'FitAR', 'Arima', 'ts' and 'numeric'.BoxCox method for 'FitAR' objects parses its input to construct  the likelihood function in Equation 10and then the MLE is determined using optimize.A plot of the likelihood function is determined to show all values of λ with plausibility greater than 1% as well as a horizontal line indicating the 95% confidence interval.In Figure 7, we compare the Box-Cox analysis for the lynx time series using our function and the boxcox function in the R package MASS (Venables and Ripley 2002).For the lynx time series we fit the ARp (1,2,4,10,11) model by least-squares using FitARp.In this case, the output from FitARp contains the design matrix for the regression as well as the column for the dependent variable.These results were then used with the MASS boxcox function to obtain the bottom plot in Figure 7.The difference in the plots is substantial.The top plot using BoxCox method for the 'FitAR' object shows that a log transformation is reasonable choice whereas the bottom plot strongly suggests otherwise and favors a transformation such as a cube-root.
The reason for the difference is that in BoxCox method for 'FitAR' objects the log-likelihood in Equation 10 is used whereas in boxcox the approximate log-likelihood corresponding to linear regression is used.
BoxCox method for 'ts' and 'numeric' objects enable one to do a Box-Cox analysis for arbitrary time series as well as simple random samples.With BoxCox method for a 'ts' object, the untransformed data is fit by a high-order AR(p) model which is then used for Box-Cox analysis.In practice, this procedure often works well even for non-stationary time series as is shown in Figure 8 for the R dataset AirlinePassengers.An example of using BoxCox method for 'numeric' objects for the R built-in rivers data set is shown in Figure 9.
The BoxCox method for 'Arima' objects implements a Box-Cox analysis for time series fit with the R arima function.An illustrative example is provided by the time series of annual production of tobacco in the U.S. for the period 1871-1984 (See Figure 10).Note the increase in variability.Since this variability is related to level, a model of the Box-Cox transformed series provides a much simpler alternative to the more complex ARIMA-GARCH approach (Wei 2005, p.379).Wei (2005, p.120) also suggested an ARIMA(0, 1, 1) to the logarithms of this series.However Box-Cox analysis indicates that a square-root transformation works much better than logarithms (Figure 11).

Diagnostic checks
A methods function is supplied so that the generic function plot can be used for a variety of diagnostic plots for 'FitAR' class objects.The default is with option terse = TRUE which produces a plot of the p-values of the Ljung-Box portmanteau test (Ljung and Box 1978) and a plot of the residual autocorrelations with simultaneous 95% limits (Hosking and Ravishanker 1993).Both these plots appear as a panel display.These are usually the most important diagnostics.The default diagnostic plots shown in Figure 12 are produced by: R> ans <-FitARp(log(lynx), c(1, 2, 4, 10, 11)) R> plot(ans) When terse = FALSE, after the plots shown in Figure 12, seven more plots are produced in the order as follows: Normal probability plot and the p-value for an ominbus normality test of the Jarque-Bera test (Jarque and Bera 1987).This is useful in detecting outliers or thick tails.Thick tails are a characteristic of ARCH/GARCH models (Tsay 2005, Chapter 3) so this test may suggest that some form of conditional heteroscedasticity is present.
Box-and-whisker plot of the residuals.This is useful in detecting skewness.Skewness may indicate a power transformation.On the other hand, if the residuals are symmetric but outliers are present, some form of conditional heteroscedasticity may be present.
Various types of ARCH/GARCH models may be used for modeling conditional heteroscedasticity.
Plots of the time series and a bootstrap version of the time series.This gives some idea of whether or not important features in the time series appear in the simulated model.
Monotone spread plot (Cleveland 1979(Cleveland , 1993)).This type of residual diagnostic plot is useful for detecting situations where the variance depends on the level or mean.In general Cleveland (1993) suggested plotting the square root of the absolute residual vs. the fitted value and then visualizing the relationship using a suitable loess smooth.The R function lowess with f = 1 is used.For time series models, the fitted value represents a conditional mean.If the model is adequate, the loess curve should be approximately horizontal.If it is monotonic up or down, a power transformation could be used to remove this dependence.
Residual-fit spread (RFS) plot (Cleveland 1993).This plot is useful for visualizing how much of the variation in the data is explained by model and so whether the model may be useful for forecasting.Alternatively, the coefficient of determination, R 2 , defined as the ratio of the variance of the one-step forecast divided by the series variance may be used.For the lynx and ARp(1, 2, 7, 10, 11), R 2 = 85%, as shown below, R> z<-log(lynx) R> var(fitted(FitARp(z,c(1, 2, 7, 10, 11))))/var(z) [1] 0.8536185 The RFS plot is shown in Figure 13.Plots of the observed autocorrelation function, the fitted theoretical autocorrelation function and the sample autocorrelation of a parametric bootstrap.In many cases a good model can be expected to match the main features of the sample autocorrelations, particularly at the lower lags.However, caution is necessary since it is well known that the sample autocorrelations may be highly autocorrelated with large variances (Box, Jenkins, and Reinsel 1994, Section 2.1.6).In some cases, such as highly autocorrelated fractional noise, the sample autocorrelations may have very large biases (Newbold and Agiakloglou 1993).In this situation, even when the model has been correctly specified and fit, the sample autocorrelations may not resemble very well the theoretical autocorrelations.This is quite a surprising, but in hindsight, not an unreasonable result.To the degree that fractional noise may be approximated by a high-order AR model, this problem can be expected to occur with some AR models as well.In summary, we can say that it may happen that even for a correctly fit model, the sample and theoretical autocorrelations might not match very well and there may be more similarity in the sample autocorrelations of the data and the bootstrap version.
Plot of the fitted spectral density function.This plot provides a frequency domain summary of the model.
Because of the large number of graphics windows, the default behavior of plot is to delete all graphics windows before producing the plots.This default can be bypassed by using the optional argument clearGraphics = FALSE.Some of the diagnostic checks are also available as separate functions as shown in Table 4.For more details, please see the online documentation.

Other useful time series functions
We give a brief overview of the functions shown in Table 5.For more details and examples, please see the online R documentation.

AcfPlot: Correlation plots
The function AcfPlot is useful for plotting other types of correlation functions, for example, the residual autocorrelations or the inverse autocorrelations.Cleveland (1971) used the inverse autocorrelations to select an ARp(1, 2, 7) model for the Series A time series.The example given in the online help documentation for AcfPlot shows how to compute and plot the inverse autocorrelations.

AR1Est: Exact MLE for AR(1)
The function AR1Est evaluates the exact MLE for a mean-zero AR(1) model using a closed form solution (Zhang 2002).This function is useful in bootstrapping and Monte-Carlo since it is fast, accurate and reliable.

ARToMA: Moving-average approximation
This function is used in computing the variances of the residual autocorrelations and in constructing plots of the residual autocorrelations.This type of computation arises also in other computations such as in the confidence limits for forecasts from AR models.

ARToPacf and PacfToAR: Reparameterization
These functions are central to working with ARz models.For many purposes it is necessary to convert from one parameterization to the other.

BackcastResidualsAR: Compute residuals
The innovation residuals for a fitted AR(p) may be computed recursively from for t = p + 1, . . ., n.For t = 1, . . ., p, the conditional expected value of a t may be computed by backforecasting the z t for t = 0, −1, − . . ., Q for Q large enough so that the backforecast value of z t is approximately µ for all t ≤ Q.Then the required residuals may be computed directly using Equation 12.This backcasting approach is described in detail in Box et al. (1994, Section 6.4.3).

Concatenation of time series
The function cts allows values to be easily concatenated to an existing time series object.See documentation for more details.

Information matrix
The large-sample Fisher information per observation can be obtained for AR, ARp and ARz models.This is used in FitAR to obtain the estimated standard errors of the estimated parameters.Standard errors can also be obtained by using the observed Fisher information that is obtained by numerically differentiating the log-likelihood function.Our preference is to use the theoretical information matrix since it is well-known that numerical differentiation may be unreliable (Fröberg 1969, Chapter 9, p.192).We have found for simulations and bootstrapping applications this unreliability does occasionally arise resulting in negative variances and other difficulties.Standard errors may also be estimated by bootstrapping (Box and Luceño 1997).for 'FitAR' class objects.It took about 18 seconds for the 100 bootstrap replications used in Table 6.

PacfDL: Partial autocorrelations
PacfDL implements the Durbin-Levinsion algorithm to compute the partial autocorrelations and the optimal linear predictor given a sequence of autocovariances.This function is not used in our package but we include it since it can be useful in many applications.For example, Hipel and McLeod (1977) and Hipel and McLeod (1994) pointed out that the inverse partial autocorrelations may be useful in identifying MA(q) time series models.After computing the inverse autocorrelations we can then directly compute the inverse partial autocorrelations using PacfDL.An example of this is given in the online documentation.

Prediction
FitAR produces a 'FitAR' class object and a predict method is implemented for this object.This function uses TrenchForecast in ltsa (McLeod et al. 2007).See documentation for illustrative examples of its usage.

Readts: Time series input
The Readts function is another function which is not used directly in our package but that we have found very helpful.Many time series data are usefully stored in ASCII files with titles and comments that serve as documentation.These ASCII datasets may then be input to other software although sometimes some editing is necessary to remove the documentation, as for example, if the R function scan is used.Our Readts function is especially convenient since it can input ASCII data files containing documentation and title information.In interactive mode, Readts can be helpful in setting up the necessary parameters for a 'ts' object.When Readts is used, a 'ts' object is created with attribute title.This attribute, if present, is used to put a title on plots created by TimeSeriesPlot.

SimulateGaussianAR: Autoregressive simulation
The built-in R arima method for 'sim' objects provides for simulation of AR and more general ARIMA models with non-Gaussian innovations.It simulates by using the model equation directly with arbitrary starting values.Strictly speaking this process is not stationary but after a burn-in period, it may closely approximate a stationary process.To avoid the problem of choosing a burn-in period, which may not be adequate or which may be more than adequate and hence computationally inefficient, we can use an exact technique for Gaussian time series suggested by McLeod (1975) for the ARMA case.In the AR(p) case, we simply generate p initial time series values z 1 , . . ., z p using the appropriate multivariate normal distribution.Then z t , t = p + 1, . . ., n may be calculated directly from the model equation.

TacvfAR and TacvfMA: Theoretical autocovariance functions
The built-in function ARMAacf is quite complicated and uses an interface to C. Our function TacvfAR is based on the exact algorithm given by McLeod (1975) and is much simpler.It is very easy to translate our TacvfAR function into other programming environments such as MATLAB.
The TacvfMA is used by GetMeanMLEAR.TacfMA uses an efficient vectorized approach to compute the theoretical autocorrelations for a moving-average process.The R code is very simple and easy to understand.For an MA(q), z t = a t −θ 1 a t−1 −. ..−θ q a t−q , the non-zero autocovariances can be written as the product of a square matrix of order q + 1 times a column vector.In R, this matrix multiplication is most efficiently computed using crossprod (Venables and Ripley 2002, Section 3.9).

Conclusion
A complete suite of functions is described for selecting, identifying and fitting autoregressive and subset autoregressive models.Two families of subset autoregressive models are fully supported.The subset autoregressive models may provide a parsimonious alternative to the more complex ARMA models.In particular, the ARz family of models is suitable for modeling complex long time series with high-order lags.For such time series, it is difficult to select and estimate the parameters when using standard ARMA models.The ARz with the UBIC model selection criterion provides a practical and perhaps better alternative to the ARMA model for such series.Our package FitAR is available from the Comprehensive R Archive Network at http://CRAN.R-project.org/package=FitAR.
R scripts to generate all figures and simulations reported in Table 3

Figure 1 :Figure 2 :
Figure1: The lynx series with default aspect-ratio with plot method for 'ts' objects.

FirFigure 3 :
Figure3: Trellis or lattice style graphics are useful for displaying long time series as well as for choosing a suitable aspect-ratio.

Figure 4 :
Figure 4: Partial autocorrelation plot for identifying an ARz model to the lynx time series.

Figure 5 :
Figure 5: Plot of the output from SelectModel with Best = 3 for ARp and ARz models using BIC and UBIC for subset selection for the logged lynx series.

Figure 6 :
Figure 6: Logged spectral density functions of the non-subset AR and subset ARz models fitted to the square-root of the sunspots series.

Figure 10 :Figure 11 :Figure 12 :
Figure 10: Time series plot for annual tobacco production series.

Table 2 :
Table 2 lists the three functions useful in initial model selection.Model identification functions.

Table 4 :
Model diagnostic checking functions.

Table 5 :
Other useful time series functions.

Table 6 :
Table 6 below compares the bootstrap estimates with the large-sample ones.The script to generate these results is included in the online documentation for Boot method Comparison of bootstrap and large-sample estimates of the standard deviations of the parameter estimates in an ARz(1, 2, 4, 7, 10, 11)fitted to the logged lynx time series using 100 bootstrap replications.
and Table 6 are available in the examples in the help for FitAR-package.Please see, R> help("FitAR-package")