Title: | Fitting and Forecasting Gegenbauer ARMA Time Series Models |
---|---|
Description: | Methods for estimating univariate long memory-seasonal/cyclical Gegenbauer time series processes. See for example (2022) <doi:10.1007/s00362-022-01290-3>. Refer to the vignette for details of fitting these processes. |
Authors: | Richard Hunt [aut, cre] |
Maintainer: | Richard Hunt <[email protected]> |
License: | GPL-3 |
Version: | 0.9.23 |
Built: | 2024-11-12 03:25:36 UTC |
Source: | https://github.com/rlph50/garma |
Approximate AIC for model.
## S3 method for class 'garma_model' AIC(object, ...)
## S3 method for class 'garma_model' AIC(object, ...)
object |
The garma_model object |
... |
Other parameters. Ignored. |
(double) Approximate AIC - uses approximation of whichever methoid is used to find model params.
The ggplot function generates a ggplot of actuals and predicted values for a "garma_model" object. This adds in sensible titles etc as best it can determine.
## S3 method for class 'garma_model' autoplot(object, h = 24, include_fitted = FALSE, ...)
## S3 method for class 'garma_model' autoplot(object, h = 24, include_fitted = FALSE, ...)
object |
(garma_model) The garma_model from which to ggplot the values. |
h |
(int) The number of time periods to predict ahead. Default: 24 |
include_fitted |
(bool) whether to include the 1-step ahead 'fitted' values in the plot. Default: FALSE |
... |
other parameters passed to ggplot. |
A ggplot2 "ggplot" object. Note that the standard ggplot2 "+" notation can be used to enhance the default output.
library(ggplot2) data(AirPassengers) ap <- as.numeric(diff(AirPassengers, 12)) mdl <- garma(ap, order = c(9, 1, 0), k = 0, method = "CSS", include.mean = FALSE) autoplot(mdl)
library(ggplot2) data(AirPassengers) ap <- as.numeric(diff(AirPassengers, 12)) mdl <- garma(ap, order = c(9, 1, 0), k = 0, method = "CSS", include.mean = FALSE) autoplot(mdl)
Model Coefficients/parameters.
## S3 method for class 'garma_model' coef(object, ...)
## S3 method for class 'garma_model' coef(object, ...)
object |
The garma_model object |
... |
Other parameters. Ignored. |
(double) array of parameter value estimates from the fitted model.
For a Gegenbauer process, transform to remove Gegenbauer long memory component to get a short memory (ARMA) process.
extract_arma(x, ggbr_factors)
extract_arma(x, ggbr_factors)
x |
(num) This should be a numeric vector representing the Gegenbauer process. |
ggbr_factors |
(class) Each element of the list represents a Gegenbauer factor and includes f, u and fd elements. |
An object of same class as x. Any time series attributes of x are copied to the returned object.
data(AirPassengers) ap <- as.numeric(diff(AirPassengers, 12)) # find semiparametric estimates of the Gegenbauer parameters. sp <- ggbr_semipara(ap) # extract the underlying short-memory ARMA process ap_arma <- extract_arma(ap, sp) summary(arima(ap_arma, order = c(1, 0, 0)))
data(AirPassengers) ap <- as.numeric(diff(AirPassengers, 12)) # find semiparametric estimates of the Gegenbauer parameters. sp <- ggbr_semipara(ap) # extract the underlying short-memory ARMA process ap_arma <- extract_arma(ap, sp) summary(arima(ap_arma, order = c(1, 0, 0)))
Fitted values are 1-step ahead predictions.
## S3 method for class 'garma_model' fitted(object, ...)
## S3 method for class 'garma_model' fitted(object, ...)
object |
The garma_model object |
... |
Other parameters. Ignored. |
(double) array of 1-step ahead fitted values for the model.
The forecast function predicts future values of a "garma_model" object, and is exactly the same as the "predict" function with slightly different parameter values.
## S3 method for class 'garma_model' forecast(object, h = 1, newdata = NULL, ...)
## S3 method for class 'garma_model' forecast(object, h = 1, newdata = NULL, ...)
object |
(garma_model) The garma_model from which to forecast the values. |
h |
(int) The number of time periods to predict ahead. Default: 1 |
newdata |
(real vector or matrix) If the original model was fitted with the 'xreg=' option then this will provide the xreg values for predictions. If this is a vector then its length should be 'h'; if it is a matrix then it should have 'h' rows. It should have columns with the same names as the original xreg matrix. |
... |
Other parameters passed to the forecast function. For "garma_model" objects, these are ignored. |
- a "ts" object containing the requested forecasts.
library(forecast) data(AirPassengers) ap <- as.numeric(diff(AirPassengers, 12)) mdl <- garma(ap, order = c(9, 1, 0), k = 0, method = "CSS", include.mean = FALSE) forecast(mdl, h = 12)
library(forecast) data(AirPassengers) ap <- as.numeric(diff(AirPassengers, 12)) mdl <- garma(ap, order = c(9, 1, 0), k = 0, method = "CSS", include.mean = FALSE) forecast(mdl, h = 12)
The GARMA package provides the main function "garma" as well as print, summary, predict, forecast and plot/ggplot options.
The garma function is the main function for the garma package. Depending on the parameters it will calculate the parameter estimates for the GARMA process, and if available the standard errors (se's) for those parameters.
garma( x, order = c(0L, 0L, 0L), periods = NULL, k = 1, include.mean = (order[2] == 0L), include.drift = FALSE, xreg = NULL, method = "Whittle", d_lim = c(0, 0.5), opt_method = c("cobyla", "solnp"), control = NULL )
garma( x, order = c(0L, 0L, 0L), periods = NULL, k = 1, include.mean = (order[2] == 0L), include.drift = FALSE, xreg = NULL, method = "Whittle", d_lim = c(0, 0.5), opt_method = c("cobyla", "solnp"), control = NULL )
x |
(num) This should be a numeric vector representing the process to estimate. A minimum length of 96 is required. |
order |
(numeric vector) This should be a vector (similar to the stats::arima order parameter) which will give the order of the process to fit. The format should be list(p,d,q) where p, d, and q are all positive integers. p represents the degree of the autoregressive process to fit, q represents the order of the moving average process to fit and d is the (integer) differencing to apply prior to any fitting. WARNING: Currently only d==0 or d==1 are allowed. |
periods |
(num) This parameter can be used to specify a fixed period or set of periods for the Gegenbauer periodicity. For instance if you have monthly data, then it might be sensible (after an examination of the periodogram) to set periods = 12. The default value is NULL. Either 'periods' or 'k' parameters must be specified but not both - 'periods' implies fixed period(s) are to be used and 'k' implies that the periods should be estimated. |
k |
(int) This parameter indicates that the algorithm should estimate the 'k' frequencies as a part of the model. An alternative is the 'periods' parameter which can be used to specify exactly which periods should be used by the model. This parameter can also be interpreted as specifying the number of (multiplicative) Gegenbauer terms to fit in the model. |
include.mean |
(bool) A boolean value indicating whether a mean should be fit. Note that no mean term is fit if the series is integer differenced. |
include.drift |
(bool) A boolean value indicating whether a 'drift' term should be fit to the predictions. The default is to fit a drift term to the predictions if the process is integer-differenced. |
xreg |
(numeric matrix) A numerical vector or matrix of external regressors, which must have the same number of rows as x. It should not have any NA values. It should not be a data frame. The default value is NULL. Note that the algorithm used here is that if any 'xreg' is supplied, then a linear regression model is fit first, and the GARMA model is then based on the residuals from that regression model. |
method |
(character) This defines the estimation method for the routine. The valid values are 'CSS', 'Whittle', and 'WLL'. The default ('Whittle') method will generally return very accurate estimates quite quickly, provided the assumption of a Gaussian distribution is even approximately correct, and is probably the method of choice for most users. For the theory behind this, refer Giraitis et. al. (2001). The 'CSS' method is a conditional 'sum-of-squares' technique and can be quite slow. Reference: Robinson (2006), Chung (1996). Note that the paper of Chung (1996) was partially critisised by Giraitis et. al. (2001), however still contains useful results. 'WLL' is a new technique, originally developed by the author of this package and which appears to work well even if the
|
d_lim |
(list) the limits for the d parameter. The default is 'c(0,0.5)', which restricts the model to be stationary. However sometimes it is desirable to understand what the unrestricted value might be. |
opt_method |
(character) This names the optimisation method used to find the parameter estimates. This may be a list of methods, in which case the methods are applied in turn, each using the results of the previous one as the starting point for the next. The default is to use c('solnp', 'cobyla'). For some data or some models, however, other methods may work well. Supported algorithms include:
Note that the algorithms are selected to be those which do not require derivatives, even numerically calculated derivatives. The function being optimised by 'garma()' has a point of discontinuity at the minimum value - the point we are trying to find. This means that standard algorithms like BFGS et al. perform very poorly here. Note further that if you specify a value of 'k' > 1, then inequality constraints are required, and this will further limit the list of supported routines. |
control |
(list) list of optimisation routine specific values. |
The GARMA model is specified as
where
represents the short-memory Autoregressive component of order p,
represents the short-memory Moving Average component of order q,
represents the long-memory Gegenbauer component (there may in
general be k of these),
represents the degree of integer differencing, where as
represents the degree of fractional
differencing. Note that
is a value supplied by the user (the second number on the 'order=' parameter - similarly
to the way that the base R 'arima' function works) whereas
is estimated by this function.
represents the observed process,
represents the random component of the model - these are assumed to be uncorrelated but
identically distributed variates. Generally the routines in this package will work best if these have an approximate
Gaussian distribution.
represents the Backshift operator, defined by
.
when k=0, then this is just a short memory model as fit by the stats "arima" function.
An S3 object of class "garma_model".
Richard Hunt
C Chung. A generalized fractionally integrated autoregressive moving-average process. Journal of Time Series Analysis, 17(2):111-140, 1996. DOI: https://doi.org/10.1111/j.1467-9892.1996.tb00268.x
L Giraitis, J Hidalgo, and P Robinson. Gaussian estimation of parametric spectral density with unknown pole. The Annals of Statistics, 29(4):987–1023, 2001. DOI: https://doi.org/10.1214/AOS/1013699989
R Hunt, S Peiris, N Webe. A General Frequency Domain Estimation Method for Gegenbauer Processes. Journal of Time Series Econometrics, 13(2):119-144, 2021. DOI: https://doi.org/10.1515/jtse-2019-0031
R Hunt, S Peiris, N Weber. Estimation methods for stationary Gegenbauer processes. Statistical Papers 63:1707-1741, 2022. DOI: https://doi.org/10.1007/s00362-022-01290-3
P. Robinson. Conditional-sum-of-squares estimation of models for stationary time series with long memory. IMS Lecture Notes Monograph Series, Time Series and Related Topics, 52:130-137, 2006. DOI: https://doi.org/10.1214/074921706000000996.
Useful links:
data(AirPassengers) ap <- as.numeric(diff(AirPassengers, 12)) print(garma(ap, order = c(9, 1, 0), k = 0, method = "CSS", include.mean = FALSE)) # Compare with the built-in arima function print(arima(ap, order = c(9, 1, 0), include.mean = FALSE))
data(AirPassengers) ap <- as.numeric(diff(AirPassengers, 12)) print(garma(ap, order = c(9, 1, 0), k = 0, method = "CSS", include.mean = FALSE)) # Compare with the built-in arima function print(arima(ap, order = c(9, 1, 0), include.mean = FALSE))
For a Gegenbauer process, use semi-parametric methods to obtain short memory version of the process, then run a ggtsdisplay().
garma_ggtsdisplay(x, k = 1, ...)
garma_ggtsdisplay(x, k = 1, ...)
x |
(num) This should be a numeric vector representing the process to estimate. |
k |
(int) The number of Gegenbauer factors |
... |
additional parameters to pass to ggtsdisplay |
The purpose of this function is to ease the process of identifying the underlying short memory process.
A ggplot object.
data(AirPassengers) ap <- as.numeric(diff(AirPassengers, 12)) garma_ggtsdisplay(ap)
data(AirPassengers) ap <- as.numeric(diff(AirPassengers, 12)) garma_ggtsdisplay(ap)
Display the raw periodogram for a time series, and not on a log scale.
gg_raw_pgram(x, k = 1)
gg_raw_pgram(x, k = 1)
x |
(num) This should be a numeric vector representing the process to estimate. |
k |
(int) The number of Gegenbauer factors |
The standard "R" functions display periodograms on a log scale which can make it more difficult to locate high peaks in the spectrum at differing frequencies. This routine will display the peaks on a raw scale.
A ggplot object representing the raw periodogram
data(AirPassengers) ap <- as.numeric(diff(AirPassengers, 12)) sp <- ggbr_semipara(ap) print(sp)
data(AirPassengers) ap <- as.numeric(diff(AirPassengers, 12)) sp <- ggbr_semipara(ap) print(sp)
For a Gegenbauer process, use semi-parametric methods to estimate the Gegenbauer frequency and fractional differencing.
ggbr_semipara(x, periods = NULL, k = 1, alpha = 0.8, method = "gsp")
ggbr_semipara(x, periods = NULL, k = 1, alpha = 0.8, method = "gsp")
x |
(num) This should be a numeric vector representing the process to estimate. |
periods |
(num) This parameter can be used to specify a fixed period or set of periods for the Gegenbauer periodicity. For instance if you have monthly data, then it might be sensible (after an examination of the periodogram) to set 'periods = 12'. The default value is NULL. Either 'periods' or 'k' parameters must be specified but not both - 'periods' implies fixed period(s) are to be used and 'k' implies that the periods should be estimated. |
k |
(int) This parameter indicates that the algorithm should estimate the 'k' frequencies semi-parametrically, before estimating the degree of fractional differencing at each period. An alternative is the 'periods' parameter which can be used to specify exactly which periods should be used by the model. |
alpha |
(num) Default = 0.8 - This is the bandwidth for the semiparametric estimate, and should be between 0 and 1. Robinson (1994) indicated optimality for a (scaled) version of 'alpha' = 0.8, at least for the "lpr" 'method'. |
method |
(char) One of "gsp" or "lpr" - lpr is the log-periodogram-regression technique, "gsp" is the Gaussian semi-parametric technique. "gsp" is the default. Refer Arteche & Robinson (1998). |
An object of class "garma_semipara".
J Arteche and P Robinson. Semiparametric inference in seasonal and cyclical long memory processes. Journal of Time Series Analysis, 21(1):1–25, 2000. DOI: https://doi.org/10.1111/1467-9892.00170
P Robinson. Rates of convergence and optimal spectral bandwidth for long range dependence. Probability Theory and Related Fields, 99:443–473, 1994. DOI: https://doi.org/10.1007/BF01199901.
data(AirPassengers) ap <- as.numeric(diff(AirPassengers, 12)) sp <- ggbr_semipara(ap) print(sp)
data(AirPassengers) ap <- as.numeric(diff(AirPassengers, 12)) sp <- ggbr_semipara(ap) print(sp)
Provides a goodness-of-fit test for a GARMA Model, using Bartletts Tp test. This has been justified for long memory and for GARMA models by Delgado, Hidalgo and Velasco (2005).
gof(object)
gof(object)
object |
(garma_model) The garma_model to test. |
Invisibly returns the array of p-values from the test.
M Delgado, J Hidalgo, and C Velasco. Distribution free goodness-of-fit tests for linear processes. The Annals of Statistics, 33(6):2568–2609, 2005. DOI: https://doi.org/10.1214/009053605000000606.
The approximate likelihood for the model.
## S3 method for class 'garma_model' logLik(object, ...)
## S3 method for class 'garma_model' logLik(object, ...)
object |
The garma_model object |
... |
Other parameters. Ignored. |
Object of class "logLik" with values for the (approx) log-likelihood for the model
The plot function generates a plot of actuals and predicted values for a "garma_model" object.
## S3 method for class 'garma_model' plot(x, ...)
## S3 method for class 'garma_model' plot(x, ...)
x |
(garma_model) The garma_model from which to plot the values. |
... |
other arguments to be passed to the "plot" function, including h (int) - the number of periods ahead to forecast. |
An R "plot" object.
data(AirPassengers) ap <- as.numeric(diff(AirPassengers, 12)) mdl <- garma(ap, order = c(9, 1, 0), k = 0, method = "CSS", include.mean = FALSE) plot(mdl)
data(AirPassengers) ap <- as.numeric(diff(AirPassengers, 12)) mdl <- garma(ap, order = c(9, 1, 0), k = 0, method = "CSS", include.mean = FALSE) plot(mdl)
Predict ahead using algorithm of Godet (2009).
## S3 method for class 'garma_model' predict(object, n.ahead = 1, newdata = NULL, ...)
## S3 method for class 'garma_model' predict(object, n.ahead = 1, newdata = NULL, ...)
object |
(garma_model) The garma_model from which to predict the values. This should have been generated by the [garma()] function. |
n.ahead |
(int) The number of time periods to predict ahead. Default: 1 |
newdata |
(real vector or matrix) If the original model was fitted with the 'xreg=' option then this will provide the xreg values for predictions. If this is a vector then its length should be 'n.ahead'; if it is a matrix then it should have 'n.ahead' rows. It should have columns with the same names as the original xreg matrix. |
... |
Other parameters. Ignored. |
A "ts" object containing the requested forecasts.
Godet, F. Linear prediction of long-range dependent time series, ESAIM: PS (2009) 13 115-134. DOI: https://doi.org/10.1051/ps:2008015
data(AirPassengers) ap <- as.numeric(diff(AirPassengers, 12)) mdl <- garma(ap, order = c(9, 1, 0), k = 0, method = "CSS", include.mean = FALSE) predict(mdl, n.ahead = 12)
data(AirPassengers) ap <- as.numeric(diff(AirPassengers, 12)) mdl <- garma(ap, order = c(9, 1, 0), k = 0, method = "CSS", include.mean = FALSE) predict(mdl, n.ahead = 12)
The print function prints a summary of a "garma_model" object, printed to the output.
## S3 method for class 'garma_model' print(x, ...)
## S3 method for class 'garma_model' print(x, ...)
x |
(garma_model) The garma_model from which to print the values. |
... |
Other arguments. Ignored. |
(null)
data(AirPassengers) ap <- as.numeric(diff(AirPassengers, 12)) mdl <- garma(ap, order = c(9, 1, 0), k = 0, method = "CSS", include.mean = FALSE) print(mdl)
data(AirPassengers) ap <- as.numeric(diff(AirPassengers, 12)) mdl <- garma(ap, order = c(9, 1, 0), k = 0, method = "CSS", include.mean = FALSE) print(mdl)
Print a 'ggbr_factors' object.
## S3 method for class 'ggbr_factors' print(x, ...)
## S3 method for class 'ggbr_factors' print(x, ...)
x |
An object of class ggbr_factors |
... |
further parameters for print function |
null
Response Residuals from the model.
## S3 method for class 'garma_model' residuals(object, type = "response", h = 1, ...)
## S3 method for class 'garma_model' residuals(object, type = "response", h = 1, ...)
object |
The garma_model object |
type |
(chr) The type of residuals. Must be 'response'. |
h |
(int) The number of periods ahead for the residuals. Must be 1. |
... |
Other parameters. Ignored. |
(double) array of resideuals from the model.
The summary function provides a summary of a "garma_model" object, printed to the output.
## S3 method for class 'garma_model' summary(object, ...)
## S3 method for class 'garma_model' summary(object, ...)
object |
(garma_model) The garma_model from which to print the values. |
... |
Other arguments. Ignored. |
(null)
data(AirPassengers) ap <- as.numeric(diff(AirPassengers, 12)) mdl <- garma(ap, order = c(9, 1, 0), k = 0, method = "CSS", include.mean = FALSE) summary(mdl)
data(AirPassengers) ap <- as.numeric(diff(AirPassengers, 12)) mdl <- garma(ap, order = c(9, 1, 0), k = 0, method = "CSS", include.mean = FALSE) summary(mdl)
Produces diagnostic plots of the model fit. This function is copied from stats::tsdiag but modifies the fit_df for the Ljung-Box test for use with garma models.
## S3 method for class 'garma_model' tsdiag(object, gof.lag = 10, ...)
## S3 method for class 'garma_model' tsdiag(object, gof.lag = 10, ...)
object |
(garma_model) The garma_model to produce the diagnostic plots for. |
gof.lag |
(int) The number of lags to examine for the Ljung-Box white noise test. |
... |
further arguments to be passed to particular methods. |
None. Diagnostics are generated.
The stats package tsdiag function: https://stat.ethz.ch/R-manual/R-patched/library/stats/html/tsdiag.html.
data(AirPassengers) ap <- as.numeric(diff(AirPassengers, 12)) mdl <- garma(ap, order = c(9, 1, 0), k = 0, method = "CSS", include.mean = FALSE) tsdiag(mdl)
data(AirPassengers) ap <- as.numeric(diff(AirPassengers, 12)) mdl <- garma(ap, order = c(9, 1, 0), k = 0, method = "CSS", include.mean = FALSE) tsdiag(mdl)
Covariance matrix of parameters if available
## S3 method for class 'garma_model' vcov(object, ...)
## S3 method for class 'garma_model' vcov(object, ...)
object |
The garma_model object |
... |
Other parameters. Ignored. |
(double) estimated variance-covariance matrix of the parameter estimates