Communications in Nonlinear Science and Numerical Simulation The stochastic θ -SEIHRD model: Adding randomness to the COVID-19 spread

In this article we mainly extend a newly introduced deterministic model for the COVID- 19 disease to a stochastic setting. More precisely, we incorporated randomness in some coefficients by assuming that they follow a prescribed stochastic dynamics. In this way, the model variables are now represented by stochastic process, that can be simulated by appropriately solving the system of stochastic differential equations. Thus, the model becomes more complete and flexible than the deterministic analogous, as it incorporates additional uncertainties which are present in more realistic situations. In particular, confidence intervals for the main variables and worst case scenarios can be computed.


Introduction
The coronavirus disease 2019, renamed as COVID-19, is an infectious disease produced by the virus SARS-CoV-2, which was declared pandemic by the World Health Organization (WHO) in March 11, 2020. This disease has posed many novel scientific challenges, ranging from contagious patterns, medical treatments or vaccine developments, to data analytics, spread modelling or evolution forecast. The research on all of these topics has been intensive in the last few months, as so will be in the near future. Thus, many disciplines from interconnected fields, like bio-sciences, mathematical and statistical modelling or artificial intelligence, will be involved and they will play an important role to rapidly overcome this exceptional emergency situation.
From the mathematical modelling point of view, epidemiological compartmental models have been often used to understand and analyse the behaviour of the contagious diseases, with the article [1] being the pioneering work. These models provide useful tools to make predictions about the future evolution of the epidemic and to control its propagation. In the literature of epidemic modelling, many examples of general-purpose models have been proposed (see [2], for a review), each of them accounting for some specific characteristics of the diseases, like the well-known examples SIR (Susceptible-Infected-Recovered), SEIR (Susceptible-Exposed-Infected-Recovered) or SIS (Susceptible-Infected-Susceptible) models.
In addition, compartment models are building blocks of more involved approaches like the Agent-Based (AB) models, see [3] for example. This type of models are potentially interesting, as they are designed to include heterogeneity in the modelling, in contrast to the homogeneity that characterizes the compartmental counterparts. Such models are then able to incorporate individual attributes and the interaction among the agents themselves and with the environment, at cost of increasing the model complexity. However, AB models exhibit certain relevant disadvantages. Firstly, they are not suitable candidates at the initial stages of diseases when the lack of high quality data, specially when considering lower (micro) levels, can cause a problematic calibration of model parameters and erroneous predictions. As stated in [4], this fact can produce catastrophic consequences due to the butterfly effect. In addition, overparameterized AB models can result in a poor balance between complexity and interpretability, a key aspect for policy makers to decide the control measures and deliver clear messages. Furthermore, AB models are computationally intensive, while the analogous compartmental ones often turn out to be much faster alternatives. Adding extra compartments to the classical SIR-like models to accounting for different groups of interest can be understood as an intermediate strategy, aiming to preserve the advantages of both kind of models, while mitigating its own disadvantages. All in all, here we focus on an advanced compartmental model that includes enough information keeping a sufficient degree of flexibility.
Typically, the compartmental models are originally formulated following a deterministic approach, i.e., in terms of a system of Ordinary Differential Equations (ODEs), although, usually, they are readily extended to a stochastic version. There are two common approaches to include stochasticity into a deterministic model, relying on either the Continuous Time Markov Chain (CTMC) or Stochastic Differential Equations (SDEs), see [5] and the references therein. The stochastic models allow to capture many kinds of circumstances including all types of uncertainty that may influence the compartments dynamics: behavioural effects, public interventions, seasonal patterns, environmental factors, etc. For example, in [6,7], the authors consider a statistical inference-based approach to include random noise arising from different sources/components. Furthermore, while the solution of a deterministic model is given by a set of functions of time uniquely dependent on the initial data, the solution of the stochastic model is a set of stochastic processes, containing much more information than the deterministic analogous. In fact, at each time instant we can exploit the information provided by the probability distribution associated to the underlying random variable. A statistical analysis can be therefore performed, producing useful quantities like expected outcomes, quantiles or worst case scenarios.
In this work we present a stochastic extension of the deterministic compartmental model in [8], that has been proposed as an ad hoc model for the COVID-19 disease. However, note that the stochastic extension proposed here is modelindependent and it can be exploited in a similar fashion under any compartmental-based approach. Unlike to the classical models (SIR, SEIR, SIS), the selected model is adapted to the specific characteristics of the COVID-19, taking into account, besides the usual factors, the undetected infectious cases, the hospitalized population or the deaths, for example. Thus, the so-called θ-SEIHRD model proposed in [8] was developed as a very general and sophisticated model (based on a previously introduced model [9]) in order to be able to study the spread of the disease worldwide. In the same work, the authors proposed a simplified version which reduces the mathematical complexity towards a more tractable model, but yet preserving the ability to capture the most relevant features. Other recent studies addressing the mathematical or statistical modelling in the context of the COVID-19 include [10][11][12][13], for example. As the literature focusing on COVID-19 modelling is increasing a lot during last months and weeks, it results quite difficult for the authors to select a list of the more relevant papers more recently arising in the topic, surely some of them being developed in parallel to this contribution.
Here, we will follow the SDE-based approach where the randomness is typically achieved by incorporating a Brownian motion, also known as Wiener process, to the ODEs of the deterministic system. There are two common ways of addressing this stochastic extension. On the one hand, an arbitrary random noise can be added to some of the equations in the ODE system, thus transforming them into SDEs. On the other hand, one (or more) of the existing model parameters can be perturbed, meaning that it (they) becomes a random variable or a stochastic process. Although both approaches have been well investigated, recent examples are [5,14,15] for the random noise approach and [16][17][18] for the parameter perturbation, here we prefer the second alternative. When adding random noise to an equation, it often comes with either an extra parameter to control the level of volatility or a covariance term (usually relying on the Cholesky decomposition as in [5]). Both methodologies might lack biological interpretation, specially when employing a mid-high number of compartments (SDEs), as it is the case of the model considered in this work. Furthermore, in practice, the uncertainty will have impact on a particular model component, typically represented by a model parameter. Therefore, a randomly perturbed parameter can be reasonably explained in terms of the variability produced by the source of the considered uncertainty.
In this work, we will employ a randomly perturbed disease contact rates. The approach proposed here is however conceptually different than the ones typically found in the literature, aiming to deal with some of the observed inconsistencies in the COVID-19 context. As a natural choice, many stochastic SIR-like models rely on a normally distributed perturbation for the disease contact rate parameters (see [16,18], for example). However, this approach allows negative values for the rates, which is biologically non-sensical, potentially appearing when the rate is close to zero (a pattern observed in practice for the COVID-19, see [19]). This issue is overcome in [17], by employing an exponential Ornstein-Uhlenbeck (OU) process to model the contact rate. Although this process indeed ensures positiveness, it presents another undesirable feature: an increasing variance in time. When analysing the disease transmission dynamics, there is not an objective evidence of more variability in the disease contact rates in the long-term [20]. Actually, one can expect more control of the disease spread patterns over the course of the epidemic. For all the reasons mentioned above, we propose to model the contact disease rates by means of the so-called CIR process, named after its authors Cox, Ingersoll and Ross [21]. The CIR process is widely employed to simulate the evolution of interest rates in the quantitative finance framework. In some sense, the interest rates in computational finance and the disease contact rates in epidemiology present a rather similar behaviour: positiveness, controlled variability and long-term stability.
The article is organized as follows. In Section 2 we introduce the proposed mathematical model. Section 3 describes the numerical methods for solving the model and compares it with a simpler alternative. Section 4 contains the obtained results and their numerical and statistical analysis. Finally, Section 5 points out some conclusions and possible future research lines.

