Modeling COVID-19 pandemic with financial markets models: The case of Jaén (Spain)

1.   Introduction

Coronavirus disease 2019 (COVID-19) is a lung disease caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2). In December 2019, the Chinese authorities reported different cases of this virus in Wuhan. This disease spread rapidly throughout the world from less than 30 cases at the end of December 2019 to more than 8, 455, 738 confirmed cases on June 20, 2021.

The first case in Spain was a German tourist, on January 31, 2021. Since then, several cases were confirmed throughout the country. In Andalusia, the first positive case was reported in Seville on February 26, 2020. Two days later the first case was confirmed in the province of Jaén, and on March 6 the first case of coronavirus in the city of Jaén.

As the days went by, due to the high number of infections in the country, the state of alarm was decreed and a national lockdown was imposed on March 14, which was extended until June 21. The application of measures such as the use of a mask, perimeter confinements, the closure of non-essential services, etc., reduced the infection rate two weeks after the declaration of quarantine. During these years, many researchers from various disciplines have used various modelling tools to analyze the impact of the pandemic at the global and local levels. In our case, we are going to focus on modelling the contagion in Jaén in two different ways. The first approach is based on the Autoregressive Integrated Moving Average with Explanatory Variable (ARIMAX) model ^[1] which, we found, provides better performances than the Autoregressive Integrated Moving Average (ARIMA) (on the same line, see also ^[2,3,4,5,6]). This is a common model in time series forecasting and is often adopted in finance ^[7,8,9,10].

The second approach is based on the Cox, Ingersoll and Ross (CIR) model. This is a model designed for interest rate pricing that we turn into a forecasting tool. The advantage of this model lies in its simplicity (being a single factor model) and analytical tractability ^[11]. The CIR model is mean reverting and is a squared Ornstein–Uhlenbeck process ^[12]. Notice that the "Ornstein-Uhlenbeck process is a natural model to consider in a biological context because it stabilizes around some equilibrium point. This corresponds to the homeostasis often observed in biology" ^[13], and the quasi-stationarity is relevant with reference to mortality plateaus ^[13]. Furthermore, we point out a connection to models for short-term interest rates in financial modelling.

We prove that the proposed transformation of the CIR model, which we have denoted CIR*, outperforms the classical ARIMAX. For this purpose both the Euler-Maruyama method and the Milstein method were used, and this is one of the main contributions of the present study. Notice that the suggested approach not only extends the models available to scholars to model pandemics but, also, paves the way for similar approaches where financial models can be converted into econometric models.

For the implementation of the models on real-world data, we use the data from the moving averages for 14 days of the daily cases of the city of Jaén, that is, each value represents the average number of people infected in each of the previous 14 days. This is because, due to the weekend effect and occasional misreporting, we found that the moving average is a more reliable target. This is due to two reasons: a) the relatively small size of the city of Jaén, which affects the number of cases, and b) the effect of the weekend, when reporting is altered. The two effects lead to highly irregular behaviour in the time series considered. Also, the 14-day moving average is the standard way of reporting data in Spain during the pandemic, both for publishing data and deciding the different measures (lockdowns, mobility restrictions, etc.). Over that time series, we estimate the parameters of the considered models, in such a way that they best fit the data. The forecasts have been validated with 95% confidence intervals and with statistical measures of goodness of fit, such as the RMSE.

This article is organized as follows. Section 2 briefly summarizes the literature. Section 3 reports the data and describes the CIR model, as well as the methodology on which it is based. That is followed by an explanation of the suggested adaptation to forecasting, its calibration and a sketch of the ARIMAX model. Section 4 shows the obtained results of the two models by comparing them. Section 5 concludes.

2.   Literature review

