# Forecasting the Effectiveness of COVID-19 Vaccination Using Vector Autoregressive with an Exogenous Variable: On the Cases of COVID-19 in Indonesia

Authors: Sukono Sukono; Riza Andrian Ibrahim; Riaman Riaman; Elis Hertini; Yuyun Hidayat; Jumadil Saputra
Journal: Discrete Dynamics in Nature and Society
 (2023)

---

## Abstract

This study aims to forecast the COVID-19 spread in Indonesia involving vaccination factors using vector autoregressive with exogenous variables (VARX). The COVID-19 spread represented by active, recovered, and death case rate indicators acts as endogenous variables, while the COVID-19 vaccination represented by second-dose vaccination rates acts as exogenous variables. Because the sum of three COVID-19 spread indicators in one day is one, only two indicators with the highest correlation rates are involved in VARX modelling. The other indicator is practically projected by subtracting one from the sum of two indicator projection results. Based on the analysis results, the active and recovered case rates are two indicators chosen in VARX modelling. Using Akaike information criterion, the most suitable VARX model to project the case and recovered case rates are VARX (7, 1). This model is expected to help the Indonesian government project the COVID-19 spread in Indonesia.


---

## Body

## 1. Introduction

COVID-19 has spread in Indonesia for over a year and a half since 11 March 2020. The spread of COVID-19 has continued to record terrible records. Based on data from the Ministry of Health of the Republic of Indonesia, on 23 June 2021, the number of daily confirmed positive cases of COVID-19 in Indonesia recorded a record high, namely 15,308 new cases. This record continued to be broken in the days that followed. The highest record recorded until 7 November 2021 was 56,757 cases, set on 15 July 2021. Then, the highest number of daily deaths due to COVID-19 recorded until 7 November 2021 was 2,069 cases, where the incident occurred on 21 July 2021. The Indonesian government has made various policies to reduce COVID-19 spread rates. One of the main policies is the National Vaccination Program [1–3]. In the first two months after the start of the policy on 13 January 2021, this policy showed positive developments toward the COVID-19 spread [4]. Active and death case rates of COVID-19 tend to decrease, while COVID-19 recovered case rates increase. As of 7 November 2021, active, recovered, and death case rates of COVID-19 are 0.2548%, 96.3662%, and 3.3790%, respectively. It further gives the Indonesian government a sense of optimism about reducing the COVID-19 spread through vaccination [5].Knowledge of models that can assist this projection process is needed to project how many vaccination rates must be achieved every day to realize the desired COVID-19 spread rates. Therefore, this study aims to design a model for projecting the COVID-19 spread rates by involving the COVID-19 vaccination in Indonesia. The model used is vector autoregressive with exogenous variables (VARX). This model is chosen because it can explain the causal relationship between endogenous variables, where exogenous variables influence endogenous variables. The COVID-19 spread rates can be measured using various indicators. This study uses three indicators: active, recovered, and death case rates. The sum of three indicators in one day is one. So, to simplify the projection, only two indicators are involved in VARX modelling. The projection of the other indicator is made by subtracting one from the projection results of the two indicators modelled by VARX. Then, the second-dose COVID-19 vaccination rates represent the COVID-19 vaccination factor against the Indonesian population. It is done because people in Indonesia have been wholly vaccinated if they have been vaccinated twice [6]. This research is expected to help the Indonesian government project vaccination rates the next day to control the COVID-19 spread in Indonesia.


## 2. Literature Review

Forecasting is closely related to time series analysis [7, 8]. Therefore, articles on modelling and forecasting the COVID-19 spread in Indonesia using the concept of time series analysis are discussed briefly in this section. The articles discussed are English articles indexed in the Scopus database. Kirana and Bhawiyuga [9] modelled and forecasted the COVID-19 spread using naive methods. The COVID-19 spread used is the number of monthly COVID-19 cases. Their model’s mean absolute percentage error (MAPE) is 15.85%. On the basis of the accurate classification of the model, Wei [10] stated that this model is suitable for forecasting the number of monthly COVID-19 cases in the next month. Then, Djakaria and Saleh [11] modelled and forecasted the COVID-19 spread in Gorontalo Province, Indonesia, using the Holt-Winters smoothing method. The cumulative number of COVID-19 cases is used as the COVID-19 spread indicator. The MAPE of the model obtained is 6.14%, so the model is perfect for forecasting the cumulative number of COVID-19 cases on the next day based on the model’s accurate classification by Wei [10].Then, Aditya Satrio et al. [12] analyzed and forecasted the COVID-19 spread in Indonesia, where the COVID-19 spread indicators used are the cumulative number of cases, recovered, and death of COVID-19. Each indicator is modelled separately using the autoregressive integrated moving average (ARIMA) model. Using the same method, Fadly and Sari [13] modelled and forecasted the number of funerals carried out in DKI Jakarta Province, Indonesia, since COVID-19 hit. The articles described in the previous paragraph use univariate time series analysis, whereas multivariate time series analysis can also be used. Fitriani et al. [14] simultaneously model and forecast the number of daily COVID-19-positive cases in Indonesia and Singapore. The multivariate time series model used is vector autoregressive integrated (VARI). The MAPE obtained in each place was 32.74% and 37%, respectively. The MAPE values are enormous, so based on the model’s accurate classification by Wei [10], the model is unsuitable for forecasting. Then, Meimela et al. [15] modelled and forecasted the COVID-19 spread in Indonesia using the vector autoregressive integrated moving average (VARIMA). The COVID-19 spread indicators used are the daily positive and death cases. VARIMA (1, 1, 1) is the best model got.Also, Akbar et al. [16] analyzed the relationship between COVID-19 active case rates assumed to be influenced by population growth rates. Bandung City and Purwakarta Regency are the two areas considered, whereas the population growth rate in West Java Province represents the population growth rate. The vector autoregressive integrated with exogenous variables (VARX) describes the relationship. The MAPE obtained in each place is 1% and 11.8%, respectively. It aligns with Wei’s accurate classification model (2006); the model is perfect for forecasting the COVID-19 spread in Bandung City, while the other is only good. As explained in the previous paragraphs, no articles in the Scopus database have explored modelling and forecasting active, recovered, and death case rates of COVID-19 in Indonesia, considering vaccination factors using the VARX model. Therefore, this is a good opportunity to perform such modelling and forecasting.


## 3. Materials and Methods

### 3.1. Materials

The data used in this study are (1) active, recovered, and death case rates of COVID-19 data and (2) second-dose COVID-19 vaccination rates data. The period considered is from 27 January 2021, when the second-dose vaccination program starts, until 7 November 2021. All data are obtained from the official website of the Science and Technology Index, National Research and Innovation Agency, Ministry of Research and Technology of the Republic of Indonesia.


### 3.2. Methods

#### 3.2.1. Vector Autoregressive with Exogenous Variables

