Autoregressive moving average (ARMA) modeling of time series has been successfully applied in different fields such as economics, finance, and engineering. Bayesian analysis of ARMA models is difficult due to the nonlinearity in the coefficients of moving average (MA) part, which complicates the likelihood function and makes the posterior density analytically intractable (Amin, 2019a). Accordingly, the exact posterior analysis of ARMA models requires the use of numerical integration, which is computationally expensive. In order to overcome this problem, three well-known analytic approximations were proposed by Newbold (1973), Zellner and Reynolds (1978), and Broemeling and Shaarawy (1984, 1988) to approximate the posterior density to be an analytically tractable closed-form distribution.

Newbold (1973) proposed the first analytic approximation by expanding unobserved errors as a linear function in the model coefficients using the first order Taylor’s expansion. In addition, Zellner and Reynolds (1978) proposed the second analytic approximation, denoted by Zellner-Reynolds, by expanding the errors sum of squares, rather than the errors, of ARMA model as a quadratic function in the model cofficients using the second order Taylor’s expansion. Moreover, Broemeling and Shaarawy (1984, 1988) proposed the third analytic approximation, denoted by Broemeling-Shaarawy, by simply replacing unobserved lagged errors with their estimates to linearize the errors as a function in the model coefficients.

Because of the simplicity of Broemeling-Shaarawy analytic approximation compared to other proposed approximations, it has been widely extended to the Bayesian analysis of many time series models such as seasonal ARMA models (Amin, 2009; Ismail and Amin, 2014), multivariate MA models (Shaarawy and Ali, 2012), double seasonal MA models (Amin, 2017b), and double seasonal ARMA models (Amin, 2017a, 2018). However, few work have been introduced to evaluate and compare the accuracy of these three analytic approximations. Ismail (1994) used a small scale of simulation study to investigate the accuracy of the three approximations in the case of MA model of order one, and Ali (1998) and Soliman (1999) extended his work to the MA model of order two and to ARMA and seasonal ARMA of order one, respectively. Therefore, there is a need in the Bayesian time series analysis to study how these analytic approximations are mathematically related and to comprehensively evaluate their accuracy, aiming to help researchers to make a trade-off between the simplicity and accuracy of these approximations. As a motivation to this work, the outcome of these analytic approximations for complicated time series models can be used as inputs to advanced methods such as Markov chain Monte Carlo (MCMC) methods (Amin and Ismail, 2015; Amin, 2020, 2022). For example, in the case of seasonal ARMA models where the posterior density is nonlinear function in the model coefficients and thus analytically intractable. So, these analytic approximations can be used to approximate the posterior density, which in turn is used to derive the full conditional posteriors as a main requirement to apply MCMC methods to introduce the Bayesian estimation and prediction for the underlying seasonal ARMA models.

In order to fill this gap, we first review these three analytic approximations and then we use the Kullback-Leibler (KL) divergence (Kullback and Leibler, 1951) to study mathematically the relation between them and to measure the distance between their derived approximate posteriors for ARMA models. In addition, we evaluate the impact of the approximate posteriors distance in the Bayesian estimates of the model coefficients and their precision by generating a large number of Monte Carlo simulations from the approximate posteriors. Therefore, we can summarize our work in this paper as follows. First, we derive the KL divergence between Newbold, Zellner-Reynolds, and Broemeling-Shaarawy approximate posteriors and compute its calibration to study mathematically the relation between them and to measure their distances. Second, we use a large number of Monte Carlo simulations from the approximate posteriors to evaluate the impact of the posteriors distance in the Bayesian estimates of mean and precision of the model coefficients. Finally, we use real-world time series datasets to illustrate the use of the KL divergence to measure the distance between the approximate posteriors and show the impact of this distance in the Bayesian estimates.

The remainder of this paper is organized as follows. In Section 2 we present the background of the ARMA models and related Bayesian concepts, and in Section 3 we summarize the analytic approximations. In Section 4 we derive the KL divergence between the approximate posteriors and its calibration. In Section 5 we present the simulation study and real-world time series datasets to measure the distance between the approximate posteriors and evaluate the impact of this distance in the Bayesian estimates. Finally, we give the conclusions in Section 6.

Time series {y_t} can be modeled by an autoregressive moving average (ARMA) model of order (p, q), simply denoted by ARMA(p, q), and written as (Box et al., 2015):