A stochastic compartmental model for the COVID-19
As mentioned in the previous introduction, we base our approach on the the so called θ-SEIHRD model first proposed in [8], which provides an advanced extension of the classical SIR-like models and adds to the usual compartments in the SEIR model the very specific ones for COVID-19: undetected infectious, hospitalized and deaths compartments. In fact, it results in a very sophisticated, well motivated and general model, with a significant number of free parameters, allowing an enormous flexibility (including several territories). In order to make this model more practically tractable, the authors also proposed a simplified version by imposing several assumptions: single territory, regular natality/mortality is neglected, and no imported/exported cases. Furthermore, some predefined forms for the open parameters of the model are adopted.
For the sake of completeness, we devote the following section to briefly describe the most relevant components of the simplified θ-SEIHRD model, before presenting its extension with randomly perturbed parameters towards the stochastic θ-SEIHRD model here proposed.

Original model
The original model consists of nine compartments, including classical ones like susceptible, exposed or infected (although some have a slightly different interpretation in comparison with the usual meaning) and other ones that are COVID-related, i.e., they account for the particular features of the disease. As it is common in this kind of models, an individual stays in a compartment a period of time and then moves to another compartment according to some transition rates. More precisely, the simplified θ-SEIHRD model in [8] is given by the following system of differential equations: where the first two compartments, S(t) and E(t), denote the persons not affected by the disease pathogen (Susceptible) and the ones in incubation after being infected without clinical signs (Exposed), respectively. The compartment I(t) (Infectious) includes the persons in the very preliminary stage of the infection, where nobody is expected to be detected yet, although they may infect other people after finishing the incubation period. After this period, they can either remain undetected and enter in the compartment I u (t) (Undetected Infectious), or be taken in charge by sanitary authorities and enter in the hospitalized/quarantined group. People in compartment I u can infect others and develop the disease, although they are not reported by sanitary authorities. Moreover, it is assumed that people in compartment I u present no (or very weak) symptoms, so that they directly pass to recovered state after an infectious period of time. Among hospitalized persons (which include those ones in quarantine at home), we distinguish between those ones who will recover, placed in compartment H R (t), and others who will die, placed in compartment H D (t). Note that people in hospital can still infect others. Among recovered people, we distinguish between persons who were previously detected, R d (t), from those ones who were undetected, R u (t). Individuals in recovered compartments are not contagious anymore and have developed immunity. Finally, D(t) denotes the persons who did not overcome the disease and died from COVID-19. As every compartment presents a dependence in time, we incorporate this time dependence by referring to each compartment's situation at instant t. More detailed explanations about this deterministic model can be found in [8].
The right hand sides of equations in system (1) contain the balances between the entries and exits of persons in each compartment. In order to properly pose the initial value problem, the system of equations is completed with the initial data S(t 0 ), and D(t 0 ), representing the compartment's composition at inception time, t 0 . Note that the last three equations of the system (1) are uncoupled with the six previous ones, so that the expression of their solution can be obtained and is given by once the first six ones have been solved in a first stage. Thus, only the first six coupled equations need to be solved.
The model includes a number of open parameters, which are determined either from the available data/literature or by calibration. In the following, a detailed explanation of each parameter is provided.
Here, only one control measure is assumed (mobility restrictions, for example), implemented at date λ 1 and represented by the following decaying function: with the parameter κ 1 ∈ [0, 0.2] entering in the calibration procedure. The generalization to more/individualized control measures is straightforward.
• The fraction of infected individuals which are detected and documented by the authorities: θ. We assume that all deaths are detected, so that θ ∈ [ω, 1], and that it changes in time. More precisely, we consider the expression with θ, θ, λ 1 , and λ 2 inferred from the data.
• The compartment transition rates: Thus, represents the decrease of the duration of d I due to the application of control measures at time t, with d g being the maximum number of days that d I can be decreased due to the control measures. The parameter δ R is included in the calibration. • The disease contact rates: β E , β I , β Iu , β H R , β H D ∈ R + . The parameter β I is assumed to be given, after calibration to the available data. Furthermore, the following relation between β I and the rest of the contact rates is considered: where β I = C u β I , with C E , C H (t) and C u ∈ [0, 1]. Parameters C E and C u are also obtained by calibration, while C H (t) is determined by the following expression: with α H being the percentage of cases (healthcare workers) infected by individuals in compartments H R or H D .
Expression (3) results from the use of the available data for infection transmissions within hospitals, which is a trustworthy and more easily gathered information, specially at the beginning of the pandemic. The term incorporates a relation with the other disease contact and transition rates (including a temporal dependency).