The general form of the vector autoregressive with exogenous variables (VARX) model with the order of the endogenous variablesp and the order of the endogenous variables q, VARX (p,q), is expressed as follows [17–21]:(1)zt=μ+∑i=1pϕizt−i+∑j=0qθjxt−j+et,where zt represents the n-dimensional vector containing the endogenous variables at time t, ϕi,i=1,2,…,p express the n×n-order coefficient matrix of zt−i, xt represents the r-dimensional vector containing the exogenous variables at time t, θj,j=0,1,…,q express the n×r-order coefficient matrix of xt−j, and et represents the n-dimensional random error vector. Equation (1) can also be reformulated as follows [22]:(2)ϕBzzt=μ+θBxxt+et,where ϕBz=In−∑i=1pϕiBzi, θBx=∑j=0pθjBxj, Bzizt=zt−i, and Bxjxt=xt−j.There are three assumptions of the VARX model presented as follows [23]:(1)
The VARX (p, q) model is stationary if the zero generators of the ϕv and θw determinants are outside the unit circle(2)
et is independent and identically multivariate-normally distributed with zero-vector mean and the constant covariance-matrix variance, and Σ. Σ is a positive semidefinite matrix expressed as follows: Σ=σek1,t,ek2,t, where k1,k2=1,2,…,n(3)
et and xt are independentMany criteria can be used to determine the order of the VARX model. Lui et al. [22] stated that the order of the VARX model is determined through the Akaike information criterion (AIC). The model with the smallest AIC is the best. The AIC value of a model can be determined by the following equation [24]:(3)AIC=lnΣ+2LT,where T is the size of observation data, and L is many model parameters.


#### 3.2.2. Augmented Dicky-Fuller Test

The unit root test equation for Augmented Dicky-Fuller (ADF) from data that have constant and trend elements is as follows [25–27]:(4)zk,t′=ηk+ψkt+βkzk,t−1+∑i=1p−1ζk,izk,t−i′+ek,t,where zk,t represents observation data k at the time t, zk,t′ represents the difference between zk,t and zk,t−1, ηk is a constant, ψk represents the trend parameter, ζk,i is the parameter of zk,t−i′, p represents the optimum lag, and ek,t is random error that is independent and identically normally distributed with zero mean and constant variance. The followings are the test hypotheses used: H0: βk=0, and H1: βk<0. The test statistic used is the t-ratio of the parameter estimator βk estimated using the least square method stated as follows:(5)τ=β^ksβ^k,where β^k is an estimator βk, and sβ^k is the standard deviation of β^k. If τ less than the critical value, then H0 rejected which means the data are stationary and vice versa. If H0 is not rejected, then data must be differentiated until it is stationary [28, 29].


#### 3.2.3. Granger Causality Test

Suppose that there are two stationary time seriesz1,t and z2,t. In certain cases, z1,t can be affected by z1,t−i and z2,t−i, and z2,t can also be affected by z1,t−i and z2,t−i. This relationship is called a bidirectional causality relationship [30, 31]. The Granger causality test is a popular statistical test used to check the relationship.The relationship betweenz1,t with z1,t−i, which is also called a restricted relationship, is stated as follows [32]:(6)z1,t=μ1+∑i=1pϕ1,iz1,t−i+e1,t.Then, the relationship betweenz1,t with z1,t−i and z2,t−i, which is also called the unrestricted relationship, is stated as follows [33–35]:(7)z1,t=μ1+∑i=1pϕ1,iz1,t−i+∑i=1pϕ2,iz2,t−i+e1,t.The test hypotheses used are as follows:H0:∀i,ϕ2,i=0, and H1:∃i,ϕ2,i≠0. Test statistics F determined by the following equation [36]:(8)F=Rur2−Rr2T−p−11−Rur2p,where Rur2 and Rr2 represent the coefficient of determination in equations (6) and (7), respectively. If the F test statistic is greater than Fα,p,T−p−1, then H0 is rejected, which means z1,t influenced by z1,t−i and z2,t−i [37]. z1,t and z2,t have a bidirectional causality if z1,t influenced by z1,t−i and z2,t−i, and z2,t influenced by z1,t−i and z2,t−i [38, 39].


#### 3.2.4. VARX Parameter Estimation Using Maximum Likelihood (ML) Method

The VARX models withn=2,r=1,p=1, and q=0 can be stated as follows:(9)yt=Xtβ+et,where(10)yt=z1,tz2,t,Xt=10z1,t−1z2,t−100x1,t00100z1,t−1z2,t−10x1,t,β=μ1μ2ϕ111ϕ121ϕ211ϕ221θ110θ210,et=e1,2e2,2.Based on equation (9), et can be stated as follows:(11)et=yt−Xtβ.Becauseet∼N0,Σ, then yt−Xtβ∼N0,Σ. Therefore, the distribution function of et can be expressed as follows [40, 41]:(12)fetβ,Σ=12πΣexp−12yt−Xtβ′∑−1yt−Xtβ.Vector estimator ofβ obtained by maximizing the likelihood function of et which is stated as follows:(13)lβ,Σe2,…,eT=12πΣT−1exp−12∑t=2Tyt−Xtβ′∑−1yt−Xtβ.Generally, equation (14) is transformed into a natural logarithm to facilitate maximization. This function is then referred to as the log-likelihood function, which is expressed as follows:(14)llnlβ,Σe2,…,eT=−T−1ln2πΣ−12∑t=2Tyt−Xtβ′∑−1yt−Xtβ.Vector estimator ofβ is the zero-generator vector of the derivative of equation (15) on β. The following is the equation used to estimate the vector β:(15)β^=∑t=2TXt′Σ−1Xt−1∑t=2TXt′Σ−1yt.


#### 3.2.5. Diagnostic Test

The diagnostic test is a feasibility test of the model for forecasting. This feasibility test includes checking the assumption of error and the model’s mean absolute percentage error (MAPE). The model errors are tested to determine whether they are independent and identically normally distributed with zero mean and constant variance. This assumption is also known as the white noise assumption. One of the statistical tests commonly used to test whether the errors are independent or not is the Ljung-Box test. The null hypothesisH0 in this test is that the errors are independent, while this test’s alternative hypothesis H1 is the opposite. The test statistic with the lag length M denoted by QM is determined using the following equation [42, 43]:(16)QM=TT+2∑m=1Mrm2T−m,where(17)rm=∑t=m+1Tek,tek,t−m∑t=1Tek,t2.RejectH0 if QM>χ1−α,dfek,t2 [44, 45].Meanwhile, the naked eye can check the assumption that the error is normally distributed with zero mean and constant variance using the quantile-quantile plot (Q-Q Plot) [46]. Suppose FNkek,t is the value of the normal cumulative distribution function with zero mean and constant variance of ek,t. If the scatter of point pairs ek,t,FNk−1ek,t is around a line with a gradient one, then ek,t is visually normally distributed with zero mean and constant variance. If the assumption of white noise in the model error is met and the MAPE of the model is less than 20%, then the model obtained is suitable for use in forecasting.


## 3.1. Materials

The data used in this study are (1) active, recovered, and death case rates of COVID-19 data and (2) second-dose COVID-19 vaccination rates data. The period considered is from 27 January 2021, when the second-dose vaccination program starts, until 7 November 2021. All data are obtained from the official website of the Science and Technology Index, National Research and Innovation Agency, Ministry of Research and Technology of the Republic of Indonesia.


## 3.2. Methods

### 3.2.1. Vector Autoregressive with Exogenous Variables

