Title:	Beta Autoregressive Moving Average Models
Version:	1.0.1
Description:	Fits Beta Autoregressive Moving Average (BARMA) models for time series data distributed in the standard unit interval (0, 1). The estimation is performed via the conditional maximum likelihood method using the Broyden-Fletcher-Goldfarb-Shanno (BFGS) quasi-Newton algorithm. The package includes tools for model fitting, diagnostic checking, and forecasting. Based on the work of Rocha and Cribari-Neto (2009) <doi:10.1007/s11749-008-0112-z> and the associated erratum Rocha and Cribari-Neto (2017) <doi:10.1007/s11749-017-0528-4>. The original code was developed by Fabio M. Bayer.
License:	MIT + file LICENSE
URL:	https://github.com/Everton-da-Costa/betaARMA
BugReports:	https://github.com/Everton-da-Costa/betaARMA/issues
Encoding:	UTF-8
RoxygenNote:	7.3.2
Language:	en-US
Imports:	forecast
Suggests:	knitr, xtable, here, moments, rmarkdown, tseries, lbfgs, ggplot2, foreach, zoo, dplyr, gridExtra
Depends:	R (≥ 3.5)
NeedsCompilation:	no
Packaged:	2026-03-24 13:22:53 UTC; cost
Author:	Everton da Costa [aut, cre], Francisco Cribari-Neto [ctb, ths] (Theoretical foundations), Vinicius Scher [ctb]
Maintainer:	Everton da Costa <everto.cost@gmail.com>
Repository:	CRAN
Date/Publication:	2026-03-29 14:50:09 UTC

Internal Helper to Calculate Initial Phi

Description

Calculates the initial value for the precision parameter phi based on the residuals of an initial 'lm.fit' object, following Ferrari & Cribari-Neto (2004).

This helper uses the *exact* mathematical logic from the original 'start_values.R' file to ensure identical numerical output.

Usage

._get_phi_start(fit, n_eff, linkinv, mu.eta)

Arguments

Value

A single numeric value for phi_start.

Fit Beta Autoregressive Moving Average (BARMA) Models via Maximum Likelihood

Description

Fits a Beta Autoregressive Moving Average (BARMA) model to time series data valued in (0, 1) using Maximum Likelihood Estimation (MLE). The function performs complete model estimation including parameter estimation, hypothesis testing infrastructure, and model diagnostics.

Usage

barma(y, ar = NA, ma = NA, link = "logit", xreg = NULL)

Arguments

Details

Fit Beta Autoregressive Moving Average (BARMA) Models

This function fits the BARMA(p,q) model as proposed by Rocha & Cribari-Neto (2009, with erratum 2017). It serves as the main wrapper for the optimization process, calling specialized helper functions for likelihood computation, gradient calculation, and Fisher Information Matrix estimation.

**Model Specification**: The BARMA model is defined as:

g(\mu_t) = \alpha + X_t\beta + \sum_{i=1}^{p} \varphi_i (g(y_{t-i}) - X_{t-i}\beta) + \sum_{j=1}^{q} \theta_j \epsilon_{t-j}

where y_t | F_{t-1} \sim Beta(\mu_t\phi, (1-\mu_t)\phi), g(\cdot) is the link function, and F_{t-1} is the information set at time t-1.

**Model Types** (specified via 'ar' and 'ma' arguments):

• **BARMA(p,q):** Both 'ar' and 'ma' are specified.

• **BAR(p):** Only 'ar' is specified; set 'ma = NA'.

• **BMA(q):** Only 'ma' is specified; set 'ar = NA'.

**External Regressors**: Covariates can be included via 'xreg'. The model becomes a regression BARMA, where the mean depends on current covariates and lagged responses/errors.

**Optimization**: The function uses the BFGS quasi-Newton algorithm via optim with analytic gradients for efficiency. Initial values are obtained from the 'start_values()' function.

**Implementation Notes**: - The conditional log-likelihood is computed, conditioning on the first m = \max(p,q) observations (burn-in period). - The 2017 Erratum corrections are implemented for correct handling of moving average components in the score vector and Fisher Information Matrix. - All computations are vectorized where possible for efficiency.

Value

An object of class '"barma"' containing:

Note

The original version of this function was developed by Fabio M. Bayer (Federal University of Santa Maria, bayer@ufsm.br). It has been substantially modified and improved by Everton da Costa, with suggestions and contributions from Francisco Cribari-Neto.

Author(s)

Everton da Costa (Federal University of Pernambuco, everto.cost@gmail.com); Francisco Cribari-Neto (Federal University of Pernambuco, francisco.cribari@ufpe.br)