Among those works that adopted the ARIMA (and the like models) to estimate the cases of the COVID-19 pandemic, we mention Ekinci ^[14], who considered data from USA, India, Brazil, France, Russia, UK, Italy, Spain and Germany. When comparing ARMA-GARCH, ARMA-TGARCH and ARMA-EGARCH models, it was found that while considering the conditional variance effect improves the forecasting power, the asymmetric effect (such as asymmetric GARCH models) has mixed results. Sahai et al. ^[15] adopted the ARIMA for analyzing the trend of COVID-19 cases in Spain, Italy, France, Germany and the US. The authors claim that their model provides considerable forecast accuracy and could be a useful tool for governments to ramp up their healthcare preparations. Subramaniam et al. ^[8] drew a parallel between forecasting stock prices and cases of the pandemic by means of the ARIMA model. Then, they explored the correlation between the predictive efficiency of the ARIMA model and variation in the data. Katoch et al. ^[16] adopted an ARIMA model to analyzed the temporal dynamics of the COVID-19 outbreak in India from 30 January 2020 to 16 September 2020. Their approach suggests "varying epidemic's inflexion point and final size for underlying states and the mainland, India". Regarding the alternative between ARIMA and ARIMAX, in the literature, it has been found that the latter may yield better forecast compared to the seasonal ARIMA (SARIMA) model and Neural Networks (e.g., see Suhartono ^[17]). This is because the ARIMAX is most suited to deal with calendar effects ^[2,3,4,5,18].

As regards the Cox, Ingersoll and Ross (CIR) model ^[19,20], as already mentioned, it has been proposed for the pricing of interest rates. At the time of its introduction, it quickly gained popularity in finance because it was perceived as "an improvement on the Vasicek model ^[21], not allowing for negative rates and introducing rate-dependent volatility, as well as for its relatively handy implementation and analytical tractability" ^[11]. The CIR process is a variant of a squared Bessel process and, in the physical community, its analogues were used both to enrich the list of modified stochastic processes used for the description of various time series ^[22,23] as well as to model some financial historical data ^[24].

Other applications of the CIR model include stochastic volatility modelling in option pricing problems ^[25,26], FX ^[27] or default intensities in credit risk ^[28,29]. As mentioned, the Ornstein–Uhlenbeck (OU) and the Cox–Ingersoll–Ross (CIR) processes are strongly linked, and the reflected OU (ROU) ^[30] is used in queuing theory ^[30,31], in population dynamics ^[32], catastrophes ^[33,34], etc. In this study, similar to what has been done by Orlando et al. ^{[11,35,36,37,38]} when developing the CIR# model, we transform the original CIR model into a forecasting tool and compare its performance with that of the well-known ARIMAX model.

3.   Materials and methods

3.1.   Data

The available data are the moving averages of 14 days of infections in the city of Jaén, from February 2, 2020 to October 8, 2021 and have been provided by the Health and Family Council of the Andalusian Regional Government (see Figure 1). Jaén is a relatively small city (around 110, 000 inhabitants), which represents a challenge since the number of COVID-19 daily cases is small and with large fluctuations, even after the moving averages are computed. Therefore, we shall have the opportunity of testing different methods in unfavourable circumstances.

A visual inspection employing the partial autocorrelation function (PACF) confirms that the data is correlated (see Figure 2).

3.2.   ARIMAX model

Leaving aside the cases where data show evidence of non-stationarity where an initial differencing removes the integrated (I) part, the ARMAX model can be described as:

with $y(t)$, dependent/output variable at time $t$; $n_a$, number of poles; $n_b$ number of zeroes plus 1; $n_c$, number of $c$ coefficients; $n_k$, dead time in the system. Moreover, $y(t-1) \ldots y(t-n_a)$ denotes the dependence between the current output and the previous outputs, $x(t-n_k) \ldots x(t-n_k-n_b+1)$ indicates the dependence between the current output and both the previous and delayed inputs, and $\epsilon(t)$ expresses a white-noise error.

The orders of the ARMAX model are given by the parameters $n_a, n_b$ and $n_c$. $n_k$ is the delay, and $q$ is the delay operator. The ARMAX in compact form can be written as

such that,