Stochastic extension
Our aim is to introduce stochasticity in the simplified θ-SEIHRD model previously described. Moreover, in order to keep the interpretability of the proposed stochastic model, we will add some randomness to a set of parameters. More precisely, we will add randomness on the disease contact rates, β's. Note that in the simplified deterministic version proposed in [8] all the β's depend on β I , as indicated in Eq. (2). Therefore, we propose to add uncertainty in β I by turning it into a stochastic process, i.e., a random variable at each time instant.
Then, we start by writing the model in terms of β I . From Eq.
(2), we have ) . By using the previous notation, the simplified θ-SEIHRD model (1) where As previously indicated, we incorporate a stochastic component into the model by replacing the parameter β I of the deterministic model with a random variable. Instead of considering a constant value for β I , the disease contact rate in compartment I follows a newly introduced stochastic processβ I (t). In order to preserve the positiveness in the parameter definition (imposed by its biological interpretation), we choose the well-known CIR process [21]. The main advantage of the CIR process is that it theoretically ensures the spacial states to be non-negative. Moreover, as it is a mean-reverting process, the dynamics of the CIR process tends to a prescribed value in the long term. This property can have the following biological interpretation: at the early stages of the disease, the contacts between people are less controlled (so more volatile), while in the mid-long term the individual contacts and the disease spread patterns are more studied and the ranges of variability are more reduced. Accordingly, the dynamics ofβ I satisfies the following stochastic differential equation: where ν β I is the mean reverting speed, µ β I is the long-term average, σ β I is the volatility, and dW (t) is a Brownian motion increment. Therefore, the stochastic θ-SEIHRD model is governed by the following system of Random Ordinary Differential Equations (RODEs): with M(t) as defined in Eq. (5). Although only one source of randomness is introduced, the solution of the whole system becomes a set of stochastic processes, due to the dependence of the remaining equations on the first equations for S and E, and their own dependence onβ I . Note also that we could add time dependency toβ I both in the deterministic model (which is not the case in [8]) as well as in the stochastic version (by considering either ν β I , µ β I or σ β I time dependent).
Remark. Note that our approach could be generalized to a setting where some of the parameters were independent, in this case we would consider each one as a Gaussian random variable with possible correlations between (some or) all of them. In this more general setting, a certain number of different (possibly correlated) Brownian motion processes would come into place.