The general form of the vector autoregressive with exogenous variables (VARX) model with the order of the endogenous variablesp and the order of the endogenous variables q, VARX (p,q), is expressed as follows [17–21]:(1)zt=μ+∑i=1pϕizt−i+∑j=0qθjxt−j+et,where zt represents the n-dimensional vector containing the endogenous variables at time t, ϕi,i=1,2,…,p express the n×n-order coefficient matrix of zt−i, xt represents the r-dimensional vector containing the exogenous variables at time t, θj,j=0,1,…,q express the n×r-order coefficient matrix of xt−j, and et represents the n-dimensional random error vector. Equation (1) can also be reformulated as follows [22]:(2)ϕBzzt=μ+θBxxt+et,where ϕBz=In−∑i=1pϕiBzi, θBx=∑j=0pθjBxj, Bzizt=zt−i, and Bxjxt=xt−j.There are three assumptions of the VARX model presented as follows [23]:(1)
The VARX (p, q) model is stationary if the zero generators of the ϕv and θw determinants are outside the unit circle(2)
et is independent and identically multivariate-normally distributed with zero-vector mean and the constant covariance-matrix variance, and Σ. Σ is a positive semidefinite matrix expressed as follows: Σ=σek1,t,ek2,t, where k1,k2=1,2,…,n(3)
et and xt are independentMany criteria can be used to determine the order of the VARX model. Lui et al. [22] stated that the order of the VARX model is determined through the Akaike information criterion (AIC). The model with the smallest AIC is the best. The AIC value of a model can be determined by the following equation [24]:(3)AIC=lnΣ+2LT,where T is the size of observation data, and L is many model parameters.


### 3.2.2. Augmented Dicky-Fuller Test

The unit root test equation for Augmented Dicky-Fuller (ADF) from data that have constant and trend elements is as follows [25–27]:(4)zk,t′=ηk+ψkt+βkzk,t−1+∑i=1p−1ζk,izk,t−i′+ek,t,where zk,t represents observation data k at the time t, zk,t′ represents the difference between zk,t and zk,t−1, ηk is a constant, ψk represents the trend parameter, ζk,i is the parameter of zk,t−i′, p represents the optimum lag, and ek,t is random error that is independent and identically normally distributed with zero mean and constant variance. The followings are the test hypotheses used: H0: βk=0, and H1: βk<0. The test statistic used is the t-ratio of the parameter estimator βk estimated using the least square method stated as follows:(5)τ=β^ksβ^k,where β^k is an estimator βk, and sβ^k is the standard deviation of β^k. If τ less than the critical value, then H0 rejected which means the data are stationary and vice versa. If H0 is not rejected, then data must be differentiated until it is stationary [28, 29].


### 3.2.3. Granger Causality Test

Suppose that there are two stationary time seriesz1,t and z2,t. In certain cases, z1,t can be affected by z1,t−i and z2,t−i, and z2,t can also be affected by z1,t−i and z2,t−i. This relationship is called a bidirectional causality relationship [30, 31]. The Granger causality test is a popular statistical test used to check the relationship.The relationship betweenz1,t with z1,t−i, which is also called a restricted relationship, is stated as follows [32]:(6)z1,t=μ1+∑i=1pϕ1,iz1,t−i+e1,t.Then, the relationship betweenz1,t with z1,t−i and z2,t−i, which is also called the unrestricted relationship, is stated as follows [33–35]:(7)z1,t=μ1+∑i=1pϕ1,iz1,t−i+∑i=1pϕ2,iz2,t−i+e1,t.The test hypotheses used are as follows:H0:∀i,ϕ2,i=0, and H1:∃i,ϕ2,i≠0. Test statistics F determined by the following equation [36]:(8)F=Rur2−Rr2T−p−11−Rur2p,where Rur2 and Rr2 represent the coefficient of determination in equations (6) and (7), respectively. If the F test statistic is greater than Fα,p,T−p−1, then H0 is rejected, which means z1,t influenced by z1,t−i and z2,t−i [37]. z1,t and z2,t have a bidirectional causality if z1,t influenced by z1,t−i and z2,t−i, and z2,t influenced by z1,t−i and z2,t−i [38, 39].


### 3.2.4. VARX Parameter Estimation Using Maximum Likelihood (ML) Method

The VARX models withn=2,r=1,p=1, and q=0 can be stated as follows:(9)yt=Xtβ+et,where(10)yt=z1,tz2,t,Xt=10z1,t−1z2,t−100x1,t00100z1,t−1z2,t−10x1,t,β=μ1μ2ϕ111ϕ121ϕ211ϕ221θ110θ210,et=e1,2e2,2.Based on equation (9), et can be stated as follows:(11)et=yt−Xtβ.Becauseet∼N0,Σ, then yt−Xtβ∼N0,Σ. Therefore, the distribution function of et can be expressed as follows [40, 41]:(12)fetβ,Σ=12πΣexp−12yt−Xtβ′∑−1yt−Xtβ.Vector estimator ofβ obtained by maximizing the likelihood function of et which is stated as follows:(13)lβ,Σe2,…,eT=12πΣT−1exp−12∑t=2Tyt−Xtβ′∑−1yt−Xtβ.Generally, equation (14) is transformed into a natural logarithm to facilitate maximization. This function is then referred to as the log-likelihood function, which is expressed as follows:(14)llnlβ,Σe2,…,eT=−T−1ln2πΣ−12∑t=2Tyt−Xtβ′∑−1yt−Xtβ.Vector estimator ofβ is the zero-generator vector of the derivative of equation (15) on β. The following is the equation used to estimate the vector β:(15)β^=∑t=2TXt′Σ−1Xt−1∑t=2TXt′Σ−1yt.


### 3.2.5. Diagnostic Test

The diagnostic test is a feasibility test of the model for forecasting. This feasibility test includes checking the assumption of error and the model’s mean absolute percentage error (MAPE). The model errors are tested to determine whether they are independent and identically normally distributed with zero mean and constant variance. This assumption is also known as the white noise assumption. One of the statistical tests commonly used to test whether the errors are independent or not is the Ljung-Box test. The null hypothesisH0 in this test is that the errors are independent, while this test’s alternative hypothesis H1 is the opposite. The test statistic with the lag length M denoted by QM is determined using the following equation [42, 43]:(16)QM=TT+2∑m=1Mrm2T−m,where(17)rm=∑t=m+1Tek,tek,t−m∑t=1Tek,t2.RejectH0 if QM>χ1−α,dfek,t2 [44, 45].Meanwhile, the naked eye can check the assumption that the error is normally distributed with zero mean and constant variance using the quantile-quantile plot (Q-Q Plot) [46]. Suppose FNkek,t is the value of the normal cumulative distribution function with zero mean and constant variance of ek,t. If the scatter of point pairs ek,t,FNk−1ek,t is around a line with a gradient one, then ek,t is visually normally distributed with zero mean and constant variance. If the assumption of white noise in the model error is met and the MAPE of the model is less than 20%, then the model obtained is suitable for use in forecasting.


## 3.2.1. Vector Autoregressive with Exogenous Variables

The general form of the vector autoregressive with exogenous variables (VARX) model with the order of the endogenous variablesp and the order of the endogenous variables q, VARX (p,q), is expressed as follows [17–21]:(1)zt=μ+∑i=1pϕizt−i+∑j=0qθjxt−j+et,where zt represents the n-dimensional vector containing the endogenous variables at time t, ϕi,i=1,2,…,p express the n×n-order coefficient matrix of zt−i, xt represents the r-dimensional vector containing the exogenous variables at time t, θj,j=0,1,…,q express the n×r-order coefficient matrix of xt−j, and et represents the n-dimensional random error vector. Equation (1) can also be reformulated as follows [22]:(2)ϕBzzt=μ+θBxxt+et,where ϕBz=In−∑i=1pϕiBzi, θBx=∑j=0pθjBxj, Bzizt=zt−i, and Bxjxt=xt−j.There are three assumptions of the VARX model presented as follows [23]:(1)
The VARX (p, q) model is stationary if the zero generators of the ϕv and θw determinants are outside the unit circle(2)
et is independent and identically multivariate-normally distributed with zero-vector mean and the constant covariance-matrix variance, and Σ. Σ is a positive semidefinite matrix expressed as follows: Σ=σek1,t,ek2,t, where k1,k2=1,2,…,n(3)
et and xt are independentMany criteria can be used to determine the order of the VARX model. Lui et al. [22] stated that the order of the VARX model is determined through the Akaike information criterion (AIC). The model with the smallest AIC is the best. The AIC value of a model can be determined by the following equation [24]:(3)AIC=lnΣ+2LT,where T is the size of observation data, and L is many model parameters.


