Package {ModalForecast}

Title:	Parametric Modal ARIMA and Seasonal ARIMA Models using the SKD Family
Version:	0.2.0
Description:	Implements parametric modal Autoregressive Integrated Moving Average (ARIMA) and seasonal ARIMA (SARIMA) models utilizing the Skewed Distribution (SKD) family, in which the conditional mode, rather than the conditional mean, follows the (seasonal) ARIMA recursion. Current distributions supported are the Skew-Normal, Skewed Student-t, and Skewed Laplace. The parameters are estimated by maximum likelihood using analytical gradients. Includes residual diagnostics, simulation envelopes, automatic order selection, joint and marginal modal forecasts, exact and parametric bootstrap prediction intervals, and classical asymptotic inference via the Fisher Information matrix. Methods are described in Galarza, C.E., Lachos, V.H., Cabral, C.R.B., & Castro, L.M. (2017) <doi:10.1002/sta4.140>.
URL:	https://github.com/chedgala/ModalForecast
BugReports:	https://github.com/chedgala/ModalForecast/issues
Depends:	R (≥ 3.5.0)
License:	GPL-3
Encoding:	UTF-8
RoxygenNote:	7.3.3
Imports:	stats, utils, graphics, forecast, ggplot2, gridExtra, scales, grid
Suggests:	rmarkdown, testthat (≥ 3.0.0), knitr
Config/testthat/edition:	3
NeedsCompilation:	no
Packaged:	2026-09-24 22:17:00 UTC; chedgala
Author:	Christian Galarza [aut, cre], Sergio Luis Mercado Londoño [ctb], Víctor Hugo Lachos [ctb]
Maintainer:	Christian Galarza <chedgala@espol.edu.ec>
Repository:	CRAN
Date/Publication:	2026-09-25 07:10:07 UTC

ModalForecast: Parametric Modal ARIMA and Seasonal ARIMA Models using the SKD Family

Description

Implements parametric modal Autoregressive Integrated Moving Average (ARIMA) and seasonal ARIMA (SARIMA) models utilizing the Skewed Distribution (SKD) family, in which the conditional mode, rather than the conditional mean, follows the (seasonal) ARIMA recursion. Current distributions supported are the Skew-Normal, Skewed Student-t, and Skewed Laplace. The parameters are estimated by maximum likelihood using analytical gradients. Includes residual diagnostics, simulation envelopes, automatic order selection, joint and marginal modal forecasts, exact and parametric bootstrap prediction intervals, and classical asymptotic inference via the Fisher Information matrix. Methods are described in Galarza, C.E., Lachos, V.H., Cabral, C.R.B., & Castro, L.M. (2017) doi:10.1002/sta4.140.

Author(s)

Maintainer: Christian Galarza chedgala@espol.edu.ec

Other contributors:

• Sergio Luis Mercado Londoño luis.mercado@udea.edu.co [contributor]

• Víctor Hugo Lachos hlachos@uconn.edu (ORCID) [contributor]

See Also

Useful links:

• https://github.com/chedgala/ModalForecast

• Report bugs at https://github.com/chedgala/ModalForecast/issues

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

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

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

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

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

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

Value

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

Objects exported from other packages

Description

These objects are imported from other packages. Follow the links below to see their documentation.

forecast

forecast