References

Rocha, A.V., & Cribari-Neto, F. (2009). Beta autoregressive moving average models. TEST, 18(3), 529-545. doi:10.1007/s11749-008-0112-z

Rocha, A.V., & Cribari-Neto, F. (2017). Erratum to: Beta autoregressive moving average models. TEST, 26, 451-459. doi:10.1007/s11749-017-0528-4

See Also

simu_barma for simulation, loglik_barma for likelihood computation, score_vector_barma for gradient computation, fim_barma for Fisher Information Matrix

Examples

  # Example 1: Fit a BAR(1) model
  set.seed(2025)
  y_sim_bar <- simu_barma(
    n = 250,
    alpha = 0.0,
    varphi = 0.6,
    phi = 25.0,
    link = "logit",
    freq = 12
  )

  # Fit the model
  fit_bar <- barma(y_sim_bar, ar = 1, link = "logit")

  # View results
  summary(fit_bar)
  coef(fit_bar)

  # Example 2: Fit a BARMA(1,1) model
  set.seed(2025)
  y_sim_barma <- simu_barma(
    n = 250,
    alpha = 0.0,
    varphi = 0.6,
    theta = 0.3,
    phi = 25.0,
    link = "logit",
    freq = 12
  )

  # Fit ARMA structure
  fit_barma <- barma(y_sim_barma, ar = 1, ma = 1, link = "logit")
  summary(fit_barma)

  # Example 3: BARMA(1,1) model with harmonic seasonal regressors
  hs <- sin(2 * pi * seq_along(y_sim_barma) / 12)
  hc <- cos(2 * pi * seq_along(y_sim_barma) / 12)

  # Create regressor matrix
  X <- cbind(hs = hs,
             hc = hc)

  fit_barma_xreg <- barma(
    y_sim_barma,
    ar = 1, ma = 1,
    link = "logit", xreg = X
  )
  summary(fit_barma_xreg)

Extract Coefficients from a barma Model

Description

S3 method for extracting the vector of estimated coefficients from a fitted model object of class '"barma"'.

Usage

## S3 method for class 'barma'
coef(object, ...)

Arguments

Value

A named numeric vector of all estimated coefficients (e.g., alpha, varphi, theta, phi).

Compute Fisher Information Matrix for Beta Autoregressive Moving Average Models

Description

Computes the observed Fisher Information Matrix (FIM) of a Beta Autoregressive Moving Average (BARMA) model. This function also efficiently returns auxiliary values like fitted values and residuals.

This function is designed for users who:

• Compute standard errors of parameter estimates

• Construct confidence intervals

• Perform hypothesis tests

• Verify theoretical properties

• Conduct simulation studies

Usage

fim_barma(y, ar, ma, alpha, varphi, theta, phi, link, xreg = NULL, beta = NULL)

Arguments

Details

Fisher Information Matrix for BARMA Models

The Fisher Information Matrix is computed from the outer product of score vectors at the MLE, which provides the observed FIM. The FIM is used to obtain standard errors via sqrt(diag(solve(FIM))).

**Non-Diagonal Structure**: Unlike generalized linear models, the FIM is not block-diagonal due to the coupling introduced by the ARMA dynamics. This is documented in the Rocha & Cribari-Neto (2009) paper.

**Important**: This function implements the corrections from the 2017 Erratum (Rocha & Cribari-Neto, 2017) for moving average components. See References section for details.

**Parameter Order**: Parameters should be supplied in the order: alpha, varphi (AR), theta (MA), phi, beta (regressors). This matches the parameter order used by barma.

Value

A list containing:

References

Rocha, A.V., & Cribari-Neto, F. (2009). Beta autoregressive moving average models. TEST, 18(3), 529-545. doi:10.1007/s11749-008-0112-z

Rocha, A.V., & Cribari-Neto, F. (2017). Erratum to: Beta autoregressive moving average models. TEST, 26, 451-459. doi:10.1007/s11749-017-0528-4

See Also

barma for model fitting, loglik_barma for log-likelihood computation, score_vector_barma for score vector (gradient)