## 3.2.2. Augmented Dicky-Fuller Test

The unit root test equation for Augmented Dicky-Fuller (ADF) from data that have constant and trend elements is as follows [25–27]:(4)zk,t′=ηk+ψkt+βkzk,t−1+∑i=1p−1ζk,izk,t−i′+ek,t,where zk,t represents observation data k at the time t, zk,t′ represents the difference between zk,t and zk,t−1, ηk is a constant, ψk represents the trend parameter, ζk,i is the parameter of zk,t−i′, p represents the optimum lag, and ek,t is random error that is independent and identically normally distributed with zero mean and constant variance. The followings are the test hypotheses used: H0: βk=0, and H1: βk<0. The test statistic used is the t-ratio of the parameter estimator βk estimated using the least square method stated as follows:(5)τ=β^ksβ^k,where β^k is an estimator βk, and sβ^k is the standard deviation of β^k. If τ less than the critical value, then H0 rejected which means the data are stationary and vice versa. If H0 is not rejected, then data must be differentiated until it is stationary [28, 29].


## 3.2.3. Granger Causality Test

Suppose that there are two stationary time seriesz1,t and z2,t. In certain cases, z1,t can be affected by z1,t−i and z2,t−i, and z2,t can also be affected by z1,t−i and z2,t−i. This relationship is called a bidirectional causality relationship [30, 31]. The Granger causality test is a popular statistical test used to check the relationship.The relationship betweenz1,t with z1,t−i, which is also called a restricted relationship, is stated as follows [32]:(6)z1,t=μ1+∑i=1pϕ1,iz1,t−i+e1,t.Then, the relationship betweenz1,t with z1,t−i and z2,t−i, which is also called the unrestricted relationship, is stated as follows [33–35]:(7)z1,t=μ1+∑i=1pϕ1,iz1,t−i+∑i=1pϕ2,iz2,t−i+e1,t.The test hypotheses used are as follows:H0:∀i,ϕ2,i=0, and H1:∃i,ϕ2,i≠0. Test statistics F determined by the following equation [36]:(8)F=Rur2−Rr2T−p−11−Rur2p,where Rur2 and Rr2 represent the coefficient of determination in equations (6) and (7), respectively. If the F test statistic is greater than Fα,p,T−p−1, then H0 is rejected, which means z1,t influenced by z1,t−i and z2,t−i [37]. z1,t and z2,t have a bidirectional causality if z1,t influenced by z1,t−i and z2,t−i, and z2,t influenced by z1,t−i and z2,t−i [38, 39].


## 3.2.4. VARX Parameter Estimation Using Maximum Likelihood (ML) Method

The VARX models withn=2,r=1,p=1, and q=0 can be stated as follows:(9)yt=Xtβ+et,where(10)yt=z1,tz2,t,Xt=10z1,t−1z2,t−100x1,t00100z1,t−1z2,t−10x1,t,β=μ1μ2ϕ111ϕ121ϕ211ϕ221θ110θ210,et=e1,2e2,2.Based on equation (9), et can be stated as follows:(11)et=yt−Xtβ.Becauseet∼N0,Σ, then yt−Xtβ∼N0,Σ. Therefore, the distribution function of et can be expressed as follows [40, 41]:(12)fetβ,Σ=12πΣexp−12yt−Xtβ′∑−1yt−Xtβ.Vector estimator ofβ obtained by maximizing the likelihood function of et which is stated as follows:(13)lβ,Σe2,…,eT=12πΣT−1exp−12∑t=2Tyt−Xtβ′∑−1yt−Xtβ.Generally, equation (14) is transformed into a natural logarithm to facilitate maximization. This function is then referred to as the log-likelihood function, which is expressed as follows:(14)llnlβ,Σe2,…,eT=−T−1ln2πΣ−12∑t=2Tyt−Xtβ′∑−1yt−Xtβ.Vector estimator ofβ is the zero-generator vector of the derivative of equation (15) on β. The following is the equation used to estimate the vector β:(15)β^=∑t=2TXt′Σ−1Xt−1∑t=2TXt′Σ−1yt.


## 3.2.5. Diagnostic Test

The diagnostic test is a feasibility test of the model for forecasting. This feasibility test includes checking the assumption of error and the model’s mean absolute percentage error (MAPE). The model errors are tested to determine whether they are independent and identically normally distributed with zero mean and constant variance. This assumption is also known as the white noise assumption. One of the statistical tests commonly used to test whether the errors are independent or not is the Ljung-Box test. The null hypothesisH0 in this test is that the errors are independent, while this test’s alternative hypothesis H1 is the opposite. The test statistic with the lag length M denoted by QM is determined using the following equation [42, 43]:(16)QM=TT+2∑m=1Mrm2T−m,where(17)rm=∑t=m+1Tek,tek,t−m∑t=1Tek,t2.RejectH0 if QM>χ1−α,dfek,t2 [44, 45].Meanwhile, the naked eye can check the assumption that the error is normally distributed with zero mean and constant variance using the quantile-quantile plot (Q-Q Plot) [46]. Suppose FNkek,t is the value of the normal cumulative distribution function with zero mean and constant variance of ek,t. If the scatter of point pairs ek,t,FNk−1ek,t is around a line with a gradient one, then ek,t is visually normally distributed with zero mean and constant variance. If the assumption of white noise in the model error is met and the MAPE of the model is less than 20%, then the model obtained is suitable for use in forecasting.


## 4. Results

The symbols used for each data are stated as follows:z1,t states daily COVID-19 active case rate data, z2,t states daily COVID-19 recovered case rate data, z3,t states daily COVID-19 death case rate data, and x1,t states second-dose COVID-19 vaccination rate data. The sum of the three indicators in one day is one, so to simplify the modelling and forecasting process, and only two indicators are used in VARX modelling. The remaining indicator is practically projected by subtracting one from the two indicators’ projected results. The two indicators selected are the two indicators that have the highest correlation rates.


### 4.1. Data Descriptive