On the stochastic θ-SEIHRD model
In this section, we provide a detailed analysis of the stochastic θ-SEIHRD model defined in (7).

Existence and uniqueness of solution
First, note that (7) is a system of RODEs driven by an Itô process, which defines the CIR model. Since the paths of the Itô process are at most Hölder continuous, the solution of the system of RODEs is at most continuously differentiable independently of the smoothness of the functions that define the right hand side of the system of ODEs [14]. Therefore, the usual arguments based on Taylor expansions to obtain the order of convergence for the classical numerical methods does not apply. For example, classical Runge-Kutta methods do not achieve their order of convergence in ODEs when applied to RODEs (see [22] and the references therein).
Secondly, in order to establish the existence of solution for the system of RODEs (7) driven by the CIR process, we introduce the notation where the super-index (·) T denotes the transpose operator. Moreover, we introduce the following notation for the coefficients of the CIR process, a(β I ) = ν β I (µ β I −β I ), b(β I ) = σ β I √β I . By using the previous notation, the system (7) can be equivalently written as where F is the function that defines the right hand side of the system of RODEs. Note that (8) can be understood as a system of SDEs where the diffusion coefficients for all equations are equal to zero except for the last equation which is equal to b(β I (t)). For a set of constant initial data S(0), E(0), I(0), I u (0), H R (0), H D (0), R d (0), R u (0), D(0), andβ I (0), the system (8) has a unique strong solution. This follows from the fact that the coefficients of the system (8) of SDEs are locally Lipschitz continuous and the initial condition is a constant value (see [23], for example). Therefore, the existence and uniqueness of solution for the θ-SEIHRD model defined by (7) follows.