The ARIMAX model can be seen as a generalization of the ARIMA because adds to the structure above described an integrator in the white noise $\varepsilon(t)$ as follows:

3.2.1.   Estimation of the ARIMAX model

To estimate the ARIMAX model the following steps have been performed.

1) Ensure stationarity of the time series by conducting Augmented Dicky Fuller (ADF) Test.

2) Model identification, i.e., specification of the autoregressive (AR) and moving average (MA) terms with the help of the autocorrelation function (ACF) and partial autocorrelation function (PACF).

3) Parameter estimation according to Ljung ^[39] and related implementation in Matlab ^[40]. The best model is selected based on Akaike information criterion (AIC) values ^[41].

From now on, we refer to the ARMAX model (and not to the ARIMAX) assuming that the integration (I) has been removed.

3.3.   CIR* model

As mentioned, the CIR model emerged in 1985 from the hand of John C. Cox, Jonathan E. Ingersoll and Stephen A. Ross ^[19,20] as an improvement of the Vasicek model to prevent negative interest rates.

The CIR model is defined by the following stochastic differential equation (SDE):

Here, $\alpha, \mu$ and $\sigma$ are positive constants, $X(t)$ is the interest rate, $t$ is time, and $W_t$ denotes the standard Wiener process.

The parameters include the following:

● $\alpha (\mu-X_t)$ is the same factor as in Vasicek's model and represents a mean reversion term.

● The standard deviation factor $\sigma \sqrt{X_t}$ removes negative rates.

● $\sqrt{X_t}$ increases the standard deviation as the short-term rate increases.

This model can only have positive solutions since when the interest rate is $0$ it ends up being positive later on. Also, when it is low or close to 0, the standard deviation is close to 0.

The only solution to (3.4) is what is known as the CIR process. Integrating Eq (3.4):

and therefore

If we call $m_t = E[X_t|X_s]$ we have

whose solution is

So

and therefore

Thus, $E[X_t|X_s]-\mu$ has the same sign as $X_s-\mu$. In addition, if $\mu > 0$ and $\alpha > 0$, starting with $X_s > 0$ we conclude that $X_t > 0$.

Similarly, the variance is:

As stated at the beginning, the fundamental advantage of this model is that the solutions are nonnegative. However, the distribution of the CIR model is not Gaussian, which makes it difficult to analyze.

The density function is given by:

where

$I_q(\cdot)$ is a Bessel function of first type and order $q$:

where $\Gamma$ is the gamma function.

Let $z_t = 2cX_t$. Then, the conditional distribution of $z_t$ given $z_s$ is a non-central $chi$-squared distribution $\chi_d^2(2u)$, with $d = \frac{4\alpha\mu}{\sigma^2}$ degrees of freedom, and the non-centrality parameter is $\lambda = 2u$.

Therefore,

where:

Since $z_t = 2cX_t$, $X_t$ conditional on $X_s$ has the same distribution as $z_t/2c$ conditional on $z_s/2c$. So,

At this point, we are going to analyze the different behaviours that the deterministic solution of the CIR model equation can present in terms of the relations among the different parameters, which will be useful later to give an interpretation of the model parameters, although we know that the stochastic part would give oscillations with respect to said behaviour. We shall allow in this analysis for negative values of $\alpha$ since in certain regions of the data the calibrated values of $\alpha$ result in negative values. According to Eq (3.7) we distinguish two cases, depending on whether $X_t$ is greater or less than $\mu$ and in each case two subcases, depending on whether $\alpha$ is positive or negative:

1) If $X_t < \mu$:

● If $\alpha > 0$, $X_t$ approaches $\mu$ from below.

● If $\alpha < 0$, $X_t$ moves away from $\mu$ downwards.

2) If $X_t > \mu$:

● If $\alpha > 0$, $X_t$ approaches $\mu$ from above.

● If $\alpha < 0$, $X_t$ moves away from $\mu$ upwards.

3.3.1.   Estimation of the parameters

