Package 'ModalForecast'

Automatic selection of Parametric Modal ARIMA and SARIMA models

Description

Automatic selection of Parametric Modal ARIMA and SARIMA models

Usage

auto.modal.arima(
  y,
  d = NA,
  D = NA,
  max.p = 5,
  max.q = 5,
  max.P = 1,
  max.Q = 1,
  seasonal = TRUE,
  period = stats::frequency(y),
  ic = c("aic", "bic"),
  dist = c("normal", "t", "laplace"),
  trace = FALSE
)

Arguments

y	numeric vector or time series of observations
d	Integer, degree of differencing. If NA, it's determined automatically.
D	Integer, degree of seasonal differencing. If NA, it's determined automatically. Ignored for non-seasonal models.
max.p	Maximum AR order
max.q	Maximum MA order
max.P	Maximum seasonal AR order
max.Q	Maximum seasonal MA order
seasonal	Logical. If TRUE (default) and the seasonal period is greater than one, seasonal models are considered.
period	Seasonal period. Defaults to frequency(y).
ic	Information criterion to be used in model selection ("aic", "bic")
dist	Character string specifying the error distribution. "normal" (default) for Skew-Normal, "t" for Skewed Student-t, "laplace" for Skewed Laplace.
trace	Logical. If TRUE, prints each model and its criterion.

Details

This function performs a grid search over AR and MA orders, and over seasonal AR and MA orders when the series is seasonal, to find the Modal (S)ARIMA model with the smallest information criterion (AIC or BIC). As in forecast::auto.arima, the seasonal differencing order D is chosen with forecast::nsdiffs and then the regular order d with forecast::ndiffs on the seasonally differenced series, unless they are supplied.

Value

An object of class modal_arima.

Note

GitHub repository: https://github.com/chedgala/ModalForecast

References

Galarza, C. E., Lachos, V. H., Cabral, C. R. B., and Castro, L. M. (2017). Robust quantile regression using a generalized class of skewed distributions. Stat, 6(1), 113-130.

See Also

fit_modal_arima

Examples

library(forecast)

# Non-seasonal
fit_auto <- auto.modal.arima(log10(lynx), d = 0, max.p = 2, max.q = 2)
summary(fit_auto)

# Seasonal, with a small search grid
fit_sauto <- auto.modal.arima(log(AirPassengers), max.p = 1, max.q = 1,
                              max.P = 1, max.Q = 1)
fit_sauto

---

Diagnostic Plots for Modal ARIMA and SARIMA Models

Description

Provides visual and statistical diagnostics for the residuals of a fitted Modal (S)ARIMA model. Produces a panel with the differenced series and fitted modes, the ACF and PACF of the randomized quantile residuals (RQR), a normal QQ-plot, a histogram and Ljung-Box p-values. For seasonal models the ACF, PACF and Ljung-Box statistics cover at least two seasonal periods.

Usage

diagnostics(object, ...)

## S3 method for class 'modal_arima'
diagnostics(object, ...)

Arguments

object	An object of class modal_arima.
...	Additional arguments (unused).

Value

A list of ggplot objects (invisibly) and draws the panel.

Examples

# Non-seasonal
fit <- fit_modal_arima(log10(lynx), order = c(2, 0, 0))
diagnostics(fit)

# Seasonal airline model
fit_air <- fit_modal_arima(log(AirPassengers), order = c(0, 1, 1),
                           seasonal = list(order = c(0, 1, 1), period = 12))
diagnostics(fit_air)

---

Simulation Envelope Diagnostics for Modal ARIMA and SARIMA Models

Description

Constructs simulation envelopes for the distances $\rho_p(\epsilon_t/\sigma)$ of the modal residuals, based on their theoretical distribution under the SKD family (Galarza et al., 2017): twice the distance follows the half-normal, half-t or standard exponential distribution for the Skew-Normal, Skewed Student-t and Skewed Laplace members, respectively.

Usage

envelope(object, ...)

## S3 method for class 'modal_arima'
envelope(object, B = 100, refit = TRUE, ...)

Arguments

object	An object of class modal_arima.
...	Additional arguments (unused).
B	Number of Monte Carlo replications for envelope construction. Default is 100.
refit	Logical. If TRUE (default), re-estimate the model on each simulated series.

Details