Numerical solution
As the system (7) is nonlinear, it is not possible to obtain a closed-form expression for the solution, as it also happens with the corresponding deterministic version. Therefore, the use of numerical methods for solving (7) becomes mandatory. Here, we adopt the following strategy. First, we perform a simulation of the dynamics ofβ I (t), in accordance with the CIR process. After, we solve the resulting ODE system for each simulated path ofβ I (t). In this way, we obtain a set of random walks for each stochastic process representing a model variable.
Remark. Although in the stochastic θ-SEIHRD model (7) the solution is a finite set of stochastic processes, we keep the same notation we used for the finite set of real valued functions representing the solution of the deterministic θ-SEIHRD model.
The CIR process is a well-studied dynamics often employed in computational finance (see [24] for example), which satisfies the SDE in Eq. (6). From the mathematical point of view, one of its relevant features is that the underlying distribution is known analytically, relying on the non-central chi-squared distribution. Thus, for s < t, the conditional distribution ofβ I defined in Eq. (6) reads with the constant parameter d = 4ν β I µ β I σ 2 β I and some initial valueβ I (t 0 ) =β I (0). By employing this scheme, we can simulate n discrete sample paths ofβ I . Once the sample paths ofβ I have been obtained, we numerically solve n ODE systems, one for each of these paths. For this purpose, we choose the explicit Runge-Kutta method of order 5(4) when applied to deterministic ODE systems, known as RK45, RKDP or Dormand-Prince method, see [25]. Its practical implementation requires that the time discretization mesh includes the time points used in the simulation ofβ I . As previously indicated, the lack of enough smoothness of the sample paths of the CIR process reduce the regularity of the solution and also the order of convergence of Runge-Kutta method obtained when solving deterministic ODE systems.
Remark. Alternative numerical methods for solving (7) could be based on the use of classical numerical stochastic methods for solving the equivalent system of SDEs (8). For example, applying the Euler-Maruyama method to (8) would be equivalent to use the explicit Euler scheme to solve the eight ODEs with random coefficients combined with Euler-Maruyama scheme to approximate the paths of the CIR process. Also note that using a stochastic Runge-Kutta scheme for the system of SDEs (8) would be equivalent to use a Runge-Kutta scheme for the ODEs in (7) combined with a stochastic Runge-Kutta method to approximate the paths of the CIR process. Both previous alternatives, by Euler-Maruyama and Runge-Kutta stochastic numerical methods for SODEs, have been numerical analysed in the classical book [23]. By using the equivalence between a RODE driven by an Itô process and the corresponding Stochastic ODE, numerical methods for this kind of RODEs are analysed in [22].
In our approach, instead of using numerical methods to approximate the sample paths of the CIR process, we employ what is known as an exact simulation scheme since we exactly sample the distribution to generate each sample path. This exact simulation is combined with a Runge-Kutta method to integrate the resulting ODE for each path. In this work we do not address the numerical analysis of the proposed numerical strategy, but only use it to simulate the proposed stochastic model.
In order to illustrate the potential of the stochastic version of θ -SEIHRD model, in Figs. 1, 2 and 3 we present the solution of the deterministic model versus a number of possible scenarios (simulations) obtained by the stochastic model. Further, we include the so-called ''Average scenario'', which is nothing else than an estimation of the expected value of the model variables in time. We can clearly observe that the outcomes produced by the stochastic model provide much more information about the future evolution of each model compartment, while the deterministic one seems to be somehow an averaged version of the stochastic formulation. This gives just an initial insight of the potential of the newly introduced model, which will be completed in Section 4 with a more detailed analysis.

Comparison with a stochastic SEIR model
In this section, we argue the reasons why we choose to extend the θ-SEIHRD model to a stochastic setting, instead of any other simpler alternative. For this purpose, let us consider the well known SEIR model. The θ-SEIHRD can be easily where all the model parameters have the same interpretation as in Section 2, while we incorporate a specific disease contact rate,β I , for compartment I. Since this newly introduced SEIR-like model does not distinguish between detected and undetected (actually it considers all the cases detected, i.e., θ = 1), the disease contact rate of compartment I needs to be compensated in order to take into account this fact and allow to compare the results provided by both models. Therefore, we adopt the following natural formulation, in which we assume that the new disease contact rate of compartment I is the addition of the two rates of the original compartments for detected and undetected infections, i.e., β * I = β I + β Iu = β I + β I , where the second equality comes from Eq. (2) and θ = 1. The advantage of this approach is that it enables the use of the previously defined parameters, which are already calibrated to the data.
Next, we compare the outcomes produced by the models defined by Eqs. (1) and (9). In Fig. 4a, we observe that the curve of infected cases resulting from the SEIR model is very similar to the one obtained with the θ-SEIHRD model (see Fig. 1e). Therefore, we can conclude that both models provide equivalent outcomes in terms of the infected cases.
Let us now consider a stochastic SEIR-like model that results from simplifying the stochastic θ-SEIHRD of (7) in a similar way as for the deterministic case, i.e., taking m Iu (t) = m H R (t) = m H D (t) = 0 and θ = 1, incorporating the ''removed cases'' compartment, R. Then, we have the dynamics , where a new stochastic process,β * I (t), following a CIR dynamics is considered. By using the stochastic SEIR model defined by Eq. (10), we can now generate a number of scenarios and compute an estimation of the expected infected cases, E[I(t)], as we have done in the previous subsection. In Fig. 4b, several of those Monte Carlo scenarios and the expected value of I(t) (labelled as ''Average scenario'') are depicted. Again, we can easily conclude that this stochastic extension of the SEIR model reproduces the observed behaviour for the θ-SEIHRD model (see Figs. 1, 2 and 3), since the deterministic version can be seen as an averaging scenario of the stochastic generalization. Next, the total number of active cases predicted by both models, θ-SEIHRD and SEIR, and their corresponding stochastic versions are considered. The obtained results are presented in Fig. 5. We can readily observe that the SEIR model tends to significantly underestimate the number of active cases, an issue that is highly undesirable. This effect is even more pronounced when the stochastic models are employed, thus obtaining a gap of around 3000 cases in average (around 8000 cases in the worst case scenario). The undesired behaviour produced by the simplified SEIR alternative can be very problematic. Although this simpler model is able to capture the infected cases (see Fig. 4), an erratic estimation of the total active cases is observed. These poor predictions of the disease evolution generated with the SEIR model are caused by the concentration of the cases in a single compartment and the absence of intermediate ones (hospitalized). In addition, the impact of this effect on the stochastic model outcomes is significantly more important, which makes the introduced SEIR model an inappropriate candidate to be stochastically extended in the context of the COVID-19 disease. Remark. In order to employ the stochastic SEIR alternative presented above to model the COVID-19 disease evolution, it would require ad hoc calibration for its free parameters, both the model related and the CIR process related ones. However, our goal here is to compare the extension of an existing and already calibrated deterministic model to its stochastic counterpart.