To approximate the data well, it is necessary to give a good adjustment of the parameters. In the case of the CIR* model, we must estimate three parameters, $\alpha$, $\mu$ and $\sigma$. We will generally refer to them as the parameter vector $\theta \equiv (\alpha, \mu, \sigma)$. The procedure that will be followed to estimate the parameters is the one shown in ^[42], and it is the maximum likelihood method (MLE), which is based on maximizing the objective function under consideration. For the maximum likelihood estimation of the parameter vector $\theta \equiv (\alpha, \mu, \sigma)$, the transition densities are required. The CIR process is one of the processes for which we know its density function explicitly. Given $X_t$ at time $t$, the density of $X_{t+\Delta t}$ at time $t+\Delta t$ is:

where

where $I_q(2\sqrt{uv})$ is a Bessel function.

The likelihood function for time series with $N$ observations is:

To simplify the calculations, it is usual to work with the log-likelihood expression, which consists of taking logarithms in the Eq (3.9).

from which the log-likelihood function of the CIR process can be easily derived.

where $u_{t_i} = cX_{t_i}e^{-\alpha\Delta t}$ y $v_{t_{i+1}} = cX_{t_{i+1}}$

To find the maximum likelihood estimate $\hat{\theta}$ of the parameter vector $\theta$, we have to maximize the function (3.11) over its parameter space.

Since the logarithm function is monotonically increasing, maximizing the log-likelihood function is equivalent to maximizing the likelihood function.

To solve the problem (3.12), we resort to numerical computation. For global optimal convergence, the initial optimization points are essential, for which the method of least squares will be used. We first write the equation of the discretized CIR*:

where $W_t$ is distributed with zero mean and variance $\Delta t$.

Dividing Eq (3.13) by $\sqrt{X_t}$, we get

Based on Eq (3.14), the initial values of $\hat{\alpha}$ and $\hat{\mu}$ are found by minimizing the function:

The exact expression of the solution is on page 3 of ^[42]. The estimate of $\hat{\sigma}$ is the standard deviation of the residuals.

To optimize the objective function (3.11), we need to evaluate the Bessel function $I_q(2\sqrt{uv})$. The function besseli is implemented in Matlab, but this usually causes problems, because the function $I_q = (2\sqrt{uv})$ approaches infinity very quickly. Fortunately, Matlab allows us to give another scaled version, which we will call $I_q^1(2\sqrt{uv})$, which solves the divergence problem in such a way that

and therefore

Rewriting the expression (3.11), we get

3.3.2.   Numerical methods

To obtain an approximation of the exact solution of the CIR SDE we need to establish a numerical scheme. In our case, we will implement the Euler-Maruyama and Milstein numerical schemes, to see later if there are notable differences between them.

The Euler-Maruyama scheme or method is an extension of Euler's method for ordinary differential equations to SDEs.

Let $\{X_t, 0 \leq t \leq T\}$ be an Itô process solution of the following SDE:

where $W(t)$ represents the Wiener process, and suppose we want to solve this SDE in the time interval $[0, T]$.

The Euler-Maruyama approximation $Y_i$ to the true solution of $X$ is defined as follows:

● Divide the interval $[0, T]$ in $N$ subintervals of size $\Delta t = T/N$ being $0 = t_0 < t_1 < ... < t_N = T$.

● Set the initial condition $Y_0 = x_0$.

● Define recursively $Y_i$ for $1\leq i\leq N$,

where $\Delta W_i = W_{t_i+1}-W_{t_i}$.

The variables $\Delta W_i$ are independent and identically distributed normal random variables, that is, $\Delta W_i \sim \sqrt{\Delta t}Z$ with $Z \sim N(0, 1)$.

The Milstein method ^[43] is used to increase the accuracy of the Euler-Maruyama method. This is achieved by introducing a term of order 2 by using the partial derivative with respect to $x$ of $g(t, x)$.