With refit = TRUE (the default), each replication simulates a series from the fitted model, re-estimates the model and computes the distances of the re-estimated residuals, so the envelope accounts for parameter estimation. This matters most for the Skewed Laplace, whose maximum likelihood fit, like least absolute deviations, pulls several residuals to zero: without refitting, the smallest distances fall below the envelope even when the model is correct. With refit = FALSE, the envelope is built from simulated innovations at the estimated parameters, which is faster.

Value

A ggplot object (invisibly) and draws the envelope plot.

Examples

fit <- fit_modal_arima(log10(lynx), order = c(2, 0, 0))
envelope(fit, B = 10)

# Seasonal model; refit = FALSE is faster, refit = TRUE accounts for estimation
fit_air <- fit_modal_arima(log(AirPassengers), order = c(0, 1, 1),
                           seasonal = list(order = c(0, 1, 1), period = 12))
envelope(fit_air, B = 10, refit = FALSE)

---

Fit a Parametric Modal ARIMA or Seasonal ARIMA Model using the SKD Family

Description

Fits a Modal ARIMA or Modal SARIMA model, where the conditional mode follows an (optionally seasonal) ARIMA recursion and the innovations follow a member of the SKD (Skewed Distribution) family: Skew-Normal, Skewed Student-t or Skewed Laplace.

Usage

fit_modal_arima(
  y,
  order = c(1, 0, 0),
  seasonal = list(order = c(0, 0, 0), period = NA),
  dist = c("normal", "t", "laplace")
)

Arguments

y	numeric vector or time series of observations.
order	A specification of the non-seasonal part of the model: the three components (p, d, q) are the AR order, the degree of differencing, and the MA order.
seasonal	A specification of the seasonal part of the model, as in arima: a list with components order (the seasonal orders (P, D, Q)) and period (the seasonal period $s$; defaults to frequency(y)). A numeric vector of length 3 is taken as the seasonal order. The default is a non-seasonal model.
dist	Character string specifying the error distribution from the SKD family. One of "normal" (default) for the Skew-Normal, "t" for the Skewed Student-t (adds degrees-of-freedom parameter nu), or "laplace" for the Skewed Laplace distribution.

Details

Let $w_t = (1-B)^d (1-B^s)^D y_t$. The model is

$\phi(B)\Phi(B^s) w_t = c + \theta(B)\Theta(B^s)\epsilon_t,$

where $\epsilon_t$ are independent $SKD(0, \sigma, \gamma)$ errors with mode zero, so that $\mu_t = w_t - \epsilon_t$ is the conditional mode of $w_t$. The non-seasonal Modal ARIMA model is the special case $P = D = Q = 0$.

All members share the parameterization of Galarza et al. (2017),

$f(y \mid \mu, \sigma, \gamma) = \frac{4p(1-p)}{\sigma} g\{2\rho_p((y-\mu)/\sigma)\},$

with $p = 1/(1+\gamma^2)$ and $\rho_p(u) = u(p - I(u<0))$, where $g$ is the standard normal, Student-t or Laplace density. The mode $\mu$ is also the $p$-th quantile, and $\gamma > 1$ gives a heavier right tail.

Parameters are estimated by maximum likelihood with BFGS and analytical gradients, started from a conditional-sum-of-squares Gaussian (S)ARIMA fit. The Skewed Laplace log-likelihood is not differentiable when a residual is zero, so for dist = "laplace" the BFGS solution is refined with Nelder-Mead. Stationarity and invertibility of the regular and seasonal polynomials are enforced during optimization.

Value

An object of class modal_arima containing:

y

The original time series.

order

The non-seasonal order (p,d,q).

seasonal

List with the seasonal order (P,D,Q) and period.

coefficients

Named vector of estimated parameters.

vcov

Asymptotic covariance matrix of the coefficients, from the observed information (delta method for sigma, gamma and nu).

loglik

The maximized log-likelihood.

nobs

Number of observations used in the likelihood, $T - d - Ds$.

hessian

The observed information on the optimization scale: the Hessian of the negative log-likelihood, or the outer product of the per-observation scores (OPG) for the Skewed Laplace, whose log-likelihood is not twice differentiable, and whenever the Hessian is not positive definite.

information

Which estimator was used for hessian: "hessian" or "opg".

convergence

Convergence code from optim.

dist

The distribution used ("normal", "t", or "laplace").

fitted.values, residuals

Fitted modes and modal residuals on the scale of y (NA for the first $d + Ds$ observations).