Descriptive statistics for each data obtained using Microsoft Excel software are presented in Table1.Table 1
Descriptive statistics of each data.
Descriptive statisticsz1,t (%)z2,t (%)z3,t (%)x1,t (%)Average7.811989.28122.90697.8703Variance0.26780.24700.00080.6484Standard deviation5.17504.97030.28078.0520Table1 shows the intensity of active, recovery, death, and vaccination rates from COVID-19 in Indonesia per day, each of which is different. The intensity values are 7.8119%, 89.2812%, 2.9069%, and 7.8703% per day, respectively. The intensity values of the active and death cases appear to be lower than that of vaccination. It indicates that the number of Indonesian citizens exposed and died from COVID-19 is lower than the number of people vaccinated daily. It shows that the Indonesian government is serious about handling COVID-19 via vaccination. Then, there are deviations from the active, recovery, death, and vaccination rates from COVID-19 against each average daily. On average, the deviation values from the active, recovery, death, and vaccination rates from COVID-19 were 5.1750%, 4.9703%, 0.2807%, and 8.0520% per day, respectively. The daily deviations of the active and death rates are lower than the deviation of the vaccination rate. It indicates that Indonesian citizens are increasingly varied of COVID-19 with a high willingness to vaccinate. It prevents extreme value deviations from occurring. Finally, the variance of z1,t, z2,t, z3,t, and x1,t is 0.2678%, 0.2470%, 0.0008%, and 0.6484%, respectively. The interpretation of the variance values is similar to the interpretation of the deviation value. Then, the correlation rates between z1,t, z2,t, and z3,t obtained using Microsoft Excel software are presented in Table 2.Table 2
The correlation rates betweenz1,t, z2,t, and z3,t.
Data pairsThe correlation ratesz1,t and x1,t−0.6599z2,t and x1,t0.6352Table2 displays the correlation rate between z1,t and z2,t, −0.9993, which is the highest correlation rate among the correlation rates of the other indicator data pairs. It is not surprising because the two are rationally the opposite of each other. If the active case rate increases, the recovered case rate decreases, and vice versa. This fact can be seen visually in Figure 1. Therefore, z1,t and z2,t are selected as endogenous variables. The high negative correlation between z1,t and z2,t is not a problem in VARX modelling. It is the same as regression modelling, where if the correlation rate between variables is high, the model formed can provide an excellent description of the data.Figure 1
Graph ofz1,t (a), z2,t (b), and x1,t (c).
(a)(b)(c)Table3 captures the correlation rates between endogenous and exogenous variables above 0.6. It indicates that endogenous and exogenous variables are strong enough.Table 3
The correlation rates between endogenous and exogenous variables.
Data pairsThe correlation ratesz1,t and z2,t−0.9993z1,t and z3,t−0.7421z2,t and z3,t0.7162


### 4.2. Data Stationary Test

In this section,z1,t, z2,t, and x1,t are checked for stationary visually and formally. Visually, the stationarity of the three data is seen from each graph, while formally, the stationary of the data is tested using the ADF test. Each data graph is observed whether it is around the average line or not. If the data graph is around the average line, then the data are stationary, and vice versa. The graphs of z1,t, z2,t, and x1,t with their respective average lines are presented in Figure 1.Figure1 displays that the graphs of the three data do not appear to be around the average line. Therefore, visually, the three data are not stationary. Then, the formal stationarity test of the three data is carried out using the ADF test. The significance level and optimal lag used are 0.05 and 6, respectively. The test statistics of each data are determined using equation (5). The formal stationary test result of the three data using the ADF test is presented in Table 4.Table 4
Data stationary testz1,t, z2,t, and x1,t.
Dataz1,tz2,tx1,tτ−2.7926−2.78152.4323Critical value−3.2000−3.2000−3.2000Based on Table4, all data have a value of not less than the critical value. Therefore, H0 on each data test is not rejected, which means that all the data are not stationary, so the differencing process using equation (3) is carried out. After the data are differentiated, the next step is the three data’s visual and formal stationarity tests. After the differencing process is carried out twice, all data are stationary. The graphs of z1,t″, z2,t″, and x1,t″ with their respective average lines are presented in Figure 2.Figure 2
Graph ofz1,t″ (a), z2,t″ (b), and x1,t″ (c).
(a)(b)(c)Figure2 shows the graphs of three data differentiated twice and appear to be around the average line. Therefore, visually, they are stationary. Next, the formal stationarity test is carried out using the ADF test. The significance level and optimal lag used are 0.05 and 6, respectively. The test statistics of each data are determined using equation (5). Furthermore, a formal summary of the three data stationary tests using the ADF test is presented in Table 5.Table 5
Data stationary testz1,t″, z2,t″, and x1,t″.
Dataz1,t″z2,t″x1,t″τ−8.2813−8.3078−13.7290Critical value−3.2000−3.2000−3.2000Table5 captures the three data differentiated twice that have a τ value less than the critical value. Therefore, H0 of each test is rejected, which means that all data are stationary.


### 4.3. Causality Relationship Test

A check of the existence of a bidirectional causal relationship betweenz1,t″ and z2,t″ in this study is carried out using the Granger causality test. The significance levels, optimal lag, and residual degrees of freedom used are 0.05, 7, and 256, respectively. The effect of z1,t−1″ and z2,t−i″ on z1,t″ is tested first. The F value obtained is 2.6931. This value is greater than the critical value, 2.0442. Therefore, H0 is rejected, which means that z1,t″ is affected by z1,t−1″ and z2,t−i″. Next, the effect of z1,t−1″ and z2,t−i″ on z2,t″ is tested. The F value obtained is 2.7514. This value is also greater than the critical value, 2.0442. Therefore, H0 is rejected which means that z2,t″ is affected by z1,t−1″ and z2,t−i″. Because z1,t″ is affected by z1,t−1″ and z2,t−i″, and z2,t″ is affected by z1,t−1″ and z2,t−i″, then z1,t″ and z2,t″ have a bidirectional causality relationship.


### 4.4. VARX Order Identification

The AIC value, as described in Section2, is used as a criterion to determine the best order of endogenous and exogenous variables of the VARX model. The model with the smallest AIC value is the most suitable. The AIC value for each model is determined by equation (3). The maximum order of the endogenous and exogenous variables considered is 10. The AIC value for each model is determined using the R software. Based on the results obtained from the software, the AIC value of the VARX (7, 1), −34.2194, is the smallest AIC value among other models. Therefore, the most suitable model chosen is the VARX (7, 1).


### 4.5. VARX (7, 1) Parameter Estimation Using ML Method

The VARX (7, 1) parameters are estimated using the ML method, equation (14), as described in Section 2. The estimation process is carried out using the help of Microsoft Excel software. The parameter estimation results are presented in Table 6.Table 6
Results of parameter estimates of VARX (7, 1).
ParameterEstimatorμ11.745E − 05μ2−1.817E – 05ϕ1,110.5495ϕ1,210.7931ϕ2,11−1.0376ϕ2,21−1.2881ϕ1,121.2891ϕ1,221.4981ϕ2,12−1.5079ϕ2,22−1.7207ϕ1,135.3873ϕ1,235.6679ϕ2,13−5.5021ϕ2,23−5.7843ϕ1,145.6900ϕ1,246.0157ϕ2,14−5.8682ϕ2,24−6.1934ϕ1,158.9929ϕ1,259.4176ϕ2,15−9.0363ϕ2,25−9.4565ϕ1,163.2196ϕ1,263.3846ϕ2,16−3.3553ϕ2,26−3.5186ϕ1,170.7175ϕ1,270.6588ϕ2,17−0.8463ϕ2,27−0.7909θ1,10−0.0117θ2,100.0104θ1,110.2203θ2,11−0.2162On the basis of the parameters obtained in Table6, the VARX (7, 1) model is stated in equation (16).(18)z^1,t″z^2,t″=1.7448E−05−1.8174E−05+0.54950.7931−1.0376−1.2881z1,t−1″z2,t−1″+1.28911.4981−1.5079−1.7207z1,t−2″z2,t−2″+5.38735.6679−5.5021−5.7843z1,t−3″z2,t−3″+5.69006.0157−5.8682−6.1934z1,t−4″z2,t−4″+8.99299.4176−9.0363−9.4565z1,t−5″z2,t−5″+3.21963.3846−3.3553−3.5186z1,t−6″z2,t−6″+0.71750.6588−0.8463−0.7909z1,t−7″z2,t−7″+−0.01170.0104x1,t″+0.2203−0.2162x1,t−1″.In order to be used to forecast the actual data, equation (18) is transformed to the actual data as follows:(19)z^1,tz^2,t=1.7448E−05−1.8174E−052.54950.7931−1.03760.7119z1,t−1z2,t−1+−0.8099−0.08820.5673−0.1446z1,t−2z2,t−2+3.35863.4649−3.5239−3.6309z1,t−3z2,t−3+−3.7954−3.82203.62813.6545z1,t−4z2,t−4+3.0003.0540−2.8021−2.8541z1,t−5z2,t−5+−9.0762−9.43498.84919.2011z1,t−6z2,t−6+3.27133.3072−3.1719−3.2103z1,t−7z2,t−7+1.78452.0670−1.6627−1.9367z1,t−8z2,t−8+0.71790.6588−0.8463−0.7909z1,t−9z2,t−9+−0.01170.0104x1,t+0.2437−0.2370x1,t−1+−0.45230.4428x1,t−2+0.2203−0.2162x1,t−3.Meanwhile,z3,t can be modelled as follows:(20)z^3,t=1−z^1,t−z^2,t.