Given an Itô process $\{X_t, 0 \leq t \leq T\}$ which is a solution of the SDE (3.15), the approximation of the Milstein method $Y_i$ to the true solution of $X$ is given by the following:

● Divide the interval $[0, T]$ into subintervals of size $\Delta t = T/N$ with $0 = t_0 < t_1 < ... < t_n = N$.

● Take as initial condition $Y_0 = x_0$.

● Recursively define $Y_i$ for $1\leq i\leq N$ by:

where $\Delta W_i = W_{t_{i+1}}-W_{t_i}$.

The variables $\Delta W_i$ are independent and identically distributed normal random variables, that is, $\Delta W_i \sim \sqrt{\Delta t}Z$ with $Z \sim N(0, 1)$.

As we can see, the expression of the scheme is the same as that of the Euler-Maruyama method, except that the following term is added:

Therefore, if $\frac{\partial g(t, x)}{\partial x}$ turns out to be $0$, this method is equivalent to the Euler-Maruyama method. When a method satisfies $E\left(|Y_i-X(t_i)|\right) \leq K (\Delta t)^\gamma$ for some $\gamma$, that method is said to be a strong approximation of order $\gamma$. Applying this, the Euler-Maruyama method is a strong approximation of order $\gamma = 1/2$, while the Milstein method is a strong approximation of order $\gamma = 1$ if $f(t, X_t)$ and $g(t, X_t)$ are $\mathcal{C}^1$ functions. The functions with which we are working in these models comply with this, so the order of convergence of the Milstein method will always be greater than that of Euler-Maruyama.

4.   Results

Since we have daily data, we set the time step $\Delta t = 1$. To compare both the Euler-Maruyama and the Milstein methods, we chose a specific time window and a number of forecasts. A window of 100 data will be taken from the first real data available, and the next 500 days will be estimated. The prediction is made for the day following the last one of the windows. Then the window is moved one unit to the right, and the process is repeated until all the forecasts are performed.

Figure 6 shows a global view of the three estimates that we want to compare. Next, Figures 7 and 8 represent, for a given time window, the real data (in grey), the ARMAX model estimates (in magenta) and the ones calculated with the CIR* model (blue) with the Milstein method. Note that since the Euler-Maruyama and Milstein methods nearly overlap, the differences between the plotted curves cannot be seen unless zoomed in further. For this reason, we refer the reader to Table 1.

It can be seen that the blue curve generally fits the real data better, so it stands to reason that the CIR* model is better than the ARMAX model. To verify this rigorously, the root mean square error (RMSE) of the CIR* model has been calculated with the Euler-Maruyama and Milstein method, and the RMSE of the ARMAX model. These errors are collected in Table 1. As shown, the CIR* model gives better results than the ARMAX model.

To verify that the ARMAX filters correctly the data, the PACF depicted in Figure 9 shows that the errors/residuals are not correlated. This is further confirmed by the Durbin-Watson (DW) test. In fact, the DW statistic is 1.9868, and the p-value of 0.8753 suggests absence of autocorrelation in the residuals.

Finally, we are going to give an interpretation of the parameters of the CIR* model at different stages of the pandemic, based on the analysis of the deterministic solution of the CIR* model given in Section 3.3.1.

4.1.   Time frame in which COVID-19 cases increase

We begin by studying a stage in which infections are increasing. For the estimated parameters to have the same trend, the window from which they are estimated must also be in the growth range, which will force us to take a small window and estimate a few values, since if we look at Figure 1, we see that the periods in which the infections grow are very short. Taking this into account, we are going to focus on the section that goes from day 210 to 229, that is, 20 forecasts, taking a window of 30 days. Based on the results obtained, we observe that both the mean of $\alpha$ and $\mu$ are negative. Since the estimated data values are greater than the mean, then $X_t$ would move away from the mean upwards.

