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 = list(0, 0, 0),
k = 1,
include.mean = (order[2] == 0),
include.drift = (order[2] > 0),
method = "Whittle",
allow_neg_d = FALSE,
opt_method = NULL,
m_trunc = 50,
min_freq = 0,
max_freq = 0.5,
control = NULL
)(num) This should be a numeric vector representing the process to estimate. A minimum length of 96 is required.
(list) This should be a list (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.
(int) This is the number of (multiplicative) Gegenbauer terms to fit. Note the 'QML' method only allows k=1.
(bool) A boolean value indicating whether a mean should be fit. Note that no mean term is fit if the series is integer differenced.
(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.
(character) This defines the estimation method for the routine. The valid values are 'CSS', 'Whittle', 'QML' and 'WLL'. The default (Whittle) 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) 'CSS' is a conditional 'sum-of-squares' technique and can be quite slow. Reference: Chung (1996). 'QML' is a Quasi-Maximum-Likelihood technique, and can also be quite slow. Reference Dissanayake (2016). (k>1 is not supported for QML) 'WLL' is a new technique which appears to work well even if the \(\epsilon_{t}\) are highly skewed and/or have heavy tails (skewed and/or lepto-kurtic). However the asymptotic theory for the WLL method is not complete and so standard errors are not available for most parameters.
(bool) A boolean value indicating if a negative value is allowed for the fractional differencing component of the Gegenbauer term is allowed. This can be set to FALSE (the default) to force the routine to find a positive value.
(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('directL', 'solnp') when k<2 and 'solnp' when k>=2. The directL algorithm is used to perform a global search for the minima, and solnp to refine the values. For some data or some models, however, other methods may work well. Supported algorithms include:
cobyla algorithm in package nloptr
directL algorithm in package nloptr
BBoptim from package BB
psoptim from package pso
hjkb from dfoptim package
nmkb from dfoptim package
solnp from Rsolnp package
gosolnp from Rsolnp package
best - this option evaluates all the above options in turn and picks the one which finds the lowest value of the objective. This can be quite time consuming to run, particularly for the 'CSS' method.
Note further that if you specify a k>1, then inequality constraints are required, and this will further limit the list of supported routines.
Used for the QML estimation method. This defines the AR-truncation point when evaluating the likelihood function. Refer to Dissanayake et. al. (2016) for details.
(num) When searching for Gegenbauer peaks, this is the minimum frequency used. Default 0. Note that when there is an AR(1) component, the peaks corresponding to the AR(1) can be higher than the Gegenbauer peaks. Setting this parameter to 0.05 or above can help.
(num) default 0.5. When searching for Gegenbauer peaks, this is the maximum frequency used.
(list) list of optimisation routine specific values.
An S3 object of class "garma_model".
The GARMA model is specified as $$\displaystyle{\phi(B)\prod_{i=1}^{k}(1-2u_{i}B+B^{2})^{d_{i}}(1-B)^{id} (X_{t}-\mu)= \theta(B) \epsilon _{t}}$$
where
\(\phi(B)\) represents the short-memory Autoregressive component of order p,
\(\theta(B)\) represents the short-memory Moving Average component of order q,
\((1-2u_{i}B+B^{2})^{d_{i}}\) represents the long-memory Gegenbauer component (there may in general be k of these),
\(id\) represents the degree of integer differencing.
\(X_{t}\) represents the observed process,
\(\epsilon_{t}\) 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.
\(B\) represents the Backshift operator, defined by \(B X_{t}=X_{t-1}\).
when k=0, then this is just a short memory model as fit by the stats "arima" function.
C Chung. A generalized fractionally integrated autoregressive moving-average process. Journal of Time Series Analysis, 17(2):111<U+2013>140, 1996.
G Dissanayake, S Peiris, and T Proietti. State space modelling of Gegenbauer processes with long memory. Computational Statistics and Data Analysis, 100:115<U+2013>130, 2016.
L Giraitis, J Hidalgo, and P Robinson. Gaussian estimation of parametric spectral density with unknown pole. The Annals of Statistics, 29(4):987<U+2013>1023, 2001.
# NOT RUN {
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))
# }
Run the code above in your browser using DataLab