### 4.6. Diagnostic Test of VARX (7, 1)

Independence of bothe1,t and e2,t to time is carried out using the Ljung-Box test described in section 2. The significance levels and the degrees of freedom used are 0.05 and 7, respectively. The test statistic is determined using equation (18). The test statistics of e1,t and e2,t obtained are 3.3610 and 3.4270, respectively. The test statistic values are smaller than the critical value, 14.0671. Therefore, H0 of each independence test of e1,t and e2,t to times is not rejected, which means that e1,t and e2,t are independent of time. Next, the normality test of the VARX (7, 1) model error is visually carried out first through the Q-Q plot. Q-Q plots for e1,t and e2,t are presented in Figure 3.Figure 3
Q-Q plots of e1,t (a) and e2,t (b).
(a)(b)Figure3 indicates the scattering of points in each plot appears to be around a line with one gradient. Therefore, e1,t and e2,t visually follow a normal distribution with zero mean and constant variances, namely 0.00108 and 0.00107. Based on this, the model equation (20) errors fulfill the assumption that they are independent and identically normally distributed with zero mean and constant variance. The next step is checking the error size of the model. This study measures the error size using the mean absolute percentage error (MAPE). The MAPE value of z1,t and z2,t is 1.5751% and 0.0833%, respectively. The MAPE value is less than 10%, even close to zero. Therefore, according to Wei [10], the VARX (7, 1) model is very feasible to forecast z1,t and z2,t in the next day based on the feasibility criteria of the MAPE-based model.The feasibility of the VARX (7, 1) model is also checked using the visual comparison between the forecast results and the actual data. The visual comparison between the forecast results and the actual data is presented in Figure4.Figure 4
The visual comparison between thez1,t and z1,t forecast (a) and z2,t and z2,t forecast (b).
(a)(b)Figure4 shows that the forecast line of z1,t and z2,t appears very close to the actual data of z1,t and z2,t. It means that the VARX (7, 1) model is visually feasible to forecast z1,t and z2,t in the next day. The VARX (7, 1) model is also practically forecast z1,t and z2,t in the next few days.


### 4.7. Forecasting

Forecastingz1,t, z2,t, and z3,t are carried out on six days, from 8 November 2021 until 13 November 2021. These forecasts are practically carried out using equation (19) for z1,t and z2,t and equation (20) for z3,t. The second-dose COVID-19 vaccination rates used for this forecast are the same as the actual on that date, 29.5200%, 29.7817%, 30.1154%, 30.4491%, 30.7828%, and 31.0465%. The z1,t, z2,t, and z3,t forecasting results obtained are presented in Table 7.Table 7
Thez1,t, z2,t, and z3,t forecasting results.
Datez1,t forecast (%)z2,t forecast (%)z3,t forecast (%)8 November 20210.227596.39303.37959 November 20210.235096.38623.378810 November 20210.224796.39603.379311 November 20210.238996.38223.378812 November 20210.222096.39883.379113 November 20210.216596.40423.3793Table7 shows that the forecasting results of z1,t, z2,t, and z3,t from 8 November 2021 until 13 November 2021 tend to decrease. To determine the accuracy of the forecast results presented in Table 7, we will look at the mean absolute percentage error (MAPE) size. The MAPE of z1,t, z2,t, and z3,t forecasting results on 8 November 2021 until 13 November 2021 is 2.2853%, 0.0051%, and 0.0082%, respectively. This MAPE of z1,t, z2,t, and z3,t forecasting results is tiny. Therefore, equations (19) and (20) are suitable for forecasting z1,t, z2,t, and z3,t for the next six days.


## 4.1. Data Descriptive

Descriptive statistics for each data obtained using Microsoft Excel software are presented in Table1.Table 1
Descriptive statistics of each data.
Descriptive statisticsz1,t (%)z2,t (%)z3,t (%)x1,t (%)Average7.811989.28122.90697.8703Variance0.26780.24700.00080.6484Standard deviation5.17504.97030.28078.0520Table1 shows the intensity of active, recovery, death, and vaccination rates from COVID-19 in Indonesia per day, each of which is different. The intensity values are 7.8119%, 89.2812%, 2.9069%, and 7.8703% per day, respectively. The intensity values of the active and death cases appear to be lower than that of vaccination. It indicates that the number of Indonesian citizens exposed and died from COVID-19 is lower than the number of people vaccinated daily. It shows that the Indonesian government is serious about handling COVID-19 via vaccination. Then, there are deviations from the active, recovery, death, and vaccination rates from COVID-19 against each average daily. On average, the deviation values from the active, recovery, death, and vaccination rates from COVID-19 were 5.1750%, 4.9703%, 0.2807%, and 8.0520% per day, respectively. The daily deviations of the active and death rates are lower than the deviation of the vaccination rate. It indicates that Indonesian citizens are increasingly varied of COVID-19 with a high willingness to vaccinate. It prevents extreme value deviations from occurring. Finally, the variance of z1,t, z2,t, z3,t, and x1,t is 0.2678%, 0.2470%, 0.0008%, and 0.6484%, respectively. The interpretation of the variance values is similar to the interpretation of the deviation value. Then, the correlation rates between z1,t, z2,t, and z3,t obtained using Microsoft Excel software are presented in Table 2.Table 2
The correlation rates betweenz1,t, z2,t, and z3,t.
Data pairsThe correlation ratesz1,t and x1,t−0.6599z2,t and x1,t0.6352Table2 displays the correlation rate between z1,t and z2,t, −0.9993, which is the highest correlation rate among the correlation rates of the other indicator data pairs. It is not surprising because the two are rationally the opposite of each other. If the active case rate increases, the recovered case rate decreases, and vice versa. This fact can be seen visually in Figure 1. Therefore, z1,t and z2,t are selected as endogenous variables. The high negative correlation between z1,t and z2,t is not a problem in VARX modelling. It is the same as regression modelling, where if the correlation rate between variables is high, the model formed can provide an excellent description of the data.Figure 1
Graph ofz1,t (a), z2,t (b), and x1,t (c).
(a)(b)(c)Table3 captures the correlation rates between endogenous and exogenous variables above 0.6. It indicates that endogenous and exogenous variables are strong enough.Table 3
The correlation rates between endogenous and exogenous variables.
Data pairsThe correlation ratesz1,t and z2,t−0.9993z1,t and z3,t−0.7421z2,t and z3,t0.7162


## 4.2. Data Stationary Test

