CC BY
59
Bayesian Analysis of ARMA and BSTS models for COVID-19 data using R and Stan
Muhammed Navas T and Athar Ali Khan
•
Department of Statistics and Operations Research Aligarh Muslim University, Aligarh-202002, India navasmku@gmail.com atharkhan1962@gmail.com
Abstract
This study compares the performance of Bayesian ARIMA and BSTS models for COVID-19 data using Bayesian approach. Many studies in the literature have compared the BSTS model and classical ARIMA models for infectious disease modelling, and the BSTS model performs well. Apart from the literature, this study is trying to prove the Bayesian ARIMA model gives a better result than the BSTS model. This study uses a different modelling and model comparison method to compare widely used autoregressive integrated moving average (ARIMA) models with their Bayesian structural time series (BSTS) models for COVID-19 data using the Bayesian approach. It is essential to find the order of the ARIMA model before doing bayesian analysis. We find the order of the ARIMA model using measurement LOO information criteria, using the Hamiltonian Montecarlo algorithm and rstan estimate the parameters of ARIMA and BSTS models for COVID-19 data. Furthermore, compare both models using Looic and Waic values; Bayesian ARIMA models outperform in this study.
Keywords: BSTS, COVID-19, ARIMA, MCMC, LOOIC, WAIC, Stan
1. Introduction
The literature contains a variety of traditional studies that use ARIMA models. The statistical models and techniques for evaluating discrete time series are discussed in [1], along with some of the methodology's most significant applications. The class of autoregressive integrated moving average (ARIMA) models and various extensions of these models are among the models taken into consideration. [2] used the ARIMA model, a version of ARMA, and fitted the same to non-seasonal data by identifying autoregressive and moving average terms with the help of PACF and ACF. [3] describes how to find ARIMA models using the extended autocorrelation function. Over the past few decades, the Bayesian method's significance in econometrics has grown significantly. In this sense, significant references include [4], [5], [6], [7], and [8], etc. [9] clarifies the autoregressive and moving average parameters are implicitly constrained to the stationary and invertible region using a straightforward reparameterisation. [10] offers a few new transformation concepts and looks at how they fit into an effective numerical integration strategy for ARIMA models. The ARMA model can be used to model a variety of data sets due to its universal structure. Typically, the theory does not specify which model should be chosen, so it must be chosen from among a variety of competing models. The choice of an ARIMA model is vital for both statistical inference and prediction.
The bayesian structural time series (BSTS) model, which is based on [11] technique, is another model we used in this work. Create numerous layers utilising this model, such as trends, seasonality, and regression components. In addition, unlike with traditional ARIMA models,
stationary time series are not required for BSTS models, since the BSTS model can deal with structural changes in time series. In contrast to classical ARIMA and machine learning models, a notable advantage of the BSTS model is its distinct interpretable structure in both observable and unobservable dynamic components.
2. Methodology
In many studies comparison between classical ARIMA and BSTS models for COVID-19 data, the BSTS model performs well in these studies [12, 13, 14]. Apart from the literature, this study is trying to prove the Bayesian ARIMA model gives a better result than the BSTS model for COVID-19 data. The study's goal is to provide a Bayesian approach to assess ARIMA models and the BSTS model, utilizing the Hamiltonian Monte Carlo (HMC) algorithm and the programming language Rstan. This is being done for the COVID-19 cumulative number of cases and the cumulative number of deaths from the date 01 March 2020, to 30 June 2021. Performing this study in two steps; first, for the analysis of Bayesian ARIMA models, it is important to find the order of the ARIMA model. We find the optimized order using the LOOIC value, and the model with the lowest LOOIC value is considered the best model [15]. In the second step, using various priors, we estimate the parameters of the BSTS model and the Bayesian ARIMA model. In the following phase, we compare the BSTS model and the Bayesian ARIMA model using various measures for the corresponding data, such as the looic and waic values [16]. The model with the lowest looic and waic values is considered the best model [17].
3. ARIMA model
ARMA model with AR coefficient < 's and MA coefficient 9j's is defined as
p q
yt = № + E <yt-i + E 6i €H + £t (1)
i=1 j=1
in equation 1, p and q are order of ARMA model, yt's are the data points collected over time, ^0 is the intercept, and €'s are the error terms distributed according to IID normal variate with mean zero and variance a2. Equation 1 is denoted by ARIMA (p,q) model. Invertibility enables us to estimate the noise (€t) recursively by
pq
€t = yt - V0 - E <yt-i - E 6i€t-j, t = 1,2,..., T (2)
i=1 j=1
As an approximation, we consider the starting values yt = €t = 0 for t < 0. Choosing p + q < 2 greatly simplifies the model, and it is frequently observed that higher order ARMA models are more difficult to justify in practise. (see, for example, [6]). A careful review of the literature also suggests that the problems of stationarity and invertibility can be easily tackled for small p and q, say p + q < 2.
Advantages of Bayesian ARIMA models:
• Flexibility: The Bayesian ARIMA model can handle a wide range of time series patterns, such as non-stationary, multi-seasonal, and multi-trend data.
• Better forecasting: The Bayesian approach to ARIMA models allows for more accurate predictions by taking into account uncertainties in the model parameters.
• Model selection: The Bayesian framework enables the use of model selection methods, such as LOOIC, which helps to determine the optimal number of AR and MA terms in the model.
• Prior information: The Bayesian ARIMA model allows for incorporating prior information or domain knowledge into the model, making it more robust and reliable.
• Model uncertainty: The Bayesian ARIMA model provides a quantification of model uncertainty, which can be useful in decision-making processes.
• Model comparison: The Bayesian ARIMA model allows for comparison with other models, and selection of the best model based on a range of criteria.
3.1. Bayesian Model Formulation of ARMA model
Let y : y1, y2,..., yT be the time series observations from a strictly stationary and invertible ARMA model 1. To write the approximate likelihood function for an ARMA(p, q) model with yt's are the time series observed data is
yt — N ((uo + yt_i + ... + fayt-p + Oiet-i + ... + dqet-q), o2) under the assumption of yt = et = 0 for t < 0. The conditional density of yt given yt-1, yt_2,..., yt_p can then be written as,
f (yt | yt_i, yt_2,..., yt_p) Uo, ©) «
i \ [if p q \2\ (3)
02) exp I _ 202 (^t _ Uo _ E fryt_i _ E 0et_j^
The likelihood function is defined as:
T
f (y 1 ^ ^ ©) « n f (yt 1 yt_v yt_2,..., yt_p; ^ ©) (4)
t=p+i
which reduces to:
f[l1 ^ ©) «
1 \(T_ p)/2 / 1 T / p q (5)
02) exp 1 _ 202 e i v _ uo _ E ^iyt_i _ E0jet_
4. BSTS model
The BSTS model is defined by a pair of equations,
yt = Zjttt + £ t, at+i = Tt at + Rt nt,
(6)
where £t — N (o, of) and nt — N (o, Qt), both error terms are independent of all other unknowns [11]. In this study, we are using a local linear trend model (a BSTS model without regression components) that assumes level and slope as random walk components.
4.1. Local linear trend model
The local linear trend model is a popular option for trend modelling because it responds rapidly to local variation, which is important when making short-term forecasts. The equation as follows:
yt = ut + et et - N (o,of) (7)
the equation of the level component is:
Ut+1 = Ut + ¿t + nt nt - N (0,o2) (8)
and the equation of the slope is:
¿t+1 = ¿t + Zt Zt - N (o, of) (9)
The local linear trend is based on the supposition that both the mean and slope components follow random walks.
5. Bayesian estimation using R and Stan
5.1. Bayesian estimation of ARIMA models
In this session, use the corresponding priors to express the following stan code for the ARIMA(1,0,1) model in R:
p0 ~ studentt(0,1,6) ~ normal (0,0.5) 0i ~ normal (0,0.5) c ~ studentt (0,1,7)
For all ARIMA models, we use the same priors for ц and c; for all AR components, we use prior as normal(0,0.5), and for all MA components, we use prior as normal(0,0.5). Stan code for the ARIMA(1,0,1) model is given below:
armamodel="data { int<lower=1> N;
real y[N]; }
parameters { real mu0; real phil; real thetal;
real<lower=0> sigma; }
model {
vector[N] nu;
vector[N] err;
psi[1] = mu0 + phil * mu0;
err[1] = y[1] - psi[1];
for ( n in 2:N) {
psi[n] = mu0 + phi1 * y[n-1] + theta1 * err[n-1];
err[n] = y[n] - psi[n]; }
mu0 ~ studentt(0, 1, 6); phi1 ~ normal(0, 0.5); theta1 ~ normal(0, 0.5); sigma ~ studentt(0, 1, 7);
err ~ normal(0, sigma); }
Only the ARIMA(1,0,1) model code is expressed here; other types of ARIMA models were used in this study but were not displayed owing to space limitations.
5.2. BSTS models
In this session, the Stan formulation of the Bayesian structural time series model is discussed. The appropriate priors for the parameters are discussed below:
Uerr ^^ N(0,1) verr ~ N(0,1) Cslope ~ N(0,0.5)
auvd - N(0,0.5) aobs - N(5,10)
Stan code for BSTS model is:
\label{stan:2} bstsmodel=" data {
int <lower=0> N;
vector[N] y; }
parameters { vector[N] a_err; vector[N] b_err; real beta;
real <lower=0> sigma_obs; real <lower=0> sigma_slope;
real <lower=0> sigma_level; }
transformed parameters {
vector[N] a;
vector[N] b;
a[1] = a_err[1];
b[1] = b_err[1];
for (n in 2:N) {
a[n] = a[n-1] + b[n-1] + s_level*u_err[n] ;
b[n] = b[n-1] + s_slope*v_err[n] ; }
}
model {
a_err ~ normal(0,1); b_err ~ normal(0,1); sigma_slope~normal(0,0.5); sigma_level~normal(0,0.5); sigma_obs~normal(5,10); for(n in 1:N){
y[n] ~ normal (a[n] ,sigma_obs); }
}
generated quantities{ vector[N] yrepg; real log_lik[N];
for(n in 1:N){log_lik[n]=normal_lpdf(y[n]|a[n],sigma_obs); }
}
In this chapter for the BSTS model, we apply the same priors and stan code for the remaining data sets; when we use the same priors for all data sets, we can more precisely compare our looic and waic values.
6. Results and discussions
During this session, we reviewed the estimated parameters obtained after executing Stan code for both Bayesian ARIMA and BSTS models.
6.1. ARIMA model Bayesian estimation 6.1.1 ARIMA model selection
Finding the AR component, MA component, and order of differencing for the model is the first step before using an ARIMA model. Table 1 lists various p and q order ARIMA models and their associated LOOIC values. We assess alternative order ARIMA models with LOOIC values for the cumulative number of cases and the cumulative number of deaths for each of the five countries. We evaluate the model with the lowest LOOIC value as the ideal model, and we employ that model in further bayesian analysis. For instance, from table 1, take the cumulative number of cases in the USA, the minimum LOOIC value is 10720.14 for the order of ARIMA (1, 2, 1), so we find our ARIMA model using the minimum LOOIC value.Similarly, for the cumulative number of deaths in the USA, the minimum LOOIC value is 7000.51 for the model ARIMA (3, 2, 2), so we use this model in further Bayesian analysis. Similarly, we ordered ARIMA for further Bayesian analysis for all other countries' cumulative cases and deaths.
Table 1: ARIMA model selection
Countries Cases model AIC value Deaths model AIC value
USA ARIMA(2„2,2) ARIMA(0,2,1) ARIMA(1,2,2) ARIMA(2,2,1) ARIMA(1,2,1) ARIMA(2,2,0) 10721.76 10730.94 10721.83 10722.18 10720.14 10751.45 ARIMA(2,2,2) ARIMA(1,2,2) ARIMA(2,2,1) ARIMA(3,2,2) ARIMA(3,2,1) ARIMA(3,2,3) 7007.206 7176.245 7131.75 7000.51 7115.77 7012.184
UK ARIMA(2,2,2) ARIMA(0,2,0) ARIMA(1,2,0) ARIMA(0,2,1) ARIMA(1,2,1) ARIMA(1,2,2) 8921.798 8935.118 8919.275 8917.05 8918.23 8920.127 ARIMA(2,2,2) ARIMA(1,2,2) ARIMA(2,2,1) ARIMA(5,2,2) ARIMA(5,2,1) ARIMA(4,2,1) 6181.108 6204.96 6237.94 5801.33 5968.22 6114.78
UAE ARIMA(0,2,0) ARIMA(1,2,0) ARIMA(0,2,1) ARIMA(1,2,1) ARIMA(0,2,2) 6906.30 6788.24 6729.45 6732.08 6731.07 ARIMA(2,2,2) ARIMA(0,2,1) ARIMA(1,2,1) ARIMA(2,2,1) ARIMA(1,2,2) 1987.16 1984.84 1983.08 1985 1973.254
Bahrain ARIMA(2,2,2) ARIMA(1,2,2) ARIMA(1,2,1) ARIMA(2,2,1) 5911.402 5908.39 5929.98 5931.204 ARIMA(2,2,2) ARIMA(0,2,1) ARIMA(1,2,1) ARIMA(1,2,2) 2096.992 2096.94 2099.94 2101.53
India ARIMA(2,2,2) ARIMA(2,2,1) ARIMA(3,2,2) ARIMA(3,2,1) 10000.67 10033.83 9941.79 9990.04 ARIMA(2,2,2) ARIMA(1,2,1) ARIMA(0,2,1) ARIMA(1,2,2) 6895.87 6895.35 6892.32 6895.8
Using the above mentioned stan code in 5.1 and using different priors to estimate unknown parameters shown in table 2. Figure 1a and 1b show the posterior density plot and trace plot of the ARIMA model for the cumulative number of UK cases and the cumulative number of Bahrain
deaths, respectively. Due to space restrctions, here is the convergence plot of two countries: the UK and Bahrain are only displayed.
-0.60 (b)
Figure 1: (a) ARIMA model for UK cases Posterior density plot, and Traceplot (b) ARIMA model for bahrain deaths Posterior density plot, and Traceplot
Table 2 displays the anticipated outcomes for the ARIMA models. For all countries, it is simple to verify that the majority of estimated parameters across the ARIMA models are significantly different from zero. Let's use the results for deaths in the UK as an example. Except for all other parameters are statistically significant in the ARIMA(5,2,2) model. This model's parameters have a 95% credible interval that excludes zero and demonstrates the statistical significance of the estimate. Additionally, the model's ar.1 and ma.2 components contribute favourably, but ar.2, ar.3, ar.4, and ar.5 components have statistically significant negative contributions to the ARIMA(5, 2, 2) model. For the USA's cumulative number of cases, the recommended model is ARIMA (1, 2, 2). There is one AR component, two MA components, and two other parameters, and a0. All five parameters are statistically significant, which means their 95% credible interval does not contain
zero, and the first MA component contributes to the model negatively. Similarly, we can interpret other countries' parameters using table 2.
Countries Items (Models) № 7 $1 $2 $3 $4 $5 01 02
USA Cases (ARIMA(1/2/1) Deaths (ARIMA(3,2,2) 0.0545 (0.0286) 0.0367 (0.0209) 15711.35 (4.6476) 360.2125 (0.1044) 0.2375 (0.0006) 0.7407 (0.0004) -0.4641 (0.0005) -0.2742 (0.0004) -0.6858 (0.0004) -0.9896 (0.0001) 0.3402 (0.0004)
UK Cases ARIMA(0,2,1) Deaths ARIMA(5,2,2) 0.0717 (0.0277) 0.0263 (0.0211) 2433.14 (0.7048) 95.9465 (0.0277) 0.1916 (0.0004) -0.8767 (0.0003) -0.1517 (0.0005) -0.4477 (0.0003) -0.5309 (0.0004) -0.2125 (0.0004) -0.6963 (0.0004) 0.6722 (0.0003)
Cases ARIMA(0,2,1) 1.2845 (0.0371) 252.93 (0.1261) -0.6439 (0.0005)
UAE
Deaths ARIMA(1,2,2) 0.0135 (0.0003) 1.8694 (0.0005) -0.2635 (0.0020) -0.5254 (0.0021) -0.1067 (0.0016)
Cases ARIMA(1,2,2) 0.0278 (0.0279) 108.622 (0.0564) 0.4209 (0.0027) -0.8067 (0.0026) 0.3154 (0.0010)
Bahrain Deaths ARIMA(0,2,1) 0.0054 (0.0003) 2.1048 (0.0006) -0.6267 (0.0003)
Table 2: ARIMA model Bayesian estimation: estimation of the posterior means.
6.2. BSTS model Bayesian estimation
From table 3, the majority of the posterior estimates for the BSTS model are statistically significant, which is consistent with the results from the BSTS models. For UK deaths, for instance, a posterior parameter estimate of aobs with a range of 142.04 to 155.34 implies that zero excludes from a credible interval of 95 %, indicating that the parameter is statistically significant. Similarly, for parameter aslope and &level, UK deaths have statistical significance. Similarly, for parameter aslope and &level, UK deaths have statistical significance. The parameter value of &slope is 0.02, and &ievel is 0.38 for the UK cumulative number of deaths. We can provide a similar justification for the parameter estimate for the other four countries. Figure 2 to 9 show the autocorrelation plot, caterpillar plot, posterior density plot, and trace plot of the BSTS model for the UK cumulative number of cases and deaths, respectively.From the trace plot, we can interpret the convergence of the Markov chain, due to space restrictions not displaying the plots of all countries.
Table 3: BSTS model Bayesian estimation: estimation of the posterior means.
Countries Items 7obs 7slope 7level
Cases 1948.43 3.25 0.40
USA (0.04) (0.001) (0.001)
Deaths 284.66 0.03 0.41
(0.07) (0.001) (0.02)
Cases 1239.84 11.51 0.39
UK (0.10) (0.01) (0.01)
Deaths 148.45 0.02 0.38
(0.04) (0.001) (0.001)
Cases 218.46 0.03 0.40
(0.07) (0.001) (0.01)
UAE
Deaths 2.41 0.02 0.02
(0.08) (0.001) (0.001)
Cases 107.37 0.03 0.41
(0.08) (0.001) (0.02)
Bahrain
Deaths 2.50 0.01 0.02
(0.08) (0.001) (0.001)
Cases 1240.02 11.53 0.40
India (0.10) (0.01) (0.01)
Deaths 233.30 0.04 0.39
(0.09) (0.001) (0.001)
Figure 2: BSTS model for UK cases Figure 3: BSTS model for UK cases Figure 4: BSTS model for UK cases
autocorrelation plot
caterpillar plot
posterior density plot
L
Figure 5: BSTS model for UK cases Figure 6: BSTS model for UK deaths Figure 7: BSTS model for UK deaths trace plot autocorrelation plot caterpillar plot
Figure 8: BSTS model for UK deaths posterior density plot
Figure 9: BSTS model for UK deaths trace plot
7. Conclusion
In this study, we used a Bayesian estimation to evaluate the widely used ARIMA and BSTS models while modelling COVID-19 cumulative number of cases and the cumulative number of deaths in five countries. From table 4, we discover that the ARIMA model LOOIC value is lower than that of their corresponding BSTS models. The model with the lowest LOOIC value is considered to be the best model. For example, Bahrain's cumulative number of cases LOOIC value for the Bayesian ARIMA model is 5945.3 and for the BSTS model is 6028.5. Here Looic value is the minimum for the ARIMA model and which is the best model. Similarly, for Bahrain cumulative number of deaths LOOIC value for the Bayesian ARIMA model is 2110.7, and for the BSTS model is 2273.2, therefore the minimum LOOIC value is for the ARIMA model. Therefore
almost all countries' data set LOOIC value is minimum for the Bayesian ARIMA model. Hence, when compared to the corresponding BSTS model, the Bayesian ARIMA model performs the best.
Table 4: ARIMA and BSTS model validation using LOOIC and WAIC.
Looic Waic Looic Waic
Cases 10771.5 10771.6 47343.5 47887.2
USA
Deaths 7221.3 7260 7665.8 7665
Cases 8958.4 8958.3 24255.4 24362.8
UK
Deaths 5831.3 5831.2 6448.6 6448.6
Cases 6786 6800.9 7131.2 7180.2
UAE
Deaths 1993.4 1993.5 2236 2236
Cases 5945.3 5945.3 6028.5 6028.5
Bahrain
Deaths 2110.7 2110.9 2273.2 2273.2
Cases 10095.5 10096.1 24219.5 24326.3
India
Deaths 6981.2 6983.5 7279.4 7361.7
References
[1] BOX, G. (1976). Time series analysis, forecasting and control.
[2] Tsay, R. S. and Tiao, G. C. (1984). Consistent estimates of autoregressive parameters and extended sample autocorrelation function for stationary and nonstationary arma models.
Journal of the American Statistical Association, 79(385):84-96.
[3] Pankratz, A. (2012). Forecasting with dynamic regression models. John Wiley & Sons.
[4] Zellner, A. (1971). An introduction to bayesian inference in econometrics. Technical report.
[5] Chib, S. and Greenberg, E. (1994). Bayes inference in regression models with arma (p, q) errors. Journal of Econometrics, 64(1-2):183-206.
[6] Marriott, J., Ravishanker, N., Gelfand, A. E., and Pai, J. (1996). Bayesian analysis of arma processes: Complete sampling-based inference under exact likelihoods. In Bayesian analysis in statistics and econometrics: essays in honor of Arnold Zellner, pages 243-256.
[7] Kleibergen, F. R. and Hoek, H. (2000). Bayesian analysis of arma models.
[8] Fan, J. and Yao, Q. (2003). Nonlinear time series: nonparametric and parametric methods,volume 20. Springer.
[9] Monahan, J. F. (1984). A note on enforcing stationarity in autoregressive-moving average models. Biometrika, 71(2):403-404.
[10] Marriott, J. and Smith, A. (1992). Reparametrization aspects of numerical bayesian methodology for autoregressive moving-average models. Journal of Time Series Analysis, 13(4):327-343.
[11] Scott, S. L. and Varian, H. R. (2014). Predicting the present with bayesian structural time series. International Journal of Mathematical Modelling and Numerical Optimisation, 5(1-2):4-23.
[12] Navas Thorakkattle, M., Farhin, S., et al. (2022). Forecasting the trends of covid-19 and causal impact of vaccines using bayesian structural time series and arima. Annals of Data Science, 9(5):1025-1047.
[13] Feroze, N. (2020). Forecasting the patterns of covid-19 and causal impacts of lockdown in top five affected countries using bayesian structural time series models. Chaos, Solitons & Fractals, 140:110196.
[14] Zhang, Y. and Fricker, J. D. (2021). Quantifying the impact of covid-19 on non-motorized transportation: A bayesian structural time series model. Transport policy, 103:11-20.
[15] Voss, M. S. and Feng, X. (2002). Arma model selection using particle swarm optimization and aic criteria. IFAC Proceedings Volumes, 35(1):349-354.
[16] Vehtari, A., Gelman, A., and Gabry, J. (2015). Efficient implementation of leave-one-out cross-validation and waic for evaluating fitted bayesian models. arXiv preprint arXiv:1507.04544.
[17] Khan, Athar Ali and others (2022). Bayesian Analysis of Type II Generalized Topp-Leone Accelerated Failure Time Models Using R and Stan. Reliability: Theory & Applications, 17(4 (71)), 477-493.
