3 A Case of Numerical Instability
To motivate the penalized estimator, we focus on a specific subsample of the Brasília dataset. In Cribari-Neto, Costa, and Fonseca (2025), this subsample was used to illustrate convergence failure under standard CMLE, and it was shown that ridge penalization successfully restores convergence.
While that paper established the numerical stability of the penalized estimator, the objective here is different: to verify the statistical adequacy of the penalized fit through residual diagnostics, and to show how this workflow can guide model reduction. To that end, we specify a \(\beta\text{ARMA}\) model with \(\text{AR} = (10, 18)\), \(\text{MA} = (1, 11, 13)\), and harmonic regressors.
3.1 Train–Test Split
y_subsample_ts <- window(
brasilia_ts,
start = c(2006, 02),
end = c(2019, 07)
)
y_train <- window(
y_subsample_ts,
end = c(2018, 07)
)
y_test <- window(
y_subsample_ts,
start = c(2018, 08)
)The figure below shows the full time series with the training period delimited by dashed vertical lines.
Relative humidity in Brasília. Dashed red lines delimit the training period used in the instability analysis.
3.2 Descriptive Statistics
| Min | Max | Median | Mean | SD | Skewness | Excess_kurtosis |
|---|---|---|---|---|---|---|
| 0.33 | 0.88 | 0.75 | 0.69 | 0.16 | -0.72 | 2.11 |
3.3 Feature Engineering
To capture annual seasonality we include two harmonic regressors, \(s_t = \sin(2\pi t / 12)\) and \(c_t = \cos(2\pi t / 12)\), which together represent a smooth annual wave.
freq <- frequency(y_train)
n_test <- length(y_test)
trend_index_train <- seq_along(y_train)
trend_index_test <- (max(trend_index_train) + 1):
(max(trend_index_train) + n_test)
hs_train <- sin(2 * pi * trend_index_train / freq)
hc_train <- cos(2 * pi * trend_index_train / freq)
hs_test <- sin(2 * pi * trend_index_test / freq)
hc_test <- cos(2 * pi * trend_index_test / freq)
X_train <- cbind(hs = hs_train, hc = hc_train)
X_test <- cbind(hs = hs_test, hc = hc_test)3.4 Standard CMLE
ar_lags <- c(10, 18)
ma_lags <- c(1, 11, 13)
fit_cmle <- barma(
y = y_train,
ar = ar_lags,
ma = ma_lags,
penalty = FALSE,
xreg = X_train
)summary(fit_cmle)
#> Beta Autoregressive Moving Average Model
#>
#> Call:
#> barma(y = y_train, ar = ar_lags, ma = ma_lags, xreg = X_train,
#> penalty = FALSE)
#>
#> Link function: logit
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(>|z|)
#> alpha 0.72525 0.11752 6.171 6.77e-10 ***
#> varphi10 -0.09790 0.06146 -1.593 0.111167
#> varphi18 0.28234 0.08558 3.299 0.000969 ***
#> theta1 0.92386 0.04349 21.242 < 2e-16 ***
#> theta11 0.16907 0.05591 3.024 0.002496 **
#> theta13 0.15547 0.05554 2.799 0.005120 **
#> hs 0.51653 0.06040 8.552 < 2e-16 ***
#> hc 0.71515 0.05372 13.313 < 2e-16 ***
#> phi 53.56433 6.54492 8.184 2.22e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Warning: Optimization algorithm [BFGS (stats::optim)] did not converge (Code: 1).
#> ---
#> Log-likelihood: 186.2
#> AIC: -354.3 BIC: -327.2 HQ: -343.3The optimizer failed to converge (conv = 1), and the
estimates are unreliable. Without a penalized alternative, an analyst
would have little recourse but to abandon this lag structure and search
for a different specification.
3.5 Penalized CMLE (PCMLE)
Applying ridge penalization to the same specification restores convergence.
summary(fit_pcmle)
#> Beta Autoregressive Moving Average Model
#>
#> Call:
#> barma(y = y_train, ar = ar_lags, ma = ma_lags, xreg = X_train,
#> penalty = TRUE)
#>
#> Link function: logit
#> Penalization : Ridge (lambda = 0.01234)
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(>|z|)
#> alpha 0.76929 0.10587 7.266 3.69e-13 ***
#> varphi10 -0.18060 0.07521 -2.401 0.01634 *
#> varphi18 0.23408 0.08074 2.899 0.00374 **
#> theta1 0.47711 0.07387 6.459 1.06e-10 ***
#> theta11 -0.05077 0.07329 -0.693 0.48843
#> theta13 0.31205 0.07846 3.977 6.97e-05 ***
#> hs 0.53171 0.04495 11.828 < 2e-16 ***
#> hc 0.77621 0.04482 17.317 < 2e-16 ***
#> phi 59.57003 7.28618 8.176 2.22e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Optimization: BFGS (stats::optim) converged after 111 objective function evaluations.
#> ---
#> Log-likelihood: 191.8
#> AIC: -365.7 BIC: -338.6 HQ: -354.7The penalized estimator successfully converges. The ridge penalty value \(\lambda_n = 1/(n - a)^{0.9}\) is chosen to balance bias and variance, as established in Cribari-Neto, Costa, and Fonseca (2025).
This illustrates a practical use case for PCMLE. When standard CMLE presents a convergence failure, an analyst might discard the entire lag structure without further investigation. By applying PCMLE, convergence is achieved and the parameter estimates become available for inspection. Furthermore, the diagnostic plots below for the PCMLE fit confirm that the residuals exhibit no serial correlation (Scher, Cribari-Neto, Pumi, and Bayer, 2020), indicating that the lag structure adequately captures the dynamics of the series.
Diagnostic plots for the PCMLE fit: AR=(10,18), MA=(1,11,13).
The observed and fitted values track each other closely. All ACF and PACF spikes lie within the confidence bands, and both the Ljung-Box and Monti p-values remain consistently above 0.05 across all evaluated lags, providing no evidence of residual serial correlation.
Upon reviewing the PCMLE summary, we note that the \(\theta_{11}\) estimate is not statistically significant (p-value = 0.488). This suggests that the MA(11) lag does not contribute meaningfully to the model dynamics. Removing it yields the reduced specification \(\text{AR} = (10, 18)\), \(\text{MA} = (1, 13)\), for which the likelihood surface is no longer flat and standard CMLE converges. Without PCMLE, the non-significance of \(\theta_{11}\) would never have been discovered, and the valid underlying structure of the series would have remained inaccessible.
3.6 Reduced Model
ar_lags_reduced <- c(10, 18)
ma_lags_reduced <- c(1, 13)
fit_cmle_reduced <- barma(
y = y_train,
ar = ar_lags_reduced,
ma = ma_lags_reduced,
penalty = FALSE,
xreg = X_train
)summary(fit_cmle_reduced)
#> Beta Autoregressive Moving Average Model
#>
#> Call:
#> barma(y = y_train, ar = ar_lags_reduced, ma = ma_lags_reduced,
#> xreg = X_train, penalty = FALSE)
#>
#> Link function: logit
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(>|z|)
#> alpha 0.79153 0.10931 7.241 4.45e-13 ***
#> varphi10 -0.18959 0.07669 -2.472 0.013435 *
#> varphi18 0.22200 0.08253 2.690 0.007147 **
#> theta1 0.47857 0.07509 6.373 1.85e-10 ***
#> theta13 0.29543 0.07676 3.849 0.000119 ***
#> hs 0.53129 0.04584 11.591 < 2e-16 ***
#> hc 0.77560 0.04571 16.969 < 2e-16 ***
#> phi 59.35648 7.26010 8.176 2.22e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Optimization: BFGS (stats::optim) converged after 97 objective function evaluations.
#> ---
#> Log-likelihood: 193.2
#> AIC: -370.4 BIC: -346.3 HQ: -360.6The reduced CMLE model converges and all coefficients are statistically significant. We proceed with residual diagnostics and forecasting using this parsimonious specification.
3.7 Residual Diagnostics
The panels below assess whether the fitted model adequately captures the dynamics of the Brasília relative humidity series.
Diagnostic plots for the reduced CMLE model AR=(10,18), MA=(1,13).
The observed and fitted values track each other closely throughout the sample period. All ACF and PACF spikes lie within the confidence bands, and both the Ljung-Box and Monti p-values remain consistently above 0.05 across all evaluated lags, providing no evidence of residual serial correlation.
Normal Q-Q plot of quantile residuals for the reduced CMLE model.
The quantile residuals fall closely along the reference line within the 95% confidence band, supporting the normality assumption. One observation in the upper tail departs slightly from the line but does not indicate a systematic departure from the assumed distribution.
The largest residual corresponds to October 2011 (quantile residual = 4.37), a month that marks the climatological transition from the dry to the rainy season in Brasília — the period of highest interannual humidity variability. This observation coincides with the onset of an extreme drought period documented across Brazil from 2011 onward (Cunha et al., 2019). The second largest residual corresponds to January 2015 (−3.44), coinciding with the strong 2015–2016 El Niño event, which suppressed rainfall across central Brazil. These isolated departures reflect genuine meteorological anomalies rather than systematic model misspecification.