Note

GitHub repository: https://github.com/chedgala/ModalForecast

References

Galarza, C. E., Lachos, V. H., Cabral, C. R. B., and Castro, L. M. (2017). Robust quantile regression using a generalized class of skewed distributions. Stat, 6(1), 113-130.

See Also

auto.modal.arima, forecast.modal_arima

Examples

library(forecast)

# Non-seasonal: Lynx data
y <- log10(lynx)
fit_n <- fit_modal_arima(y, order = c(2, 0, 0), dist = "normal")
fit_t <- fit_modal_arima(y, order = c(2, 0, 0), dist = "t")
fit_l <- fit_modal_arima(y, order = c(2, 0, 0), dist = "laplace")
c(Normal = AIC(fit_n), Student = AIC(fit_t), Laplace = AIC(fit_l))
summary(fit_n)

# Seasonal: the airline model for monthly air passengers
fit_air <- fit_modal_arima(log(AirPassengers), order = c(0, 1, 1),
                           seasonal = list(order = c(0, 1, 1), period = 12))
summary(fit_air)
fc <- forecast(fit_air, h = 24)
autoplot(fc)

---

Forecasting with Modal (S)ARIMA Models

Description

Produces modal point forecasts and prediction intervals for a fitted modal_arima object, with or without a seasonal component.

Usage

## S3 method for class 'modal_arima'
forecast(
  object,
  h = 10,
  level = c(80, 95),
  interval = c("asymptotic", "bootstrap"),
  npaths = 1000,
  point = c("joint", "marginal"),
  ...
)

Arguments

object	A modal_arima object.
h	The forecast horizon.
level	Confidence level for prediction intervals.
interval	Method for computing prediction intervals ("asymptotic" or "bootstrap").
npaths	Number of simulated paths for bootstrap intervals. Defaults to 1000.
point	Modal point forecast: "joint" (most probable trajectory, default) or "marginal" (conditional mode at each horizon).
...	Additional arguments.

Details

Because the mode is not a linear operator, two modal point forecasts are available. point = "joint" (the default) returns the joint modal trajectory: the most probable future path, obtained by iterating the model with future innovations set to zero. point = "marginal" returns, for each horizon $h$, the conditional mode of $y_{T+h}$ itself. The two coincide at $h = 1$ and for symmetric errors ($\gamma = 1$); for skewed errors and $h \ge 2$ the marginal mode equals the joint path plus the mode $\delta_h$ of the forecast error $\sum_{j=0}^{h-1}\psi_j\epsilon_{T+h-j}$, which can be substantial for integrated series.

With interval = "asymptotic", the prediction intervals use the exact distribution of the forecast error at the estimated parameters, computed by numerical convolution of the SKD densities. With interval = "bootstrap", they use empirical quantiles of npaths simulated future paths.

Value

An object of class forecast. Besides the usual components, it contains mode_path (the joint modal trajectory) and mode_shift (the correction $\delta_h$ from the joint path to the marginal mode).

Examples

library(forecast)

# Non-seasonal
fit <- fit_modal_arima(log10(lynx), order = c(2, 0, 0))
pred <- forecast(fit, h = 5, level = c(80, 95))
autoplot(pred)
accuracy(pred)

# Seasonal airline model: joint trajectory versus marginal mode
fit_air <- fit_modal_arima(log(AirPassengers), order = c(0, 1, 1),
                           seasonal = list(order = c(0, 1, 1), period = 12))
fc_joint <- forecast(fit_air, h = 24)
fc_marg <- forecast(fit_air, h = 24, point = "marginal")
round(cbind(joint = fc_joint$mean, marginal = fc_marg$mean)[c(1, 12, 24), ], 4)
fc_boot <- forecast(fit_air, h = 12, interval = "bootstrap", npaths = 200)

---

Modal Trajectory Prediction for Modal (S)ARIMA Models

Description

Returns the joint modal trajectory of the future values: the recursion is iterated with future innovations set to their mode, zero. For the one-step horizon this is the conditional mode of $y_{T+1}$; for longer horizons it is the most probable path, see forecast.modal_arima.

Usage

## S3 method for class 'modal_arima'
predict(object, n.ahead = 10, ...)

Arguments

object	A modal_arima object.
n.ahead	The forecast horizon.
...	Additional arguments (unused).

Value

A numeric vector with the modal trajectory on the scale of y.

---