Lesson 2: Simulation Of Stochastic Dynamic Models
Lesson 2:Simulation of stochastic dynamic modelsAaron A. King, Edward L. Ionides, and Kidus Asfaw1 / 54
Outline1Compartment modelsExample: the SIR modelNotationA deterministic interpretationA stochastic interpretation2Euler’s methodNumerical solution of deterministic dynamicsNumerical solution of stochastic dynamics3Compartment models in pompA basic pomp model for measlesC snippetsChoosing parameters4Exercises2 / 54
ObjectivesThis tutorial develops some classes of dynamic models relevant tobiological systems, especially for epidemiology.1Dynamic systems can often be represented in terms of flows betweencompartments.2We develop the concept of a compartmental model for which wespecify rates for the flows between compartments.3We show how deterministic and stochastic versions of acompartmental model are derived and related.4We introduce Euler’s method to simulate from dynamic models.5We specify deterministic and stochastic compartmental models inpomp using Euler method simulation.3 / 54
Compartment modelsExample: the SIR modelA basic compartment model: The SIR modelWe develop deterministic and stochastic representations of asusceptible-infected-recovered (SIR) system, a fundamental class ofmodels for disease transmission dynamics.We set up notation applicable to general compartment models (Bretóet al., 2009).4 / 54
Compartment modelsExample: the SIR modelA basic compartment model: The SIR model II C 6ρ-µ SSµ? S -µSIS : susceptibleR : recovered and/or removedIµ? I -µIRRµ? R I : infected and infectiousC : reported cases5 / 54
Compartment modelsExample: the SIR modelA basic compartment model: The SIR model IIIWe suppose that each arrow has an associated rate, so here there is arate µSI (t) at which individuals in S transition to I, and µIR atwhich individuals in I transition to R.To account for demography (births/deaths/migration) we allow thepossibility of a source and sink compartment, which is not usuallyrepresented on the flow diagram. We write µ S for a rate of birthsinto S, and denote mortality rates by µS , µI , µR .6 / 54
Compartment modelsExample: the SIR modelA basic compartment model: The SIR model IVThe rates may be either constant or varying. In particular, for asimple SIR model, the recovery rate µIR is a constant but theinfection rate has the time-varying formµSI (t) β I(t),with β being the contact rate. For the simplest SIR model, ignoringdemography, we setµ S µS µI µR 0.7 / 54
Compartment modelsNotationGeneral notation for compartment modelsTo develop a systematic notation, it turns out to be convenient tokeep track of the flows between compartments as well as the numberof individuals in each compartment. LetNSI (t)count the number of individuals who have transitioned from S to I bytime t. We say that NSI (t) is a counting process. A similarlyconstructed processNIR (t)8 / 54
Compartment modelsNotationGeneral notation for compartment models IIcounts individuals transitioning from I to R. To include demography,we could keep track of birth and death events by the countingprocesses N S (t), NS (t), NI (t), NR (t).For discrete population compartment models, the flow countingprocesses are non-decreasing and integer valued.For continuous population compartment models, the flow countingprocesses are non-decreasing and real valued.9 / 54
Compartment modelsNotationCompartment processes from counting processesThe numbers of people in each compartment can be computed viathese counting processes. Ignoring demography, we have:S(t) S(0) NSI (t)I(t) I(0) NSI (t) NIR (t)R(t) R(0) NIR (t)These equations represent conservation of individuals or what goes inmust come out.10 / 54
Compartment modelsA deterministic interpretationOrdinary differential equation interpretationTogether with initial conditions specifying S(0), I(0) and R(0), we justneed to write down ordinary differential equations (ODEs) for the flowcounting processes. These are:dNSI µSI (t) S(t)dtdNIR µIR I(t)dt11 / 54
Compartment modelsA stochastic interpretationContinuous-time Markov chain interpretationContinuous-time Markov chains are the basic tool for building discretepopulation epidemic models.The Markov property lets us specify a model by the transitionprobabilities on small intervals (together with the initial conditions).For the SIR model, we have P NSI (t δ) NSI (t) 1 µSI (t) S(t) δ o(δ) P NSI (t δ) NSI (t) 1 µSI (t) S(t) δ o(δ) P NIR (t δ) NIR (t) 1 µIR I(t) δ o(δ) P NIR (t δ) NIR (t) 1 µIR (t) I(t) δ o(δ)Here, we are using little o notation We write h(δ) o(δ) to meanlimδ 0 h(δ)δ 0.12 / 54
Compartment modelsA stochastic interpretationExercise 2.1What is the link between little o notation and the derivative? Explain whyf (x δ) f (x) δg(x) o(δ)is the same statement asdf g(x).dxWhat considerations might help you choose which of these notations touse?Worked solution to the Exercise13 / 54
Compartment modelsA stochastic interpretationSimple counting processesA simple counting process is one which cannot count more thanone event at a time.Technically, the SIR Markov chain model we have written is simple.One may want to model the extra randomness resulting from multiplesimultaneous events: someone sneezing in a bus; large gatherings atfootball matches; etc. This extra randomness may even be critical tomatch the variability in data.Later in the course, we may see situations where this extrarandomness plays an important role. Setting up the model usingcounting processes, as we have done here, turns out to be useful forthis.14 / 54
Euler’s methodOutline1Compartment modelsExample: the SIR modelNotationA deterministic interpretationA stochastic interpretation2Euler’s methodNumerical solution of deterministic dynamicsNumerical solution of stochastic dynamics3Compartment models in pompA basic pomp model for measlesC snippetsChoosing parameters4Exercises15 / 54
Euler’s methodNumerical solution of deterministic dynamicsEuler’s method for ordinary differential equationsEuler (1707–1783) wanted a numeric solution of an ordinarydifferential equation (ODE) dx/dt h(x) with an initial conditionx(0).He supposed this ODE has some true solution x(t) which could notbe worked out analytically. He wanted an approximation x̃(t) of x(t).He initialized the numerical solution at the known starting value,x̃(0) x(0).16 / 54
Euler’s methodNumerical solution of deterministic dynamicsEuler’s method for ordinary differential equations IIFor k 1, 2, . . . , he supposed that the gradient dx/dt isapproximately constant over the small time intervalkδ t (k 1)δ. Therefore, he defined x̃ (k 1)δ x̃(kδ) δ h x̃(kδ) .This only defines x̃(t) when t is a multiple of δ, but suppose x̃(t) isconstant between these discrete times.We now have a numerical scheme, stepping forwards in timeincrements of size δ, that can be readily evaluated by computer.17 / 54
Euler’s methodNumerical solution of deterministic dynamicsEuler’s method versus other numerical methodsMathematical analysis of Euler’s method says that, as long as thefunction h(x) is not too exotic, then x(t) is well approximated by x̃(t)when the discretization time-step, δ, is sufficiently small.Euler’s method is not the only numerical scheme to solve ODEs.More advanced schemes have better convergence properties, meaningthat the numerical approximation is closer to x(t). However, there are3 reasons we choose to lean heavily on Euler’s method:123Euler’s method is the simplest (cf. the KISS principle).Euler’s method extends naturally to stochastic models, bothcontinuous-time Markov chains models and stochastic differentialequation (SDE) models.Close approximation of the numerical solutions to a continuous-timemodel is less important than it may at first appear, a topic to bediscussed.18 / 54
Euler’s methodNumerical solution of deterministic dynamicsContinuous-time models and discretized approximationsIn some physical and engineering situations, a system follows an ODEmodel closely. For example, Newton’s laws provide a very goodapproximation to the motions of celestial bodies.In many biological situations, ODE models only become closemathematical approximations to reality at reasonably large scale. Onsmall temporal scales, models cannot usually capture the full scope ofbiological variation and biological complexity.If we are going to expect substantial error in using x(t) to model abiological system, maybe the numerical solution x̃(t) represents thesystem being modeled as well as x(t) does.19 / 54
Euler’s methodNumerical solution of deterministic dynamicsContinuous-time models and discretized approximations IIIf our model fitting, model investigation, and final conclusions are allbased on our numerical solution x̃(t) (i.e., we are sticking entirely tosimulation-based methods) then we are most immediately concernedwith how well x̃(t) describes the system of interest. x̃(t) becomesmore important than the original model, x(t).20 / 54
Euler’s methodNumerical solution of deterministic dynamicsNumerical solutions as scientific modelsIt is important that a scientist fully describe the numerical modelx̃(t). Arguably, the main purpose of the original model x(t) is to givea succinct description of how x̃(t) was constructed.All numerical methods are, ultimately, discretizations.Epidemiologically, setting δ to be a day, or an hour, can be quitedifferent from setting δ to be two weeks or a month. Forcontinuous-time modeling, we still require that δ is small compared tothe timescale of the process being modeled, so the choice of δ shouldnot play an explicit role in the interpretation of the model.Putting more emphasis on the scientific role of the numerical solutionitself reminds you that the numerical solution has to do more thanapproximate a target model in some asymptotic sense: the numericalsolution should be a sensible model in its own right.21 / 54
Euler’s methodNumerical solution of deterministic dynamicsEuler’s method for a discrete SIR modelRecall the simple continuous-time Markov chain interpretation of theSIR model without demography: P NSI (t δ) NSI (t) 1 µSI (t) S(t)δ o(δ), P NIR (t δ) NIR (t) 1 µIR I(t)δ o(δ). We want a numerical solution with state variables S̃(kδ), I(kδ),R̃(kδ).22 / 54
Euler’s methodNumerical solution of deterministic dynamicsEuler’s method for a discrete SIR model IIThe counting processes for the flows between compartments areÑSI (t) and ÑIR (t). The counting processes are related to thenumbers of individuals in the compartments by the same flowequations we had before: S(0) ÑSI (kδ) I(0) ÑSI (kδ) ÑIR (kδ)R̃(kδ) R(0) ÑIR (kδ)S̃(kδ) I(kδ)We focus on a numerical solution to NSI (t), since the same methodscan also be applied to NIR (t).23 / 54
Euler’s methodNumerical solution of stochastic dynamicsThree different stochastic Euler solutions(1) A Poisson approximation. ÑSI (t δ) ÑSI (t) Poisson µSI I(t)S̃(t) δ ,where Poisson(µ) is a Poisson random variable with mean µ and µSI I(t) β I(t).(2) A binomial approximation, ÑSI (t δ) ÑSI (t) Binomial S̃(t), µSI I(t)δ ,where Binomial(n, p) is a binomial random variable with mean np and variance np(1 p). Here, p µSI I(t)δ.(3) A binomial approximation with exponential transition probabilities. δ .ÑSI (t δ) ÑSI (t) Binomial S̃(t), 1 exp µSI I(t)Analytically, it is usually easiest to reason using (1) or (2).Practically, it is usually preferable to work with (3).24 / 54
Euler’s methodNumerical solution of stochastic dynamicsCompartment models as stochastic differential equationsThe Euler method extends naturally to stochastic differentialequations (SDEs).A natural way to add stochastic variation to an ODE dx/dt h(x) isdBdX h(X) σdtdtwhere {B(t)} is Brownian motion and so dB/dt is Brownian noise.An Euler approximation X̃(t) is X̃ (k 1)δ X̃(kδ) δ h X̃kδ) σ δ Zkwhere Z1 , Z2 , . . . is a sequence of independent standard normalrandom variables, i.e., Zk N [0, 1].25 / 54
Euler’s methodNumerical solution of stochastic dynamicsCompartment models as stochastic differential equations IIAlthough SDEs are often considered an advanced topic in probability,the Euler approximation doesn’t demand much more than familiaritywith the normal distribution.26 / 54
Euler’s methodNumerical solution of stochastic dynamicsExercise 2.2. Euler’s method vs Gillespie’s algorithmA widely used, exact simulation method for continuous time Markovchains is Gillespie’s algorithm. We do not put much emphasis onGillespie’s algorithm here. Why? When would you prefer animplementation of Gillespie’s algorithm to an Euler solution?Worked solution to the ExerciseNumerically, Gillespie’s algorithm is often approximated using so-calledtau-leaping methods. These are closely related to Euler’s approach. In thiscontext, the Euler method has sometimes been called tau-leaping.27 / 54
Compartment models in pompOutline1Compartment modelsExample: the SIR modelNotationA deterministic interpretationA stochastic interpretation2Euler’s methodNumerical solution of deterministic dynamicsNumerical solution of stochastic dynamics3Compartment models in pompA basic pomp model for measlesC snippetsChoosing parameters4Exercises28 / 54
Compartment models in pompA basic pomp model for measlesThe Consett measles outbreakAs an example that we can probe in some depth, let’s look at outbreak ofmeasles that occurred in the small town of Consett in England in 1948.The town had population of 38820, with 737 births over the course of theyear.29 / 54
Compartment models in pompA basic pomp model for measlesThe Consett measles outbreak IIWe download the data and examine them:library(tidyverse)read m/","Measles Consett 1948.csv")) % %select(week,reports cases) - measmeas % % as.data.frame() % % head()week reports10203240536030 / 54
Compartment models in pompA basic pomp model for measlesThe Consett measles outbreak III31 / 54
Compartment models in pompA basic pomp model for measlesA simple POMP model for measlesThese are incidence data: The reports variable counts the numberof reports of new measles cases each week.Let us model the outbreak using the simple SIR model.Our tasks will be, first, to estimate the parameters of the SIR and,second, to decide whether or not the SIR model is an adequatedescription of these data.The rate at which individuals move from S to I is the force ofinfection, µSI β I/N , while that at which individuals move intothe R class is µIR .32 / 54
Compartment models in pompA basic pomp model for measlesFraming the SIR as a POMP modelThe unobserved state variables, in this case, are the numbers ofindividuals, S(t), I(t), R(t) in the S, I, and R compartments,respectively.It’s reasonable in this case to view the population sizeN S(t) I(t) R(t), as fixed at the known population size of38,000.The numbers that actually move from one compartment to anotherover any particular time interval are modeled as stochastic processes.In this case, we’ll assume that the stochasticity is purelydemographic, i.e., that each individual in a compartment at any giventime faces the same risk of exiting the compartment.Demographic stochasticity is the unavoidable randomness thatarises from chance events occurring in a discrete and finite population.33 / 54
Compartment models in pompA basic pomp model for measlesImplementing the SIR model in pompTo implement the model in pomp, the first thing we need is astochastic simulator for the unobserved state process.We follow method 3 above, modeling the number, NSI , movingfrom S to I over interval t as I NSI Binomial S, 1 e β N t ,and the number moving from I to R as NIR Binomial I, 1 e µIR t .34 / 54
Compartment models in pompA basic pomp model for measlesImplementing the SIR model in pomp IIsir step - function (S, I, R, N, Beta, mu IR, delta.t, .) {dN SI - rbinom(n 1,size S,prob 1-exp(-Beta*I/N*delta.t))dN IR - rbinom(n 1,size I,prob 1-exp(-mu IR*delta.t))S - S - dN SII - I dN SI - dN IRR - R dN IRc(S S, I I, R R)}At day zero, we’ll assume that I 1 but we don’t know how manypeople are susceptible, so we’ll treat this fraction, η, as a parameterto be estimated.sir rinit - function (N, eta, .) {c(S round(N*eta), I 1, R round(N*(1-eta)))}35 / 54
Compartment models in pompA basic pomp model for measlesImplementing the SIR model in pomp IIIWe fold these basic model components, with the data, into a pompobject thus:library(pomp)meas % %pomp(times "week",t0 0,rprocess euler(sir step,delta.t 1/7),rinit sir rinit) - measSIRNow assume the case reports result from a process by which newinfections are diagnosed and reported with probability ρ, which wecan think of as the probability that a child’s parents take the child tothe doctor, who recognizes measles and reports it to the authorities.36 / 54
Compartment models in pompA basic pomp model for measlesImplementing the SIR model in pomp IVMeasles symptoms tend to be quite recognizable, and children withmeasles tend to be confined to bed. Therefore diagnosed cases have,presumably, a much lower transmission rate. Accordingly, let’s treateach week’s reports as being related to the number of individualswho have moved from I to R over the course of that week.We need a variable to track these daily counts. We modify ourrprocess function above, adding a variable H to tally the trueincidence.37 / 54
Compartment models in pompA basic pomp model for measlesImplementing the SIR model in pomp Vsir step - function (S, I, R, H, N, Beta, mu IR, delta.t, .){dN SI - rbinom(n 1,size S,prob 1-exp(-Beta*I/N*delta.t))dN IR - rbinom(n 1,size I,prob 1-exp(-mu IR*delta.t))S - S - dN SII - I dN SI - dN IRR - R dN IRH - H dN IR;c(S S, I I, R R, H H)}sir rinit - function (N, eta, .) {c(S round(N*eta), I 1, R round(N*(1-eta)), H 0)}38 / 54
Compartment models in pompA basic pomp model for measlesImplementing the SIR model in pomp VIIn pomp terminology, H is an accumulator variable. Since we wantH to tally only the incidence over the week, we’ll need to reset it tozero at the beginning of each week. We accomplish this using theaccumvars argument to pomp:measSIR % %pomp(rprocess euler(sir step,delta.t 1/7),rinit sir rinit,accumvars "H") - measSIRNow, we’ll model the data as a binomial process,reportst Binomial (H(t), ρ) .39 / 54
Compartment models in pompA basic pomp model for measlesImplementing the SIR model in pomp VIINow, to include the observations in the model, we must write either admeasure or an rmeasure component, or both:sir dmeas - function (reports, H, rho, log, .) {dbinom(x reports, size H, prob rho, log log)}sir rmeas - function (H, rho, .) {c(reports rbinom(n 1, size H, prob rho))}40 / 54
Compartment models in pompA basic pomp model for measlesImplementing the SIR model in pomp VIIIWe then put these into our pomp object:measSIR % %pomp(rmeasure sir rmeas,dmeasure sir dmeas) - measSIR41 / 54
Compartment models in pompC snippetsSpecifying model components using C snippetsAlthough we can always specify basic model components using Rfunctions, as above, we’ll typically want the computational speed-upthat we can obtain only by using compiled native code.pomp provides a facility for doing so with ease, using C snippets.C snippets are small pieces of C code used to specify basic modelcomponents.For example, a C snippet encoding the rprocess for an sir model isas follows.42 / 54
Compartment models in pompC snippetsSpecifying model components using C snippets IIsir step - Csnippet("double dN SI rbinom(S,1-exp(-Beta*I/N*dt));double dN IR rbinom(I,1-exp(-mu IR*dt));S - dN SI;I dN SI - dN IR;R dN IR;H dN IR;")43 / 54
Compartment models in pompC sn
A deterministic interpretation A stochastic interpretation 2 Euler’s method Numerical solution of deterministic dynamics Numerical solution of stochastic dynamics 3 Compartment models in pomp A basic pomp model for measles C snippets Choosing parameters 4 Exercises 15/54
4 Step Phonics Quiz Scores Step 1 Step 2 Step 3 Step 4 Lesson 1 Lesson 2 Lesson 3 Lesson 4 Lesson 5 Lesson 6 Lesson 7 Lesson 8 Lesson 9 Lesson 10 Lesson 11 Lesson 12 Lesson 13 Lesson 14 Lesson 15 . Zoo zoo Zoo zoo Yoyo yoyo Yoyo yoyo You you You you
Jul 09, 2010 · Stochastic Calculus of Heston’s Stochastic–Volatility Model Floyd B. Hanson Abstract—The Heston (1993) stochastic–volatility model is a square–root diffusion model for the stochastic–variance. It gives rise to a singular diffusion for the distribution according to Fell
are times when the fast stochastic lines either cross above 80 or below 20, while the slow stochastic lines do not. By slowing the lines, the slow stochastic generates fewer trading signals. INTERPRETATION You can see in the figures that the stochastic oscillator fluctuates between zero and 100. A stochastic value of 50 indicates that the closing
Participant's Workbook Financial Management for Managers Institute of Child Nutrition iii Table of Contents Introduction Intro—1 Lesson 1: Financial Management Lesson 1—1 Lesson 2: Production Records Lesson 2—1 Lesson 3: Forecasting Lesson 3—1 Lesson 4: Menu Item Costs Lesson 4—1 Lesson 5: Product Screening Lesson 5—1 Lesson 6: Inventory Control Lesson 6—1
sion analysis on discrete-time stochastic processes. We now turn our focus to the study of continuous-time stochastic pro-cesses. In most cases, it is di cult to exactly describe the probability dis-tribution for continuous-time stochastic processes. This was also di cult for discrete time stochastic processes, but for them, we described the .
(e.g. bu er stocks, schedule slack). This approach has been criticized for its use of a deterministic approximation of a stochastic problem, which is the major motivation for stochastic programming. This dissertation recasts this debate by identifying both deterministic and stochastic approaches as policies for solving a stochastic base model,
Stochastic Programming Stochastic Dynamic Programming Conclusion : which approach should I use ? Objective and constraints Evaluating a solution Presentation Outline 1 Dealing with Uncertainty Objective and constraints Evaluating a solution 2 Stochastic Programming Stochastic Programming Approach Information Framework Toward multistage program
STOCHASTIC CALCULUS AND STOCHASTIC DIFFERENTIAL EQUATIONS 5 In discrete stochastic processes, there are many random times similar to (2.3). They are non-anticipating, i.e., at any time n, we can determine whether the cri-terion for such a random time is met or not solely by the “history” up to time n.