Examples

  # Example 1: Fisher Information Matrix for a BAR(1) model
  set.seed(2025)
  y_sim_bar <- simu_barma(
    n = 250,
    alpha = 0.0,
    varphi = 0.6,
    phi = 25.0,
    link = "logit",
    freq = 12
  )

  result_bar <- fim_barma(
    y = y_sim_bar,
    ar = 1,
    ma = NA,
    alpha = 0.0,
    varphi = 0.6,
    theta = numeric(0),
    phi = 25.0,
    link = "logit"
  )

  # Check positive definiteness
  fim <- result_bar$fisher_info_mat
  all(eigen(fim)$values > 0)

  # Standard errors from inverse of FIM
  sqrt(diag(solve(fim)))

  # Example 2: Fisher Information Matrix for a BARMA(1,1) model
  set.seed(2025)
  y_sim_barma <- simu_barma(
    n = 250,
    alpha = 0.0,
    varphi = 0.6,
    theta = 0.3,
    phi = 25.0,
    link = "logit",
    freq = 12
  )

  result_barma <- fim_barma(
    y = y_sim_barma,
    ar = 1,
    ma = 1,
    alpha = 0.0,
    varphi = 0.6,
    theta = 0.3,
    phi = 25.0,
    link = "logit"
  )

  # Standard errors
  sqrt(diag(solve(result_barma$fisher_info_mat)))

Extract Fitted Values from a barma Model

Description

S3 method for extracting the fitted mean values (mu-hat) from a fitted model object of class '"barma"'.

Usage

## S3 method for class 'barma'
fitted(object, ...)

Arguments

Details

The fitted values are returned as a time series ('ts') object, matching the time properties of the original input 'y'. The first 'max_lag' observations are 'NA', as they cannot be fitted.

Value

A 'ts' object of the fitted mean values.

Forecast a barma Model

Description

S3 method for producing forecasts from a fitted 'barma' model object.

Usage

## S3 method for class 'barma'
forecast(object, h = 6, xreg = NULL, ...)

Arguments

Details

This function computes dynamic, multi-step-ahead point forecasts. It implements a "Regression with ARMA errors" logic, where the AR components are applied to the deviations from the regression line: AR(g(y_{t-k}) - x_{t-k}^\top \beta).

Value

A 'ts' object containing the point forecasts for h steps ahead.

Compute Conditional Log-Likelihood for Beta Autoregressive Moving Average Models

Description

Computes the conditional log-likelihood of a Beta Autoregressive Moving Average (BARMA) model. This function is designed for users who:

• Implement custom optimization algorithms

• Verify theoretical properties

• Conduct simulation studies

• Debug model fitting

• Integrate BARMA models into their own workflows

Usage

loglik_barma(
  y,
  ar,
  ma,
  alpha,
  varphi,
  theta,
  phi,
  link,
  xreg = NULL,
  beta = NULL
)

Arguments

Details

Log-Likelihood for BARMA Models

The log-likelihood is computed as:

\ell = \sum_{t=m+1}^{n} \log f(y_t | \mu_t, \phi)

where f is the Beta density with shape parameters shape1 = \mu_t * \phi and shape2 = (1-\mu_t)*\phi.

The linear predictor is constructed as:

\eta_t = \alpha + X_t \beta + \sum_{i=1}^{p} \varphi_i (y_{t-i} - X_{t-i}\beta) + \sum_{j=1}^{q} \theta_j \epsilon_{t-j}

**Important**: This function implements the corrections from the 2017 Erratum (Rocha & Cribari-Neto, 2017) for moving average components. See References section for details.

**Parameter Order**: Parameters should be supplied in the order: alpha, varphi (AR), theta (MA), phi, beta (regressors). This matches the parameter order used by barma.

Value

A numeric scalar representing the conditional log-likelihood value. Returns -Inf if:

• phi is non-positive or non-finite

• Insufficient observations for specified lag structure

• Fitted values are outside (0, 1)

• Any numerical issues occur

References

Rocha, A.V., & Cribari-Neto, F. (2009). Beta autoregressive moving average models. TEST, 18(3), 529-545. doi:10.1007/s11749-008-0112-z

Rocha, A.V., & Cribari-Neto, F. (2017). Erratum to: Beta autoregressive moving average models. TEST, 26, 451-459. doi:10.1007/s11749-017-0528-4

See Also

barma for model fitting, score_vector_barma for gradient computation, fim_barma for Fisher Information Matrix

Examples

  # Example 1: Log-likelihood for a BAR(1) model
  set.seed(2025)
  y_sim_bar <- simu_barma(
    n = 250,
    alpha = 0.0,
    varphi = 0.6,
    phi = 25.0,
    link = "logit",
    freq = 12
  )

  loglik_barma(
    y = y_sim_bar,
    ar = 1,
    ma = NA,
    alpha = 0.0,
    varphi = 0.6,
    theta = numeric(0),
    phi = 25.0,
    link = "logit"
  )

  # Example 2: Log-likelihood for a BARMA(1,1) model
  set.seed(2025)
  y_sim_barma <- simu_barma(
    n = 250,
    alpha = 0.0,
    varphi = 0.6,
    theta = 0.3,
    phi = 25.0,
    link = "logit",
    freq = 12
  )

  loglik_barma(
    y = y_sim_barma,
    ar = 1,
    ma = 1,
    alpha = 0.0,
    varphi = 0.6,
    theta = 0.3,
    phi = 25.0,
    link = "logit"
  )