To corroborate that this is true, another stage of increase in cases of the pandemic has been taken, specifically from day 310 to 319, that is, 10 days, and a window of 30 days. The number of forecasts had to be reduced because, as previously mentioned, the window of days must be in the growth range for good analysis. In Figure 10 one can see both sections and in Table 2, the exact values of the mean of the parameters in the two stages.

4.2.   Tme frame in which COVID-19 cases decrease

Keeping the ideas we used from the previous case, we now take a period in which the cases decrease. We estimate the interval that goes from day 265 to 284. In this case, we again obtain a negative mean of $\alpha$, which makes sense, since in periods of strong growth or decrease the values move away from the mean. On the contrary, now the mean of the parameter $\mu$ is very large and exceeds the mean of infections in that section, therefore, the values $X$ are far from the mean, but this time below.

As in the previous case, another period has been taken to verify the results. The estimates of both stages can be seen in Figure 11, and the comparison of the mean values of the parameters shown in Table 3.

4.3.   Time frame in which COVID-19 cases are stable

Finally, we are going to interpret the parameters in a section where the values are relatively constant. As can be seen in Figure 1, there are only a few sections where this occurs. We are going to take days between 110 and 139, in which it is observed that the COVID cases are close to the 0 value. This may seem surprising, but it makes sense because it precisely coincides with the end of the state of alarm. In fact, on the right side of the graph, from day 150 the cases begin to increase, coinciding with the de-escalation process and the summer of 2020 when restrictions were relaxed.

At the aforementioned stage, the mean of the parameters $\alpha$ is 0.1397, that is, positive, unlike the previous cases. The mean of the parameters $\mu$ is 0.3303, and that of $\sigma$ is 0.1448, so we have a smaller deviation than the previously analysed stages.

According to the analysis of the deterministic solution, $X_t$ will oscillate around the mean, which makes sense, since the real mean of those days is around $0.33$, and the value of the data varies between $0.25$ and $0.55$ cases.

5.   Conclusions

The main objective of this work is the study and development of some stochastic models typically used in financial markets applied to the COVID-19 pandemic in the city of Jaén.

For solving stochastic equations, both the Euler-Maruyama method and the Milstein method were used with reference to the CIR* stochastic process. Over the reported COVID-19 daily cases of the pandemic in the city of Jaén (Spain), the maximum likelihood method was used for parameter calibration. The forecasts given with this model have been validated with 95% confidence intervals and with statistical measures of goodness of fit, such as the RMSE (root mean square error). The results obtained are consistent with the observations and quite accurate. For comparison, the classical ARIMAX model has been used, resulting in more accurate predictions for the suggested CIR* model. The reason could be the relatively small size of the city of Jaén, causing large fluctuations in the number of cases that are not sufficiently softened by the moving averages, resulting in a worse behaviour of ARIMAX in comparison with CIR*. The importance of the suggested approach is twofold because it not only extends the models available to scholars to model pandemics but, also, paves the way for similar approaches in which financial models can be converted into econometric models.

Future research could be aimed at enlarging the scope to the whole province of Jaén and to other cities and/or provinces of Andalusia. In such a case, we could expect that a greater number of data could imply a longer time for the trend to change. In addition, the number of healed and deceased could also be studied, although the latter, being much smaller, will present the aforementioned problems. In terms of considered models, future research could include a comparison with the more advanced CIR# by Orlando et al. ^{[11,36,37,38]}.

In addition, although our work has only been done for one equation, it could also be generalized to systems of equations to discover the interrelation between different cities or health districts. Finally, note that stochastic differential equations are not only a very powerful tool for modelling economic-financial variables, but also in the epidemiological field, being proven from this practical point of view ^[32].

Acknowledgments

G.O. is a member of the research group of GNAMPA - INdAM (Italy) and acknowledges the project on anomalous diffusion and its applications to fractal domains. J.G. acknowledges Junta de Andalucía through the project FEDER-UJA-1381026. The Consejería de Salud y Familia of Junta de Andalucía is acknowledged for providing the COVID-19 data.

Conflict of interest

The authors declare there is no conflict of interest.