In this section,z1,t, z2,t, and x1,t are checked for stationary visually and formally. Visually, the stationarity of the three data is seen from each graph, while formally, the stationary of the data is tested using the ADF test. Each data graph is observed whether it is around the average line or not. If the data graph is around the average line, then the data are stationary, and vice versa. The graphs of z1,t, z2,t, and x1,t with their respective average lines are presented in Figure 1.Figure1 displays that the graphs of the three data do not appear to be around the average line. Therefore, visually, the three data are not stationary. Then, the formal stationarity test of the three data is carried out using the ADF test. The significance level and optimal lag used are 0.05 and 6, respectively. The test statistics of each data are determined using equation (5). The formal stationary test result of the three data using the ADF test is presented in Table 4.Table 4
Data stationary testz1,t, z2,t, and x1,t.
Dataz1,tz2,tx1,tτ−2.7926−2.78152.4323Critical value−3.2000−3.2000−3.2000Based on Table4, all data have a value of not less than the critical value. Therefore, H0 on each data test is not rejected, which means that all the data are not stationary, so the differencing process using equation (3) is carried out. After the data are differentiated, the next step is the three data’s visual and formal stationarity tests. After the differencing process is carried out twice, all data are stationary. The graphs of z1,t″, z2,t″, and x1,t″ with their respective average lines are presented in Figure 2.Figure 2
Graph ofz1,t″ (a), z2,t″ (b), and x1,t″ (c).
(a)(b)(c)Figure2 shows the graphs of three data differentiated twice and appear to be around the average line. Therefore, visually, they are stationary. Next, the formal stationarity test is carried out using the ADF test. The significance level and optimal lag used are 0.05 and 6, respectively. The test statistics of each data are determined using equation (5). Furthermore, a formal summary of the three data stationary tests using the ADF test is presented in Table 5.Table 5
Data stationary testz1,t″, z2,t″, and x1,t″.
Dataz1,t″z2,t″x1,t″τ−8.2813−8.3078−13.7290Critical value−3.2000−3.2000−3.2000Table5 captures the three data differentiated twice that have a τ value less than the critical value. Therefore, H0 of each test is rejected, which means that all data are stationary.


## 4.3. Causality Relationship Test

A check of the existence of a bidirectional causal relationship betweenz1,t″ and z2,t″ in this study is carried out using the Granger causality test. The significance levels, optimal lag, and residual degrees of freedom used are 0.05, 7, and 256, respectively. The effect of z1,t−1″ and z2,t−i″ on z1,t″ is tested first. The F value obtained is 2.6931. This value is greater than the critical value, 2.0442. Therefore, H0 is rejected, which means that z1,t″ is affected by z1,t−1″ and z2,t−i″. Next, the effect of z1,t−1″ and z2,t−i″ on z2,t″ is tested. The F value obtained is 2.7514. This value is also greater than the critical value, 2.0442. Therefore, H0 is rejected which means that z2,t″ is affected by z1,t−1″ and z2,t−i″. Because z1,t″ is affected by z1,t−1″ and z2,t−i″, and z2,t″ is affected by z1,t−1″ and z2,t−i″, then z1,t″ and z2,t″ have a bidirectional causality relationship.


## 4.4. VARX Order Identification

The AIC value, as described in Section2, is used as a criterion to determine the best order of endogenous and exogenous variables of the VARX model. The model with the smallest AIC value is the most suitable. The AIC value for each model is determined by equation (3). The maximum order of the endogenous and exogenous variables considered is 10. The AIC value for each model is determined using the R software. Based on the results obtained from the software, the AIC value of the VARX (7, 1), −34.2194, is the smallest AIC value among other models. Therefore, the most suitable model chosen is the VARX (7, 1).


## 4.5. VARX (7, 1) Parameter Estimation Using ML Method

The VARX (7, 1) parameters are estimated using the ML method, equation (14), as described in Section 2. The estimation process is carried out using the help of Microsoft Excel software. The parameter estimation results are presented in Table 6.Table 6
Results of parameter estimates of VARX (7, 1).
ParameterEstimatorμ11.745E − 05μ2−1.817E – 05ϕ1,110.5495ϕ1,210.7931ϕ2,11−1.0376ϕ2,21−1.2881ϕ1,121.2891ϕ1,221.4981ϕ2,12−1.5079ϕ2,22−1.7207ϕ1,135.3873ϕ1,235.6679ϕ2,13−5.5021ϕ2,23−5.7843ϕ1,145.6900ϕ1,246.0157ϕ2,14−5.8682ϕ2,24−6.1934ϕ1,158.9929ϕ1,259.4176ϕ2,15−9.0363ϕ2,25−9.4565ϕ1,163.2196ϕ1,263.3846ϕ2,16−3.3553ϕ2,26−3.5186ϕ1,170.7175ϕ1,270.6588ϕ2,17−0.8463ϕ2,27−0.7909θ1,10−0.0117θ2,100.0104θ1,110.2203θ2,11−0.2162On the basis of the parameters obtained in Table6, the VARX (7, 1) model is stated in equation (16).(18)z^1,t″z^2,t″=1.7448E−05−1.8174E−05+0.54950.7931−1.0376−1.2881z1,t−1″z2,t−1″+1.28911.4981−1.5079−1.7207z1,t−2″z2,t−2″+5.38735.6679−5.5021−5.7843z1,t−3″z2,t−3″+5.69006.0157−5.8682−6.1934z1,t−4″z2,t−4″+8.99299.4176−9.0363−9.4565z1,t−5″z2,t−5″+3.21963.3846−3.3553−3.5186z1,t−6″z2,t−6″+0.71750.6588−0.8463−0.7909z1,t−7″z2,t−7″+−0.01170.0104x1,t″+0.2203−0.2162x1,t−1″.In order to be used to forecast the actual data, equation (18) is transformed to the actual data as follows:(19)z^1,tz^2,t=1.7448E−05−1.8174E−052.54950.7931−1.03760.7119z1,t−1z2,t−1+−0.8099−0.08820.5673−0.1446z1,t−2z2,t−2+3.35863.4649−3.5239−3.6309z1,t−3z2,t−3+−3.7954−3.82203.62813.6545z1,t−4z2,t−4+3.0003.0540−2.8021−2.8541z1,t−5z2,t−5+−9.0762−9.43498.84919.2011z1,t−6z2,t−6+3.27133.3072−3.1719−3.2103z1,t−7z2,t−7+1.78452.0670−1.6627−1.9367z1,t−8z2,t−8+0.71790.6588−0.8463−0.7909z1,t−9z2,t−9+−0.01170.0104x1,t+0.2437−0.2370x1,t−1+−0.45230.4428x1,t−2+0.2203−0.2162x1,t−3.Meanwhile,z3,t can be modelled as follows:(20)z^3,t=1−z^1,t−z^2,t.


## 4.6. Diagnostic Test of VARX (7, 1)