Create Link Function Structure for BARMA Models

Description

A helper function that constructs a list containing the link function, its inverse, and the derivative of the mean function. It extends the standard 'make.link' by adding support for the "loglog" link.

Usage

make_link_structure(link = link)

Arguments

Details

This function is primarily used internally by the 'barma()' function to process the 'link' argument. It standardizes the way link functions are handled within the model fitting process.

For the "logit", "probit", and "cloglog" links, the function acts as a wrapper around the base R 'make.link'.

For the "loglog" link, which is not available in 'make.link', the necessary components are defined explicitly:

• Link function: g(\mu) = -\log(-\log(\mu))

• Inverse link function: g^{-1}(\eta) = \exp(-\exp(-\eta))

• Derivative \frac{d\mu}{d\eta}: \exp(-\exp(-\eta)) \times \exp(-\eta)

If an unsupported link is provided, the function will stop and return an error message listing the available options.

Value

A list with three components:

Author(s)

Original R code by: Fabio M. Bayer (Federal University of Santa Maria, <bayer@ufsm.br>) Modified and improved by: Everton da Costa (Federal University of Pernambuco, <everto.cost@gmail.com>)

See Also

barma, make.link

Examples

# --- Create a logit link structure ---
logit_link <- make_link_structure(link = "logit")

# Apply the link function
mu <- 0.5
eta <- logit_link$linkfun(mu)
print(eta) # Should be 0

# Apply the inverse link function
mu_restored <- logit_link$linkinv(eta)
print(mu_restored) # Should be 0.5

# --- Create a loglog link structure ---
loglog_link <- make_link_structure(link = "loglog")

# Apply the loglog link function
mu_loglog <- 0.8
eta_loglog <- loglog_link$linkfun(mu_loglog)
print(eta_loglog)

# Apply the inverse loglog link function
mu_restored_loglog <- loglog_link$linkinv(eta_loglog)
print(mu_restored_loglog) # Should be ~0.8

Print Method for a barma Model

Description

S3 method for printing a summary of a fitted '"barma"' model object.

Usage

## S3 method for class 'barma'
print(x, ...)

Arguments

Details

This function provides a concise summary, showing the function call that generated the model, the link function used, and the final estimated coefficients.

Value

Invisibly returns the original object 'x'.

Print Method for a barma Summary

Description

S3 method for printing the detailed summary of a '"barma"' model object.

Usage

## S3 method for class 'summary.barma'
print(x, ...)

Arguments

Details

This function formats and displays the summary list created by 'summary.barma()', including the call, coefficient table, and information criteria.

Value

Invisibly returns the original object 'x'.

Calculate Residuals for a barma Model Object

Description

Computes various types of residuals for a fitted Beta Autoregressive Moving Average (BARMA) model object of class '"barma"'.

Usage

## S3 method for class 'barma'
residuals(object, type = "standardized", ...)

Arguments

Details

This function is an S3 method for the generic 'residuals' function, tailored for objects returned by the 'barma()' function.

By default ('type = "standardized"'), it calculates standardized residuals on the predictor scale, defined as:

r_t = \frac{g(y_t) - \hat{\eta}_t}{\sqrt{Var(g(y_t) - \hat{\eta}_t)}}