Numerical and statistical analysis
The experiments have been conducted in a computer system with the following characteristics: CPU Intel Core i7-4720HQ 2.6 GHz, 16 GB RAM memory and GPU GeForce GTX 970M. The numerical codes have been implemented in Python programming language. We consider a uniform time grid, i.e., ∆t := t i+1 − t i , ∀i, with time step ∆t = 1 6 (around 4 h).
In this section, we perform a numerical and statistical study of the proposed stochastic model. As mentioned, the solution of the system of SDEs in (7) is a set of stochastic processes, meaning that we can ''extract'' a random variable at each time point. In this way, not only a single value (like for the deterministic case), but also some statistics can be provided for a prescribed time. In this work, we consider the mean, the interquantile interval, [Q 1 , Q 3 ] (with Q 1 and Q 3 being the first and the third quantiles, respectively), and a worst case scenario (WS), applied to both the evolution of the model variables and the possible model outputs.
Here we take advantage of the experiments conducted in [8], referred to the case of China. In Tables 1 and 2

Model variables
We firstly test the evolution of the model variables. Thus, we extract some simulation-based statistics at a couple of time instants, the 8th February (inflection point) and the 29th March (final point). Furthermore, the impact of different levels of uncertainty is also reported by considering several representative values for σ β I , thus reflecting situations of no    uncertainty (σ β I = 0), low uncertainty (σ β I = 0.1) and high uncertainty (σ β I = 0.5). In all cases, the long-term mean of the perturbation is set to the calibrated value of β I for the deterministic model, i.e., µ β I = β I . The mean reverting speed parameter is chosen as ν β I = 1, thus representing a regular (not too high, not too low) reversion speed. In Table 3, the obtained results are presented. We can clearly observe the significant impact of the uncertainty in the disease evolution. In the case of higher uncertainty, the number of infections, in average, is almost doubled and the worst case scenario multiplies this value by six. Even when a lower volatility is considered, the increment of cases and deaths in the worst scenario becomes important, up to 50%. Next to the previous experiment, in Fig. 6 we present the histograms for the model variable I(t) to give an insight of the impact of the uncertainty in the infection evolution. First, we clearly observe that the produced distribution is skewed with a fatter right tail. Secondly, the bigger the volatility of processβ I (t) (denoted by σ β I ), the fatter the right tail becomes, thus indicating more probability of extreme events.

Outputs
In [8] the authors proposed a set of possible outputs that can be useful for the authorities to plan the resources allocation (like the number of clinical beds, among other indicators). These outputs are: • Hospitalized people, Hos(t), at time t: where p(t) is the fraction of people in compartment H R that are hospitalized and C o is the period of convalescence.
• Maximum number of hospitalized people in the interval [t 0 , t]: Hos(τ ).  We therefore perform a similar statistical analysis as before, although now reporting the model outputs. The results are shown in Table 4. Again, it is clear that the uncertainty can significantly affect the disease evolution, justifying why it should be taken into consideration. For example, the people requiring hospitalization may vary from around 4500 with no stochasticity included, to around 30000 including randomness. As previously pointed out, one of the major advantages of the θ-SEIHRD model is that it accounts for important aspects of the COVID-19 pandemic, directly affecting the population or the healthcare systems. Particularly, the evolution of some curves like infected, hospitalized, and deaths is typically reported by the authorities in both cumulative and daily fashion. In Figs. 7 and 8, we show the stochastic model outcomes for these specific curves, considering two uncertainty levels,