Independence of bothe1,t and e2,t to time is carried out using the Ljung-Box test described in section 2. The significance levels and the degrees of freedom used are 0.05 and 7, respectively. The test statistic is determined using equation (18). The test statistics of e1,t and e2,t obtained are 3.3610 and 3.4270, respectively. The test statistic values are smaller than the critical value, 14.0671. Therefore, H0 of each independence test of e1,t and e2,t to times is not rejected, which means that e1,t and e2,t are independent of time. Next, the normality test of the VARX (7, 1) model error is visually carried out first through the Q-Q plot. Q-Q plots for e1,t and e2,t are presented in Figure 3.Figure 3
Q-Q plots of e1,t (a) and e2,t (b).
(a)(b)Figure3 indicates the scattering of points in each plot appears to be around a line with one gradient. Therefore, e1,t and e2,t visually follow a normal distribution with zero mean and constant variances, namely 0.00108 and 0.00107. Based on this, the model equation (20) errors fulfill the assumption that they are independent and identically normally distributed with zero mean and constant variance. The next step is checking the error size of the model. This study measures the error size using the mean absolute percentage error (MAPE). The MAPE value of z1,t and z2,t is 1.5751% and 0.0833%, respectively. The MAPE value is less than 10%, even close to zero. Therefore, according to Wei [10], the VARX (7, 1) model is very feasible to forecast z1,t and z2,t in the next day based on the feasibility criteria of the MAPE-based model.The feasibility of the VARX (7, 1) model is also checked using the visual comparison between the forecast results and the actual data. The visual comparison between the forecast results and the actual data is presented in Figure4.Figure 4
The visual comparison between thez1,t and z1,t forecast (a) and z2,t and z2,t forecast (b).
(a)(b)Figure4 shows that the forecast line of z1,t and z2,t appears very close to the actual data of z1,t and z2,t. It means that the VARX (7, 1) model is visually feasible to forecast z1,t and z2,t in the next day. The VARX (7, 1) model is also practically forecast z1,t and z2,t in the next few days.


## 4.7. Forecasting

Forecastingz1,t, z2,t, and z3,t are carried out on six days, from 8 November 2021 until 13 November 2021. These forecasts are practically carried out using equation (19) for z1,t and z2,t and equation (20) for z3,t. The second-dose COVID-19 vaccination rates used for this forecast are the same as the actual on that date, 29.5200%, 29.7817%, 30.1154%, 30.4491%, 30.7828%, and 31.0465%. The z1,t, z2,t, and z3,t forecasting results obtained are presented in Table 7.Table 7
Thez1,t, z2,t, and z3,t forecasting results.
Datez1,t forecast (%)z2,t forecast (%)z3,t forecast (%)8 November 20210.227596.39303.37959 November 20210.235096.38623.378810 November 20210.224796.39603.379311 November 20210.238996.38223.378812 November 20210.222096.39883.379113 November 20210.216596.40423.3793Table7 shows that the forecasting results of z1,t, z2,t, and z3,t from 8 November 2021 until 13 November 2021 tend to decrease. To determine the accuracy of the forecast results presented in Table 7, we will look at the mean absolute percentage error (MAPE) size. The MAPE of z1,t, z2,t, and z3,t forecasting results on 8 November 2021 until 13 November 2021 is 2.2853%, 0.0051%, and 0.0082%, respectively. This MAPE of z1,t, z2,t, and z3,t forecasting results is tiny. Therefore, equations (19) and (20) are suitable for forecasting z1,t, z2,t, and z3,t for the next six days.


## 5. Discussions

First, an analysis of the influence of the second-dose COVID-19 vaccination rates on active and recovered case rates of COVID-19 is discussed. Generally, equation (19) shows that the second-dose COVID-19 vaccination rate in the previous three days influences active and recovered case rates of COVID-19 today. It does not always decrease active COVID-19 case rates but increases active COVID-19 case rates. It is also an increase and decrease in the COVID-19 recovered case rate. The active case rates of COVID-19 today are −1.17% of second-dose COVID-19 vaccination rate today, 24.37% of second-dose COVID-19 vaccination rate yesterday, −45.23% of second-dose COVID-19 vaccination rates the previous two days, and 22.03% of second-dose COVID-19 vaccination rate previous three days. A negative percentage means that the second-dose COVID-19 vaccination rate will decrease today’s COVID-19 active case rate. In contrast, a positive percentage means that the second-dose COVID-19 vaccination rate will increase today’s COVID-19 active case rate.Meanwhile, the recovered case rates of COVID-19 today are 1.04% of the second-dose COVID-19 vaccination rate today, −23.70% of the second-dose COVID-19 vaccination rate yesterday, 44.28% of the second-dose COVID-19 vaccination rate previous two days, and −21.62% of second-dose COVID-19 vaccination rate previous three days. A positive percentage means that the second-dose COVID-19 vaccination rate will increase today’s COVID-19 recovered case rate. In contrast, a negative percentage means that the second-dose COVID-19 vaccination rate will reduce COVID-19 recovered case rate today. These facts do not mean that vaccination in Indonesia is not adequate. Suppose the percentage of the effect of today’s vaccination up to the previous three days on the active and recovered case rates is added. In that case, the results have negative and positive values, respectively. Negative values indicate that vaccination reduces the daily active case rate, and positive values indicate that vaccination generally increases the daily recovery case rate. Next is an analysis of the influence of active and recovered case rates of COVID-19 on itself. In general, equation (19) shows that the active and recovered rates of COVID-19 today are influenced by themselves on the previous three days. The percentage of each day can be seen in equation (19).The last is an analysis of the effect of second-dose COVID-19 vaccination rates on 8 November 2021 on active and recovered case rates of COVID-19. Suppose the second-dose COVID-19 vaccination rates on 8 November 2021 are assumed not to increase from the previous day. Thus, the forecasted active and recovered case rates of COVID-19 on 8 November 2021 are 0.2321% and 96.3886%, respectively. The forecast value of the COVID-19 active case rates is greater than the forecast value obtained in the forecasting section. The forecast value of the COVID-19 recovered case rates is also smaller than that obtained in the forecasting section. It shows that if the second-dose COVID-19 vaccination rates on 8 November 2021 are not increased, the decline in the COVID-19 spread rate in Indonesia will be slower. Therefore, although the second-dose COVID-19 vaccination rates do not always decrease the COVID-19 active case rates and increase the COVID-19 recovered case rates every day, the second-dose COVID-19 vaccination rates on 8 November 2021 must still be increased. It also applies from 9 November 2021 until 13 November 2021.


## 6. Conclusion

The VARX model (7, 1) is the best VARX model to describe a causal relationship between active and recovered case rates of COVID-19 in Indonesia influenced by the exogenous variable, second-dose COVID-19 vaccination rates. The assumption of white noise is fulfilled, and the MAPE of this model is tiny, so this model is very suitable to be used in the forecasting process the next day. Forecasting results in the six days also give very close results to the actual data (even if the absolute percentage error is almost zero). There is an oddity in this VARX (7.1) model. The second-dose COVID-19 vaccination rates carried out daily should always reduce COVID-19 active case rates and increase COVID-19 recovered case rates.However, in this model, the impact of the second-dose COVID-19 vaccination rates carried out in the previous three days on COVID-19 active case rates can decrease and increase simultaneously. It also happened to the impact of the second-dose COVID-19 vaccination rates on the COVID-19 recovered case rates. The second-dose COVID-19 vaccination rates carried out in the previous three days on COVID-19 recovered case rates can increase and decrease simultaneously. However, it does not indicate that vaccination does not affect daily activities and recovery case rates. The sum of the percentage effect of today’s vaccination to the previous three days on the active and recovered case rates has negative and positive values, respectively. These negative and positive values indicate that vaccination generally decreases the daily active case rate and increases the daily recovery case rate.Another analysis result obtained is that if the second-dose COVID-19 vaccination rates do not increase from the previous day, COVID-19 active case rates will only decrease slightly, and COVID-19 recovered case rates will only increase slightly. It has caused the COVID-19 spread in Indonesia to decline more slowly. Based on this, the vaccination rate must still be increased. With this model, the Indonesian government can project the second-dose COVID-19 vaccination rates that must be carried out today so that active and recovered case rates of COVID-19 can be achieved by that date.



---

*Source: 6285328-2023-03-14.xml*