where Var(g(y_t) - \hat{\eta}_t) \approx (g'(\hat{\mu}_t))^2 \cdot

\frac{\hat{\mu}_t(1-\hat{\mu}_t)}{1+\hat{\phi}}

. These residuals are useful for diagnostic checking, as they are approximately standard normal if the model is correctly specified.

Value

A numeric vector or 'ts' object containing the requested residuals. For standardized residuals, the first 'max_lag' values will be 'NA'.

Score Vector for the BARMA Model

Description

Computes the score vector (gradient of the log-likelihood) for the Beta Autoregressive Moving Average (BARMA) model at a given parameter vector. Used internally by barma during optimization via the BFGS algorithm.

Usage

score_vector_barma(
  y,
  ar,
  ma,
  alpha,
  varphi,
  theta,
  phi,
  link,
  xreg = NULL,
  beta = NULL
)

Arguments

Value

A numeric vector of the same length as the parameter vector, giving the partial derivatives of the log-likelihood with respect to each parameter.

See Also

barma, loglik_barma, fim_barma

Simulate a Beta Autoregressive Moving Average (BARMA) Time Series

Description

Generates a random time series from a Beta Autoregressive Moving Average (BARMA) model. The function can simulate BARMA(p,q), BAR(p), or BMA(q) processes.

Usage

simu_barma(
  n,
  alpha = 0,
  varphi = NA,
  theta = NA,
  phi = 20,
  freq = 12,
  link = "logit"
)

Arguments

Details

The model type is determined by the 'varphi' and 'theta' parameters:

• **BARMA(p,q):** Both 'varphi' and 'theta' are provided.

• **BAR(p):** Only 'varphi' is provided ('theta' is 'NA').

• **BMA(q):** Only 'theta' is provided ('varphi' is 'NA').

Value

A 'ts' object of length 'n' containing the simulated Beta-distributed time series.

Author(s)

Original R code by: Fabio M. Bayer (Federal University of Santa Maria, <bayer@ufsm.br>) Modified and improved by: Everton da Costa (Federal University of Pernambuco, <everto.cost@gmail.com>)

See Also

barma, make_link_structure

Examples

# Set seed for reproducibility
# --- Example 1: Simulate a BAR(1) process ---
# y_t depends on y_{t-1}
set.seed(2025)
bar1_series <- simu_barma(n = 250, alpha = 0.0, varphi = 0.5, phi = 20,
link = "logit")
plot(bar1_series, main = "Simulated BAR(1) Process", ylab = "Value")

# --- Example 2: Simulate a BMA(1) process ---
# y_t depends on the previous error term
set.seed(2025)
bma1_series <- simu_barma(n = 250, alpha = 0.0, theta = -0.2, phi = 20)
plot(bma1_series, main = "Simulated BMA(1) Process", ylab = "Value")

# --- Example 3: Simulate a BARMA(2,1) process with a cloglog link ---
set.seed(2025)
barma21_series <- simu_barma(
  n = 200,
  alpha = 0.0,
  varphi = c(0.4, 0.2), # AR(2) components
  theta = -0.3,         # MA(1) component
  phi = 20,
  link = "cloglog"
)
plot(barma21_series, main = "Simulated BARMA(2,1) Process", ylab = "Value")

Generate Initial Values for BARMA Model Estimation

Description

This function calculates reasonable starting values for the parameters of various Beta Autoregressive Moving Average (BARMA) models. The method is based on the approach proposed by Ferrari & Cribari-Neto (2004) for beta regression, adapted here for the time series context.

Usage

start_values(y, link, ar = NA, ma = NA, X = NA)

Arguments

Details

The function computes initial values by fitting a linear model ('lm.fit') to the link-transformed response variable 'g(y)'. This provides a computationally cheap and stable way to initialize the main optimization algorithm.

The specific procedure depends on the model's structure:

• For models with autoregressive (AR) terms, the lagged values of the transformed series 'g(y)' are used as predictors to estimate the initial 'alpha' (intercept) and 'varphi' (AR) coefficients.

• If exogenous variables 'X' are included, they are added as predictors in the linear model to obtain initial 'beta' values.

• Moving average (MA) coefficients ('theta') are initialized to zero, a standard practice in ARMA model estimation.

• The initial value for the precision parameter 'phi' is derived from the variance of the residuals of the initial linear fit, following the methodology from Ferrari & Cribari-Neto (2004).

• For a pure BMA model (no AR or X components), a simpler method is used where 'alpha' is the mean of 'g(y)' and 'phi' is based on the unconditional variance of 'y'.

This function is called internally by the main model fitting function but is exported for standalone use and inspection.

Value

A named numeric vector containing the initial values for the model parameters ('alpha', 'varphi', 'theta', 'phi', 'beta'), ready to be used by an optimization routine. Returns 'NULL' if the model specification is not recognized.

Author(s)

Original code by Fabio M. Bayer (bayer@ufsm.br). Substantially modified and improved by Everton da Costa (everto.cost@gmail.com).

References

Ferrari, S. L. P., & Cribari-Neto, F. (2004). Beta regression for modelling rates and proportions. *Journal of Applied Statistics*, 31(7), 799-818. <doi:10.1080/0266476042000214501>

Summarize a barma Model Fit

Description

(Full documentation block...)

Usage

## S3 method for class 'barma'
summary(object, ...)

Arguments

Value

A list object of class '"summary.barma"'...