key: cord-0542030-emdb1vvp authors: Baltazar-Larios, Fernando; Delgado-Vences, Francisco; Diaz-Infante, Saul title: Maximum likelihood estimation for a stochastic SEIR system for COVID-19 date: 2021-11-24 journal: nan DOI: nan sha: 902de02b865b340d6bc11170c5af520df24b613f doc_id: 542030 cord_uid: emdb1vvp The parameter estimation of epidemic data-driven models is a crucial task. In some cases, we can formulate a better model by describing uncertainty with appropriate noise terms. However, because of the limited extent and partial information, (in general) this kind of model leads to intractable likelihoods. Here, we illustrate how a stochastic extension of the SEIR model improves the uncertainty quantification of an overestimated MCMC scheme based on its deterministic model to count reported-confirmed COVID-19 cases in Mexico City. Using a particular mechanism to manage missing data, we developed MLE for some parameters of the stochastic model, which improves the description of variance of the actual data. Parameter calibration of epidemiological models with data is a demanding but crucial task. The current COVID-19 pandemic has dramatically revealed the complexities and subtleties of developing epidemic data-driven models. Although computational cost would grow when we develop a model with stochastic differential equations, in some particular cases this would lead to tractable estimators. By quantifying uncertainty by Brownian motion and applying methods based on maximum likelihood, we can improve the calibration of a deterministic base with conventional methods. When we calibrate parameters of an epidemic model, we argue that we can sometimes improve the quality and computing effort by adding noise and proper information management. Models based on stochastic differential equations (SDEs) have the advantage of characterising the central tendencies-usually a regarding deterministic model to describe the mean. This stochastic structure could also capture sources of variability (see e.g. Allen, 2007; Han and Kloeden, 2017; Gardiner, 2009 ). These models represent an alternative to analyse the skewed data that is commonly sampled in the observations. Stochastic models, and particularly SDEs, attempt to capture in their formulation the random variability of the data, and therefore quantify the uncertainty that is naturally related to parameters or phenomena Allen (2007) . One of the crucial challenges in likelihood-based inference for scalar SDEs is to compute of the underlying transition density, and hence the likelihood function. Thus, maximum likelihood estimation on the true likelihood of the solution of SDEs is a rare case. Moreover, because the available sampled solution is a sequence of discrete time observations, it is usually incomplete and limited, and we must instead appeal to approximations Iacus (2008) . Furthermore, this situation worsens for a coupled system of SDEs. Regarding stochastic models for epidemiology, and in particular for COVID-19, in previous studies there are two types of works related to the one presented here. On one hand, some manuscripts report similar ideas but for simpler structures and with numeric experiments of synthetic nature. For example, Pan et al. reports the estimation of a SIS structure that only needs a scalar version of the model to approximate the likelihood Pan et al. (2014) . The authors use a pseudo-likelihood and least square method to estimate the parameters in the model. They also calculate confidence intervals by applying least squares. In contrast, the results presented here deal with an epidemic model of five nonlinear coupled SDEs and with real data for COVID-19 from Mexico City, we also have a natural discretised version of the MLE to estimate the parameters in the model. On the other hand, there are complex stochastic models similar to the one considered in this manuscript that are only studied on a theoretical basis and/or with numerical simulations. However, this does not address estimations methods. For example, in Faranda and Alberti (2020) , the authors propose a SEIR model, with four classes of individuals, given by a random differential equation. Indeed, they consider a model where two parameters are Brownian motion with drift and one parameter is the solution of a SDE with a log-normal perturbation. Nevertheless, the authors do not estimate parameters from actual data to fit the model, the parameters are instead taken from previous works. Additionally, they used date information about confinement in Italy and France to validate the model. In Djević et al. (2021) , the authors study a version of the stochastic SEIR model for the . Their model is a system of stochastic differential equations driven by two independent standard Brownian motions. More precisely, their model depends on several parameters, including the constant transmission coefficients related to symptomatic infected individuals (β), hospitalised individuals and to superspreaders (β ). In their study, the authors formulate a stochastic model via a stochastic perturbation of the constant transmission rates β, β . They show the existence and uniqueness of the global positive solution. They also study conditions under which extinction and persistence in mean hold. Finally, they provide numerical simulations of the stochastic model. Nevertheless, estimation methods to fit the model is an open problem for their model. This manuscript focus on the calibration of a stochastic SEIR (Susceptible-Exposed-Infective-Recovered) model with data of symptomatic reported and confirmed cases in Mexico City. According to a specific time window -where the outbreak follows an increasing profile-, our approach relies on the maximum likelihood estimation method. Here, we formulate our stochastic SEIR model through a stochastic perturbation of the natural death rate of the population by a Brownian motion. We then deduce a nonlinear and coupled Itô stochastic differential equation (SDE), driven by an affine Gaussian multiplicative noise. Our main objective pursuit the development of consistent Maximum Likelihood Estimator (MLE) for a stochastic SEIR structure. Parameter estimation for SDEs is an active area in the statistical literature. For example, in Young (1981) Young gives an overview of parameter estimation methods with continuous-time observation. Sørensen treat the case where the data correspond to the solution of SDEs observed in a sequence of discrete time points. We refer the reader to Sørensen (2004) for an introduction of several inference methods as estimating functions, analytical and numerical approximations of the likelihood function, Bayesian analysis Multi Chain Monte Carlo (MCMC) methods, indirect inference, and Expectation-Maximisation (EM) algorithm. Due to the COVID-19 pandemic, the estimation of parameters and development of epidemic models are under special attention. The most recent advances combine several models and data estimation. Among others, we see relevant methods based on Kalman filters, data aggregation Fintzi et al. (2017) , and techniques from data science such as Artificial Intelligence and machine learning Bragazzi et al. (2018) . Our work is related to Liu ( Our contribution consists of the development estimators for some crucial parameters of the stochastic SEIR model. Here, we apply the quadratic variation of all of thee processes to estimate the volatility parameter and calculate the MLE for parameters representing the symptomatic infection rate, asymptomatic infection rate, and proportion of symptomatic individuals. Then, from a theoretical result, we deduce a target parameters likelihood without assuming a prior data distribution. In short, this likelihood and a mechanism to manage incomplete data form the core mechanism that allows us to estimate the parameters. Data-driven epidemic models usually treat observations as counters. Assuming that data follows a given distribution, this could be a Poisson, negative binomial, and so on. O'Neill (2010) . Then, by applying some statistical inference tools to estimate some parameters, the model is fitted to the incidence data. Here, we avoid this assumption. A description of the procedure follows. We use the Radon-Nikodym derivative of two equivalent measures-which, in our case, are the respective solutions of the stochastic SEIR system with two different sets of values of the parameters-to obtain the likelihood of the parameters of interest. This means that we use the Cameron-Martin-Girsanov Theorem to set up the likelihood. Note that the method to establish the likelihood does not assume a particular distribution of the solution, therefore,-we assume that the data have an abstract distribution. We pursue estimators that satisfy strong consistency and numerical schemes of high order to its computing. With these ideas, we also illustrate their efficacy in the time of computing and performance with real-data by simulation. Our simulations suggest that the methodology developed here improves a particular outcome of a fitting with MCMC. The scheme of prior distributions returns a bias fitting with this method. In contrast, the fitting of the SDE model with MLE improves the capture of the data variance and the time of computing. According to our Theorem about consistency, by computing residuals we confirm that the stochastic model follows the profile of the regarding data and is coherent with our modelling assumptions. The results presented here point to a methodology that would boost the uncertainty quantification by SDE epidemic models with incomplete information and with dramatically less computing time. Here, we report the first pieces of evidence. After this brief introduction, in Section 2 we treat a deterministic SEIR model for COVID-19 and its calibration with MCMC. Section 4 forms the main content of this article: the construction of estimators for important parameters and the analysis of its consistency. We illustrate our methodology via simulation in Sections 5 and 6 by synthetic generated data and data of confirmed symptomatic cases from Mexico City. We consider that susceptible individuals become infected when they are in contact with asymptomatic individuals or individuals with symptoms. We propose that a proportion of asymptomatic individuals have a way to get relief and not die. A proportion of individuals infected with symptoms may die of the disease or may be relieved. We will introduce two types of control, vaccination of asymptomatic infected and susceptible individuals and treatment of individuals infected with symptoms. We use a common SEIR deterministic structure that has been applied for COVID-19. Our formulation splits the entire population N according to the following compartments: Suceptible (S) This compartment denotes the populations' member under risk of contagious. Exposed (E) Members of this class come from the susceptible class and represent the population fraction that hosts but can not transmit the SARS-CoV-2. Then, after a period of κ −1 days, they develop the faculty to spread this virus and became infected. Asymptomatic Infected I a This class denotes the population fraction that spreads the SARS-CoV-2 but does not develop symptoms or just develop mild symptoms. This class also enclose all non-reported COVID-19 cases. Symptomatic Infected I s Represents the COVID-19 reported and confirmed cases that develop symptoms. Recovered After a time period of length α −1 s , a member of the infected class enters this compartment to become recovered and immune. Death Here we count the cases that lamentably dies due to COVID-19. In symbols, our deterministic base reads (1) Table 1 enclose a description of the underlying parameters. To calibrate parameters of this deterministic base (1), we deploy a MCMC by counting the cumulative incidence of new infected and reported symptomatic cases. We denote by Y t this cumulative incidence at time t and suppose that its profile follows a Poisson distribution with data mean λ t = EY t . Thus, following ideas from Acuña-Zegarra et al. (2020) we postulate priors for p, β s , β a and count the cumulative reported-confirmed cases in the CDMX-Valle de Mexico database Gobierno de México (2021 (Accessed January 04, 2021). µ Natural death rate β s Symptomatic infection rate β a Asymptomatic infection rate κ Transfer rate from the exposed class to the infected stage, thus κ −1 is the average incubation time. p Proportion of asymptomatic individuals θ Proportion of symptomatic individuals who die due to the disease α −1 s Symptomatic recovery period α −1 a Asymptomatic recovery average time γ −1 Disease-immunity period Table 1 : Parameters definition of the model in Equation (1). State Prior 26 446 062.901 024 N Table 3 : Initial conditions of the system of equations (1). In this section, we consider a stochastic version of the model given by (1). We aim to quantify uncertainty. To this end, we conveniently perturb the natural mortality rate µ by a standard Wiener process W (t). That is, we add a stochastic infinitesimal fluctuation of the parameter µ in the interval (t, t + dt) by adding the corresponding Wiener increment dW (t). In symbols, we have This perturbation infinitesimally quantifies the environmental fluctuations of the population mortality, due, for example, to differences in the health of individual conditions and other-regarding environmental factors. Furthermore, this particular form allows us to deduce a closed form for the underlying likelihood. Because our model couples the spread and other processes, this perturbation also implicitly captures random fluctuations from all involved compartments. In other words, we describe this effect by a stochastic process with mean zero and dispersion proportional to σ √ t. According to Wiener statistical properties, we have E [µdt + σdW (t)] = µdt and Var [µdt + σdW (t)] = σ 2 dt. This approach has been employed in the emergent research of mathematical epidemiology, see for example Chang et al. (2019) ; Acuña Zegarra and Díaz-Infante (2018); Liu et al. (2018) . Substituting perturbation (2) in equation (1) results, Existence and uniqueness of the solution in a compact time interval [0, T ] follows from the boundedness, local Lipschitz and linear growth conditions on the coefficients of each SDE component (seeØksendal (2003) Theorem 5.2.1). Furthermore, Burkholder-Davis-Gundy and Holder inequalities imply that, for with constant C > 0 depending on the model parameters. We rewrite the SDEs in Equation (3) by applying Lamperti's transformation to each component (see Iacus (2008) ). Thus, the susceptible component S becomes (4) Again, by Lamperti's transformation for the equation of process E, we obtain Analogously, we obtain the Lamperti transformation, for the components processes I a , I s , R: Then, we can write these SDEs as a five dimensions system of SDEs: where, We apply the version of the Girsanov formula presented in (Särkkä and Solin, 2019, see Theorem 7.4) for the system in Equation (9). We assume that the true value of the parameters β s,0 , β a,0 , p 0 are unknown. We denote by P β,p the law for the solution of (9). It is known that the measures P β,p are equivalent for different values of β, p (see for instance Iacus (2008) , Bishwal (2008) , Särkkä and Solin (2019) . Then, by Theorem 7.4 in Särkkä and Solin (2019), we have that the likelihood ratio, or Radon-Nikodyn derivative has the form where Q, in our case, is the five-identity matrix and Note that the MLE for β a and β s is coupled, but independent of p. Therefore, we calculate the MLE for β a and β s together and for p separately. Thus, We denote by = a or s, thus β will mean the variable β s or β a and similarly I means I a or I s . In addition, when, for instance, = a then c = s and vice versa. Then, by maximising the log likelihood with respect to β we have that then, From this expression we obtain that Therefore, (11) Meanwhile, deriving the log likelihood with respect to p yields This implies that with , I a (t 0 ), I s (t 0 ))} initial state where all populations classes are strictly positive. Denote by ϕ := {µ, β s , β a , κ, p, θ, α s , α a , γ} , a model parameter configuration. Accordingly, to the van den Driessche definition, the reproductive number for the deterministic version of our dynamics results Thus given the initial condition X + 0 and a parameter configuration ϕ, we define the set Ω * = Ω(X + 0 , ϕ, T * ) by To prove consistency and asymptotic normality of the parameter, we will make the following hypotheses. Hypotheses 1. There exists T 0 > 0 such that for all t ∈ [0, T 0 ] 1. R D 0 > 1 2. The initial condition X + 0 and parameter configuration ϕ are such that Ω * = ∅ This implies that our results are only valid in a sufficiently large time window. Remark 1. In other words, Hypothesis 1 assures conditions to estimate parameters in a growth phase of the outbreak. Under these conditions, the susceptible population decreases while the Exposed and Infected classes are in the growth phase. Furthermore, according to the definition of set Ω * , this growth phase occurs in a time window of size (T * − t 0 ). Hypotheses 1 implies that Now we prove the consistency of our estimators. Theorem 2. Assume that Hypotheses 1 are satisfied. The estimators (β s,M L ,β a,M L ,p M L ) are strongly consistent; that is, The proof of this theorem is deferred to supplementary material. To estimate the parameter of diffusion σ, we can use the quadratic variation over [0, T ] of the processes of system of SDEs denoted by < * , * > T . The estimator is given byσ wherê Here, we use the Milstein scheme Kloeden and Platen (1992) to generate paths of processes of system (3) in a time interval [0, T ]. Then we apply (14) to estimate σ and expressions in eqs. (20) to (22) to compute the MLE of β s , β a , and p. Using the Milstein scheme, we divide the interval [0, T ] into n sub-intervals of length ∆ = T n and we obtain observations of each process at the times 0 = t 0 < t 1 < . . . < t n = T where t i − t i−1 = ∆, i = 1, . . . , n. Based on these observations the diffusion parameter σ can be estimated using the expression (14) wherê we deduce thatσ 2 is an unbiased estimator of σ 2 (see Miao, 2006) . Meanwhile, to obtain the MLE of β a , β s , and p, we use the discrete time version of expressions (20), (21) and (22). Now, we present two examples with different parameters to be estimated and initial conditions. In both examples, we use the parameters of Table 2 for µ, κ, α a , α s and γ. Example 1. We simulate a 1000 datasets in the time interval [0, 1]. To generate the datasets we use β a = 0.251521, β s = 0.45639, p = 0.1213, σ = 1/5000 and ∆ = 1/850 for the approximation in Milstein method and the initial conditions from Table 4 . Using the observation at times t 0 = 0, t 1 = ∆, t 2 = 2∆, . . . , t 850 = 1, the values of Table 2 for µ, κ, α a , α s , γ and ∆ = 1/850, and the discrete version of expressions (14), (20), (21) and (22) we find the corresponding estimators. The average of the 1, 000 datasets and the corresponding standard deviations are presented in Table 5 . (3). Combining I s and Milstein scheme, we generate the others four processes of system (3) using Algorithm 1. With the parameters of Table 2 , σ = 1/100 and ∆ = 1/1000 Algorithm 1 works as follows. Algorithm 1 Construction of dataset. 1: Choose the initial values of Table 6 for E(0), I a (0) and R(0), . 7: If n < 47 make n = n + 1 and go to 2. In other case stop. With this procedure, we have 47 records in the same observation period for each of the five processes that should fit to the stochastic model given by (3). In other words, we have observations of the processes of system (3) at the times t 0 = 0, t 1 = 0.001, t 2 = 0.002, . . ., t 46 = 0.046. We assume the following parameters as given: γ = 1/365, κ = 0.196 078, p = 0.585 504 9, α a = 0.167 504 α s = 0.092 506 9. Thus, using the described methodology in Section 6.1 to generate data and applying (14), we estimate σ and using the expressions (20), (21), (22) we estimate β s , β a , and p. Next, we generate 1000 datasets using the data of Mexico City. Table 8 gives the average of these estimators and the corresponding 95% confidence interval. We refer the reader to https://chart-studio.plotly.com/~sauld/47 for a plot visualisation. closer way than the corresponding Bayesian calibration. Further, all observations drop in the corresponding confidence interval. So, this suggests that our proposed methodology and the stochastic extension improves the Bayesian fit of this dataset. To validate our stochastic model, we compute the residuals corresponding to data of reported and confirmed COVID-19 cases of Mexico City. Our target is to show that the residuals are Gaussian. If this is the case, then the noise W in our stochastic model would be consistent with Brownian motion. (16) and a theoretical Gaussian distribution distribution with zero mean and variance 1. Using equation (4), we define the increments where ∆ = t k − t k−1 and f β (t k ) is the evaluation of the function f β in time t k ; that is, If the right-hand side of (16) is Gaussian N (0, ∆) for k = 1, . . . , 46 , then this will imply that the increments are Gaussian and also implies that the stochastic process W is a Brownian motion, and therefore the stochastic model is wellfitted. We note that in expression (16) there are four components of the stochastic processes, namely S, R, I a , I s , but we only observed the component I s and the other three processes had to be approximated using the methodology from the previous sections. We use the parameter values of Section 6.2, estimators of Table 8 and ∆ = 0.001 into (16) and we get the residual ∆ W (t k ). Then, we prove that ∆ W (t k ) are Gaussian with zero mean and finite variance, which leaves us to conclude that the model fits properly to the actual data. Indeed, Figure 3 shows a QQplot of standardised residual of 1000 paths of the real data versus a sample of same size of standard normal random variable. Since an important part of the quantile frequency in Figure 3 follows the reference profile, we can conclude that the residuals are Gaussian. This suggests that our stochastic model is well-fitted to the data. We studied the MLE for some key parameters of a stochastic SEIR model. These parameters are important to study and quantify uncertainty related to the stochastic nature of data. By using synthetic data, we verified that our approximation of the MLE is a good estimator for the infection transmission, symptomatic ratio and noise intensity parameters. In the supplementary material, we provide a detailed proof for the Gaussian consistency of our estimators and other important properties. We then fitted our stochastic model to a dataset of COVID-19 confirmed symptomatic cases from México City. Because this information is limited and incomplete in the context of the regarding model, the database lacks of reliable information of the susceptible, exposed, asymptomatic and recover classes, we implemented the ideas as was explained in Section 6.1. Thus, via simulation, we obtain confidence intervals that capture more observations than the given MCMC outcome with its respective deterministic version, see Figure 2 . One of our goals was to improve the uncertainty quantification of COVID-19 data by a stochastic SEIR structure based on Itô SDEs. Our simulations suggest that our MLE-estimators improve the variance description of a prescribed overfitting MCMC outcome. Previous results confirmed this behavior in other applications, see for example Chatzilena et al. (2019) . Although our findings do not implicate a general methodology, most of the presented techniques would be applicable to problems of the same sort. The estimators developed here capture the noisy nature of the data with more precision. Furthermore, our simulations run in much less time with respect to the MCMC run with the STAN implementation. However, our conclusions rely on a very particular input for the underlying outcome, and we depend on the algebraic form of a given model to construct an estimator. We have confirmed the well-fitness of our stochastic model to real data, via error residual computing. However, it is necessary to explore our ideas in other epidemic structures, and with other datasets and target parameters. We should also develop more and stronger theories to construct consistent and computable versions of our estimators. As was discussed in the introduction, Pan et al. reports the estimation of a simple stochastic SIS structure that needs only a scalar version of the model to approximate the likelihood. There are more complex stochastic models similar to the one studied here, but they have not addressed estimation methods, see for instance Faranda and Alberti (2020) or Djević et al. (2021) . In contrast, our results deal with an epidemic model of five non-linear coupled SDEs. Furthermore, we combine the stochastic model with real data for COVID-19 from Mexico City and, using a natural discretised version of the MLE, we estimate the parameters. Therefore we fit the model to the actual data, and we obtain a good behavior of the stochastic model. Otunuga also derives and calibrates a similar SEIR stochastic model for COVID-19 in Otunuga (2021) . Their methodology relies on in the Stratonovich stochastic calculus and new estimation scheme , which they named 'lagged adapted generalized method of moments' (LLGMM). The LLGMM is a dynamic non-parametric method and seems to be a well-behaved method; however, the properties of the estimated parameters (e.g. consistency, bias, efficiency, etc.) are (up to our knowledge) still an open question. Since we base our proposal in the Itô framework, and the MLE method, the results of our contribution are complementary. One limitation of our work is that we have only one noise in all of the system. We then we pretend to extend the model in such manner that each process has an independent noise. One open question is if such model could be susceptible to be studied with the methodology developed in this work, which means applying the MLE method. In addition, it will be necessary to study if the model is well-fitted to the actual data. We also want to investigate if the model studied in this manuscript (and a possible extension) could be applied to other type of diseases or phenomena. Another question is to determine, via simulation, the minimum number of observations that ensure a nice approximation to the true parameters. Furthermore, this work missed a proof of the asymptotic normality of the estimators, which we guess that could be achieved by using the so-called Malliavin-Stein method. In addition, we want to study rates of convergence of the discretised version of the estimator to the continuous one. This section is devoted to present the proof of Theorem 2. For the sake of completeness, we first present a version of a classical Strong Law of Large Numbers that we use in getting our result. The proof of this particular version can be found, for instance, in Shiryaev (Shiryaev, 1996, Theorem IV.3 .2)). Theorem 3 (Strong Law of Large Numbers). Let ξ k , k ≥ 1, be independent random variables with the following properties: • Eξ k = 0, Eξ 2 k > 0, • There exist real numbers c > 0 and α ≥ −1 such that Then, with probability one, If, in addition, Eξ 4 k ≤ c 1 Eξ 2 k 2 for all k ≥ 1, with c 1 > 0 independent of k, then, also with probability one, We now present the proof of the consistency of the MLE parameters. Proof. (Of Theorem 2) We first will get some estimates for J (T ) and J sa (T ). Observe that using Hypotheses 1 we obtain In a very similar manner 2T. which implies that At other hand, since (1 − S(t)) 2 < 1 and E 2 (t) < 1 then 1 (1 − S(t)) 2 > 1, and For J sa we obtain a similar bound From this estimations we get We now focus on the quantityβ s − β s,0 , using (11) we can writê We have that We note that I 1 (T ) and I 2 (T ) are stochastic processes with zero mean and finite variance. Indeed, by the Itô isometry we get, for instance for I 1 (T ) Define the sequence of real numbers t 0 = 0 and t k = T (1 − 1 k ) for k ≥ 1. Define the random variables Thus, as before E ξ k = 0 and Then, we have that the following limit is true This means that there exists C > 0 such that lim k→∞ k −α E ξ 2 k = C, with α = −1. Then, by using the Theorem 3 we have that, with probability one, We observe that (1 − S(0))E(0) Then, from (17) we deduce that t N +1 0 (1−S(t)) + S(t)I (t) with probability one. This implies that Λ 1,1 + Λ 1,2 → 0 which consequently proves thatβ s → β s,0 with probability one. A similar procedure proves thatβ a → β a,0 with probability one. We now focus on the proof of the consistency forp. We rewrite (12) aŝ with J 2 as defined before and We will show that the first term in the right side of (18) goes to zero with probability one and the second is bounded with probability one. As before, for the sequence {t k } k≥0 we define the random variables and we have that E(ξ k ) = 0 and E(ξ 2 k ) > 0. Moreover, we can have the following estimate for E(ξ 2 k ): E(ξ 2 k ) < κ I s (0) + κ I a (0) and thus the following limit is satisfied lim k→∞ k E(ξ 2 k ) < κ I s (0) + κ I a (0) 2 T, Therefore, there exists a constant C > 0 such that lim k→∞ k E(ξ 2 k ) = C. Then, by the Theorem 3 we have that, with probability one, This implies that meaning that the first term in (18) goes to zero with probability one. For the second term in (18), observe that It is possible to prove (see supplementary material) that (1 − S(t)) µ + γR (1 − S) + 1 2 σ 2 dt Moreover, Modeling behavioral change and COVID-19 containment in Mexico: A trade-off between lockdown and compliance COVID-19 optimal vaccination policies: A modeling study on efficacy, natural and vaccine-induced immunity responses Modeling with Itô stochastic differential equations A conceptual introduction to hamiltonian monte carlo Parameter estimation in stochastic differential equations vaccines Meet Big Data: State-ofthe-Art and Future Prospects. From the Classical 3Is ('Isolate-Inactivate-Inject') vaccinology 1.0 to vaccinology 3.0, vaccinomics, and Beyond: A Historical Overview A new way of investigating the asymptotic behaviour of a stochastic sis system with multiplicative noise Contemporary statistical inference for infectious disease models using stan A two diffusion stochastic model for the spread of the new corona virus sars-cov-2 Modeling the second wave of covid-19 infections in france and italy via a stochastic seir model Efficient data augmentation for fitting stochastic epidemic models to prevalence data Stochastic methods. Springer Series in Synergetics Datos abiertos Random ordinary differential equations and their numerical solution The no-U-turn sampler: adaptively setting path lengths in Hamiltonian Monte Carlo Bayesian melding estimation of a stochastic seir model Simulation and inference for stochastic differential equations Numerical solution of stochastic differential equations Statistical Inference for Multivariate Stochastic Differential Equations Stationary distribution and extinction of a stochastic dengue epidemic model Estimation of diffusion parameters in diffusion processes and their asymptotic normality Analysis of sdes applied to seir epidemic models by extended kalman filter method Stochastic differential equations Introduction and snapshot review: Relating infectious disease transmission models to data Estimation of epidemiological parameters for COVID-19 cases using a stochastic SEIRS epidemic model with vital dynamics Parameter estimation for the stochastic SIS epidemic model Stochastic anlysis and statistical inference for seir models of infectious diseases Applied stochastic differential equations Probability Parametric inference for diffusion processes observed at discrete points in time: a survey An investigation of transmission control measures during the first 50 days of the COVID-19 epidemic in China Parameter estimation for continuous-time models: a survey Stochastic asymptotic analysis of a multi-host model with vector transmission