Selection of model

The ODEs in the SIR model are discretised using the finite difference scheme. There are numerous works with mainly explicit schemes regarding time-discrete SIR models in literature [Reference Bohner, Streipert and Torres21, Reference Brauer24]; however Allen [Reference Allen20], and Wacker and Schluter [Reference Wacker and Schlüter22] proposed an implicit time-discrete edition of this classical SIR model and showed that this time-discrete variant maintained various time-continuous properties. In this study, we follow the implicit numerical algorithm in [Reference Wacker and Schlüter22] and discretise the SIR model as

(3)$$\left\{\matrix{\displaystyle{{S_{i + 1}-S_i} \over {t_{i + 1}-t_i}} = \displaystyle{{-\alpha_{i + 1}} \over N}S_{i + 1}I_{i + 1} \hfill \cr \displaystyle{{I_{i + 1}-I_i} \over {t_{i + 1}-t_i}} = \displaystyle{{\alpha_i + 1} \over N}S_{i + 1}I_{i + 1}-\beta_{i + 1}I_{i + 1} \hfill \cr \displaystyle{{R_{i + 1}-R_i} \over {t_{i + 1}-t_i}} = \beta_{i + 1}I_{i + 1}. \hfill} \right.$$

Assuming that our time interval [0, T] can be divided by a strictly increasing sequence $\{ {t_i} \} _{i = 1}^M$ with M ∈ ℕ subintervals, Eq. 3 gives a fully implicit structure of the time-continuous SIR model for all i ∈ (1,  2,  … M − 1) and N = S _i+1 + I _i+1 + R _i+1 = S _i + I _i + R _i [Reference Wacker and Schlüter22].

Selection and pre-treatment of COVID-19 data

The data for COVID-19 was acquired from the Fiji's Ministry of Health and Medical Services (MOHFiji) website (www.health.gov.fj) [6]. The data includes the cumulative number of infected cases ($\hat{I}$), the cumulative number of recovered cases ($\hat{R}$) and the cumulative number of death cases ($\hat{D}$) in Fiji. Following [Reference Wacker and Schlüter22], we define $R_i = \widehat{{R_i}} + \widehat{{D_i}}$ and $I_i = \widehat{{I_i}}-R_i.$ We considered the COVID-19 data of Fiji for the second wave of the pandemic from 4 May 2021 (t = 1) to 9 December 2021 (t = 220) of which the data corresponding to time $\{ {t_i} \} _{i = 1}^{160}$ is used for estimating the parameters of the discrete-time SIR model and the data for time $\{ {t_i} \} _{i = 161}^{220}$ is used to validate the estimated model. Figure 3a shows the cumulative cases of reported infected people and the cumulative cases of reported recovered people in Fiji. Apparently, there are two jumps in the plot of total recoveries around days 152 and 154. The corresponding jumps are subsequently visible in Figure 3b showing the number of the active cases. Possible explanations for these sudden jumps would be late testing, identification and recording of recovered patient numbers.

Fig. 3. Unprocessed (observed) COVID-19 data for Fiji from t ₁ (4 May 2021) to t ₁₆₀ (10 October 2021). (a) shows the cumulative infection and cumulative recovered cases (R). (b) shows the daily active cases (I) and cumulative recovered cases (R).

A 3-day moving average (3MA) and a 5-day moving average (5MA) filter was employed by Law et al. [Reference Law9] to smoothen daily infected data for Malaysia. Cartocci et al. [Reference Cartocci, Cevenini and Barbini25] in their study on analysing gender and age-grouped data for Italy used the 7-day moving average (7MA) method to remove noise and excess variability from data. The authors in [Reference Law9] and [Reference Cartocci, Cevenini and Barbini25] also highlighted the importance of choosing a reasonable smoothing moving average window as a wider window would allow for a higher degree of noise filter, but meaningful rapid variations of the pandemic data would not be obtained. A 7-day window for the moving average filter is a fair compromise, and widely acceptable for smoothening noisy and erratic for the COVID-19 data as the mean incubation period of the virus is also 7 days [Reference McAloon26]. For this study, a 7MA was applied to the recovered data, $\hat{R}$, to eliminate sudden jumps and thus to obtain a smooth sequential data. The processed data using 7MA is shown in Figure 4, showing the improvement in smoothness and continuity in the recovered and hence the active cases.

Fig. 4. (a) Smoothening the recovery (R) curve with a 7MA filter, and (b) subsequently improving the curve for active cases (I).

Calculation of time-varying α and β from observed data

We follow Wacker and Schlüter [Reference Wacker and Schlüter22] and summarise the algorithm used to estimate the time-varying model parameters and the procedure adopted to forecast the COVID-19 cases in Fiji. The time-varying transmission and recovery rates are determined from Eq. (3) using discrete values of I and R observed at time t _i for i = 1,  …,  M − 1 where M = 160. Assuming S _i+1 > 0 and I _i+1 > 0, the time-varying parameters are computed as

(4)$$\alpha _{i + 1} = \displaystyle{{N( S_i-S_{i + 1}) } \over {I_{i + 1}\cdot S_{i + 1}\cdot \Delta _{i + 1}}}$$

and

(5)$$\beta _{i + 1} = \displaystyle{{R_{i + 1}-R_i} \over {I_{i + 1}\cdot \Delta _{i + 1}}}, \;$$

where Δ_i+1 = t _i+1 − t _i = 1 for i = 1,  …. M − 1.

After the time-varying model parameters are computed, the implicit time-discrete solution of the SIR model is estimated. Assuming that 0 < α _i < 1 and 0 < β _i < 1 are known for all i ∈ (1,  2,  … M − 1) and that the initial values of S ₁ > 0, I ₁ > 0 and R ₁ ≥ 0 are known, an implicit solution scheme of Eq. (3) reads

(6)$$\left\{\matrix{S_{i + 1} = \displaystyle{{S_i} \over {1 + \alpha_{i + 1}\cdot \Delta_{i + 1}\cdot \displaystyle{{I_{i + 1}} \over N}}} \hfill \cr I_{i + 1} = \displaystyle{{I_i} \over {1 + \beta_{i + 1}\cdot \Delta_{i + 1}-\alpha_{i + 1}\cdot \Delta_{i + 1}\cdot \displaystyle{{S_{i + 1}} \over N}}} \hfill \cr R_{i + 1} = R_j + \beta_{i + 1}\cdot \Delta_{i + 1}\cdot I_{i + 1}. \hfill} \right.$$

Assuming that S _i > 0, I _i > 0 and R _i ≥ 0, Eq. 6 is uniquely solvable for all i ∈ (1,  2,  … M − 1) [Reference Wacker and Schlüter22]. Here, I _i+1 for i = 1,  …,  M − 1 is first calculated using

(7)$$I_{i + 1} = {-}\displaystyle{{\Phi _{i + 1}} \over {\Theta _{i + 1}}} + \sqrt {\displaystyle{{\Phi _{i + 1}^2 } \over {\Theta _{i + 1}^2 }} + \displaystyle{{N\cdot I_i} \over {\Theta _{i + 1}}}} , \;$$

where Θ and Φ are defined as

(8)$$\Theta _{i + 1}: = ( 1 + \beta _{i + 1}\cdot \Delta _{i + 1}) \cdot ( \alpha _{i + 1}\cdot \Delta _{i + 1}) $$

and

(9)$$\Phi _{i + 1}: = \displaystyle{{( 1 + \beta _{i + 1}\cdot \Delta _{i + 1}) \cdot N-\alpha _{i + 1}\cdot \Delta _{i + 1}\cdot ( S_i + I_i) } \over 2}.$$

The solution scheme for S _i+1 and R _i+1 follows after computing I _i+1. Note that α _i+1 ≠ 0, which implies that Δ_S(t) = S _i+1 − S _i ≠ 0 from (4).

The resulting implicit time-discrete scheme, i.e. Eq. 6, is the forward difference approximation of the SIR model, i.e. Eq. 3. The global existence, global uniqueness, non-negativity and boundedness of the solution, monotonicity properties, and error analysis of the implicit scheme have been well established. We refer the readers to [Reference Wacker and Schlüter22] for a detailed analysis of the above properties. Additionally, the numerical solution implicit scheme is uniquely solvable for all time steps [Reference Wacker and Schlüter22]. Moreover, the implicit time-discrete scheme is a rewritten version of the well-known implicit Eulerian time-stepping scheme, which is known to be unconditionally stable [Reference Wanner and Hairer27]. This method was implemented in [Reference Wacker and Schlüter22] to successfully model the spread of COVID-19 in Germany and Iran.

Model validation and short term forecast

To validate the model for $\{ {t_i} \} _{i = 161}^{220}$ and make a short-term forecast of COVID-19 cases in Fiji, time-varying model parameters, for t _i (i > 220), are estimated following the approach in [Reference Wacker and Schlüter22]. Through inspection of the time-varying transmission rate as shown in Figure 5a, it can be inferred that the model follows an exponential decay function. Similarly, except for a few areas of instability, as explained later in Section ‘Results’, the recovery rates are mostly constant for the duration of the study as seen in Figure 5b as well. Similar observations were made for COVID-19 cases in Germany and Iran [Reference Wacker and Schlüter22], and Bulgaria [Reference Margenov, Pasheva, Popivanov and Venkov28]. Hence, the time-varying transmission and recovery rate is assumed to take the following form:

(10)$$\alpha ( t) : = \alpha _1\cdot e^{-\alpha _2.t}$$

and

(11)$$\beta ( t) : = \beta , \;$$

for t ≥ 1. The real constants α ₁, α ₂ and β which is are determined using the time-varying parameters $\{ {\alpha_i} \} _{i = 2}^M$and $\{ {\beta_i} \} _{i = 2}^M$ as defined by Eqs. (4) and (5).

Fig. 5. Time-varying transmission and recovery rates from processed data for Fiji. Superimposed are the estimator functions. (a) The estimated parameters are α ₁ ≈ 0.2144 and α ₂ ≈ 0.0210. (b) The estimated recovery rate is β ≈ 0.0403 (the mean value on the full interval).

Since both $\{ {\alpha_i} \} _{i = 2}^M > 0$ and $\{ {\beta_i} \} _{i = 2}^M > 0$ for i = 2,  … M, it is assumed that constants α ₁ > 0 and β > 0. The non linear relation in Eq. (10) is linearised using the parametric transformation

(12)$$\ln ( \alpha ( t) ) = \ln ( \alpha _1) -\alpha _2\cdot t = \delta _1 + \delta _2\cdot t, \;$$

where α ₁: = ln(δ ₁) and α ₂: = −δ ₂ as in the case of maximum log-likelihood estimation (MLE). MLEs draw conclusions about the population most likely to have generated a sample, especially the joint probability distribution of the random variables {y ₁,  y ₂,  y ₃,  …,  y _n}. This method of parameter estimation's many optimal properties and working algorithm have been thoroughly discussed in [Reference Pan, Fang, Bühlmann, Diggle, Gather and Zeger29]. To find suitable estimators $\widehat{{\alpha _1}}$, $\widehat{{\alpha _2}}$ and $\hat{\beta }$ for α ₁, α ₂ and β, respectively, a cost function ${\rm \jmath }\colon {\opf R}^3\to [ 0, \;\,\infty )$, is defined as:

(13)$${\rm \jmath }( \delta _1, \;\,\delta _2, \;\,\beta ) : = \sum\limits_{i = 2}^M {[ {\delta_1 + \delta_2\cdot t_i-\ln {( \alpha_i) }^2} ] } + \sum\limits_{i = 2}^M {{( \beta -\beta _i) }^2} .$$

The solution rests in showing that the function ${\rm \jmath }$ possesses a unique local minimiser $\widehat{{\alpha _1}}$, $\widehat{{\alpha _2}}$ and $\hat{\beta }$. The local minimisers are obtained by setting all the partial derivatives of the cost function equal to zero. Hence, setting ${{\partial {\rm \jmath }} \over {\partial \beta }}( \widehat{{\delta _1}}, \;\,\widehat{{\delta _2}}, \;\,\hat{\beta }) ={{\partial {\rm \jmath }} \over {\partial \delta _1}}( \widehat{{\delta _1}}, \;\,\widehat{{\delta _2}}, \;\,\hat{\beta }) =$ $ {{\partial {\rm \jmath }} \over {\partial \delta _2}}( \widehat{{\delta _1}}, \;\,\widehat{{\delta _2}}, \;\,\hat{\beta }) = 0$ yields [Reference Wacker and Schlüter22]

(14)$$\hat{\beta } = \displaystyle{1 \over {M-1}}\sum\limits_{i = 2}^M {\beta _i} , \;$$

(15)$$\widehat{{\delta _2}} = \displaystyle{{\sum\nolimits_{i = 2}^M {t_i\cdot \ln ( \alpha _i) } -\displaystyle{1 \over {M-1}}\left[{\sum\nolimits_{i = 2}^M {\ln ( \alpha_i) } } \right]\cdot \sum\nolimits_{i = 2}^M {t_i} .} \over {\sum\nolimits_{i = 2}^M {t_i^2 } -\displaystyle{1 \over {M-1}}\cdot {\left({\sum\nolimits_{i = 2}^M {t_i} } \right)}^2}}, \;$$

and

(16)$$\widehat{{\delta _1}} = \displaystyle{1 \over {M-1}}\sum\limits_{i = 2}^M {( \ln ( } \alpha _i) -t_i\cdot \widehat{{\delta _2}}) .$$

Using the local minimisers (Eqs. (14)–(16)) and the estimated time-varying model parameters (Eqs. (10) and (11)), the estimated model is validated for $\{ {t_i} \} _{i = 161}^{220}$. In the absence of an absolute standard method in forecast verification, the model's performance has been done in line with Lal et al. [Reference Lal, Huang and Li13] and Salgotra et al. [Reference Salgotra, Gandomi and Gandomi30]. The model's accuracy was evaluated by computing the Relative Mean Absolute Error (RMAE) of the simulated model state as

(17)$${\rm RMAE} = \displaystyle{1 \over N}\mathop {\mathop \sum \nolimits^ }\limits_N^{i = 1} \displaystyle{{\vert {y( i ) -s( i ) } \vert } \over {y( i ) }}, \;$$

where y(i) is the observed state, s(i) is the simulated state using the modelled parameters, and N is the size of the observed data. Willmott and Matsuura [Reference Willmott and Matsuura31] indicate that RMAE is a more natural and accurate measure of average error, and (unlike RMSE or related measures) is unambiguous. It can be further stated that for a model to be reliable and accurate, the correlation coefficient between the desired and projected values must be strong [Reference Chicco, Warrens and Jurman32]. Hence the coefficient of determination, r ², values are calculated, using [Reference Lal, Huang and Li13, Reference Chicco, Warrens and Jurman32]

(18)$$r^2 = 1-\left[{\displaystyle{{\sum\nolimits_i^N {{( y( i) -s( i) ) }^2} } \over {\sum\nolimits_i^N {{( y( i) -\bar{y}) }^2} }}} \right].$$

Furthermore, short term (30 days) forecast is made for the cumulative number of infected and recovered COVID-19 cases in Fiji. In this study, all computations and simulations were performed using MATLAB R2016a software.

Time-varying basic reproduction number

The basic reproduction number $\Re _0$ of an infectious disease is the number of people who contract the disease from an infected person assuming the whole population is susceptible [Reference Allen20, Reference Arroyo-Marioli33]. In the SIR model, the time-varying basic reproduction number is computed as

(19)$$\Re _0( t_c) : = \displaystyle{{\alpha ( t_c) } \over {\beta ( t_c) }}$$

for arbitrary c ∈ (1,  2,  … M) and assuming β(t _c) > 0 [Reference Arroyo-Marioli33].