where φ_p(B) = (1 – φ₁B – φ₂B² – · · · – φ_pB^p) and θ_q(B) = (1 + θ₁B + θ₂B² + · · · + θ_qB^q) are the autoregressive and moving average polynomials with orders p and q respectively, and B is a backshift operator defined as B^kz_t = z_t–k. {ɛ_t} is a sequence of independent normal errors with zero mean and variance σ² < ∞. The ARMA model (2.1) is stationary if all the roots of φ_p(B) = 0 lie outside the unit circle, and if all the roots of θ_q(B) = 0 lie outside the unit circle this model is invertibile. For more details about the stationarity and invertibility properties of time series models see (Box et al., 2015). The model (2.1) can be simplified as:

where X_t = (y_t–1, . . ., y_t–p, ɛ_t–1, . . ., ɛ_t–q) is the t^th row of the design matrix X of order n × m, i.e. m = p + q, and β = (φ₁, . . ., φ_p, θ₁, . . ., θ_q)^T is the model coefficients.

For ARMA model with normally distributed errors, the likelihood function is given by

where Q_n = ɛ^Tɛ is the errors sum of squares, and the natural conjugate prior is normal-gamma distribution. So, suppose β ~ N_m(μ_β, σ²∑_β) and σ² ~ IG(ν/2, λ/2), the joint natural conjugate prior distribution of β and σ² can be written as:

where μ_β, ∑_β, ν and λ are hyperparameters need to be estimated.

In case of little or no information is available about the model parameters, Jeffreys’ prior can be employed, and it is given as

Multiplying the likelihood function (2.3) by each one of these two joint prior distributions results in the following joint posteriors. For the natural conjugate prior, the joint posterior of β and σ² is:

For Jeffreys’ prior, the joint posterior of β and σ² is:

It is worth observing that the unknown lagged errors are part of the design matrices X which complicates the likelihood function and makes the joint posteriors of β and σ² analytically intractable. As a result, analytic approximations were introduced to approximate these joint posteriors as we review in the following section.

3.1. Newbold approximation

Newbold (1973) expanded the errors as a linear function in the model coefficients β around their non-linear least squares estimates (NLSE) by using the first order Taylor’s expansion. So, the approximate errors for Newbold (ɛtN) can be calculated using the following formula:

where ε̂_t are the errors estimated recursively using NLSE β̂ of the model coefficients as:

and X̂_t is the vector of the explanatory variables (y_t–1, y_t–2, . . ., y_t–p, ε̂_t–1, ε̂_t–2, . . ., ε̂_t–q) after replacing the error terms by their NLSE’s. U_t is a derivative vector that can be defined as:

Accordingly, the Newbold approximate errors sum of squares QnN can be obtained as:

where Q̂_n = ε̂^Tε̂ and ε̂_t is computed as in Equation (3.2), and U is a matrix with the t^th row is U_t defined in Equation (3.3). Therefore, the Newbold approximate likelihood function is given by

and the Newbold approximate joint posterior of β and σ², in case of natural conjugate prior, can be obtained as:

Using some mathematical manipulations, we can prove that the marginal approximate posterior of the model coefficients β is a multivariate t distribution with degrees of freedom v_n = (n + ν) and location vector and dispersion matrix are respectively:

where An-1=(UTU+Σβ-1)-1 and Bn=(UTUβ^+Σβ-1μβ). By setting λ = 0, Σβ-1=0, and ν = −m, we can obtain the Newbold approximate marginal posterior of β in the case of Jeffreys’ prior as a multivariate t distribution with degrees of freedom v_j = (n – m) and location vector and dispersion matrix are respectively:

3.2. Zellner-Reynolds approximation

Zellner and Reynolds (1978) expanded the errors sum of squares Q_n, instead of the errors, as a quadratic function in the model coefficients using the second order of Taylor’s expansion. The Zellner- Reynolds approximation of the errors sum of squares can be calculated as:

where R is an m × m matrix of the second derivatives that can be defined as:

Accordingly, the Zellner-Reynolds approximate likelihood function is given by

and by employing the natural conjugate prior we can obtain the approximate joint posterior of β and σ² as:

With some manipulations we can simply derive the marginal approximate posterior of the model coefficients β as a multivariate t distribution with degrees of freedom v_n = (n + ν) and location vector and dispersion matrix are respectively:

where Az-1=((1/2)R+Σβ-1)-1 and Bz=((1/2)Rβ^+Σβ-1μβ). By setting λ = 0, Σβ-1=0, and ν = −m, we can obtain the approximate marginal posterior of β in the case of Jeffreys’ prior as a multivariate t distribution with degrees of freedom v_j = (n – m) and location vector and dispersion matrix are respectively:

3.3. Broemeling-Shaarawy approximation

Broemeling and Shaarawy (1984, 1988) approximated the unobserved errors by simply replacing the unobserved lagged errors with their NLSE’s to be a linear function in the model coefficients as follows:

where X̂_t is t^th row of the design matrix after replacing the unobserved lagged error terms by their NLSE’s. Therefore, the Broemeling-Shaarawy approximate errors sum of squares QnBS can be obtained as:

and the approximate likelihood function is given by

and the Broemeling-Shaarawy approximate joint posterior of β and σ², in case of natural conjugate prior, can be obtained as:

Using some manipulations, we can derive the marginal approximate posterior of the model coefficients β to be a multivariate t distribution with degrees of freedom v_n = (n + ν) and location vector and dispersion matrix are respectively:

where Ab-1=(X^TX^+Σβ-1)-1 and Bb=(X^T y+Σβ-1μβ). By setting λ = 0, Σβ-1=0, and ν = −m, the approximate marginal posterior of β in the case of Jeffreys’ prior can be obtained as a multivariate t distribution with degrees of freedom v_j = (n – m) and location vector and dispersion matrix are respectively:

The Kullback-Leibler (KL) divergence (Kullback and Leibler, 1951) can be used to measure the distance between any two of the approximate posteriors of the ARMA model coefficients introduced in Section 3. To explain the idea, suppose we have two multivariate t posteriors for β: the first is ζ₁(β | y) with parameters mean μ₁, dispersion matrix V₁, and degrees of freedom v₁, and the second is ζ₂(β | y) with parameters mean μ₂, dispersion matrix V₂, and degrees of freedom v₂. The KL divergence between these two posteriors can be written as:

This KL divergence is not symmetric and always non-negative, and it equals to zero when the two posteriors are the same. This implies that the smaller its value, the closer are the two posteriors. We can simplify the KL divergence (4.1) as:

where CH[ζ₁(β | y), (ζ₂(β | y)] and H [ζ₁(β | y)] are known as the cross-entropy and Shannon entropy respectively (Amin, 2017c). In our previous work (Amin, 2019b), we derived the approximate KL divergence between the two multivariate t posteriors ζ₁(β | y) and ζ₂(β | y) as:

This derived approximate KL divergence is not symmetric; however, we can compute the symmetric distance KL^* [ζ₁(β | y), ζ₂(β | y)] as the average of two KL divergences as follows:

where V1*=((1-2/v2)/(1+m/v1))V1, V2*=((1-2/v1)/(1+m/v2))V2, and V-1=(1+m/v1)V1-1+(1+m/v2)V2-1.

Now, we can use this derived approximate KL divergence in Equation (4.4) to study the relation between the analytic approximations summarized in Section 3 and to measure the distance between their approximate posteriors. To do that we have to observe two things. First, the derived approximate KL divergence depends on the difference between the approximate posteriors means and the ratio of the posteriors dispersions. Second, all the marginal approximate posteriors of β derived in Section 3 are multivariate t distribution with the same degrees of freedom but with different location vector and dispersion matrix. Aiming to simplify the comparison, we compute the KL divergence between the marginal approximate posteriors of β in the case of employong Jeffreys’ prior.

Using the location vectors and dispersion matrices given in Equations (3.8) and (3.14), we can simplify the KL divergence between the Newbold and Zellner-Reynolds approximate posteriors of β as:

which indicates that the distance between the Newbold and Zellner-Reynolds approximate posteriors depends mainly on the relation between the two matrices U^TU and (1/2)R. To investigate this relation, we can redefine the R matrix in Equation (3.10) by letting its i j^th element to be defined as:

which means (1/2)R = D + U^TU and implies that U^TU is a part of (1/2)R. Therefore, the distance between the Newbold and Zellner-Reynolds approximate posteriors is expected to be very small, i.e. they are very close to each other, especially when the matrix D closes to 0, which is evaluated by the simulation study in the next section. Similarly, we use the location vectors and dispersion matrices given in Equations (3.8) and (3.20) to simplify the KL divergence between the Newbold and Broemeling-Shaarawy approximate posteriors of β as:

which depends on the difference between the two approximate posteriors means and the ratio of the two matrices X̂^TX̂ and U^TU. Since the matrix U contains more information about the unobserved errors compared to the matrix X̂^T, it is expected that the distance between the Newbold and Broemeling- Shaarawy approximate posteriors is very large, i.e. they are strongly different, which we evaluate by the simulation study in the next section. We can obtain similar conclusion about the distance between the Zellner-Reynolds and Broemeling-Shaarawy approximate posteriors, since their KL divergence can be simplified as:

It has to be noted that possible values of the approximate KL divergence are non-negative with no maximum limit, however, we can calibrate the values to be in the interval (0.5,1.0) to be able to decide about the distance between the approximate posteriors. When the calibration value closes to 0.5, it implies the two approximate posteriors are almost the same; and when its value closes to 1.0, the two approximate posteriors are strongly different. Suppose KL[ζ₁(β | y), ζ₂(β | y)] = k, using the idea of calibration proposed by McCulloch (1989) we can compute the calibration of KL divergence (KLC) between the two approximate posteriors as:

All the derived approximate KL divergences of the approximate posteriors and their calibration will be used in the next section to measure the distance between these posteriors and the impact of this distance in the Bayesian estimates using simulated and real-world time series datasets.

We have two parts of work in this section. First, we present a simulation study to evaluate the KL divergence (and its calibration) between the considered approximate posteriors of the ARMA model coefficients and evaluate the impact of posteriors divergence in the Bayesian estimates. Second, we present two applications of our work to real world time series datasets.

5.1. Simulation study

In this simulation study, we have two objectives: (1) evaluating the KL divergence (and its calibration) between the considered approximate posteriors of the ARMA model coefficients, and (2) evaluating the impact of this posteriors divergence in the Bayesian estimates for several simulated time series data with different sample sizes, different ARMA model orders, and different values of the ARMA model coefficients. In the following we discuss the simulation study design and results for each objective.

5.2. Application to real-world time series

We introduce two real-world time series to demonstrate how the KL divergence can be used in reality to measure the distance between the derived approximate posteriors and to show the impact of this distance in the Bayesian estimates of the model coefficients. These two real-world time series datasets represent viscosity and concentration outputs from two different chemical processes (Box et al., 2015). These datasets were collected on full-scale processes to show the effect on these outputs of uncontrolled disturbances such as variations in feedstock. These two time series datasets are presented in Figure 1.

It can observe from Figure 1 that these two time series are not stationary, so we use the first difference to stationarize them; and we then identify the best suitable order of the ARMA model for the differenced time series datasets by using the corrected Akaike’s information criterion (AIC_c) with assuming the maximum order is five for both autoregressive and moving average parts. We compute the AIC_c for all the combinations of ARMA models and select the best model with smallest AIC_c value, and accordingly we identify the best model for both the differenced datasets is ARMA(1, 1). Using the identified models for these two time series datasets and employing the Jeffreys’ priors, we first obtain the Newbold, Zellner-Reynolds and Broemeling-Shaarawy approximate posteriors of the model coefficients. Second, we compute the KL divergence and its calibration between these approximate posteriors and present their results in Table 12. Finally, we compute the Bayesian estimates for the coefficients of these models and present results in Table 13. It is worth noting that these results are consistent with our conclusions in the simulation study in previous subsection.

In this paper we first reviewed the Newbold, Zellner-Reynolds, and Broemeling-Shaarawy analytic approximations proposed in the literature to approximate the posterior density of ARMA models to be analytically tractable. We then derived the KL divergence between these approximate posteriors and compute its calibration to measure the distance between these posteriors. We used a large number of Monte Carlo simulations from the approximate posteriors to evaluate the impact of the posteriors distance in the Bayesian estimates of mean and precision of the model coefficients. Simulation results confirmed that the Newbold and Zellner-Reynolds approximate posteriors are very close to each other, and they are strongly different from the Broemeling-Shaarawy approximate posterior. In particular, the results demonstrated that the Newbold and Zellner-Reynolds approximations show higher precision of the model coefficients compared to the Broemeling-Shaarawy approximation, which is reflected in the coefficients credible intervals. We applied our work to real-world time series datasets and their results are consistent with those of the simulation study.