From 0b7ed6ecd609a235a5de7ba4d99802aa5eb28286 Mon Sep 17 00:00:00 2001 From: Jesse Perla Date: Mon, 22 Jun 2020 23:17:45 -0700 Subject: [PATCH 01/16] First stab at the ODEs --- source/rst/tools_and_techniques/index.rst | 1 + .../tools_and_techniques/seir_model_sde.rst | 469 ++++++++++++++++++ 2 files changed, 470 insertions(+) create mode 100644 source/rst/tools_and_techniques/seir_model_sde.rst diff --git a/source/rst/tools_and_techniques/index.rst b/source/rst/tools_and_techniques/index.rst index cefec116..7019bd82 100644 --- a/source/rst/tools_and_techniques/index.rst +++ b/source/rst/tools_and_techniques/index.rst @@ -20,6 +20,7 @@ tools and techniques. :maxdepth: 2 linear_algebra + seir_model_sde orth_proj lln_clt linear_models diff --git a/source/rst/tools_and_techniques/seir_model_sde.rst b/source/rst/tools_and_techniques/seir_model_sde.rst new file mode 100644 index 00000000..3e1b5e11 --- /dev/null +++ b/source/rst/tools_and_techniques/seir_model_sde.rst @@ -0,0 +1,469 @@ +.. include:: /_static/includes/header.raw + +.. highlight:: julia + +***************************************************************** +:index:`Modeling COVID 19 with (Stochastic) Differential Equations` +***************************************************************** + +.. contents:: :depth: 2 + + + +Overview +============= + +This is a Julia version of the code for analyzing the COVID-19 pandemic. + +The purpose of these notes is to introduce economists to quantitative modeling +of infectious disease dynamics, and to modeling with ordinary and stochastic differential +equations. + +The main objective is to study the impact of suppression through social +distancing on the spread of the infection. + +Dynamics are modeled using a standard SEIR (Susceptible-Exposed-Infected-Removed) model +of disease spread, represented as a system of ordinary differential +equations when the number of agents is large and there are no exogenous stochastic shocks. + +The focus is on US outcomes but the parameters can be adjusted to study +other countries. + + +The first part of the model follows the notes from +provided by `Andrew Atkeson `__ + +* `NBER Working Paper No. 26867 `__ +* `COVID-19 Working papers and code `__ + +See further variations on the classic SIR model in Julia `here `__. + + + +Setup +------------------ + +.. literalinclude:: /_static/includes/deps_generic.jl + :class: hide-output + +.. code-block:: julia + + using LinearAlgebra, Statistics, Random, SparseArrays + +.. code-block:: julia + :class: Test + + using Test # Put this before any code in the lecture. + +In addition, we will be exploring packages within the `SciML ecosystem `__. + +.. code-block:: julia + + using OrdinaryDiffEq, StochasticDiffEq + + +The SEIR Model +============= + +In the version of the SEIR model we will analyze there are four states. + +All individuals in the population are assumed to be in one of these four states. + +The states are: susceptible (S), exposed (E), infected (I) and removed (R). + +Comments: + +* Those in state R have been infected and either recovered or died. + +* Those who have recovered are assumed to have acquired immunity. + +* Those in the exposed group are not yet infectious. + +Time Path +---------- + +The flow across states follows the path :math:`S \to E \to I \to R`. + + +All individuals in the population are eventually infected when +the transmission rate is positive and :math:`i(0) > 0`. + +The interest is primarily in + +* the number of infections at a given time (which determines whether or not the health care system is overwhelmed) and +* how long the caseload can be deferred (hopefully until a vaccine arrives) + +Using lower case letters for the fraction of the population in each state where we maintain a constant population throughout, the +dynamics are + +.. math:: + \begin{aligned} + \frac{d s}{d t} & = - \beta \, s \, i + \\ + \frac{d e}{d t} & = \beta \, s \, i - \sigma e + \\ + \frac{d i}{d t} & = \sigma e - \gamma i + \\ + \frac{d r}{d t} & = \gamma i + \end{aligned} + :label: seir_system + +In these equations, + +* :math:`\beta(t)` is called the *transmission rate* (the rate at which individuals bump into others and expose them to the virus). +* :math:`\sigma` is called the *infection rate* (the rate at which those who are exposed become infected) +* :math:`\gamma` is called the *recovery rate* (the rate at which infected people recover or die). +* :math:`dy/dt` represents the time derivative for the particular variable + + +Note that the states form a partition, so we can reconstruct the "removed" fraction of the population is :math:`r = 1 - s - e - i`. + +However, it may be convenient to leave the :math:`r(t)` in the system as well as :math:`c = i + r`, which is the cumulative caseload +(i.e., all those who have or have had the infection). Differentiating that expression and substituing from the time-derivatives of :math:`i(t), r(t)` yields :math:`\frac{d c}{d t} = \sigma e` + +We will assume that the transmission rate follows a process with a reversion to a mean :math:`b` which will remain constant for now, but could conceivably be a policy parameter. + +.. math:: + \begin{aligned} + \frac{d \beta}{d t} &= \eta (b - \beta) + \end{aligned} + +Finally, let :math:`v(t)` be the mortality rate, which we will leave constant for now. The cumulative deaths can be integrated through the flow :math:`\gamma i` entering the "Removed" state and define the cumulative number of deaths as :math:`m(t)`. The differential equation +follows, + +.. math:: + + \begin{aligned} + \frac{d m}{d t} &= v \gamma i + \end{aligned} + +While we could conveivably integate the total deaths given the solution to the model, it is convenient to use the integrator built into the ODE solver. + + +The system :eq:`seir_system` and the supplemental equations can be written in vector form in terms of the vector :math:`x := (s, e, i, r, \beta, c, m)` with parameter vector :math:`p := (\sigma, \gamma, b, eta)` + +.. math:: + \frac{d x}{d t} = \begin{bmatrix} + - \beta \, s \, i + \\ + \beta \, s \, i - \sigma e + \\ + \sigma e - \gamma i + \\ + \gamma i + \\ + \eta (b - \beta) + \\ + \sigma e + \\ + v \gamma i + \end{bmatrix} =: F(x; p) + + :label: dfcv + +Here we have maintained the time independence of the :math:`F(x)` function, but we could also have time-varying terms. + +Parameters +---------- + +Both :math:`\sigma` and :math:`\gamma` are thought of as fixed, biologically determined parameters. + +As in Atkeson's note, we set + +* :math:`\sigma = 1/5.2` to reflect an average incubation period of 5.2 days. +* :math:`\gamma = 1/18` to match an average illness duration of 18 days. +* :math:`b = 1.6 \gamma` to match an **effective reproduction rate** of 1.6 + +In addition, the transmission rate can be interpreted as + +* :math:`\beta(t) := R(t) \gamma` where :math:`R(t)` is the *effective reproduction number* at time :math:`t`. + +(The notation is standard in the epidemiology literature - though slightly confusing, since :math:`R(t)` is different to +:math:`R`, the symbol that represents the removed state.) + +Rather than set :math:`\eta`, we will begin by looking at the case where :math:`\beta(0) = \bar{\beta}`, and hence it remains constant. + + +Implementation +============== + +First we set the population size to match the US. + +.. code-block:: julia + + pop_size = 3.3e8 + +Next we fix parameters as described above. + +.. code-block:: julia + + γ = 1 / 18 + σ = 1 / 5.2 + +Now we construct a function that represents :math:`F` in :eq:`dfcv` + +.. code-block:: julia + + #def F(x, t, R0=1.6): + # """ + # Time derivative of the state vector. + # + # * x is the state vector (array_like) + # * t is time (scalar) + # * R0 is the effective transmission rate, defaulting to a constant + # + # """ + # s, e, i = x + # + # # New exposure of susceptibles + # β = R0(t) * γ if callable(R0) else R0 * γ + # ne = β * s * i + # + # # Time derivatives + # ds = - ne + # de = ne - σ * e + # di = σ * e - γ * i + # + # return ds, de, di + +Note that ``R0`` can be either constant or a given function of time. + +The initial conditions are set to + +.. code-block:: julia + + # initial conditions of s, e, i + i_0 = 1e-7 + e_0 = 4.0 * i_0 + s_0 = 1.0 - i_0 - e_0 + +In vector form the initial condition is + +.. code-block:: julia + + x_0 = s_0, e_0, i_0 + +We solve for the time path numerically using `odeint`, at a sequence of dates +``t_vec``. + +.. code-block:: julia + + #def solve_path(R0, t_vec, x_init=x_0): + # """ + # Solve for i(t) and c(t) via numerical integration, + # given the time path for R0. + # + # """ + # G = lambda x, t: F(x, t, R0) + # s_path, e_path, i_path = odeint(G, x_init, t_vec).transpose() + # + # c_path = 1 - s_path - e_path # cumulative cases + # return i_path, c_path + + + +Experiments +=========== + +Let's run some experiments using this code. + +The time period we investigate will be 550 days, or around 18 months: + +.. code-block:: julia + + t_length = 550 + grid_size = 1000 + #t_vec = np.linspace(0, t_length, grid_size) + + + +Experiment 1: Constant R0 Case +------------------------------ + + +Let's start with the case where ``R0`` is constant. + +We calculate the time path of infected people under different assumptions for ``R0``: + +.. code-block:: julia + + #R0_vals = np.linspace(1.6, 3.0, 6) + #labels = [f'$R0 = {r:.2f}$' for r in R0_vals] + #i_paths, c_paths = [], [] + # + #for r in R0_vals: + # i_path, c_path = solve_path(r, t_vec) + # i_paths.append(i_path) + # c_paths.append(c_path) + +Here's some code to plot the time paths. + +.. code-block:: julia + + #def plot_paths(paths, labels, times=t_vec): + # + # fig, ax = plt.subplots() + # + # for path, label in zip(paths, labels): + # ax.plot(times, path, label=label) + # + # ax.legend(loc='upper left') + # + # plt.show() + +Let's plot current cases as a fraction of the population. + +.. code-block:: julia + + #plot_paths(i_paths, labels) + +As expected, lower effective transmission rates defer the peak of infections. + +They also lead to a lower peak in current cases. + +Here is cumulative cases, as a fraction of population: + +.. code-block:: julia + + #plot_paths(c_paths, labels) + + + +Experiment 2: Changing Mitigation +--------------------------------- + +Let's look at a scenario where mitigation (e.g., social distancing) is +successively imposed. + +Here's a specification for ``R0`` as a function of time. + +.. code-block:: julia + + #def R0_mitigating(t, r0=3, η=1, r_bar=1.6): + # R0 = r0 * exp(- η * t) + (1 - exp(- η * t)) * r_bar + # return R0 + +The idea is that ``R0`` starts off at 3 and falls to 1.6. + +This is due to progressive adoption of stricter mitigation measures. + +The parameter ``η`` controls the rate, or the speed at which restrictions are +imposed. + +We consider several different rates: + +.. code-block:: julia + + #η_vals = 1/5, 1/10, 1/20, 1/50, 1/100 + #labels = [fr'$\eta = {η:.2f}$' for η in η_vals] + +This is what the time path of ``R0`` looks like at these alternative rates: + +.. code-block:: julia + + #fig, ax = plt.subplots() + # + #for η, label in zip(η_vals, labels): + # ax.plot(t_vec, R0_mitigating(t_vec, η=η), label=label) + # + #ax.legend() + #plt.show() + +Let's calculate the time path of infected people: + +.. code-block:: julia + + #i_paths, c_paths = [], [] + # + #for η in η_vals: + # R0 = lambda t: R0_mitigating(t, η=η) + # i_path, c_path = solve_path(R0, t_vec) + # i_paths.append(i_path) + # c_paths.append(c_path) + + +This is current cases under the different scenarios: + +.. code-block:: julia + + #plot_paths(i_paths, labels) + +Here are cumulative cases, as a fraction of population: + +.. code-block:: julia + + #plot_paths(c_paths, labels) + + + +Ending Lockdown +=============== + + +The following replicates `additional results `__ by Andrew Atkeson on the timing of lifting lockdown. + +Consider these two mitigation scenarios: + +1. :math:`R_t = 0.5` for 30 days and then :math:`R_t = 2` for the remaining 17 months. This corresponds to lifting lockdown in 30 days. + +2. :math:`R_t = 0.5` for 120 days and then :math:`R_t = 2` for the remaining 14 months. This corresponds to lifting lockdown in 4 months. + +The parameters considered here start the model with 25,000 active infections +and 75,000 agents already exposed to the virus and thus soon to be contagious. + +.. code-block:: julia + + # initial conditions + i_0 = 25_000 / pop_size + e_0 = 75_000 / pop_size + s_0 = 1.0 - i_0 - e_0 + x_0 = s_0, e_0, i_0 + +Let's calculate the paths: + +.. code-block:: julia + + #R0_paths = (lambda t: 0.5 if t < 30 else 2, + # lambda t: 0.5 if t < 120 else 2) + # + #labels = [f'scenario {i}' for i in (1, 2)] + # + #i_paths, c_paths = [], [] + # + #for R0 in R0_paths: + # i_path, c_path = solve_path(R0, t_vec, x_init=x_0) + # i_paths.append(i_path) + # c_paths.append(c_path) + + +Here is the number of active infections: + +.. code-block:: julia + + #plot_paths(i_paths, labels) + +What kind of mortality can we expect under these scenarios? + +Suppose that 1\% of cases result in death + +.. code-block:: julia + + #ν = 0.01 + +This is the cumulative number of deaths: + +.. code-block:: julia + + #paths = [path * ν * pop_size for path in c_paths] + #plot_paths(paths, labels) + +This is the daily death rate: + +.. code-block:: julia + + #paths = [path * ν * γ * pop_size for path in i_paths] + #plot_paths(paths, labels) + +Pushing the peak of curve further into the future may reduce cumulative deaths +if a vaccine is found. + From 078e525082b230a13f570669b76985a1fc213fdc Mon Sep 17 00:00:00 2001 From: Jesse Perla Date: Mon, 22 Jun 2020 23:35:16 -0700 Subject: [PATCH 02/16] System of ODEs for the simple case --- .../tools_and_techniques/seir_model_sde.rst | 19 ++++++++++++++++++- 1 file changed, 18 insertions(+), 1 deletion(-) diff --git a/source/rst/tools_and_techniques/seir_model_sde.rst b/source/rst/tools_and_techniques/seir_model_sde.rst index 3e1b5e11..4d0d86f5 100644 --- a/source/rst/tools_and_techniques/seir_model_sde.rst +++ b/source/rst/tools_and_techniques/seir_model_sde.rst @@ -204,7 +204,24 @@ Now we construct a function that represents :math:`F` in :eq:`dfcv` .. code-block:: julia - #def F(x, t, R0=1.6): + function F(t, x, p) + s, e, i, r, β, c, m = x + σ, γ, b, η = p + return [-β * s * i; # for d s(t) + β * s * i - σ * e; + σ * e - γ * i; + γ * i; + η (b - β) + v * γ * i + σ * e] + + + \frac{d m}{d t} &= v \gamma i + + end + + + # """ # Time derivative of the state vector. # From 83abbebda817054489a2a4472a370a7d9774785a Mon Sep 17 00:00:00 2001 From: Jesse Perla Date: Tue, 23 Jun 2020 15:27:47 -0700 Subject: [PATCH 03/16] Converted to the b, specification --- .../tools_and_techniques/seir_model_sde.rst | 250 ++++++++---------- 1 file changed, 112 insertions(+), 138 deletions(-) diff --git a/source/rst/tools_and_techniques/seir_model_sde.rst b/source/rst/tools_and_techniques/seir_model_sde.rst index 4d0d86f5..0b7fea4d 100644 --- a/source/rst/tools_and_techniques/seir_model_sde.rst +++ b/source/rst/tools_and_techniques/seir_model_sde.rst @@ -2,9 +2,9 @@ .. highlight:: julia -***************************************************************** +****************************************************************** :index:`Modeling COVID 19 with (Stochastic) Differential Equations` -***************************************************************** +****************************************************************** .. contents:: :depth: 2 @@ -55,11 +55,13 @@ Setup using Test # Put this before any code in the lecture. -In addition, we will be exploring packages within the `SciML ecosystem `__. +In addition, we will be exploring packages within the `SciML ecosystem `__ and +others covered in previous lectures .. code-block:: julia using OrdinaryDiffEq, StochasticDiffEq + using Parameters, StaticArrays, Plots The SEIR Model @@ -79,13 +81,13 @@ Comments: * Those in the exposed group are not yet infectious. -Time Path ----------- +Changes in the Infected State +------------------------------- The flow across states follows the path :math:`S \to E \to I \to R`. -All individuals in the population are eventually infected when +Since we are using a continuum approximation, all individuals in the population are eventually infected when the transmission rate is positive and :math:`i(0) > 0`. The interest is primarily in @@ -93,8 +95,7 @@ The interest is primarily in * the number of infections at a given time (which determines whether or not the health care system is overwhelmed) and * how long the caseload can be deferred (hopefully until a vaccine arrives) -Using lower case letters for the fraction of the population in each state where we maintain a constant population throughout, the -dynamics are +Assume that there is a constant population of size :math:`N` throughout, then define the proportion of people in each state as :math:`s := S/N` etc. With this, the SEIR model can be written as .. math:: \begin{aligned} @@ -114,54 +115,63 @@ In these equations, * :math:`\sigma` is called the *infection rate* (the rate at which those who are exposed become infected) * :math:`\gamma` is called the *recovery rate* (the rate at which infected people recover or die). * :math:`dy/dt` represents the time derivative for the particular variable +* Since the states form a partition, so we can reconstruct the "removed" fraction of the population as :math:`r = 1 - s - e - i`. However, it convenient to :math:`r(t)` in the system for graphing + +In addition, we are interested in calculatin the cumulative caseload (i.e., all those who have or have had the infection) as :math:`c = i + r`. Differentiating that expression and substituing from the time-derivatives of :math:`i(t), r(t)` yields :math:`\frac{d c}{d t} = \sigma e` -Note that the states form a partition, so we can reconstruct the "removed" fraction of the population is :math:`r = 1 - s - e - i`. -However, it may be convenient to leave the :math:`r(t)` in the system as well as :math:`c = i + r`, which is the cumulative caseload -(i.e., all those who have or have had the infection). Differentiating that expression and substituing from the time-derivatives of :math:`i(t), r(t)` yields :math:`\frac{d c}{d t} = \sigma e` +Implementing the system of ODEs in :math:`s, e, i` would be enough to implement the model, but we will extend the basic model to enable some policy experiments. -We will assume that the transmission rate follows a process with a reversion to a mean :math:`b` which will remain constant for now, but could conceivably be a policy parameter. +Evolution of Parameters +----------------------- + +We will assume that the transmission rate follows a process with a reversion to a value :math:`b` which could conceivably be a policy parameter. The intuition is that even if the targetted :math:`b(t)` was changed, lags in behavior and implementation would smooth out the transition, where :math:`\eta` governs the speed of :math:`\beta(t)` moves towards :math:`b(t)`. .. math:: \begin{aligned} \frac{d \beta}{d t} &= \eta (b - \beta) \end{aligned} -Finally, let :math:`v(t)` be the mortality rate, which we will leave constant for now. The cumulative deaths can be integrated through the flow :math:`\gamma i` entering the "Removed" state and define the cumulative number of deaths as :math:`m(t)`. The differential equation -follows, +Finally, let :math:`v(t)` be the mortality rate, which we will leave constant for now, i.e. :math:`\frac{d v}{d t} = 0`. The cumulative deaths can be integrated through the flow :math:`\gamma i` entering the "Removed" state and define the cumulative number of deaths as :math:`m(t)`. The differential equations then +follow, .. math:: - \begin{aligned} + \begin{aligned}\\ + \frac{d v}{d t} &= 0\\ \frac{d m}{d t} &= v \gamma i \end{aligned} -While we could conveivably integate the total deaths given the solution to the model, it is convenient to use the integrator built into the ODE solver. +While we could conveivably integate the total deaths given the solution to the model, it is more convenient to use the integrator built into the ODE solver. That is, we added :math:`d m(t)/dt` rather than calculating :math:`m(t) = \int_0^t \gamma v(\tau) i(\tau) d \tau` after generating the full :math:`i(t)` path. +This is a common trick when solving systems of ODEs. While equivalent in principle if you used an appropriate quadrature scheme, this trick becomes especially important and convenient when adaptive time-stepping algorithms are used to solve the ODEs (i.e. there is no fixed time grid). -The system :eq:`seir_system` and the supplemental equations can be written in vector form in terms of the vector :math:`x := (s, e, i, r, \beta, c, m)` with parameter vector :math:`p := (\sigma, \gamma, b, eta)` +The system :eq:`seir_system` and the supplemental equations can be written in vector form in terms of the vector :math:`x := (s, e, i, r, c, m, \beta, v)` with parameter vector :math:`p := (\sigma, \gamma, b, \eta)` .. math:: - \frac{d x}{d t} = \begin{bmatrix} + \begin{aligned} + \frac{d x}{d t} = F(x,t;p) := \begin{bmatrix} - \beta \, s \, i \\ \beta \, s \, i - \sigma e \\ - \sigma e - \gamma i + \sigma \, e - \gamma i \\ - \gamma i - \\ - \eta (b - \beta) + \gamma \, i \\ \sigma e \\ - v \gamma i - \end{bmatrix} =: F(x; p) - + v \, \gamma \, i + \\ + \eta (b(t) - \beta) + \\ + 0 + \end{bmatrix} + \end{aligned} :label: dfcv -Here we have maintained the time independence of the :math:`F(x)` function, but we could also have time-varying terms. +Here note that if :math:`b(t)` is time-invariant, then :math:`F(x)` is time-invariant as well. Parameters ---------- @@ -172,137 +182,123 @@ As in Atkeson's note, we set * :math:`\sigma = 1/5.2` to reflect an average incubation period of 5.2 days. * :math:`\gamma = 1/18` to match an average illness duration of 18 days. -* :math:`b = 1.6 \gamma` to match an **effective reproduction rate** of 1.6 +* :math:`\bar{b} / \gamma = 1.6` to match an **effective reproduction rate** of 1.6, and initially time-invariant +* :math:`v = 0.01` for a one-percent mortality rate In addition, the transmission rate can be interpreted as -* :math:`\beta(t) := R(t) \gamma` where :math:`R(t)` is the *effective reproduction number* at time :math:`t`. +* :math:`R(t) := \beta(t) / \gamma` where :math:`R(t)` is the *effective reproduction number* at time :math:`t`. (The notation is standard in the epidemiology literature - though slightly confusing, since :math:`R(t)` is different to -:math:`R`, the symbol that represents the removed state.) - -Rather than set :math:`\eta`, we will begin by looking at the case where :math:`\beta(0) = \bar{\beta}`, and hence it remains constant. +:math:`R`, the symbol that represents the removed state. Throughout the rest of the lecture, we will always use :math:`R` to represent the reproduction number) +As we will initially consider the case where :math:`\beta(0) = \bar{b}`, the value of :math:`\eta` will drop out of this first experiment. Implementation ============== -First we set the population size to match the US. - -.. code-block:: julia - - pop_size = 3.3e8 - -Next we fix parameters as described above. +First we set the population size to match the US and the parameters as described .. code-block:: julia + N = 3.3e8 # US Population γ = 1 / 18 σ = 1 / 5.2 + η = 1 / 20 # a placeholder, drops out of firs texperiments. Now we construct a function that represents :math:`F` in :eq:`dfcv` .. code-block:: julia - function F(t, x, p) - s, e, i, r, β, c, m = x - σ, γ, b, η = p - return [-β * s * i; # for d s(t) - β * s * i - σ * e; - σ * e - γ * i; - γ * i; - η (b - β) - v * γ * i - σ * e] - - - \frac{d m}{d t} &= v \gamma i - + # Reminder: could solve dynamics of SEIR states with just first 3 equations + function F(u, p, t) + s, e, i, r, c, m, β, v = u + @unpack σ, γ, b, η = p + + return [-β * s * i; # ds/dt = -βsi + β * s * i - σ * e; # de/dt = βsi - σe + σ * e - γ * i; # di/dt = σe -γi + γ * i; # dr/dt = γi + σ * e; # dc/dt = σe + v * γ * i; # dm/dt = vγi + η * (b(t, p) - β); # dβ/dt = η(b(t) - β) + 0.0 # dv/dt = 0 + ] end +The baseline parameters are put into a named tuple generator (see previous lectures using ``Parameters.jl``) with default values discussed above. - # """ - # Time derivative of the state vector. - # - # * x is the state vector (array_like) - # * t is time (scalar) - # * R0 is the effective transmission rate, defaulting to a constant - # - # """ - # s, e, i = x - # - # # New exposure of susceptibles - # β = R0(t) * γ if callable(R0) else R0 * γ - # ne = β * s * i - # - # # Time derivatives - # ds = - ne - # de = ne - σ * e - # di = σ * e - γ * i - # - # return ds, de, di +.. code-block:: julia -Note that ``R0`` can be either constant or a given function of time. + (t,p) = p.b̄ + p_gen = @with_kw (T = 550.0, γ = 1.0 / 18, σ = 1 / 5.2, η = 1.0 / 20, + b̄ = 1.6 * γ, b = (t, p) -> p.b̄) -The initial conditions are set to +Note that the default :math:`b(t)` function is simply the constant function :math:`\bar{b}` + + +Setting initial conditions, we will assume a fixed :math:`i, e` along with +assuming :math:`r = m = c = 0`, and that :math:`\beta(0) = \bar{b}` and :math:`v(0) = 0.01` .. code-block:: julia - # initial conditions of s, e, i + p = p_gen() # use all default parameters + i_0 = 1e-7 e_0 = 4.0 * i_0 s_0 = 1.0 - i_0 - e_0 -In vector form the initial condition is + u_0 = [s_0, e_0, i_0, 0.0, 0.0, 0.0, p.b̄, 0.01] + tspan = (0.0, p.T) + prob = ODEProblem(F, u_0, tspan, p) -.. code-block:: julia - x_0 = s_0, e_0, i_0 +The ``tspan`` determines that the :math:`t` used by the sovler, where the scale needs to be consistent with the arrival +rate of the transition probabilities (i.e. the :math:`\gamma, \sigma` were chosen based on daily data). +The time period we investigate will be 550 days, or around 18 months: -We solve for the time path numerically using `odeint`, at a sequence of dates -``t_vec``. +Experiments +=========== -.. code-block:: julia +Let's run some experiments using this code. - #def solve_path(R0, t_vec, x_init=x_0): - # """ - # Solve for i(t) and c(t) via numerical integration, - # given the time path for R0. - # - # """ - # G = lambda x, t: F(x, t, R0) - # s_path, e_path, i_path = odeint(G, x_init, t_vec).transpose() - # - # c_path = 1 - s_path - e_path # cumulative cases - # return i_path, c_path +First, we can solve the ODE using an appropriate algorthm (e.g. a good default for non-stiff ODEs might be ``Tsit5()``, which is the Tsitouras 5/4 Runge-Kutta method). +Most high-performance ODE solvers appropriate for this class of problems will have adaptive time-stepping, so you +will not specify any sort of grid +.. code-block:: julia -Experiments -=========== + sol = solve(prob, Tsit5()) # TODO: change the accuracy? + @show length(sol.t); -Let's run some experiments using this code. +We see that the adaptive time-stepping used approximately 45 time-steps to solve this problem to the desires accuracy. Evaluating the solver at points outside of those time-steps uses the an interpolator consistent with the +solution to the ODE. -The time period we investigate will be 550 days, or around 18 months: +See `here `__ for details on analyzing the solution, and `here `__ for plotting tools. The built-in plots for the solutions provide all of the `attributes `__. -.. code-block:: julia +.. code_block:: julia - t_length = 550 - grid_size = 1000 - #t_vec = np.linspace(0, t_length, grid_size) + # TODO: Chris, We could plot something else? Also labels broken + plot(sol, vars = [1, 2, 3, 4], label = ["s", "i", "e", "r"], title = "SIER Proportions") -Experiment 1: Constant R0 Case +Experiment 1: Constant Reproduction Case ------------------------------ +Let's start with the case where :math:`R = \beta / b` is constant. -Let's start with the case where ``R0`` is constant. - -We calculate the time path of infected people under different assumptions for ``R0``: +We calculate the time path of infected people under different assumptions. .. code-block:: julia + γ_base = 1.0/18.0 + R_vals = range(1.6, 3.0, length = 6) + b_vals = R_vals / γ_base + sols = [solve(ODEProblem(F, u_0, tspan, p_gen(b = b_vals)), Tsit5()) for b in b_vals] + + # TODO: Probably clean ways to plot this #R0_vals = np.linspace(1.6, 3.0, 6) #labels = [f'$R0 = {r:.2f}$' for r in R0_vals] @@ -313,9 +309,7 @@ We calculate the time path of infected people under different assumptions for `` # i_paths.append(i_path) # c_paths.append(c_path) -Here's some code to plot the time paths. - -.. code-block:: julia + # Here's some code to plot the time paths. #def plot_paths(paths, labels, times=t_vec): # @@ -345,24 +339,25 @@ Here is cumulative cases, as a fraction of population: #plot_paths(c_paths, labels) - Experiment 2: Changing Mitigation --------------------------------- Let's look at a scenario where mitigation (e.g., social distancing) is successively imposed. -Here's a specification for ``R0`` as a function of time. - -.. code-block:: julia +To do this, we will have :math:`\beta(0) \neq b` and examine the dynamics using the :math:`\frac{d \beta}{d t} &= \eta (b - \beta)` differential equation. - #def R0_mitigating(t, r0=3, η=1, r_bar=1.6): - # R0 = r0 * exp(- η * t) + (1 - exp(- η * t)) * r_bar - # return R0 +.. Mathematica Verification +.. (\[Beta][t] /. +.. First@DSolve[{\[Beta]'[t] == \[Eta] (b - \[Beta][t]), \[Beta][ +.. 0] == \[Beta]0}, \[Beta][t], +.. t] ) == \[Beta]0 E^(-t \[Eta]) + (1 - +.. E^(-t \[Eta])) b // FullSimplify -The idea is that ``R0`` starts off at 3 and falls to 1.6. + +Note that in the simple case, where :math:`b` is independent of the state, the solution to the ODE with :math:`\beta(0) = \beta_0` is :math:`\beta(t) = \beta_0 e^{-\eta t} + b(1 - e^{-\eta t})` -This is due to progressive adoption of stricter mitigation measures. +We will examine the case where :math:`R(t)` starts off at 3 and falls to 1.6 due to the progressive adoption of stricter mitigation measures. The parameter ``η`` controls the rate, or the speed at which restrictions are imposed. @@ -374,19 +369,7 @@ We consider several different rates: #η_vals = 1/5, 1/10, 1/20, 1/50, 1/100 #labels = [fr'$\eta = {η:.2f}$' for η in η_vals] -This is what the time path of ``R0`` looks like at these alternative rates: - -.. code-block:: julia - - #fig, ax = plt.subplots() - # - #for η, label in zip(η_vals, labels): - # ax.plot(t_vec, R0_mitigating(t_vec, η=η), label=label) - # - #ax.legend() - #plt.show() - -Let's calculate the time path of infected people: +Let's calculate the time path of infected people, current cases, and mortality .. code-block:: julia @@ -398,17 +381,8 @@ Let's calculate the time path of infected people: # i_paths.append(i_path) # c_paths.append(c_path) - -This is current cases under the different scenarios: - -.. code-block:: julia - #plot_paths(i_paths, labels) -Here are cumulative cases, as a fraction of population: - -.. code-block:: julia - #plot_paths(c_paths, labels) From 1108bbadf882594411a8b2051062db9a94c1c144 Mon Sep 17 00:00:00 2001 From: Jesse Perla Date: Tue, 23 Jun 2020 16:28:59 -0700 Subject: [PATCH 04/16] Converted to R instead of beta --- .../tools_and_techniques/seir_model_sde.rst | 151 ++++++++++-------- 1 file changed, 81 insertions(+), 70 deletions(-) diff --git a/source/rst/tools_and_techniques/seir_model_sde.rst b/source/rst/tools_and_techniques/seir_model_sde.rst index 0b7fea4d..d96588c7 100644 --- a/source/rst/tools_and_techniques/seir_model_sde.rst +++ b/source/rst/tools_and_techniques/seir_model_sde.rst @@ -81,27 +81,44 @@ Comments: * Those in the exposed group are not yet infectious. + + +The interest is primarily in + +* the number of infections at a given time (which determines whether or not the health care system is overwhelmed) and +* how long the caseload can be deferred (hopefully until a vaccine arrives) + + + Changes in the Infected State ------------------------------- -The flow across states follows the path :math:`S \to E \to I \to R`. +Within the SIER model, the flow across states follows the path :math:`S \to E \to I \to R`. Extensions of this model, such as SIERS, relax lifetime immunity and allow transitions from :math:` R \to S`. +The transitions between those states are governed by the following rates -Since we are using a continuum approximation, all individuals in the population are eventually infected when -the transmission rate is positive and :math:`i(0) > 0`. +* :math:`\beta(t)` is called the *transmission rate* (the rate at which individuals bump into others and expose them to the virus). +* :math:`\sigma` is called the *infection rate* (the rate at which those who are exposed become infected) +* :math:`\gamma` is called the *recovery rate* (the rate at which infected people recover or die). -The interest is primarily in -* the number of infections at a given time (which determines whether or not the health care system is overwhelmed) and -* how long the caseload can be deferred (hopefully until a vaccine arrives) +In addition, the transmission rate can be interpreted as: math:`R(t) := \beta(t) / \gamma` where :math:`R(t)` is the *effective reproduction number* at time :math:`t`. For this reason, we will work with the :math:`R(t)` reparameterization. + +(The notation is standard in the epidemiology literature - though slightly confusing, since :math:`R(t)` is different to +:math:`R`, the symbol that represents the removed state. Throughout the rest of the lecture, we will always use :math:`R` to represent the effective reproduction number, unless stated otherwise) + Assume that there is a constant population of size :math:`N` throughout, then define the proportion of people in each state as :math:`s := S/N` etc. With this, the SEIR model can be written as + +Since we are using a continuum approximation, all individuals in the population are eventually infected when +the transmission rate is positive and :math:`i(0) > 0`. + .. math:: \begin{aligned} - \frac{d s}{d t} & = - \beta \, s \, i + \frac{d s}{d t} & = - \gamma \, R \, s \, i \\ - \frac{d e}{d t} & = \beta \, s \, i - \sigma e + \frac{d e}{d t} & = \gamma \, R \, s \, i - \sigma e \\ \frac{d i}{d t} & = \sigma e - \gamma i \\ @@ -109,28 +126,23 @@ Assume that there is a constant population of size :math:`N` throughout, then de \end{aligned} :label: seir_system -In these equations, - -* :math:`\beta(t)` is called the *transmission rate* (the rate at which individuals bump into others and expose them to the virus). -* :math:`\sigma` is called the *infection rate* (the rate at which those who are exposed become infected) -* :math:`\gamma` is called the *recovery rate* (the rate at which infected people recover or die). -* :math:`dy/dt` represents the time derivative for the particular variable -* Since the states form a partition, so we can reconstruct the "removed" fraction of the population as :math:`r = 1 - s - e - i`. However, it convenient to :math:`r(t)` in the system for graphing +Here, :math:`dy/dt` represents the time derivative for the particular variable. +Since the states form a partition, we can reconstruct the "removed" fraction of the population as :math:`r = 1 - s - e - i`, and will drop it from the system -In addition, we are interested in calculatin the cumulative caseload (i.e., all those who have or have had the infection) as :math:`c = i + r`. Differentiating that expression and substituing from the time-derivatives of :math:`i(t), r(t)` yields :math:`\frac{d c}{d t} = \sigma e` +In addition, we are interested in calculating the cumulative caseload (i.e., all those who have or have had the infection) as :math:`c = i + r`. Differentiating that expression and substituing from the time-derivatives of :math:`i(t), r(t)` yields :math:`\frac{d c}{d t} = \sigma e` -Implementing the system of ODEs in :math:`s, e, i` would be enough to implement the model, but we will extend the basic model to enable some policy experiments. +Implementing the system of ODEs in :math:`s, e, i` would be enough to implement the model given fixed parameters, but we will extend the basic model to enable some policy experiments. Evolution of Parameters ----------------------- -We will assume that the transmission rate follows a process with a reversion to a value :math:`b` which could conceivably be a policy parameter. The intuition is that even if the targetted :math:`b(t)` was changed, lags in behavior and implementation would smooth out the transition, where :math:`\eta` governs the speed of :math:`\beta(t)` moves towards :math:`b(t)`. +We will assume that the transmission rate follows a process with a reversion to a value :math:`B(t)` which could conceivably be a policy parameter. The intuition is that even if the targetted :math:`B(t)` was changed, lags in behavior and implementation would smooth out the transition, where :math:`\eta` governs the speed of :math:`R(t)` moves towards :math:`B(t)`. .. math:: \begin{aligned} - \frac{d \beta}{d t} &= \eta (b - \beta) + \frac{d R}{d t} &= \eta (B - R) \end{aligned} Finally, let :math:`v(t)` be the mortality rate, which we will leave constant for now, i.e. :math:`\frac{d v}{d t} = 0`. The cumulative deaths can be integrated through the flow :math:`\gamma i` entering the "Removed" state and define the cumulative number of deaths as :math:`m(t)`. The differential equations then @@ -147,31 +159,29 @@ While we could conveivably integate the total deaths given the solution to the m This is a common trick when solving systems of ODEs. While equivalent in principle if you used an appropriate quadrature scheme, this trick becomes especially important and convenient when adaptive time-stepping algorithms are used to solve the ODEs (i.e. there is no fixed time grid). -The system :eq:`seir_system` and the supplemental equations can be written in vector form in terms of the vector :math:`x := (s, e, i, r, c, m, \beta, v)` with parameter vector :math:`p := (\sigma, \gamma, b, \eta)` +The system :eq:`seir_system` and the supplemental equations can be written in vector form in terms of the vector :math:`x := (s, e, i, c, m, R, v)` with parameter vector :math:`p := (\sigma, \gamma, B, \eta)` .. math:: \begin{aligned} \frac{d x}{d t} = F(x,t;p) := \begin{bmatrix} - - \beta \, s \, i + - \gamma \, R \, s \, i \\ - \beta \, s \, i - \sigma e + \gamma \, R \, s \, i - \sigma e \\ \sigma \, e - \gamma i \\ - \gamma \, i - \\ \sigma e \\ v \, \gamma \, i \\ - \eta (b(t) - \beta) + \eta (B(t) - R) \\ 0 \end{bmatrix} \end{aligned} :label: dfcv -Here note that if :math:`b(t)` is time-invariant, then :math:`F(x)` is time-invariant as well. +Here note that if :math:`B(t)` is time-invariant, then :math:`F(x)` is time-invariant as well. Parameters ---------- @@ -182,17 +192,11 @@ As in Atkeson's note, we set * :math:`\sigma = 1/5.2` to reflect an average incubation period of 5.2 days. * :math:`\gamma = 1/18` to match an average illness duration of 18 days. -* :math:`\bar{b} / \gamma = 1.6` to match an **effective reproduction rate** of 1.6, and initially time-invariant +* :math:`B = R = 1.6` to match an **effective reproduction rate** of 1.6, and initially time-invariant * :math:`v = 0.01` for a one-percent mortality rate -In addition, the transmission rate can be interpreted as -* :math:`R(t) := \beta(t) / \gamma` where :math:`R(t)` is the *effective reproduction number* at time :math:`t`. - -(The notation is standard in the epidemiology literature - though slightly confusing, since :math:`R(t)` is different to -:math:`R`, the symbol that represents the removed state. Throughout the rest of the lecture, we will always use :math:`R` to represent the reproduction number) - -As we will initially consider the case where :math:`\beta(0) = \bar{b}`, the value of :math:`\eta` will drop out of this first experiment. +As we will initially consider the case where :math:`R(0) = B`, the value of :math:`\eta` will drop out of this first experiment. Implementation ============== @@ -212,17 +216,17 @@ Now we construct a function that represents :math:`F` in :eq:`dfcv` # Reminder: could solve dynamics of SEIR states with just first 3 equations function F(u, p, t) - s, e, i, r, c, m, β, v = u - @unpack σ, γ, b, η = p - - return [-β * s * i; # ds/dt = -βsi - β * s * i - σ * e; # de/dt = βsi - σe - σ * e - γ * i; # di/dt = σe -γi - γ * i; # dr/dt = γi - σ * e; # dc/dt = σe - v * γ * i; # dm/dt = vγi - η * (b(t, p) - β); # dβ/dt = η(b(t) - β) - 0.0 # dv/dt = 0 + + s, e, i, c, m, R, v = u + @unpack σ, γ, B, η = p + + return [-γ*R*s*i; # ds/dt = -γRsi + γ*R*s*i - σ*e; # de/dt = γRsi - σe + σ*e - γ*i; # di/dt = σe -γi + σ*e; # dc/dt = σe + v*γ*i; # dm/dt = vγi + η*(B(t, p) - R);# dβ/dt = η(B(t) - R) + 0.0 # dv/dt = 0 ] end @@ -231,15 +235,13 @@ The baseline parameters are put into a named tuple generator (see previous lectu .. code-block:: julia - (t,p) = p.b̄ p_gen = @with_kw (T = 550.0, γ = 1.0 / 18, σ = 1 / 5.2, η = 1.0 / 20, - b̄ = 1.6 * γ, b = (t, p) -> p.b̄) + R̄ = 1.6, B = (t, p) -> p.R̄) -Note that the default :math:`b(t)` function is simply the constant function :math:`\bar{b}` +Note that the default :math:`B(t)` function is simply the constant function :math:`\bar{R}` -Setting initial conditions, we will assume a fixed :math:`i, e` along with -assuming :math:`r = m = c = 0`, and that :math:`\beta(0) = \bar{b}` and :math:`v(0) = 0.01` +Setting initial conditions, we will assume a fixed :math:`i, e`, :math:`r=0`, :math:`R(0) = \bar{B}`, and :math:`v(0) = 0.01` .. code-block:: julia @@ -249,7 +251,7 @@ assuming :math:`r = m = c = 0`, and that :math:`\beta(0) = \bar{b}` and :math:`v e_0 = 4.0 * i_0 s_0 = 1.0 - i_0 - e_0 - u_0 = [s_0, e_0, i_0, 0.0, 0.0, 0.0, p.b̄, 0.01] + u_0 = [s_0, e_0, i_0, 0.0, 0.0, p.R̄, 0.01] tspan = (0.0, p.T) prob = ODEProblem(F, u_0, tspan, p) @@ -280,23 +282,20 @@ See `here `__ for details on anal .. code_block:: julia - # TODO: Chris, We could plot something else? Also labels broken - plot(sol, vars = [1, 2, 3, 4], label = ["s", "i", "e", "r"], title = "SIER Proportions") - + # TODO: Chris, We could plot something else as well? Deaths, etc?? Also labels broken + plot(sol, vars = [3, 5], label = ["i(t)" "c(t)"], lw = 2, title = "Current and Cumulated Infected") Experiment 1: Constant Reproduction Case ------------------------------ -Let's start with the case where :math:`R = \beta / b` is constant. +Let's start with the case where :math:`B(t) = R = R̄` is constant. We calculate the time path of infected people under different assumptions. .. code-block:: julia - γ_base = 1.0/18.0 - R_vals = range(1.6, 3.0, length = 6) - b_vals = R_vals / γ_base - sols = [solve(ODEProblem(F, u_0, tspan, p_gen(b = b_vals)), Tsit5()) for b in b_vals] + R̄_vals = range(1.6, 3.0, length = 6) + sols = [solve(ODEProblem(F, u_0, tspan, p_gen(R̄ = R̄_val)), Tsit5()) for R̄ in R̄_vals] # TODO: Probably clean ways to plot this @@ -343,9 +342,9 @@ Experiment 2: Changing Mitigation --------------------------------- Let's look at a scenario where mitigation (e.g., social distancing) is -successively imposed. +successively imposed, but the target (:math:`R̄` is fixed) -To do this, we will have :math:`\beta(0) \neq b` and examine the dynamics using the :math:`\frac{d \beta}{d t} &= \eta (b - \beta)` differential equation. +To do this, we will have :math:`R(0) \neq R̄` and examine the dynamics using the :math:`\frac{d R}{d t} &= \eta (R̄ - R)` ODE. .. Mathematica Verification .. (\[Beta][t] /. @@ -355,9 +354,9 @@ To do this, we will have :math:`\beta(0) \neq b` and examine the dynamics using .. E^(-t \[Eta])) b // FullSimplify -Note that in the simple case, where :math:`b` is independent of the state, the solution to the ODE with :math:`\beta(0) = \beta_0` is :math:`\beta(t) = \beta_0 e^{-\eta t} + b(1 - e^{-\eta t})` +Note that in the simple case, where :math:`B(t) = \bar{R}` is independent of the state, the solution to the ODE with :math:`R(0) = R_0` is :math:`R(t) = R_0 e^{-\eta t} + \bar{R}(1 - e^{-\eta t})` -We will examine the case where :math:`R(t)` starts off at 3 and falls to 1.6 due to the progressive adoption of stricter mitigation measures. +We will examine the case where :math:`R(0) = 3` and then it falls to to :math:`\bar{R} = 1.6` due to the progressive adoption of stricter mitigation measures. The parameter ``η`` controls the rate, or the speed at which restrictions are imposed. @@ -391,24 +390,36 @@ Ending Lockdown =============== -The following replicates `additional results `__ by Andrew Atkeson on the timing of lifting lockdown. +The following is inspired by replicates `additional results `__ by Andrew Atkeson on the timing of lifting lockdown. Consider these two mitigation scenarios: -1. :math:`R_t = 0.5` for 30 days and then :math:`R_t = 2` for the remaining 17 months. This corresponds to lifting lockdown in 30 days. +1. choose :math:`B(t)` to target :math:`\bar{R} = 0.5` for 30 days and then :math:`\bar{R} = 2` for the remaining 17 months. This corresponds to lifting lockdown in 30 days. -2. :math:`R_t = 0.5` for 120 days and then :math:`R_t = 2` for the remaining 14 months. This corresponds to lifting lockdown in 4 months. +2. :math:`\bar{R} = 0.5` for 120 days and then :math:`\bar{R} = 2` for the remaining 14 months. This corresponds to lifting lockdown in 4 months. + +For both of these, we will choose a large :math:`\eta` to focus on the case where rapid changes in the lockdown policy remain feasible. The parameters considered here start the model with 25,000 active infections and 75,000 agents already exposed to the virus and thus soon to be contagious. .. code-block:: julia - # initial conditions - i_0 = 25_000 / pop_size - e_0 = 75_000 / pop_size + B_lift_early(t, p) = t < 30.0 ? 0.5 : 2.0 + B_lift_late(t, p) = t < 120.0 ? 0.5 : 2.0 + + + # initial conditions + i_0 = 25000 / N + e_0 = 75000 / N s_0 = 1.0 - i_0 - e_0 - x_0 = s_0, e_0, i_0 + + u_0 = [s_0, e_0, i_0, 0.0, 0.0, 0.5, 0.01] # starting in lockdown, R=0.5 + tspan = [0.0, 30.0, 120.0, p.T] # add discontinuities to tspan for adaptive solver + + # create two problems, with rapid movement of R towards B(t) + prob_early = ODEProblem(F, u_0, tspan, p_gen(B = B_lift_early, η = 10.0)) + prob_late = ODEProblem(F, u_0, tspan, p_gen(B = B_lift_late, η = 10.0) Let's calculate the paths: From 71c7225d0308fbf7f8cb3cec75e6aef088423b78 Mon Sep 17 00:00:00 2001 From: Jesse Perla Date: Tue, 23 Jun 2020 22:35:46 -0700 Subject: [PATCH 05/16] Added material on the SDE --- .../tools_and_techniques/seir_model_sde.rst | 136 ++++++++++++++---- 1 file changed, 105 insertions(+), 31 deletions(-) diff --git a/source/rst/tools_and_techniques/seir_model_sde.rst b/source/rst/tools_and_techniques/seir_model_sde.rst index d96588c7..55314161 100644 --- a/source/rst/tools_and_techniques/seir_model_sde.rst +++ b/source/rst/tools_and_techniques/seir_model_sde.rst @@ -93,7 +93,8 @@ The interest is primarily in Changes in the Infected State ------------------------------- -Within the SIER model, the flow across states follows the path :math:`S \to E \to I \to R`. Extensions of this model, such as SIERS, relax lifetime immunity and allow transitions from :math:` R \to S`. +Within the SIER model, the flow across states follows the path :math:`S \to E \to I \to R`. Extensions of this model, such as SIERS, relax lifetime immunity and allow +transitions from :math:`R \to S`. The transitions between those states are governed by the following rates @@ -102,7 +103,7 @@ The transitions between those states are governed by the following rates * :math:`\gamma` is called the *recovery rate* (the rate at which infected people recover or die). -In addition, the transmission rate can be interpreted as: math:`R(t) := \beta(t) / \gamma` where :math:`R(t)` is the *effective reproduction number* at time :math:`t`. For this reason, we will work with the :math:`R(t)` reparameterization. +In addition, the transmission rate can be interpreted as: :math:`R(t) := \beta(t) / \gamma` where :math:`R(t)` is the *effective reproduction number* at time :math:`t`. For this reason, we will work with the :math:`R(t)` reparameterization. (The notation is standard in the epidemiology literature - though slightly confusing, since :math:`R(t)` is different to :math:`R`, the symbol that represents the removed state. Throughout the rest of the lecture, we will always use :math:`R` to represent the effective reproduction number, unless stated otherwise) @@ -128,7 +129,7 @@ the transmission rate is positive and :math:`i(0) > 0`. Here, :math:`dy/dt` represents the time derivative for the particular variable. -Since the states form a partition, we can reconstruct the "removed" fraction of the population as :math:`r = 1 - s - e - i`, and will drop it from the system +Since the states form a partition, we could reconstruct the "removed" fraction of the population as :math:`r = 1 - s - e - i`. However, for further experiments and plotting it is harmless to keep it in the system. In addition, we are interested in calculating the cumulative caseload (i.e., all those who have or have had the infection) as :math:`c = i + r`. Differentiating that expression and substituing from the time-derivatives of :math:`i(t), r(t)` yields :math:`\frac{d c}{d t} = \sigma e` @@ -144,6 +145,7 @@ We will assume that the transmission rate follows a process with a reversion to \begin{aligned} \frac{d R}{d t} &= \eta (B - R) \end{aligned} + :label: Rode Finally, let :math:`v(t)` be the mortality rate, which we will leave constant for now, i.e. :math:`\frac{d v}{d t} = 0`. The cumulative deaths can be integrated through the flow :math:`\gamma i` entering the "Removed" state and define the cumulative number of deaths as :math:`m(t)`. The differential equations then follow, @@ -159,17 +161,20 @@ While we could conveivably integate the total deaths given the solution to the m This is a common trick when solving systems of ODEs. While equivalent in principle if you used an appropriate quadrature scheme, this trick becomes especially important and convenient when adaptive time-stepping algorithms are used to solve the ODEs (i.e. there is no fixed time grid). -The system :eq:`seir_system` and the supplemental equations can be written in vector form in terms of the vector :math:`x := (s, e, i, c, m, R, v)` with parameter vector :math:`p := (\sigma, \gamma, B, \eta)` +The system :eq:`seir_system` and the supplemental equations can be written in vector form :math:`x := [s, e, i, r, c, m, R, v]` with parameter tuple :math:`p := (\sigma, \gamma, B, \eta)` .. math:: \begin{aligned} - \frac{d x}{d t} = F(x,t;p) := \begin{bmatrix} - - \gamma \, R \, s \, i + \frac{d x}{d t} &= F(x,t;p) + &:= \begin{bmatrix} + - \gamma \, R \, s \, i \\ \gamma \, R \, s \, i - \sigma e \\ \sigma \, e - \gamma i \\ + \gamma i + \\ \sigma e \\ v \, \gamma \, i @@ -217,20 +222,31 @@ Now we construct a function that represents :math:`F` in :eq:`dfcv` # Reminder: could solve dynamics of SEIR states with just first 3 equations function F(u, p, t) - s, e, i, c, m, R, v = u + s, e, i, r, c, m, R, v = u @unpack σ, γ, B, η = p - return [-γ*R*s*i; # ds/dt = -γRsi - γ*R*s*i - σ*e; # de/dt = γRsi - σe - σ*e - γ*i; # di/dt = σe -γi - σ*e; # dc/dt = σe - v*γ*i; # dm/dt = vγi - η*(B(t, p) - R);# dβ/dt = η(B(t) - R) - 0.0 # dv/dt = 0 - ] + return [-γ*R*s*i; # ds/dt = -γRsi + γ*R*s*i - σ*e; # de/dt = γRsi - σe + σ*e - γ*i; # di/dt = σe -γi + γ*i; # dr/dt = γi + σ*e; # dc/dt = σe + v*γ*i; # dm/dt = vγi + η*(B(t, p) - R);# dβ/dt = η(B(t) - R) + 0.0 # dv/dt = 0 + ] end +Written this way, we see that the first four rows represent the one-directional transition from the susceptible to removed state, where the negative terms are outflows, and the positive ones inflows. + +As there is no flow leaving the :math:`dr/dt` and all parameters are positive, unless we start with a degenerate initial condition (e.g. :math:`e(0) = i(0) = 0`) the "Removed" state is asymptoically absorbing, and :math:`\lim_{t\to \infty} r(t) = 1`. Crucial to this rsult is that individuals are arbitrarily divisible, and any arbitrarily small :math:`i > 0` leads to a strictly positive flow into the exposed state. + +We will discuss this topic further in the lecture on continuous-time +markov-chains, as well as the limitations of these approximations when the discretness becomes essential + +Parameters +------------- + The baseline parameters are put into a named tuple generator (see previous lectures using ``Parameters.jl``) with default values discussed above. .. code-block:: julia @@ -238,7 +254,7 @@ The baseline parameters are put into a named tuple generator (see previous lectu p_gen = @with_kw (T = 550.0, γ = 1.0 / 18, σ = 1 / 5.2, η = 1.0 / 20, R̄ = 1.6, B = (t, p) -> p.R̄) -Note that the default :math:`B(t)` function is simply the constant function :math:`\bar{R}` +Note that the default :math:`B(t)` function always equals :math:`\bar{R}` Setting initial conditions, we will assume a fixed :math:`i, e`, :math:`r=0`, :math:`R(0) = \bar{B}`, and :math:`v(0) = 0.01` @@ -251,7 +267,7 @@ Setting initial conditions, we will assume a fixed :math:`i, e`, :math:`r=0`, :m e_0 = 4.0 * i_0 s_0 = 1.0 - i_0 - e_0 - u_0 = [s_0, e_0, i_0, 0.0, 0.0, p.R̄, 0.01] + u_0 = [s_0, e_0, i_0, 0.0, 0.0, 0.0, p.R̄, 0.01] tspan = (0.0, p.T) prob = ODEProblem(F, u_0, tspan, p) @@ -267,12 +283,12 @@ Let's run some experiments using this code. First, we can solve the ODE using an appropriate algorthm (e.g. a good default for non-stiff ODEs might be ``Tsit5()``, which is the Tsitouras 5/4 Runge-Kutta method). -Most high-performance ODE solvers appropriate for this class of problems will have adaptive time-stepping, so you +Most accurate and high-performance ODE solvers appropriate for this class of problems will have adaptive time-stepping, so you will not specify any sort of grid .. code-block:: julia - sol = solve(prob, Tsit5()) # TODO: change the accuracy? + sol = solve(prob, Tsit5()) @show length(sol.t); We see that the adaptive time-stepping used approximately 45 time-steps to solve this problem to the desires accuracy. Evaluating the solver at points outside of those time-steps uses the an interpolator consistent with the @@ -287,13 +303,14 @@ See `here `__ for details on anal Experiment 1: Constant Reproduction Case ------------------------------- +---------------------------------------- -Let's start with the case where :math:`B(t) = R = R̄` is constant. +Let's start with the case where :math:`B(t) = R = \bar{R}` is constant. We calculate the time path of infected people under different assumptions. .. code-block:: julia + R̄_vals = range(1.6, 3.0, length = 6) sols = [solve(ODEProblem(F, u_0, tspan, p_gen(R̄ = R̄_val)), Tsit5()) for R̄ in R̄_vals] @@ -342,9 +359,9 @@ Experiment 2: Changing Mitigation --------------------------------- Let's look at a scenario where mitigation (e.g., social distancing) is -successively imposed, but the target (:math:`R̄` is fixed) +successively imposed, but the target (:math:`\bar{R}` is fixed) -To do this, we will have :math:`R(0) \neq R̄` and examine the dynamics using the :math:`\frac{d R}{d t} &= \eta (R̄ - R)` ODE. +To do this, we will have :math:`R(0) \neq \bar{R}` and examine the dynamics using the :math:`\frac{d R}{d t} &= \eta (\bar{R} - R)` ODE. .. Mathematica Verification .. (\[Beta][t] /. @@ -354,7 +371,7 @@ To do this, we will have :math:`R(0) \neq R̄` and examine the dynamics using th .. E^(-t \[Eta])) b // FullSimplify -Note that in the simple case, where :math:`B(t) = \bar{R}` is independent of the state, the solution to the ODE with :math:`R(0) = R_0` is :math:`R(t) = R_0 e^{-\eta t} + \bar{R}(1 - e^{-\eta t})` +In the simple case, where :math:`B(t) = \bar{R}` is independent of the state, the solution to the ODE with :math:`R(0) = R_0` is :math:`R(t) = R_0 e^{-\eta t} + \bar{R}(1 - e^{-\eta t})` We will examine the case where :math:`R(0) = 3` and then it falls to to :math:`\bar{R} = 1.6` due to the progressive adoption of stricter mitigation measures. @@ -415,7 +432,7 @@ and 75,000 agents already exposed to the virus and thus soon to be contagious. s_0 = 1.0 - i_0 - e_0 u_0 = [s_0, e_0, i_0, 0.0, 0.0, 0.5, 0.01] # starting in lockdown, R=0.5 - tspan = [0.0, 30.0, 120.0, p.T] # add discontinuities to tspan for adaptive solver + tstops = [0.0, 30.0, 120.0, p.T] # ensure discontinuites used with adaptive timesteps # create two problems, with rapid movement of R towards B(t) prob_early = ODEProblem(F, u_0, tspan, p_gen(B = B_lift_early, η = 10.0)) @@ -438,7 +455,7 @@ Let's calculate the paths: # c_paths.append(c_path) -Here is the number of active infections: +Here is the number of active infections and mortality with the :math:`v(t) = 0.01` baseline. .. code-block:: julia @@ -446,11 +463,6 @@ Here is the number of active infections: What kind of mortality can we expect under these scenarios? -Suppose that 1\% of cases result in death - -.. code-block:: julia - - #ν = 0.01 This is the cumulative number of deaths: @@ -469,3 +481,65 @@ This is the daily death rate: Pushing the peak of curve further into the future may reduce cumulative deaths if a vaccine is found. + + +SIER with Aggregate Shocks +=========================== + +The model above is fully deterministic: given a steady state and exogenous :math:`B(t)` and :math:`v(t)`. We will extend our model to include a particular set of shocks which make the system a Stochastic Differential Equation (SDE). + +One source of randomness which would enter the model is considreing the discretness of individuals, and modeling a Markov Jump-process. These topics include such as jump-processes and the connection between SDEs and Langevin equations typically used in the approximation of chemical reactions in well-mixed. + +Instead, here we will concentrate on randomness that comes from aggregate changes in behavior or policy. + +Shocks to Transmission Rates +------------------------------ + +First, we consider that the effective transmission rate :math:`R(t)` will depend on degrees of randomness in behavior and implementation. For example, + +* Misinformation on facebook spreading non-uniformly +* Large political rallies or protests +* Deviations in the implementation and timing of lockdown policy within demographics, locations, or businesses within the system. + +To implement this, we will add on a diffusion term to :eq:`Rode` with an instantaneous volatility of :math:`\zeta \sqrt{R}`. The scaling by R ensures that the process can never go negative since the variance converges to zero as R goes to zero. + +The notation for this `SDE `__ is then + +.. math:: + \begin{aligned} + d R&= \eta (B - R) dt + \zeta \sqrt{R} dW + \end{aligned} + :label: Rsde + +where :math:`W` is standard Brownian motion (i.e a `Weiner Process `__. This equation is used in the `Cox-Ingersoll-Ross `__ and `Heston `__ models of interest rates and stochastic volatility. + +Heuristically, if :math:`\zeta = 0`, we can divide by :math:`dt` and nest the original ODE in :eq:`Rode`. + +Migration and Transporation +---------------------------- + +A second source of shocks are associated with policies where new individuals in the Exposed state enter the geography We will maintain a constant population size and assume (without specifying) compensating outflows to match the others, and assume that Infected individuals are effectively barred from entry. + +As it is the main consideration, lets add the diffusive term to the :math:`de` dynamics, + + +.. math:: + \begin{aligned} + d e & = \left(\gamma \, R \, s \, i - \sigma e\right)dt + \theta \sqrt{e} d W + \end{aligned} + :label: esde + +With these, we can define a variance term, mostly zeros since we only have two independent shocks, we can combined them in diagonal noise term :math:`\Omega(u, t)`. Extending + +.. math:: + \begin{aligned} + diag(\Omega) &:= \begin{bmatrix} 0 & \theta \sqrt{e} & 0 & 0 & 0 & 0 & \zeta \sqrt{r} & 0 + \end{alignd} + +Finally, we can extend our existing :ref:`dfcv` definition to in + + + \begin{aligned} + d u &= F(u,t;p)dt + \Omega dW + \end{aligned} + :label: dfcvsde \ No newline at end of file From e77a8945384d98fc322de525de6ed430b079e365 Mon Sep 17 00:00:00 2001 From: Jesse Perla Date: Thu, 25 Jun 2020 00:14:44 -0700 Subject: [PATCH 06/16] Added simple example as well as sketch on stochastics. --- .../tools_and_techniques/seir_model_sde.rst | 140 +++++++++++++----- 1 file changed, 104 insertions(+), 36 deletions(-) diff --git a/source/rst/tools_and_techniques/seir_model_sde.rst b/source/rst/tools_and_techniques/seir_model_sde.rst index 55314161..638fc5b8 100644 --- a/source/rst/tools_and_techniques/seir_model_sde.rst +++ b/source/rst/tools_and_techniques/seir_model_sde.rst @@ -87,6 +87,7 @@ The interest is primarily in * the number of infections at a given time (which determines whether or not the health care system is overwhelmed) and * how long the caseload can be deferred (hopefully until a vaccine arrives) +* how to model aggregate shocks related to policy, behavior, and medical innovations @@ -105,15 +106,10 @@ The transitions between those states are governed by the following rates In addition, the transmission rate can be interpreted as: :math:`R(t) := \beta(t) / \gamma` where :math:`R(t)` is the *effective reproduction number* at time :math:`t`. For this reason, we will work with the :math:`R(t)` reparameterization. -(The notation is standard in the epidemiology literature - though slightly confusing, since :math:`R(t)` is different to -:math:`R`, the symbol that represents the removed state. Throughout the rest of the lecture, we will always use :math:`R` to represent the effective reproduction number, unless stated otherwise) +The notation is standard in the epidemiology literature - though slightly confusing, since :math:`R(t)` is different to +:math:`R`, the symbol that represents the removed state. Throughout the rest of the lecture, we will always use :math:`R` to represent the effective reproduction number, unless stated otherwise. - -Assume that there is a constant population of size :math:`N` throughout, then define the proportion of people in each state as :math:`s := S/N` etc. With this, the SEIR model can be written as - - -Since we are using a continuum approximation, all individuals in the population are eventually infected when -the transmission rate is positive and :math:`i(0) > 0`. +Assume that there is a constant population of size :math:`N` throughout, then define the proportion of people in each state as :math:`s := S/N` etc. With this, and assuming a continuum approximation, the SEIR model can be written as .. math:: \begin{aligned} @@ -132,13 +128,48 @@ Here, :math:`dy/dt` represents the time derivative for the particular variable. Since the states form a partition, we could reconstruct the "removed" fraction of the population as :math:`r = 1 - s - e - i`. However, for further experiments and plotting it is harmless to keep it in the system. -In addition, we are interested in calculating the cumulative caseload (i.e., all those who have or have had the infection) as :math:`c = i + r`. Differentiating that expression and substituing from the time-derivatives of :math:`i(t), r(t)` yields :math:`\frac{d c}{d t} = \sigma e` +Since we are using a continuum approximation, all individuals in the population are eventually infected when +the transmission rate is positive and :math:`i(0) > 0`. + +We can begin writing the minimal code to solve the dynamics from a particular ``x_0 = [s_0, e_0, i_0]`` initial condition and parameter values. First, by expressing the system +, + +.. code-block:: julia + + function F(x, p, t; γ = 1/18, R = 3.0, σ = 1/5.2) + s, e, i, r = x + + return [-γ*R*s*i; # ds/dt = -γRsi + γ*R*s*i - σ*e; # de/dt = γRsi -σe + σ*e - γ*i; # di/dt = σe -γi + γ*i; # dr/dt = γi + ] + end + i_0 = 1e-7 + e_0 = 4.0 * i_0 + s_0 = 1.0 - i_0 - e_0 + r_0 = 0.0 + x_0 = [s_0, e_0, i_0, r_0] # initial condition + + tspan = (0.0, 350.0) # ≈ 18 months + prob = ODEProblem(F, x_0, tspan) # create problem + sol = solve(prob, Tsit5()) # solve the model with a particular algorithm + plot(sol, labels = ["s" "e" "i" "r"], title = "SEIR Dynamics", lw = 2) + +Or as an alternative visualization, See the proportions over time + +.. code-block:: julia + + areaplot(sol', labels = ["s" "e" "i" "r"], title = "SIER Proportions") + -Implementing the system of ODEs in :math:`s, e, i` would be enough to implement the model given fixed parameters, but we will extend the basic model to enable some policy experiments. +While implementing the system of ODEs in :math:`s, e, i` is enough to understand basic dynamics, , but we will extend the basic model to enable some policy experiments. -Evolution of Parameters +Extending the Model ----------------------- +First, we can consider some additional calculations such as the cumulative caseload (i.e., all those who have or have had the infection) as :math:`c = i + r`. Differentiating that expression and substituing from the time-derivatives of :math:`i(t), r(t)` yields :math:`\frac{d c}{d t} = \sigma e` + We will assume that the transmission rate follows a process with a reversion to a value :math:`B(t)` which could conceivably be a policy parameter. The intuition is that even if the targetted :math:`B(t)` was changed, lags in behavior and implementation would smooth out the transition, where :math:`\eta` governs the speed of :math:`R(t)` moves towards :math:`B(t)`. .. math:: @@ -147,21 +178,21 @@ We will assume that the transmission rate follows a process with a reversion to \end{aligned} :label: Rode -Finally, let :math:`v(t)` be the mortality rate, which we will leave constant for now, i.e. :math:`\frac{d v}{d t} = 0`. The cumulative deaths can be integrated through the flow :math:`\gamma i` entering the "Removed" state and define the cumulative number of deaths as :math:`m(t)`. The differential equations then +Finally, let :math:`m(t)` be the mortality rate, which we will leave constant for now, i.e. :math:`\frac{d m}{d t} = 0`. The cumulative deaths can be integrated through the flow :math:`\gamma i` entering the "Removed" state and define the cumulative number of deaths as :math:`D(t)`. The differential equations then follow, .. math:: \begin{aligned}\\ - \frac{d v}{d t} &= 0\\ - \frac{d m}{d t} &= v \gamma i + \frac{d m}{d t} &= 0\\ + \frac{d D}{d t} &= m \gamma i \end{aligned} -While we could conveivably integate the total deaths given the solution to the model, it is more convenient to use the integrator built into the ODE solver. That is, we added :math:`d m(t)/dt` rather than calculating :math:`m(t) = \int_0^t \gamma v(\tau) i(\tau) d \tau` after generating the full :math:`i(t)` path. +While we could conveivably integate the total deaths given the solution to the model, it is more convenient to use the integrator built into the ODE solver. That is, we added :math:`d D(t)/dt` rather than calculating :math:`D(t) = \int_0^t \gamma m(\tau) i(\tau) d \tau` after generating the full :math:`i(t)` path. This is a common trick when solving systems of ODEs. While equivalent in principle if you used an appropriate quadrature scheme, this trick becomes especially important and convenient when adaptive time-stepping algorithms are used to solve the ODEs (i.e. there is no fixed time grid). -The system :eq:`seir_system` and the supplemental equations can be written in vector form :math:`x := [s, e, i, r, c, m, R, v]` with parameter tuple :math:`p := (\sigma, \gamma, B, \eta)` +The system :eq:`seir_system` and the supplemental equations can be written in vector form :math:`x := [s, e, i, r, c, D, R, m]` with parameter tuple :math:`p := (\sigma, \gamma, B, \eta)` .. math:: \begin{aligned} @@ -220,9 +251,9 @@ Now we construct a function that represents :math:`F` in :eq:`dfcv` .. code-block:: julia # Reminder: could solve dynamics of SEIR states with just first 3 equations - function F(u, p, t) + function F(x, p, t) - s, e, i, r, c, m, R, v = u + s, e, i, r, c, B, R, m = x @unpack σ, γ, B, η = p return [-γ*R*s*i; # ds/dt = -γRsi @@ -230,9 +261,9 @@ Now we construct a function that represents :math:`F` in :eq:`dfcv` σ*e - γ*i; # di/dt = σe -γi γ*i; # dr/dt = γi σ*e; # dc/dt = σe - v*γ*i; # dm/dt = vγi + v*γ*i; # dD/dt = vγi η*(B(t, p) - R);# dβ/dt = η(B(t) - R) - 0.0 # dv/dt = 0 + 0.0 # dm/dt = 0 ] end @@ -267,9 +298,9 @@ Setting initial conditions, we will assume a fixed :math:`i, e`, :math:`r=0`, :m e_0 = 4.0 * i_0 s_0 = 1.0 - i_0 - e_0 - u_0 = [s_0, e_0, i_0, 0.0, 0.0, 0.0, p.R̄, 0.01] + x_0 = [s_0, e_0, i_0, 0.0, 0.0, 0.0, p.R̄, 0.01] tspan = (0.0, p.T) - prob = ODEProblem(F, u_0, tspan, p) + prob = ODEProblem(F, x_0, tspan, p) The ``tspan`` determines that the :math:`t` used by the sovler, where the scale needs to be consistent with the arrival @@ -312,7 +343,7 @@ We calculate the time path of infected people under different assumptions. .. code-block:: julia R̄_vals = range(1.6, 3.0, length = 6) - sols = [solve(ODEProblem(F, u_0, tspan, p_gen(R̄ = R̄_val)), Tsit5()) for R̄ in R̄_vals] + sols = [solve(ODEProblem(F, x_0, tspan, p_gen(R̄ = R̄_val)), Tsit5()) for R̄ in R̄_vals] # TODO: Probably clean ways to plot this @@ -431,12 +462,12 @@ and 75,000 agents already exposed to the virus and thus soon to be contagious. e_0 = 75000 / N s_0 = 1.0 - i_0 - e_0 - u_0 = [s_0, e_0, i_0, 0.0, 0.0, 0.5, 0.01] # starting in lockdown, R=0.5 + x_0 = [s_0, e_0, i_0, 0.0, 0.0, 0.5, 0.01] # starting in lockdown, R=0.5 tstops = [0.0, 30.0, 120.0, p.T] # ensure discontinuites used with adaptive timesteps # create two problems, with rapid movement of R towards B(t) - prob_early = ODEProblem(F, u_0, tspan, p_gen(B = B_lift_early, η = 10.0)) - prob_late = ODEProblem(F, u_0, tspan, p_gen(B = B_lift_late, η = 10.0) + prob_early = ODEProblem(F, x_0, tspan, p_gen(B = B_lift_early, η = 10.0)) + prob_late = ODEProblem(F, x_0, tspan, p_gen(B = B_lift_late, η = 10.0) Let's calculate the paths: @@ -455,7 +486,7 @@ Let's calculate the paths: # c_paths.append(c_path) -Here is the number of active infections and mortality with the :math:`v(t) = 0.01` baseline. +Here is the number of active infections and mortality with the :math:`m(t) = 0.01` baseline. .. code-block:: julia @@ -486,7 +517,7 @@ if a vaccine is found. SIER with Aggregate Shocks =========================== -The model above is fully deterministic: given a steady state and exogenous :math:`B(t)` and :math:`v(t)`. We will extend our model to include a particular set of shocks which make the system a Stochastic Differential Equation (SDE). +The model above is fully deterministic: given a steady state and exogenous :math:`B(t)` and :math:`m(t)`. We will extend our model to include a particular set of shocks which make the system a Stochastic Differential Equation (SDE). One source of randomness which would enter the model is considreing the discretness of individuals, and modeling a Markov Jump-process. These topics include such as jump-processes and the connection between SDEs and Langevin equations typically used in the approximation of chemical reactions in well-mixed. @@ -507,13 +538,13 @@ The notation for this `SDE `__. This equation is used in the `Cox-Ingersoll-Ross `__ and `Heston `__ models of interest rates and stochastic volatility. -Heuristically, if :math:`\zeta = 0`, we can divide by :math:`dt` and nest the original ODE in :eq:`Rode`. +Heuristically, if :math:`\zeta_R = 0`, we can divide by :math:`dt` and nest the original ODE in :eq:`Rode`. Migration and Transporation ---------------------------- @@ -525,21 +556,58 @@ As it is the main consideration, lets add the diffusive term to the :math:`de` d .. math:: \begin{aligned} - d e & = \left(\gamma \, R \, s \, i - \sigma e\right)dt + \theta \sqrt{e} d W + d e & = \left(\gamma \, R \, s \, i - \sigma e\right)dt + \zeta_e \sqrt{e} d W \end{aligned} :label: esde -With these, we can define a variance term, mostly zeros since we only have two independent shocks, we can combined them in diagonal noise term :math:`\Omega(u, t)`. Extending + +Mortality Fluctuations +----------------------------------- + +There may be a variety of medical interventions that are short of a vaccine, but still effect the :math:`m(t)` path. In addition, imperfect mixing of different demographic groups could lead to aggregate shocks in mortality (e.g. if a retirement home is afflicted vs. an elementary school) + +We will begin by adding in sorts of random shocks, and leaving out dift or any mean-reversion for simplicity + +.. math:: + \begin{aligned} + d m & = \zeta_m \sqrt{m} d W + \end{aligned} + + + + +Combining the Shocks +--------------------- + +With these, we can define a variance term, mostly zeros since we only have two independent shocks, we can combined them in diagonal noise term :math:`\Omega(x, t)`. Extending .. math:: \begin{aligned} - diag(\Omega) &:= \begin{bmatrix} 0 & \theta \sqrt{e} & 0 & 0 & 0 & 0 & \zeta \sqrt{r} & 0 + diag(\Omega) &:= \begin{bmatrix} 0\\ + \zeta_e \sqrt{e}\\ + 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ \zeta_R \sqrt{r} \\ + & \zeta_m \sqrt{m} \end{alignd} Finally, we can extend our existing :ref:`dfcv` definition to in - +.. math:: \begin{aligned} - d u &= F(u,t;p)dt + \Omega dW + d x &= F(x,t;p)dt + \Omega dW \end{aligned} - :label: dfcvsde \ No newline at end of file + +TODO: IMPLEMENT THIS WITH SDE + +Jump Processes +================== + +We will extend the above with 2 variations: +1. Assume that major medical advancements can arrival, which drop mortality, :math:`m(t)` in half. As there remains a diffusion term, the resulting :math:`d x`, becomes a Jump diffusionAlso vaccination +.. math:: + \begin{aligned} + d m & = \zeta_m \sqrt{m_t} d W + \theta N_t + \end{aligned} + +2. Vaccines arrival, :math:`v(t)` which being the process of directly enabling a "suceptible" to "removed jump". More vaccine arrivals speed up the process. + +Will consider that vaccine arrival rates are time varying. :math:`\alpha(t)`.... As there is no diffusion, this is called a Pure jump process. From 1a2543e174f060e338b75a1fba553c8bf9a784a3 Mon Sep 17 00:00:00 2001 From: Jesse Perla Date: Thu, 25 Jun 2020 16:28:24 -0700 Subject: [PATCH 07/16] Up to stochastics... --- .../tools_and_techniques/seir_model_sde.rst | 301 ++++++++++-------- 1 file changed, 169 insertions(+), 132 deletions(-) diff --git a/source/rst/tools_and_techniques/seir_model_sde.rst b/source/rst/tools_and_techniques/seir_model_sde.rst index 638fc5b8..c35c6a13 100644 --- a/source/rst/tools_and_techniques/seir_model_sde.rst +++ b/source/rst/tools_and_techniques/seir_model_sde.rst @@ -22,13 +22,14 @@ equations. The main objective is to study the impact of suppression through social distancing on the spread of the infection. -Dynamics are modeled using a standard SEIR (Susceptible-Exposed-Infected-Removed) model -of disease spread, represented as a system of ordinary differential -equations when the number of agents is large and there are no exogenous stochastic shocks. - The focus is on US outcomes but the parameters can be adjusted to study other countries. +In the first part, dynamics are modeled using a standard SEIR (Susceptible-Exposed-Infected-Removed) model +of disease spread, represented as a system of ordinary differential +equations where the number of agents is large and there are no exogenous stochastic shocks. + + The first part of the model follows the notes from provided by `Andrew Atkeson `__ @@ -39,6 +40,11 @@ provided by `Andrew Atkeson `__ See further variations on the classic SIR model in Julia `here `__. +We then look at extending the model to include policy-relevant aggregate shocks, and +examine the three main techniques for including stochasticity to continuous-time models: +* Brownian Motion: A diffusion process with stochastic, continous paths. The prototypical Stochastic Differential Equation (SDE) with additive noise. +* Pure-Jump Processes: A variable that jumps between a discrete number of values, typically with a Poisson arrival rate. +* Jump-Diffusion Process: A stochastic process that contains both a diffusion term and arrival rates of discrete jumps. Setup ------------------ @@ -104,7 +110,7 @@ The transitions between those states are governed by the following rates * :math:`\gamma` is called the *recovery rate* (the rate at which infected people recover or die). -In addition, the transmission rate can be interpreted as: :math:`R(t) := \beta(t) / \gamma` where :math:`R(t)` is the *effective reproduction number* at time :math:`t`. For this reason, we will work with the :math:`R(t)` reparameterization. +In addition, the transmission rate can be rep-parameterized as: :math:`R(t) := \beta(t) / \gamma` where :math:`R(t)` has the interpretation as the *effective reproduction number* at time :math:`t`. For this reason, we will work with the :math:`R(t)` reparameterization. The notation is standard in the epidemiology literature - though slightly confusing, since :math:`R(t)` is different to :math:`R`, the symbol that represents the removed state. Throughout the rest of the lecture, we will always use :math:`R` to represent the effective reproduction number, unless stated otherwise. @@ -131,46 +137,68 @@ Since the states form a partition, we could reconstruct the "removed" fraction o Since we are using a continuum approximation, all individuals in the population are eventually infected when the transmission rate is positive and :math:`i(0) > 0`. -We can begin writing the minimal code to solve the dynamics from a particular ``x_0 = [s_0, e_0, i_0]`` initial condition and parameter values. First, by expressing the system -, +We can begin writing the minimal code to solve the dynamics from a particular ``x_0 = [s_0, e_0, i_0, r_0]`` initial condition and parameter values. + +First, define the system of equations .. code-block:: julia - function F(x, p, t; γ = 1/18, R = 3.0, σ = 1/5.2) + function F_simple(x, p, t; γ = 1/18, R = 3.0, σ = 1/5.2) s, e, i, r = x return [-γ*R*s*i; # ds/dt = -γRsi - γ*R*s*i - σ*e; # de/dt = γRsi -σe - σ*e - γ*i; # di/dt = σe -γi + γ*R*s*i - σ*e;# de/dt = γRsi -σe + σ*e - γ*i; # di/dt = σe -γi γ*i; # dr/dt = γi ] end - i_0 = 1e-7 + +Written this way, we see that the four equations represent the one-directional transition from the susceptible to removed state, where the negative terms are outflows, and the positive ones inflows. + +As there is no flow leaving the :math:`dr/dt` and all parameters are positive, unless we start with a degenerate initial condition (e.g. :math:`e(0) = i(0) = 0`) the "Removed" state is asymptoically absorbing, and :math:`\lim_{t\to \infty} r(t) = 1`. Crucial to this result is that individuals are perfectly divisible, and any arbitrarily small :math:`i > 0` leads to a strictly positive flow into the exposed state. + +We will discuss this topic further in the lecture on continuous-time +markov-chains, as well as the limitations of these approximations when the discretness becomes essential (e.g. continuum approximations are incapable of modeling extinguishing of an outbreak). + +Given this system, we choose an initial condition and a timespan, and create a ``ODEProblem`` encapsulating the system. + +.. code-block:: julia + + i_0 = 1E-7 e_0 = 4.0 * i_0 s_0 = 1.0 - i_0 - e_0 r_0 = 0.0 x_0 = [s_0, e_0, i_0, r_0] # initial condition - tspan = (0.0, 350.0) # ≈ 18 months - prob = ODEProblem(F, x_0, tspan) # create problem - sol = solve(prob, Tsit5()) # solve the model with a particular algorithm + tspan = (0.0, 350.0) # ≈ 350 days + prob = ODEProblem(F_simple, x_0, tspan) + +With this, we can choose an ODE algorithm (e.g. a good default for non-stiff ODEs of this sort might be ``Tsit5()``, which is the Tsitouras 5/4 Runge-Kutta method). + +.. code-block:: julia + + sol = solve(prob, Tsit5()) plot(sol, labels = ["s" "e" "i" "r"], title = "SEIR Dynamics", lw = 2) - -Or as an alternative visualization, See the proportions over time + + +We did not provide either a set of timesteps or a ``dt`` time stepsize to the ``solve``. The reason is that most accurate and high-performance ODE solvers appropriate use adaptive time-stepping, changing the stepsize based the degree of curvature in the derivatives. + + +Or, as an alternative visualization, the proportions in each state over time .. code-block:: julia - areaplot(sol', labels = ["s" "e" "i" "r"], title = "SIER Proportions") + areaplot(sol.t, sol', labels = ["s" "e" "i" "r"], title = "SIER Proportions") -While implementing the system of ODEs in :math:`s, e, i` is enough to understand basic dynamics, , but we will extend the basic model to enable some policy experiments. +While implementing the system of ODEs in :math:`(s, e, i)`, we will extend the basic model to enable some policy experiments and calculations of aggregate values. Extending the Model ----------------------- First, we can consider some additional calculations such as the cumulative caseload (i.e., all those who have or have had the infection) as :math:`c = i + r`. Differentiating that expression and substituing from the time-derivatives of :math:`i(t), r(t)` yields :math:`\frac{d c}{d t} = \sigma e` -We will assume that the transmission rate follows a process with a reversion to a value :math:`B(t)` which could conceivably be a policy parameter. The intuition is that even if the targetted :math:`B(t)` was changed, lags in behavior and implementation would smooth out the transition, where :math:`\eta` governs the speed of :math:`R(t)` moves towards :math:`B(t)`. +We will assume that the transmission rate follows a process with a reversion to a value :math:`B(t)` which could conceivably be influenced by policy. The intuition is that even if the targetted :math:`B(t)` was changed through social distancing/etc., lags in behavior and implementation would smooth out the transition, where :math:`\eta` governs the speed of :math:`R(t)` moves towards :math:`B(t)`. .. math:: \begin{aligned} @@ -178,21 +206,20 @@ We will assume that the transmission rate follows a process with a reversion to \end{aligned} :label: Rode -Finally, let :math:`m(t)` be the mortality rate, which we will leave constant for now, i.e. :math:`\frac{d m}{d t} = 0`. The cumulative deaths can be integrated through the flow :math:`\gamma i` entering the "Removed" state and define the cumulative number of deaths as :math:`D(t)`. The differential equations then -follow, +Finally, let :math:`m(t)` be the mortality rate, which we will leave constant for now, i.e. :math:`\frac{d m}{d t} = 0`. The cumulative deaths can be integrated through the flow :math:`\gamma i` entering the "Removed" state and define the cumulative number of deaths as :math:`M(t)`. .. math:: \begin{aligned}\\ \frac{d m}{d t} &= 0\\ - \frac{d D}{d t} &= m \gamma i + \frac{d M}{d t} &= m \gamma i \end{aligned} -While we could conveivably integate the total deaths given the solution to the model, it is more convenient to use the integrator built into the ODE solver. That is, we added :math:`d D(t)/dt` rather than calculating :math:`D(t) = \int_0^t \gamma m(\tau) i(\tau) d \tau` after generating the full :math:`i(t)` path. +While we could conveivably integate the total deaths given the solution to the model, it is more convenient to use the integrator built into the ODE solver. That is, we add :math:`d M(t)/dt` rather than calculating :math:`M(t) = \int_0^t \gamma m(\tau) i(\tau) d \tau` ex-post. -This is a common trick when solving systems of ODEs. While equivalent in principle if you used an appropriate quadrature scheme, this trick becomes especially important and convenient when adaptive time-stepping algorithms are used to solve the ODEs (i.e. there is no fixed time grid). +This is a common trick when solving systems of ODEs. While equivalent in principle to using the appropriate quadrature scheme, this becomes especially important and convenient when adaptive time-stepping algorithms are used to solve the ODEs (i.e. there is no fixed time grid). -The system :eq:`seir_system` and the supplemental equations can be written in vector form :math:`x := [s, e, i, r, c, D, R, m]` with parameter tuple :math:`p := (\sigma, \gamma, B, \eta)` +The system :eq:`seir_system` and the supplemental equations can be written in vector form :math:`x := [s, e, i, r, R, m, c, M]` with parameter tuple :math:`p := (\sigma, \gamma, B, \eta)` .. math:: \begin{aligned} @@ -205,14 +232,14 @@ The system :eq:`seir_system` and the supplemental equations can be written in ve \sigma \, e - \gamma i \\ \gamma i - \\ - \sigma e - \\ - v \, \gamma \, i \\ \eta (B(t) - R) \\ 0 + \\ + \sigma e + \\ + m \, \gamma \, i \end{bmatrix} \end{aligned} :label: dfcv @@ -229,7 +256,7 @@ As in Atkeson's note, we set * :math:`\sigma = 1/5.2` to reflect an average incubation period of 5.2 days. * :math:`\gamma = 1/18` to match an average illness duration of 18 days. * :math:`B = R = 1.6` to match an **effective reproduction rate** of 1.6, and initially time-invariant -* :math:`v = 0.01` for a one-percent mortality rate +* :math:`m(t) = m_0 = 0.01` for a one-percent mortality rate As we will initially consider the case where :math:`R(0) = B`, the value of :math:`\eta` will drop out of this first experiment. @@ -237,44 +264,36 @@ As we will initially consider the case where :math:`R(0) = B`, the value of :mat Implementation ============== -First we set the population size to match the US and the parameters as described +# First we set the population size to match the US and the parameters as described -.. code-block:: julia +# .. code-block:: julia - N = 3.3e8 # US Population - γ = 1 / 18 - σ = 1 / 5.2 - η = 1 / 20 # a placeholder, drops out of firs texperiments. +# N = 3.3E8 # US Population +# γ = 1 / 18 +# σ = 1 / 5.2 +# η = 1 / 20 # a placeholder, inactive in first experiments -Now we construct a function that represents :math:`F` in :eq:`dfcv` +First, construct our :math:`F` from :eq:`dfcv` .. code-block:: julia - # Reminder: could solve dynamics of SEIR states with just first 3 equations function F(x, p, t) - s, e, i, r, c, B, R, m = x + s, e, i, r, R, m, c, M = x @unpack σ, γ, B, η = p - return [-γ*R*s*i; # ds/dt = -γRsi - γ*R*s*i - σ*e; # de/dt = γRsi - σe - σ*e - γ*i; # di/dt = σe -γi - γ*i; # dr/dt = γi - σ*e; # dc/dt = σe - v*γ*i; # dD/dt = vγi - η*(B(t, p) - R);# dβ/dt = η(B(t) - R) - 0.0 # dm/dt = 0 + return [-γ*R*s*i; # ds/dt + γ*R*s*i - σ*e; # de/dt + σ*e - γ*i; # di/dt + γ*i; # dr/dt + η*(B(t, p) - R);# dR/dt + 0.0; # dm/dt + σ*e; # dc/dt + m*γ*i; # dM/dt ] end -Written this way, we see that the first four rows represent the one-directional transition from the susceptible to removed state, where the negative terms are outflows, and the positive ones inflows. - -As there is no flow leaving the :math:`dr/dt` and all parameters are positive, unless we start with a degenerate initial condition (e.g. :math:`e(0) = i(0) = 0`) the "Removed" state is asymptoically absorbing, and :math:`\lim_{t\to \infty} r(t) = 1`. Crucial to this rsult is that individuals are arbitrarily divisible, and any arbitrarily small :math:`i > 0` leads to a strictly positive flow into the exposed state. - -We will discuss this topic further in the lecture on continuous-time -markov-chains, as well as the limitations of these approximations when the discretness becomes essential - Parameters ------------- @@ -282,8 +301,8 @@ The baseline parameters are put into a named tuple generator (see previous lectu .. code-block:: julia - p_gen = @with_kw (T = 550.0, γ = 1.0 / 18, σ = 1 / 5.2, η = 1.0 / 20, - R̄ = 1.6, B = (t, p) -> p.R̄) + p_gen = @with_kw ( T = 550.0, γ = 1.0 / 18, σ = 1 / 5.2, η = 1.0 / 20, + R̄ = 1.6, B = (t, p) -> p.R̄, m_0 = 0.01) Note that the default :math:`B(t)` function always equals :math:`\bar{R}` @@ -294,13 +313,12 @@ Setting initial conditions, we will assume a fixed :math:`i, e`, :math:`r=0`, :m p = p_gen() # use all default parameters - i_0 = 1e-7 + i_0 = 1E-7 e_0 = 4.0 * i_0 s_0 = 1.0 - i_0 - e_0 - x_0 = [s_0, e_0, i_0, 0.0, 0.0, 0.0, p.R̄, 0.01] - tspan = (0.0, p.T) - prob = ODEProblem(F, x_0, tspan, p) + x_0 = [s_0, e_0, i_0, 0.0, p.R̄, p.m_0, 0.0, 0.0] + prob = ODEProblem(F, x_0, (0.0, p.T), p) The ``tspan`` determines that the :math:`t` used by the sovler, where the scale needs to be consistent with the arrival @@ -312,25 +330,20 @@ Experiments Let's run some experiments using this code. -First, we can solve the ODE using an appropriate algorthm (e.g. a good default for non-stiff ODEs might be ``Tsit5()``, which is the Tsitouras 5/4 Runge-Kutta method). - -Most accurate and high-performance ODE solvers appropriate for this class of problems will have adaptive time-stepping, so you -will not specify any sort of grid - .. code-block:: julia sol = solve(prob, Tsit5()) @show length(sol.t); -We see that the adaptive time-stepping used approximately 45 time-steps to solve this problem to the desires accuracy. Evaluating the solver at points outside of those time-steps uses the an interpolator consistent with the +We see that the adaptive time-stepping used approximately 50 time-steps to solve this problem to the desires accuracy. Evaluating the solver at points outside of those time-steps uses the an interpolator consistent with the solution to the ODE. See `here `__ for details on analyzing the solution, and `here `__ for plotting tools. The built-in plots for the solutions provide all of the `attributes `__. .. code_block:: julia - # TODO: Chris, We could plot something else as well? Deaths, etc?? Also labels broken - plot(sol, vars = [3, 5], label = ["i(t)" "c(t)"], lw = 2, title = "Current and Cumulated Infected") + # TODO: Chris, nice ways to resvale things or use two axis? + plot(sol, vars = [7, 8], label = ["c(t)" "M(t)"], lw = 2, title = "Cumulative Infected and Total Mortality") Experiment 1: Constant Reproduction Case @@ -343,9 +356,9 @@ We calculate the time path of infected people under different assumptions. .. code-block:: julia R̄_vals = range(1.6, 3.0, length = 6) - sols = [solve(ODEProblem(F, x_0, tspan, p_gen(R̄ = R̄_val)), Tsit5()) for R̄ in R̄_vals] + sols = [solve(ODEProblem(F, x_0, tspan, p_gen(R̄ = R̄)), Tsit5()) for R̄ in R̄_vals] - # TODO: Probably clean ways to plot this + # TODO: Probably clean ways to plot this , but don't know them! #R0_vals = np.linspace(1.6, 3.0, 6) #labels = [f'$R0 = {r:.2f}$' for r in R0_vals] @@ -455,114 +468,143 @@ and 75,000 agents already exposed to the virus and thus soon to be contagious. B_lift_early(t, p) = t < 30.0 ? 0.5 : 2.0 B_lift_late(t, p) = t < 120.0 ? 0.5 : 2.0 - # initial conditions + N = 3.3E8 # US Population i_0 = 25000 / N e_0 = 75000 / N s_0 = 1.0 - i_0 - e_0 + R_0 = 0.5 - x_0 = [s_0, e_0, i_0, 0.0, 0.0, 0.5, 0.01] # starting in lockdown, R=0.5 - tstops = [0.0, 30.0, 120.0, p.T] # ensure discontinuites used with adaptive timesteps + x_0 = [s_0, e_0, i_0, 0.0, R_0, p.m_0, 0.0, 0.0] # start in lockdown # create two problems, with rapid movement of R towards B(t) - prob_early = ODEProblem(F, x_0, tspan, p_gen(B = B_lift_early, η = 10.0)) - prob_late = ODEProblem(F, x_0, tspan, p_gen(B = B_lift_late, η = 10.0) + p_early = p_gen(B = B_lift_early, η = 10.0) + p_late = p_gen(B = B_lift_late, η = 10.0) + prob_early = ODEProblem(F, x_0, tspan, p_early) + prob_late = ODEProblem(F, x_0, tspan, p_late) -Let's calculate the paths: - -.. code-block:: julia - - #R0_paths = (lambda t: 0.5 if t < 30 else 2, - # lambda t: 0.5 if t < 120 else 2) - # - #labels = [f'scenario {i}' for i in (1, 2)] - # - #i_paths, c_paths = [], [] - # - #for R0 in R0_paths: - # i_path, c_path = solve_path(R0, t_vec, x_init=x_0) - # i_paths.append(i_path) - # c_paths.append(c_path) +Unlike the previous examples, the :math:`B(t)` functions have discontinuties which might occur. We can tell the adaptive time-stepping methods to ensure they include those points using ``tstops`` -Here is the number of active infections and mortality with the :math:`m(t) = 0.01` baseline. +Let's calculate the paths: .. code-block:: julia - #plot_paths(i_paths, labels) - -What kind of mortality can we expect under these scenarios? + sol_early = solve(prob_early, Tsit5(), tstops = [30.0, 120.0]) + sol_late = solve(prob_late, Tsit5(), tstops = [30.0, 120.0]) + plot(sol_early, vars =[8], title = "Total Mortality", label = "Lift Early") + plot!(sol_late, vars =[8], label = "Lift Late") - -This is the cumulative number of deaths: +To calculate the daily death, calculate the :math:`\gamma i(t) m(t)`. .. code-block:: julia - #paths = [path * ν * pop_size for path in c_paths] - #plot_paths(paths, labels) - -This is the daily death rate: + daily_early = sol_early[3,:] .* sol_early[6,:] * p_early.γ + daily_late = sol_late[3,:] .* sol_late[6,:] * p_late.γ -.. code-block:: julia - - #paths = [path * ν * γ * pop_size for path in i_paths] - #plot_paths(paths, labels) + plot(sol_early.t, daily_early, title = "Flow Deaths", label = "Lift Early") + plot!(sol_late.t, daily_late, label = "Lift Late") -Pushing the peak of curve further into the future may reduce cumulative deaths +Pushing the peak of curve further into the future may reduce cumulative deaths if a vaccine is found. +Despite its richness, the model above is fully deterministic. The policy :math:`B(t)` could change over time, but only in predictable ways. -SIER with Aggregate Shocks -=========================== +One source of randomness which would enter the model is considering the discretness of individuals. This topic, the connection to between SDEs and the Langevin equations typically used in the approximation of chemical reactions in well-mixed media are explored in our lecture on continuous time markov chains. -The model above is fully deterministic: given a steady state and exogenous :math:`B(t)` and :math:`m(t)`. We will extend our model to include a particular set of shocks which make the system a Stochastic Differential Equation (SDE). +But rather than examining how granularity leads to aggregate fluctuations, we will concentrate on randomness that comes from aggregate changes in behavior or policy. -One source of randomness which would enter the model is considreing the discretness of individuals, and modeling a Markov Jump-process. These topics include such as jump-processes and the connection between SDEs and Langevin equations typically used in the approximation of chemical reactions in well-mixed. +Aggregate Shocks to Transmission Rates +======================================= -Instead, here we will concentrate on randomness that comes from aggregate changes in behavior or policy. +We will start by extending our model to include randomness in :math:`R(t)`, which makes it a system of Stochastic Differential Equations (SDEs). Shocks to Transmission Rates ------------------------------ -First, we consider that the effective transmission rate :math:`R(t)` will depend on degrees of randomness in behavior and implementation. For example, +Consider that the effective transmission rate :math:`R(t)` could depend on degrees of randomness in behavior and implementation. For example, * Misinformation on facebook spreading non-uniformly -* Large political rallies or protests +* Large political rallies, elections, or protests * Deviations in the implementation and timing of lockdown policy within demographics, locations, or businesses within the system. +* Aggregate shocks in opening/closing industries -To implement this, we will add on a diffusion term to :eq:`Rode` with an instantaneous volatility of :math:`\zeta \sqrt{R}`. The scaling by R ensures that the process can never go negative since the variance converges to zero as R goes to zero. +To implement this, we will add on a diffusion term to :eq:`Rode` with an instantaneous volatility of :math:`\zeta \sqrt{R}`. The scaling by the :math:`\sqrt{R}` ensures that the process can never go negative since the variance converges to zero as R goes to zero. The notation for this `SDE `__ is then .. math:: \begin{aligned} - d R&= \eta (B - R) dt + \zeta_R \sqrt{R} dW + d R&= \eta (B - R) dt + \zeta \sqrt{R} dW \end{aligned} :label: Rsde where :math:`W` is standard Brownian motion (i.e a `Weiner Process `__. This equation is used in the `Cox-Ingersoll-Ross `__ and `Heston `__ models of interest rates and stochastic volatility. -Heuristically, if :math:`\zeta_R = 0`, we can divide by :math:`dt` and nest the original ODE in :eq:`Rode`. +Heuristically, if :math:`\zeta = 0`, we can divide by :math:`dt` and nest the original ODE in :eq:`Rode`. -Migration and Transporation ----------------------------- +The general form of the SDE with these sorts of continuous-shocks is an extension of our :ref:`dfcv` definition . -A second source of shocks are associated with policies where new individuals in the Exposed state enter the geography We will maintain a constant population size and assume (without specifying) compensating outflows to match the others, and assume that Infected individuals are effectively barred from entry. +.. math:: + \begin{aligned} + d x &= F(x,t;p)dt + \G(x,t;p) dW + \end{aligned} -As it is the main consideration, lets add the diffusive term to the :math:`de` dynamics, +Here, it is convenient to :math:`d W` with the same dimension as :math:`x` where we can use the matrix :math:`G(x,t;p)` to associate the shocks with the appropriate :math:`x`. +As the shock only effects :math:`dR`, which is the 5th equation, define the matrix as .. math:: - \begin{aligned} - d e & = \left(\gamma \, R \, s \, i - \sigma e\right)dt + \zeta_e \sqrt{e} d W + \begin{aligned} + diag(G) &:= \begin{bmatrix} 0 & 0 & 0 & 0 & \zeta \sqrt{R} & 0 & 0 & 0 \end{bmatrix} \end{aligned} - :label: esde -Mortality Fluctuations ------------------------------------ +Since these are additive shocks, we will not need to modify the :math:`F` from our equation. + +First create a new settings generator, and and then define a `SDEProblem`__ with Diagonal Noise. + +We solve the problem with the `SRIW1 `__ algorithm, an Adaptive strong order 1.5 and weak order 2.0 for diagonal/scalar Ito SDEs) + +.. code-block:: julia + + p_sde_gen = @with_kw ( T = 550.0, γ = 1.0 / 18, σ = 1 / 5.2, η = 1.0 / 20, + R̄ = 1.6, B = (t, p) -> p.R̄, m_0 = 0.01, ζ = 1E-3) + + # TODO: CHRIS Easy way to have scalar noise or stick with this for now? + function G(x,p,t) + R = x[5] + return [0, 0, 0, 0, p.ζ * sqrt(R), 0, 0, 0] + end + prob = SDEProblem(F, G, x_0, tspan, p_sde_gen()) + sol = solve(prob, SRIW1()) + +.. Migration and Transporation +.. ---------------------------- +.. +.. A second source of shocks are associated with policies where new individuals in the .. Exposed state enter the geography We will maintain a constant population size and .. assume (without specifying) compensating outflows to match the others, and assume that .. Infected individuals are effectively barred from entry. +.. +.. As it is the main consideration, lets add the diffusive term to the :math:`de` .. dynamics, +.. +.. +.. .. math:: +.. \begin{aligned} +.. d e & = \left(\gamma \, R \, s \, i - \sigma e\right)dt + \zeta_e \sqrt{e} d W +.. \end{aligned} +.. :label: esde +.. + + + +Vacinations and Shocks to Mortality +==================================== + +The next step of randomness that we will consider involves uncertainty in technology. + +Mortality Shocks +------------------ There may be a variety of medical interventions that are short of a vaccine, but still effect the :math:`m(t)` path. In addition, imperfect mixing of different demographic groups could lead to aggregate shocks in mortality (e.g. if a retirement home is afflicted vs. an elementary school) @@ -579,7 +621,7 @@ We will begin by adding in sorts of random shocks, and leaving out dift or any m Combining the Shocks --------------------- -With these, we can define a variance term, mostly zeros since we only have two independent shocks, we can combined them in diagonal noise term :math:`\Omega(x, t)`. Extending +With these, we can define a variance term, mostly zeros since we only have two independent shocks, we can combined them in diagonal noise term :math:`G(x, t)`. Extending .. math:: \begin{aligned} @@ -587,14 +629,9 @@ With these, we can define a variance term, mostly zeros since we only have two i \zeta_e \sqrt{e}\\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ \zeta_R \sqrt{r} \\ & \zeta_m \sqrt{m} - \end{alignd} - -Finally, we can extend our existing :ref:`dfcv` definition to in + \end{aligned} -.. math:: - \begin{aligned} - d x &= F(x,t;p)dt + \Omega dW - \end{aligned} + TODO: IMPLEMENT THIS WITH SDE From cc0140ffd55d0520b440b3bfb9d6694662f13cf2 Mon Sep 17 00:00:00 2001 From: Jesse Perla Date: Thu, 25 Jun 2020 21:17:22 -0700 Subject: [PATCH 08/16] Done with the SDE. About to do the Jumps --- .../tools_and_techniques/seir_model_sde.rst | 134 ++++++++++++++---- 1 file changed, 110 insertions(+), 24 deletions(-) diff --git a/source/rst/tools_and_techniques/seir_model_sde.rst b/source/rst/tools_and_techniques/seir_model_sde.rst index c35c6a13..3b4a5ff7 100644 --- a/source/rst/tools_and_techniques/seir_model_sde.rst +++ b/source/rst/tools_and_techniques/seir_model_sde.rst @@ -71,7 +71,7 @@ others covered in previous lectures The SEIR Model -============= +============== In the version of the SEIR model we will analyze there are four states. @@ -302,7 +302,7 @@ The baseline parameters are put into a named tuple generator (see previous lectu .. code-block:: julia p_gen = @with_kw ( T = 550.0, γ = 1.0 / 18, σ = 1 / 5.2, η = 1.0 / 20, - R̄ = 1.6, B = (t, p) -> p.R̄, m_0 = 0.01) + R̄ = 1.6, B = (t, p) -> p.R̄, m_0 = 0.01, N = 3.3E8) Note that the default :math:`B(t)` function always equals :math:`\bar{R}` @@ -468,19 +468,19 @@ and 75,000 agents already exposed to the virus and thus soon to be contagious. B_lift_early(t, p) = t < 30.0 ? 0.5 : 2.0 B_lift_late(t, p) = t < 120.0 ? 0.5 : 2.0 + p_early = p_gen(B = B_lift_early, η = 10.0) + p_late = p_gen(B = B_lift_late, η = 10.0) + # initial conditions - N = 3.3E8 # US Population - i_0 = 25000 / N - e_0 = 75000 / N + i_0 = 25000 / p_early.N + e_0 = 75000 / p_early.N s_0 = 1.0 - i_0 - e_0 R_0 = 0.5 - x_0 = [s_0, e_0, i_0, 0.0, R_0, p.m_0, 0.0, 0.0] # start in lockdown + x_0 = [s_0, e_0, i_0, 0.0, R_0, p_early.m_0, 0.0, 0.0] # start in lockdown # create two problems, with rapid movement of R towards B(t) - p_early = p_gen(B = B_lift_early, η = 10.0) - p_late = p_gen(B = B_lift_late, η = 10.0) prob_early = ODEProblem(F, x_0, tspan, p_early) prob_late = ODEProblem(F, x_0, tspan, p_late) @@ -500,11 +500,13 @@ To calculate the daily death, calculate the :math:`\gamma i(t) m(t)`. .. code-block:: julia + flow_deaths(sol, p) = p.N * sol[3,:] .* sol[6,:] * p.γ + daily_early = sol_early[3,:] .* sol_early[6,:] * p_early.γ daily_late = sol_late[3,:] .* sol_late[6,:] * p_late.γ - plot(sol_early.t, daily_early, title = "Flow Deaths", label = "Lift Early") - plot!(sol_late.t, daily_late, label = "Lift Late") + plot(sol_early.t, flow_deaths(sol_early, p_early), title = "Flow Deaths", label = "Lift Early") + plot!(sol_late.t, flow_deaths(sol_late, p_late), label = "Lift Late") Pushing the peak of curve further into the future may reduce cumulative deaths if a vaccine is found. @@ -549,7 +551,7 @@ The general form of the SDE with these sorts of continuous-shocks is an extensio .. math:: \begin{aligned} - d x &= F(x,t;p)dt + \G(x,t;p) dW + d x &= F(x,t;p)dt + G(x,t;p) dW \end{aligned} Here, it is convenient to :math:`d W` with the same dimension as :math:`x` where we can use the matrix :math:`G(x,t;p)` to associate the shocks with the appropriate :math:`x`. @@ -566,20 +568,107 @@ Since these are additive shocks, we will not need to modify the :math:`F` from o First create a new settings generator, and and then define a `SDEProblem`__ with Diagonal Noise. -We solve the problem with the `SRIW1 `__ algorithm, an Adaptive strong order 1.5 and weak order 2.0 for diagonal/scalar Ito SDEs) +We solve the problem with the `SRA `__ algorithm (Adaptive strong order 1.5 methods for additive Ito and Stratonovich SDEs) .. code-block:: julia p_sde_gen = @with_kw ( T = 550.0, γ = 1.0 / 18, σ = 1 / 5.2, η = 1.0 / 20, - R̄ = 1.6, B = (t, p) -> p.R̄, m_0 = 0.01, ζ = 1E-3) + R̄ = 1.6, B = (t, p) -> p.R̄, m_0 = 0.01, ζ = 0.03, N = 3.3E8) + + p = p_sde_gen() + i_0 = 25000 / p.N + e_0 = 75000 / p.N + s_0 = 1.0 - i_0 - e_0 + R_0 = 1.5 * p.R̄ + x_0 = [s_0, e_0, i_0, 0.0, R_0, p.m_0, 0.0, 0.0] # start in lockdown + + G(x, p, t) = [0, 0, 0, 0, p.ζ*sqrt(x[5]), 0, 0, 0] + + prob = SDEProblem(F, G, x_0, (0, p.T), p) + sol_1 = solve(prob, SRA()) + +With stochastic differential equations, a "solution" is akin to a simulation for a particular realizaiton of the noise process. + +Plotting the number of infections for these two realizaitons of the shock process + +.. code-block:: julia + + # TODO: PICK BETTER PLOT + plot(sol_1, vars = [2, 5] , title = "Stochastic R(t)", label = ["e(t)" "R(t)") + +If we solve this model a second time, and plot the flow of deaths, we can see differences over time + +.. code-block:: julia + + sol_2 = solve(prob, SRA()) + plot(sol_1.t, flow_deaths(sol_1, p), title = "Daily Deaths", label = "sim 1") + plot!(sol_2.t, flow_deaths(sol_2, p), label = "sim 2") + + + +While individual simulatations are useful, you often want to look at an ensemble of multiple trajectories of the SDE + +.. code-block:: julia + + ensembleprob = EnsembleProblem(prob) + sol = solve(ensembleprob, SRA(), EnsembleSerial(),trajectories = 10) + plot(sol, vars = [3], title = "Infection Simulations", label = "i(t)") + + +Or, more frequently, you may want to run many trajectories and plot quantiles, which can be automatically run in `parallel Date: Fri, 26 Jun 2020 14:33:08 -0700 Subject: [PATCH 09/16] Finished the Jump Diffusion --- .../tools_and_techniques/seir_model_sde.rst | 322 ++++++++++++++---- 1 file changed, 257 insertions(+), 65 deletions(-) diff --git a/source/rst/tools_and_techniques/seir_model_sde.rst b/source/rst/tools_and_techniques/seir_model_sde.rst index 3b4a5ff7..34b4bc0e 100644 --- a/source/rst/tools_and_techniques/seir_model_sde.rst +++ b/source/rst/tools_and_techniques/seir_model_sde.rst @@ -209,11 +209,11 @@ We will assume that the transmission rate follows a process with a reversion to Finally, let :math:`m(t)` be the mortality rate, which we will leave constant for now, i.e. :math:`\frac{d m}{d t} = 0`. The cumulative deaths can be integrated through the flow :math:`\gamma i` entering the "Removed" state and define the cumulative number of deaths as :math:`M(t)`. .. math:: - \begin{aligned}\\ \frac{d m}{d t} &= 0\\ \frac{d M}{d t} &= m \gamma i \end{aligned} + :label: Mode While we could conveivably integate the total deaths given the solution to the model, it is more convenient to use the integrator built into the ODE solver. That is, we add :math:`d M(t)/dt` rather than calculating :math:`M(t) = \int_0^t \gamma m(\tau) i(\tau) d \tau` ex-post. @@ -264,14 +264,6 @@ As we will initially consider the case where :math:`R(0) = B`, the value of :mat Implementation ============== -# First we set the population size to match the US and the parameters as described - -# .. code-block:: julia - -# N = 3.3E8 # US Population -# γ = 1 / 18 -# σ = 1 / 5.2 -# η = 1 / 20 # a placeholder, inactive in first experiments First, construct our :math:`F` from :eq:`dfcv` @@ -518,11 +510,30 @@ One source of randomness which would enter the model is considering the discretn But rather than examining how granularity leads to aggregate fluctuations, we will concentrate on randomness that comes from aggregate changes in behavior or policy. -Aggregate Shocks to Transmission Rates -======================================= +Introduction to SDEs: Aggregate Shocks to Transmission Rates +============================================================== We will start by extending our model to include randomness in :math:`R(t)`, which makes it a system of Stochastic Differential Equations (SDEs). +Continous-Time Stochastic Processes +----------------------------------- + +In continuous-time, there is an important distinction between randomness that leads to continuous paths vs. those which may have jumps (which are almost surely right-continous). The most tractable of these includes the theory of `Levy Processes `__ (as known as Brownian Motion) which leads to a diffusion equations, and is the only continous-time Levy process with continuous paths +#. `Poisson Processes `__ with an arrival rate of jumps in the variable. + +Every other Levy Process you will typically work with is composed of those building blocks (e.g. a `Diffusion Process `__ is a diffusion process with a Poisson arrival of jumps). + +In this section, we will look at example of a diffusion process. + + Shocks to Transmission Rates ------------------------------ @@ -539,19 +550,19 @@ The notation for this `SDE `__. This equation is used in the `Cox-Ingersoll-Ross `__ and `Heston `__ models of interest rates and stochastic volatility. -Heuristically, if :math:`\zeta = 0`, we can divide by :math:`dt` and nest the original ODE in :eq:`Rode`. +Heuristically, if :math:`\zeta = 0`, divide this equation by :math:`dt` and it nests the original ODE in :eq:`Rode`. The general form of the SDE with these sorts of continuous-shocks is an extension of our :ref:`dfcv` definition . .. math:: \begin{aligned} - d x &= F(x,t;p)dt + G(x,t;p) dW + d x_t &= F(x_t,t_t;p)dt + G(x_t,t;p) dW \end{aligned} Here, it is convenient to :math:`d W` with the same dimension as :math:`x` where we can use the matrix :math:`G(x,t;p)` to associate the shocks with the appropriate :math:`x`. @@ -632,32 +643,28 @@ While ensembles, you may want to perform transformations, such as calculating ou function save_mortality(sol, ind) total = p.N * sol[8, :] flow = flow_deaths(sol, p) -# Ls = [sol.u[x][2] for x in 1:length(sol.u)] - # Zs = [sol.u[x][1] for x in 1:length(sol.u)] + #TODO Add both total and flow return (DiffEqBase.build_solution(sol.prob, sol.alg, sol.t, total), false) end sol = solve(ensembleprob, SRA(), EnsembleThreads(), output_func = save_mortality, trajectories = 1000) -For large-scale stochastic simulations, we can also use the ``output_func`` to reduce theamount of data stored for each simulation. +For large-scale stochastic simulations, we can also use the ``output_func`` to reduce the amount of data stored for each simulation. .. code-block:: julia - # CHRIS TODO: Can't Figure out plotting these two (in separate panes) + # CHRIS TODO: Can't Figure out plotting. Maybe both of them in separate panes summ = EnsembleSummary(sol) # defaults to saving 0.05, 0.5, and 0.95 quantiles summ2 = EnsembleSummary(sol, quantiles = [0.25, 0.75]) - plot(summ, title = "Total Deaths") - plot!(summ2, idxs = (3,), labels = "Middle 50%") + plot(summ, title = "Total and Daily Deaths (TBD)") + plot!(summ2, labels = "Middle 50%") Performance of these tends to be high, for example, rerunning out 1000 trajectories is measured in seconds on most computers with multithreading enabled. -.. code-block:: julia - - @time solve(ensembleprob, SRA(), EnsembleThreads(), output_func = save_mortality, trajectories = 1000); - .. -.. CHRIS: THIS ENDED UP SLOWER. CAN ADD BACK IF WE CAN FIGURE OUT WHY +.. CHRIS: THIS ENDED UP SLOWER and I am not sure why. CAN ADD BACK IF WE CAN FIGURE OUT WHY +.. .. Furthermore, we can exploit Julia's generic programming to use a static array (or GPU, .. if available) .. .. .. code-block:: julia @@ -667,70 +674,255 @@ Performance of these tends to be high, for example, rerunning out 1000 trajector .. ensembleprob = EnsembleProblem(prob) .. sol = solve(ensembleprob, SRA(), EnsembleThreads(),trajectories = 1000) .. @time solve(ensembleprob, SRA(), EnsembleThreads(),trajectories = 1000); -.. +.. +.. TODO: link to the generic programming lecture, and explain why care was needed in the F definition. -.. Migration and Transporation -.. ---------------------------- -.. -.. A second source of shocks are associated with policies where new individuals in the .. Exposed state enter the geography We will maintain a constant population size and .. assume (without specifying) compensating outflows to match the others, and assume that .. Infected individuals are effectively barred from entry. -.. -.. As it is the main consideration, lets add the diffusive term to the :math:`de` .. dynamics, -.. -.. -.. .. math:: -.. \begin{aligned} -.. d e & = \left(\gamma \, R \, s \, i - \sigma e\right)dt + \zeta_e \sqrt{e} d W -.. \end{aligned} -.. :label: esde -.. +Technological Progress and Policy Tradeoffs +============================================== +While the randomness inherent in the :math:`R(t)` can explain some of the sources of uncertainty that come from behavioral shocks, we have been ignoring two other considerations. +First, technology, both in treatment and vaccination, is evolving and in an inherently non-detemrinistic way. We will consider that the mortality rate :math:`m(t)` may evolve over time, as well as considering how a probability of vaccination development adds a path to the Removed state that does not require infection. -Vacinations and Shocks to Mortality -==================================== -The next step of randomness that we will consider involves uncertainty in technology. +In order to add one more degree of realism to the tradeoffs, we will consider that the actual death rate is driven by the mortality :math:`m(t)` but also capacity constraints in the economy with respect to medical resources for the infected state. In particular, we assume that if :math:`i(t) > \bar{i}`, then the available medical resources are exhuasted, leading to quadratically increased death rate. -Mortality Shocks ------------------- +Second, the only social objective measure we can explore with the current framework is minimizing the total deaths. That ignores the possible policy tradeoffs between minimizing deaths and the impact on the general economy. -There may be a variety of medical interventions that are short of a vaccine, but still effect the :math:`m(t)` path. In addition, imperfect mixing of different demographic groups could lead to aggregate shocks in mortality (e.g. if a retirement home is afflicted vs. an elementary school) +While a particular planner may decide that the only relevant welfare criteria is aggregate mortality, that leads to implausibly extreme policy (e.g. set :math:`B(t) = 0` forever). Hence, we will add a proxy for economic impact of Covid and the shotdown policy, summarized by :math:`u(t)` for excess covid-related "unemployment". A particular planner can then decide the weighting of the tradeoffs. -We will begin by adding in sorts of random shocks, and leaving out dift or any mean-reversion for simplicity +The policy choice :math:`B(t)` is then made a markov-process rather than current exogeous and deterministic one. + +The inherent discretness of medical innovations and policy changes provides us an opportunity to explore the use of Poisson and jump diffusion. + +Mortality +--------- + +Imperfect mixing of different demographic groups could lead to aggregate shocks in mortality (e.g. if a retirement home is afflicted vs. an elementary school). These sorts of relatively small changes might be best models as a continous path process, so we will add a diffusion term with volatility :math:`\xi\sqrt{m}` to the :math:`m` process + +In addition, there may be a variety of smaller medical interventions that are short of a vaccine, but still effect the :math:`m(t)` path. For simplicity, we will assume that each innovation cuts the inherent mortality rate in half and arrives with intensity :math:`\alpha \geq 0`. Combining these two leads to a jump diffusion process .. math:: \begin{aligned} - d m & = \zeta_m \sqrt{m} d W + d m_t & = \xi \sqrt{m_t} d W_t - \frac{m}{2} d N_t^{\alpha} \end{aligned} + :label: dmt -Combining the Shocks ---------------------- +In that notation, the :math:`d W_t` is still the increment of the brownian motion process while the :math:`d N_t^{\alpha}` is a Poisson `counting process `__ with rate :math:`\alpha`. -With these, we can define a variance term, mostly zeros since we only have two independent shocks, we can combined them in diagonal noise term :math:`G(x, t)`. Extending + +In previous versions of model, :ref:`Mode` had deaths as the integration of the death probability :math:`m` of the flow leaving the :math:`I` state. i.e. :math:`\frac{d M}{d t} = m \gamma i`. Instead, we will add an additional term that the dealh probabilty is :math:`\pi(m, i) := m + \psi \max(0, i - \bar{i})` increasing as :math:`i < \bar{i}`. The ``min`` is used to ensure that the mortality rate never goes above 1. The differential equation for cumulative deaths is then .. math:: - \begin{aligned} - diag(G) &:= \begin{bmatrix} 0\\ - \zeta_e \sqrt{e}\\ - 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ \zeta_R \sqrt{r} \\ - & \zeta_m \sqrt{m} + \begin{aligned}\\ + \pi(m, i) &:= m + \psi \max(0, i - \bar{i})\\ + d M_t &= \pi(m_t, i_t)\gamma i_t dt \end{aligned} + :label: Mode_nl - +Vaccination +------------- -TODO: IMPLEMENT THIS WITH SDE +The develeopmment of a vaccine is a classic discrete event. Define the availability of a vaccine as :math:`V_t` where the initial condition is :math:`V(0) = 0`. We assume that with an arrival rate :math:`\theta` a vaccine will be developed, leading to the a jump to :math:`V = 1`. -Jump Processes -================== +The process for this leads to the following Jump process -We will extend the above with 2 variations: -1. Assume that major medical advancements can arrival, which drop mortality, :math:`m(t)` in half. As there remains a diffusion term, the resulting :math:`d x`, becomes a Jump diffusionAlso vaccination .. math:: \begin{aligned} - d m & = \zeta_m \sqrt{m_t} d W + \theta N_t + d V_t & = (1 - V_t) d N_t^{\theta} \end{aligned} + :label: dV + +where the increment :math:`(1 - V_t)` is 1 if the vaccine hasn't been developed, and zero if it has. While the constant arrival rate :math:`\theta` is convenient for a simple model, there is no reason it couldn't be time or state dependent in the code or model. + +With a vaccine, the model ceases to be a simple SEIR model since it is possible move from :math:`S \to R` without passing through :math:`E` and :math:`I`. + +As vaccines are not instantaneously delivered to all, we can let :math:`\nu` be the rate at which the vaccine is administered if available. This leads to the following changes to :ref:`seir_system` + +.. math:: + \begin{aligned} + d s_t & = \left(- \gamma \, R_t \, s_t \, i _t - \nu \, V_t \, s_t\right)dt + \\ + d r_t & = (\gamma i_t + \nu \, V_t \, s_t) dt + \end{aligned} + :label: seir_system_vacc + + +Unemployment and Lockdown Policy +---------------------------------- + +As a simple summary of the economic and social distortions due to the covid, consider that an aggregate state :math:`u(t)` (indicating "excess" unemployment due to covid) increases as the :math:`B(t)`policy deviates from the natural :math:`\bar{R}` level, but can also increase due to flow of deaths from covid. + +The second part of this is an important consideration: if the death rate is extremely high, then opening the economy may not help much as individuals are reluctant to return to normal economic activities. We will assume that weights on these are :math:`\mu_B` and :math:`mu_M` respetivly. + +To represent the slow process of coming back to normal, we assume that the stock value :math:`u(t)` depreciates at rate :math:`\delta`. Put together gives the ODE, + +.. math:: + \begin{aligned} + d u_t & = (\mu_B(B_t - \bar{R})^2 + \mu_M \pi(s_t,i_t) -\delta) u_t dt + \end{aligned} + :label: du + +The policy choice :math:`B(t)` is markov, and will not consider implementation of forward looking behavior in this lecture. For a simple example, consider a policy driven by myopic political incentives and driven entirely by the death rates. If :math:`\pi_t > \bar{\pi}` then set :math:`B(t) = R_0` and leave it at :math:`B(t) = \bar{R}` otherwise. + +Without getting into the measure theory, this sort of bang-bang policy doesn't work as the process ceases to be a Levy Process (and much of the intuitive of considering measurability of stochastic processes and expected present discounted value fail). + +To avoid these concerns, a standard trick is to make this a Levy Process by having a Poisson "policy evaluation" rate (which can be put arbitrarily high) where the policy is adjusted according to the rule. Let that nuisance parameter by :math:`\Xi`, which gives the following pure-jump process for :math:`B(t)` + +.. math:: + \begin{aligned} + d B_t & = J(m, i, B) \pi( dN^{\Xi}_t\\ + J(m, i, B) &:= -(\pi(i,m) > bar{\pi})\max{B - R_0, 0} + (\pi(i,m) \leq \bar{\pi})\max{B - \bar{R}, 0} + \end{aligned} + + +Implementing the Complete System +--------------------------------- + +To our model, we have added in three new variables (:math:`u, V,` and :math:`B`) and a number of new parameters :math:`\bar{i}, \xi, \alpha, \theta, \nu, \mu_B, \mu_M, R_0, \Xi, \delta, \psi` + + +Stacking the :math:`u, V, B` at the end of the existing :math:`x` vector, we add or modify using :ref:`seir_system_vacc`, :ref:`Mode_nl`, and :ref:`du`. + +Here, we will move from an "out of place" to an inplace differential equation. + +.. code-block:: julia + + function F!(dx, x, p, t) + + s, e, i, r, R, m, c, M, u, V, B = x + @unpack σ, γ, η, ī, ξ, α, θ, ν, μ_B, μ_M, R_0, Ξ, δ, π_bar, ψ, R̄, δ = p + + π = min(1.0, m + ψ * max(0.0, i > ī)) + + dx[1] = -γ*R*s*i - ν*V*s; # ds/dt + dx[2] = γ*R*s*i - σ*e; # de/dt + dx[3] = σ*e - γ*i; # di/dt + dx[4] = γ*i + ν*V*s; # dr/dt + dx[5] = η*(B(t, p) - R); # dR/dt + dx[6] = 0.0; # dm/dt + dx[7] = σ*e; # dc/dt + dx[8] = π*γ*i; # dM/dt + dx[9] = (μ_B*(B - R̄)^2 + μ_M*π - δ)*u;# du/dt + dx[10] = 0.0; # dV/dt + dx[11] = 0.0; # dB/dt + end + +Note that the ``V, B`` terms do not have a drift as they are pure jump processes. + +Next, we need to consider how the variance term of the diffusion changes. With the exception of the new brownian +motion associated with the jump diffusion in :ref:`dmt`, everything else remains unchanged or zero. + +.. code-block:: julia + + function G!(dx, x, p, t) + @unpack ξ, ζ = p + dx .= 0.0 + R = x[5] + m = x[6] + dx[5] = ζ*sqrt(R) + dx[6] = ζ*sqrt(m) + end + + +Finally, we need to add in 3 jump processes which modify ``V, B, m`` separately. The connection between a jump and a variable is not necsesary, however. + + +Implementing the vaccination process, + + +.. code-block:: julia + + rate_V(x, p, t) = p.θ # could be a function of time or state + function affect_V!(x, p, t) + integrate.u[10] = 1.0 # set the vacination state = 1 + end + jump_V = ConstantRateJump(rate_V,affect_V!) + + +If the solver simulates an arrival rate of the jump, it calls the ``affect!`` function which takes in the current values through the ``integrator`` object and modifies the values directly. + +The other two jump processes are, + +.. code-block:: julia + + rate_B(x, p, t) = p.Ξ # constant + + function affect_B!(integrator) + @unpack σ, γ, η, ī, ξ, α, θ, ν, μ_B, μ_M, R_0, Ξ, δ, π_bar, ψ, R̄, δ = integrator.p + + m = integrator.u[6] + i = integrator.u[3] + π = min(1.0, m + ψ * max(0.0, i > ī)) + + if π > π_bar + integrate.u[11] = R_0 + else + integrate.u[11] = R̄ + end + end + jump_B = ConstantRateJump(rate_B, affect_B!) + + rate_m(x, p, t) = p.α + function affect_m!(x, p, t) + x[6] = x[6]/2 # half the inherent mortatilty + end + jump_m = ConstantRateJump(rate_m, affect_m!) + + +Collecting the new parameters and providing an extension of the initial condition with no vaccine or accumulated :math:`u(t)` + +.. code-block:: julia + + p_full_gen = @with_kw ( T = 550.0, γ = 1.0 / 18, σ = 1 / 5.2, η = 1.0 / 20, + R̄ = 1.6, R_0 = 0.5, m_0 = 0.01, ζ = 0.03, N = 3.3E8, + ī = 0.05, + ξ = 0.05, + α = 0.05, + θ = 0.05, + ν = 0.05, + μ_B = 0.05, + μ_M = 0.1, + Ξ = 5.0, + δ = 0.07, + π_bar = 0.05, + ψ = 2.0 + ) + + p = p_full_gen() + i_0 = 25000 / p.N + e_0 = 75000 / p.N + s_0 = 1.0 - i_0 - e_0 + x_0 = [s_0, e_0, i_0, 0.0, p.R̄, p.m_0, 0.0, 0.0, 0.0, 0.0, R̄] + + prob = SDEProblem(F!, G!, x_0, (0, p.T), p) + jump_prob = JumpProblem(prob, Direct(), jump_V, jump_B, jump_m) + + +Solving and simulating, + +.. code-block:: julia + + sol = solve(jump_prob, SRIW1()) + + + +TODO: Chris, lets make sure this is right, pick some parameters, and solve it? + + + +ROUGH NOTES ON THE SCIML LECTURE +===================================== + +The hope is that we have all or almost all of the model elements developed here, so that we can +focus on the machine learning/etc. in the following lecture and not implement any new model features (i.e. just shut things down as required). + -2. Vaccines arrival, :math:`v(t)` which being the process of directly enabling a "suceptible" to "removed jump". More vaccine arrivals speed up the process. +Some thoughts on the SciML lecture which would be separate: -Will consider that vaccine arrival rates are time varying. :math:`\alpha(t)`.... As there is no diffusion, this is called a Pure jump process. +1. Bayesian estimation and model uncertainty. Have uncertainty on the parameter such as the :math:`m_0` parameter. Show uncertainty that comes after the estimation of the parameter. +2. Explain how bayesian priors are a form of regularization. With rare events like this, you will always have +more parameters than data. +3. Solve the time-0 optimal policy problem using Flux of choosing the :math:`B` policy in some form. For the objective, consider the uncertainty and the assymetry of being wrong. +4. If we wanted to add in a model element, I think that putting in an asymptomatic state would be very helpful for explaining the role of uncertainty in evaluating the counterfactuals. \ No newline at end of file From 9f7790ecfc67ebed1991a99126dd340f078ca53a Mon Sep 17 00:00:00 2001 From: Jesse Perla Date: Fri, 26 Jun 2020 17:19:18 -0700 Subject: [PATCH 10/16] JuMP Diffusion compiling and running --- source/rst/Manifest.toml | 26 ++-- source/rst/Project.toml | 1 + .../tools_and_techniques/seir_model_sde.rst | 115 +++++++++--------- 3 files changed, 72 insertions(+), 70 deletions(-) diff --git a/source/rst/Manifest.toml b/source/rst/Manifest.toml index 5c3a5994..36b4585e 100644 --- a/source/rst/Manifest.toml +++ b/source/rst/Manifest.toml @@ -365,10 +365,10 @@ uuid = "2b5f629d-d688-5b77-993f-72d75c75574e" version = "6.36.3" [[DiffEqJump]] -deps = ["Compat", "DataStructures", "DiffEqBase", "FunctionWrappers", "LinearAlgebra", "Parameters", "PoissonRandom", "Random", "RandomNumbers", "RecursiveArrayTools", "StaticArrays", "Statistics", "TreeViews"] -git-tree-sha1 = "4c960b9cf209426b72de274334f5f77b4d9a4dad" +deps = ["ArrayInterface", "Compat", "DataStructures", "DiffEqBase", "FunctionWrappers", "LinearAlgebra", "Parameters", "PoissonRandom", "Random", "RandomNumbers", "RecursiveArrayTools", "StaticArrays", "Statistics", "TreeViews"] +git-tree-sha1 = "0dfb3aa84f687c78f697173665e304d28651dd55" uuid = "c894b116-72e5-5b58-be3c-e6d8d4ac2b12" -version = "6.7.6" +version = "6.9.2" [[DiffEqNoiseProcess]] deps = ["DataStructures", "DiffEqBase", "LinearAlgebra", "Random", "RandomNumbers", "RecipesBase", "RecursiveArrayTools", "Requires", "ResettableStacks", "StaticArrays", "Statistics"] @@ -865,9 +865,9 @@ uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb" [[Libiconv_jll]] deps = ["Libdl", "Pkg"] -git-tree-sha1 = "e5256a3b0ebc710dbd6da0c0b212164a3681037f" +git-tree-sha1 = "c9d4035d7481bcdff2babf5a55525a818ef8ed8f" uuid = "94ce4f54-9a6c-5748-9c1c-f9c7231a4531" -version = "1.16.0+2" +version = "1.16.0+5" [[LightGraphs]] deps = ["ArnoldiMethod", "DataStructures", "Distributed", "Inflate", "LinearAlgebra", "Random", "SharedArrays", "SimpleTraits", "SparseArrays", "Statistics"] @@ -1123,9 +1123,9 @@ version = "0.3.9+4" [[OpenSSL_jll]] deps = ["Libdl", "Pkg"] -git-tree-sha1 = "d2a6f25262d568b5a7e454cf7ff5066a79d16c7d" +git-tree-sha1 = "7aaaded15bf393b5f34c2aad5b765c18d26cb495" uuid = "458c3c95-2e84-50aa-8efc-19380b2a3a95" -version = "1.1.1+2" +version = "1.1.1+4" [[OpenSpecFun_jll]] deps = ["CompilerSupportLibraries_jll", "Libdl", "Pkg"] @@ -1379,9 +1379,9 @@ version = "0.6.1" [[Rmath_jll]] deps = ["Libdl", "Pkg"] -git-tree-sha1 = "1660f8fefbf5ab9c67560513131d4e933012fc4b" +git-tree-sha1 = "d76185aa1f421306dec73c057aa384bad74188f0" uuid = "f50d1b31-88e8-58de-be2c-1cc44531875f" -version = "0.2.2+0" +version = "0.2.2+1" [[Roots]] deps = ["Printf"] @@ -1673,9 +1673,9 @@ version = "0.5.2" [[XML2_jll]] deps = ["Libdl", "Libiconv_jll", "Pkg", "Zlib_jll"] -git-tree-sha1 = "987c02a43fa10a491a5f0f7c46a6d3559ed6a8e2" +git-tree-sha1 = "6b2ffe6728bdba991da4fc1aa5980a53db46a23f" uuid = "02c8fc9c-b97f-50b9-bbe4-9be30ff0a78a" -version = "2.9.9+4" +version = "2.9.9+5" [[ZipFile]] deps = ["Libdl", "Printf", "Zlib_jll"] @@ -1685,9 +1685,9 @@ version = "0.9.2" [[Zlib_jll]] deps = ["Libdl", "Pkg"] -git-tree-sha1 = "a2e0d558f6031002e380a90613b199e37a8565bf" +git-tree-sha1 = "622d8b6dc0c7e8029f17127703de9819134d1b71" uuid = "83775a58-1f1d-513f-b197-d71354ab007a" -version = "1.2.11+10" +version = "1.2.11+14" [[Zygote]] deps = ["AbstractFFTs", "ArrayLayouts", "DiffRules", "FillArrays", "ForwardDiff", "Future", "IRTools", "InteractiveUtils", "LinearAlgebra", "MacroTools", "NNlib", "NaNMath", "Random", "Requires", "SpecialFunctions", "Statistics", "ZygoteRules"] diff --git a/source/rst/Project.toml b/source/rst/Project.toml index f6f0ba35..b6f534dd 100644 --- a/source/rst/Project.toml +++ b/source/rst/Project.toml @@ -15,6 +15,7 @@ Compat = "34da2185-b29b-5c13-b0c7-acf172513d20" DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" DataFramesMeta = "1313f7d8-7da2-5740-9ea0-a2ca25f37964" Dierckx = "39dd38d3-220a-591b-8e3c-4c3a8c710a94" +DiffEqJump = "c894b116-72e5-5b58-be3c-e6d8d4ac2b12" DiffEqOperators = "9fdde737-9c7f-55bf-ade8-46b3f136cc48" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" Expectations = "2fe49d83-0758-5602-8f54-1f90ad0d522b" diff --git a/source/rst/tools_and_techniques/seir_model_sde.rst b/source/rst/tools_and_techniques/seir_model_sde.rst index 34b4bc0e..aa4436ed 100644 --- a/source/rst/tools_and_techniques/seir_model_sde.rst +++ b/source/rst/tools_and_techniques/seir_model_sde.rst @@ -8,12 +8,13 @@ .. contents:: :depth: 2 - - Overview ============= -This is a Julia version of the code for analyzing the COVID-19 pandemic. +Coauthored with Chris Rackauckas + + +This is a Julia version of code for analyzing the COVID-19 pandemic. The purpose of these notes is to introduce economists to quantitative modeling of infectious disease dynamics, and to modeling with ordinary and stochastic differential @@ -29,8 +30,6 @@ In the first part, dynamics are modeled using a standard SEIR (Susceptible-Expos of disease spread, represented as a system of ordinary differential equations where the number of agents is large and there are no exogenous stochastic shocks. - - The first part of the model follows the notes from provided by `Andrew Atkeson `__ @@ -66,7 +65,7 @@ others covered in previous lectures .. code-block:: julia - using OrdinaryDiffEq, StochasticDiffEq + using OrdinaryDiffEq, StochasticDiffEq, DiffEqJump using Parameters, StaticArrays, Plots @@ -153,11 +152,11 @@ First, define the system of equations ] end -Written this way, we see that the four equations represent the one-directional transition from the susceptible to removed state, where the negative terms are outflows, and the positive ones inflows. +Written this way, we see that the four equations represent the one-directional transition from the susceptible to removed state, where the negative terms are outflows, and the positive ones inflows. The equation associated with the state :math:`R` has inflows and no outflows, so any agents which end up in that state are "trapped". -As there is no flow leaving the :math:`dr/dt` and all parameters are positive, unless we start with a degenerate initial condition (e.g. :math:`e(0) = i(0) = 0`) the "Removed" state is asymptoically absorbing, and :math:`\lim_{t\to \infty} r(t) = 1`. Crucial to this result is that individuals are perfectly divisible, and any arbitrarily small :math:`i > 0` leads to a strictly positive flow into the exposed state. +Consequently, with no flow leaving the :math:`dr/dt` and strictly positive parameter values the ``R`` state is absorbing, and :math:`\lim_{t\to \infty} r(t) = 1`. Crucial to this result is that individuals are perfectly divisible, and any arbitrarily small :math:`i > 0` leads to a strictly positive flow into the exposed state. -We will discuss this topic further in the lecture on continuous-time +We will discuss this topic (and related topics of irreducibility) in the lecture on continuous-time markov-chains, as well as the limitations of these approximations when the discretness becomes essential (e.g. continuum approximations are incapable of modeling extinguishing of an outbreak). Given this system, we choose an initial condition and a timespan, and create a ``ODEProblem`` encapsulating the system. @@ -181,7 +180,7 @@ With this, we can choose an ODE algorithm (e.g. a good default for non-stiff ODE plot(sol, labels = ["s" "e" "i" "r"], title = "SEIR Dynamics", lw = 2) -We did not provide either a set of timesteps or a ``dt`` time stepsize to the ``solve``. The reason is that most accurate and high-performance ODE solvers appropriate use adaptive time-stepping, changing the stepsize based the degree of curvature in the derivatives. +We did not provide either a set of timesteps or a ``dt`` time stepsize to the ``solve``. Most accurate and high-performance ODE solvers appropriate for this problem use adaptive time-stepping, changing the stepsize based the degree of curvature in the derivatives. Or, as an alternative visualization, the proportions in each state over time @@ -201,11 +200,12 @@ First, we can consider some additional calculations such as the cumulative casel We will assume that the transmission rate follows a process with a reversion to a value :math:`B(t)` which could conceivably be influenced by policy. The intuition is that even if the targetted :math:`B(t)` was changed through social distancing/etc., lags in behavior and implementation would smooth out the transition, where :math:`\eta` governs the speed of :math:`R(t)` moves towards :math:`B(t)`. .. math:: - \begin{aligned} - \frac{d R}{d t} &= \eta (B - R) - \end{aligned} + \begin{aligned} + \frac{d R}{d t} &= \eta (B - R)\\ + \end{aligned} :label: Rode + Finally, let :math:`m(t)` be the mortality rate, which we will leave constant for now, i.e. :math:`\frac{d m}{d t} = 0`. The cumulative deaths can be integrated through the flow :math:`\gamma i` entering the "Removed" state and define the cumulative number of deaths as :math:`M(t)`. .. math:: @@ -223,7 +223,7 @@ The system :eq:`seir_system` and the supplemental equations can be written in ve .. math:: \begin{aligned} - \frac{d x}{d t} &= F(x,t;p) + \frac{d x}{d t} &= F(x,t;p)\\ &:= \begin{bmatrix} - \gamma \, R \, s \, i \\ @@ -313,7 +313,7 @@ Setting initial conditions, we will assume a fixed :math:`i, e`, :math:`r=0`, :m prob = ODEProblem(F, x_0, (0.0, p.T), p) -The ``tspan`` determines that the :math:`t` used by the sovler, where the scale needs to be consistent with the arrival +The ``tspan`` of ``(0.0, p.T)`` determines that the :math:`t` used by the sovler. The time scale needs to be consistent with the arrival rate of the transition probabilities (i.e. the :math:`\gamma, \sigma` were chosen based on daily data). The time period we investigate will be 550 days, or around 18 months: @@ -330,7 +330,7 @@ Let's run some experiments using this code. We see that the adaptive time-stepping used approximately 50 time-steps to solve this problem to the desires accuracy. Evaluating the solver at points outside of those time-steps uses the an interpolator consistent with the solution to the ODE. -See `here `__ for details on analyzing the solution, and `here `__ for plotting tools. The built-in plots for the solutions provide all of the `attributes `__. +See `here `__ for details on analyzing the solution, and `here `__ for plotting tools. The built-in plots for the solutions provide all of the `attributes `__ in `Plots.jl `__. .. code_block:: julia @@ -348,7 +348,7 @@ We calculate the time path of infected people under different assumptions. .. code-block:: julia R̄_vals = range(1.6, 3.0, length = 6) - sols = [solve(ODEProblem(F, x_0, tspan, p_gen(R̄ = R̄)), Tsit5()) for R̄ in R̄_vals] + sols = [solve(ODEProblem(F, x_0, tspan, p_gen(R̄ = R̄)), Tsit5()) for R̄ in R̄_vals]; # TODO: Probably clean ways to plot this , but don't know them! @@ -397,7 +397,7 @@ Experiment 2: Changing Mitigation Let's look at a scenario where mitigation (e.g., social distancing) is successively imposed, but the target (:math:`\bar{R}` is fixed) -To do this, we will have :math:`R(0) \neq \bar{R}` and examine the dynamics using the :math:`\frac{d R}{d t} &= \eta (\bar{R} - R)` ODE. +To do this, we will have :math:`R(0) \neq \bar{R}` and examine the dynamics using the :math:`\frac{d R}{d t} = \eta (\bar{R} - R)` ODE. .. Mathematica Verification .. (\[Beta][t] /. @@ -518,7 +518,7 @@ We will start by extending our model to include randomness in :math:`R(t)`, whic Continous-Time Stochastic Processes ----------------------------------- -In continuous-time, there is an important distinction between randomness that leads to continuous paths vs. those which may have jumps (which are almost surely right-continous). The most tractable of these includes the theory of `Levy Processes `_. **TBD:** Add definition of levy processes and the intuitive connection between statoionary increments and independence of increments. Also, is there intuition why additive processes in continuous time lead to additive SDEs/etc.? Might be useful. @@ -529,7 +529,7 @@ Unlike in discrete-time, where a modeller has license to be creative, the rules #. `Weiner Processes `__ (as known as Brownian Motion) which leads to a diffusion equations, and is the only continous-time Levy process with continuous paths #. `Poisson Processes `__ with an arrival rate of jumps in the variable. -Every other Levy Process you will typically work with is composed of those building blocks (e.g. a `Diffusion Process `__ is a diffusion process with a Poisson arrival of jumps). +Every other Levy Process you will typically work with is composed of those building blocks (e.g. a `Diffusion Process `__ such as Geometric Brownian Motion is a transformation of a Weiner process, and a `jump diffusion `__ is a diffusion process with a Poisson arrival of jumps). In this section, we will look at example of a diffusion process. @@ -549,20 +549,20 @@ To implement this, we will add on a diffusion term to :eq:`Rode` with an instant The notation for this `SDE `__ is then .. math:: - \begin{aligned} - d R_t &= \eta (B_t - R_t) dt + \zeta \sqrt{R_t} dW_t + \begin{aligned} + d R_t &= \eta (B_t - R_t) dt + \zeta \sqrt{R_t} dW_t\\ \end{aligned} :label: Rsde where :math:`W` is standard Brownian motion (i.e a `Weiner Process `__. This equation is used in the `Cox-Ingersoll-Ross `__ and `Heston `__ models of interest rates and stochastic volatility. -Heuristically, if :math:`\zeta = 0`, divide this equation by :math:`dt` and it nests the original ODE in :eq:`Rode`. +Heuristically, if :math:`\zeta = 0`, divide this equation by :math:`dt` and it nests the original ODE in :eq:`Rode`. -The general form of the SDE with these sorts of continuous-shocks is an extension of our :ref:`dfcv` definition . +The general form of the SDE with these sorts of continuous-shocks is an extension of our :eq:`dfcv` definition . .. math:: \begin{aligned} - d x_t &= F(x_t,t_t;p)dt + G(x_t,t;p) dW + d x_t &= F(x_t,t;p)dt + G(x_t,t;p) dW \end{aligned} Here, it is convenient to :math:`d W` with the same dimension as :math:`x` where we can use the matrix :math:`G(x,t;p)` to associate the shocks with the appropriate :math:`x`. @@ -577,14 +577,14 @@ As the shock only effects :math:`dR`, which is the 5th equation, define the matr Since these are additive shocks, we will not need to modify the :math:`F` from our equation. -First create a new settings generator, and and then define a `SDEProblem`__ with Diagonal Noise. +First create a new settings generator, and and then define a `SDEProblem `__ with Diagonal Noise. We solve the problem with the `SRA `__ algorithm (Adaptive strong order 1.5 methods for additive Ito and Stratonovich SDEs) .. code-block:: julia p_sde_gen = @with_kw ( T = 550.0, γ = 1.0 / 18, σ = 1 / 5.2, η = 1.0 / 20, - R̄ = 1.6, B = (t, p) -> p.R̄, m_0 = 0.01, ζ = 0.03, N = 3.3E8) + R̄ = 1.6, B = (t, p) -> p.R̄, m_0 = 0.01, ζ = 0.03, N = 3.3E8) p = p_sde_gen() i_0 = 25000 / p.N @@ -596,7 +596,8 @@ We solve the problem with the `SRA `_ using multiple threads, processes, or GPUs. .. code-block:: julia sol = solve(ensembleprob, SRA(), EnsembleThreads(),trajectories = 1000) - summ2 = EnsembleSummary(sol) # defaults to saving 0.05, 0.5, and 0.95 quantiles + summ = EnsembleSummary(sol) # defaults to saving 0.05, 0.5, and 0.95 quantiles plot(summ, idxs = (3,), title = "Quantiles of Infections Ensemble", labels = "Middle 95% Quantile") @@ -705,17 +706,17 @@ In addition, there may be a variety of smaller medical interventions that are sh .. math:: \begin{aligned} - d m_t & = \xi \sqrt{m_t} d W_t - \frac{m}{2} d N_t^{\alpha} + d m_t & = \xi \sqrt{m_t} d W_t - \frac{m}{2} d N_t^{\alpha}\\ \end{aligned} :label: dmt In that notation, the :math:`d W_t` is still the increment of the brownian motion process while the :math:`d N_t^{\alpha}` is a Poisson `counting process `__ with rate :math:`\alpha`. -In previous versions of model, :ref:`Mode` had deaths as the integration of the death probability :math:`m` of the flow leaving the :math:`I` state. i.e. :math:`\frac{d M}{d t} = m \gamma i`. Instead, we will add an additional term that the dealh probabilty is :math:`\pi(m, i) := m + \psi \max(0, i - \bar{i})` increasing as :math:`i < \bar{i}`. The ``min`` is used to ensure that the mortality rate never goes above 1. The differential equation for cumulative deaths is then +In previous versions of model, :eq:`Mode` had deaths as the integration of the death probability :math:`m` of the flow leaving the :math:`I` state. i.e. :math:`\frac{d M}{d t} = m \gamma i`. Instead, we will add an additional term that the dealh probabilty is :math:`\pi(m, i) := m + \psi \max(0, i - \bar{i})` increasing as :math:`i < \bar{i}`. The ``min`` is used to ensure that the mortality rate never goes above 1. The differential equation for cumulative deaths is then .. math:: - \begin{aligned}\\ + \begin{aligned} \pi(m, i) &:= m + \psi \max(0, i - \bar{i})\\ d M_t &= \pi(m_t, i_t)\gamma i_t dt \end{aligned} @@ -729,7 +730,7 @@ The develeopmment of a vaccine is a classic discrete event. Define the availabi The process for this leads to the following Jump process .. math:: - \begin{aligned} + \begin{aligned} d V_t & = (1 - V_t) d N_t^{\theta} \end{aligned} :label: dV @@ -738,7 +739,7 @@ where the increment :math:`(1 - V_t)` is 1 if the vaccine hasn't been developed, With a vaccine, the model ceases to be a simple SEIR model since it is possible move from :math:`S \to R` without passing through :math:`E` and :math:`I`. -As vaccines are not instantaneously delivered to all, we can let :math:`\nu` be the rate at which the vaccine is administered if available. This leads to the following changes to :ref:`seir_system` +As vaccines are not instantaneously delivered to all, we can let :math:`\nu` be the rate at which the vaccine is administered if available. This leads to the following changes to :eq:`seir_system` .. math:: \begin{aligned} @@ -760,7 +761,7 @@ To represent the slow process of coming back to normal, we assume that the stock .. math:: \begin{aligned} - d u_t & = (\mu_B(B_t - \bar{R})^2 + \mu_M \pi(s_t,i_t) -\delta) u_t dt + d u_t & = (\mu_B(B_t - \bar{R})^2 + \mu_M \pi(s_t,i_t) -\delta) u_t dt\\ \end{aligned} :label: du @@ -772,8 +773,8 @@ To avoid these concerns, a standard trick is to make this a Levy Process by havi .. math:: \begin{aligned} - d B_t & = J(m, i, B) \pi( dN^{\Xi}_t\\ - J(m, i, B) &:= -(\pi(i,m) > bar{\pi})\max{B - R_0, 0} + (\pi(i,m) \leq \bar{\pi})\max{B - \bar{R}, 0} + d B_t & = J(m, i, B) dN^{\Xi}_t\\ + J(m, i, B) &:= -(\pi(i,m) > \bar{\pi})\max(B - R_0, 0) + (\pi(i,m) \leq \bar{\pi})\max(B - \bar{R}, 0) \\ \end{aligned} @@ -783,7 +784,7 @@ Implementing the Complete System To our model, we have added in three new variables (:math:`u, V,` and :math:`B`) and a number of new parameters :math:`\bar{i}, \xi, \alpha, \theta, \nu, \mu_B, \mu_M, R_0, \Xi, \delta, \psi` -Stacking the :math:`u, V, B` at the end of the existing :math:`x` vector, we add or modify using :ref:`seir_system_vacc`, :ref:`Mode_nl`, and :ref:`du`. +Stacking the :math:`u, V, B` at the end of the existing :math:`x` vector, we add or modify using :eq:`seir_system_vacc`, :eq:`Mode_nl`, and :eq:`du`. Here, we will move from an "out of place" to an inplace differential equation. @@ -800,7 +801,7 @@ Here, we will move from an "out of place" to an inplace differential equation. dx[2] = γ*R*s*i - σ*e; # de/dt dx[3] = σ*e - γ*i; # di/dt dx[4] = γ*i + ν*V*s; # dr/dt - dx[5] = η*(B(t, p) - R); # dR/dt + dx[5] = η*(B - R); # dR/dt dx[6] = 0.0; # dm/dt dx[7] = σ*e; # dc/dt dx[8] = π*γ*i; # dM/dt @@ -812,17 +813,17 @@ Here, we will move from an "out of place" to an inplace differential equation. Note that the ``V, B`` terms do not have a drift as they are pure jump processes. Next, we need to consider how the variance term of the diffusion changes. With the exception of the new brownian -motion associated with the jump diffusion in :ref:`dmt`, everything else remains unchanged or zero. +motion associated with the jump diffusion in :eq:`dmt`, everything else remains unchanged or zero. .. code-block:: julia function G!(dx, x, p, t) - @unpack ξ, ζ = p - dx .= 0.0 + @unpack ξ, ζ = p + dx .= 0.0 R = x[5] m = x[6] - dx[5] = ζ*sqrt(R) - dx[6] = ζ*sqrt(m) + dx[5] = R <= 0.0 ? 0.0 : ζ*sqrt(R) + dx[6] = m <= 0.0 ? 0.0 : ζ*sqrt(m) end @@ -835,8 +836,8 @@ Implementing the vaccination process, .. code-block:: julia rate_V(x, p, t) = p.θ # could be a function of time or state - function affect_V!(x, p, t) - integrate.u[10] = 1.0 # set the vacination state = 1 + function affect_V!(integrator) + integrator.u[10] = 1.0 # set the vacination state = 1 end jump_V = ConstantRateJump(rate_V,affect_V!) @@ -857,16 +858,16 @@ The other two jump processes are, π = min(1.0, m + ψ * max(0.0, i > ī)) if π > π_bar - integrate.u[11] = R_0 + integrator.u[11] = R_0 else - integrate.u[11] = R̄ + integrator.u[11] = R̄ end end jump_B = ConstantRateJump(rate_B, affect_B!) rate_m(x, p, t) = p.α - function affect_m!(x, p, t) - x[6] = x[6]/2 # half the inherent mortatilty + function affect_m!(integrator) + integrator.u[6] = integrator.u[6]/2 # half the inherent mortatilty end jump_m = ConstantRateJump(rate_m, affect_m!) @@ -879,8 +880,8 @@ Collecting the new parameters and providing an extension of the initial conditio R̄ = 1.6, R_0 = 0.5, m_0 = 0.01, ζ = 0.03, N = 3.3E8, ī = 0.05, ξ = 0.05, - α = 0.05, - θ = 0.05, + α = 1/100, # mean 100 days to innovation + θ = 1/300, # mean 200 days to vaccine ν = 0.05, μ_B = 0.05, μ_M = 0.1, @@ -894,7 +895,7 @@ Collecting the new parameters and providing an extension of the initial conditio i_0 = 25000 / p.N e_0 = 75000 / p.N s_0 = 1.0 - i_0 - e_0 - x_0 = [s_0, e_0, i_0, 0.0, p.R̄, p.m_0, 0.0, 0.0, 0.0, 0.0, R̄] + x_0 = [s_0, e_0, i_0, 0.0, p.R̄, p.m_0, 0.0, 0.0, 0.0, 0.0, p.R̄] prob = SDEProblem(F!, G!, x_0, (0, p.T), p) jump_prob = JumpProblem(prob, Direct(), jump_V, jump_B, jump_m) @@ -905,11 +906,11 @@ Solving and simulating, .. code-block:: julia sol = solve(jump_prob, SRIW1()) - + plot(sol, vars = [8, 10]) TODO: Chris, lets make sure this is right, pick some parameters, and solve it? - +TODO: What sorts of experiments afte rthis? ROUGH NOTES ON THE SCIML LECTURE From f43b6816f0862aa559b01367f8be8de906c234be Mon Sep 17 00:00:00 2001 From: Jesse Perla Date: Sat, 27 Jun 2020 08:13:49 -0700 Subject: [PATCH 11/16] A few more visualization ideas --- source/rst/tools_and_techniques/seir_model_sde.rst | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/source/rst/tools_and_techniques/seir_model_sde.rst b/source/rst/tools_and_techniques/seir_model_sde.rst index aa4436ed..44380d7a 100644 --- a/source/rst/tools_and_techniques/seir_model_sde.rst +++ b/source/rst/tools_and_techniques/seir_model_sde.rst @@ -908,10 +908,11 @@ Solving and simulating, sol = solve(jump_prob, SRIW1()) plot(sol, vars = [8, 10]) - -TODO: Chris, lets make sure this is right, pick some parameters, and solve it? -TODO: What sorts of experiments afte rthis? - +TODO: Chris, lets make sure this is right, pick some parameters. +TODO: Lets show the same SIER decomposition we had before as a proportion, but with a vaccination arrival. +TODO: Maybe show a slice at time T = 550 or whatever of the distribution M(t) and u(t) as a histogram with 5, 50, and 95 percent confidence intervals? +TODO: Then show how those two change as the myopic B(t) policy changes? +TODO: Show how M(t) and u(t) change as the eta changes for a given policy? ROUGH NOTES ON THE SCIML LECTURE ===================================== @@ -926,4 +927,5 @@ Some thoughts on the SciML lecture which would be separate: 2. Explain how bayesian priors are a form of regularization. With rare events like this, you will always have more parameters than data. 3. Solve the time-0 optimal policy problem using Flux of choosing the :math:`B` policy in some form. For the objective, consider the uncertainty and the assymetry of being wrong. -4. If we wanted to add in a model element, I think that putting in an asymptomatic state would be very helpful for explaining the role of uncertainty in evaluating the counterfactuals. \ No newline at end of file +4. If we wanted to add in a model element, I think that putting in an asymptomatic state would be very helpful for explaining the role of uncertainty in evaluating the counterfactuals. +5. Solve for the optimal markovian B(t) policy through simulation and using Flux. From ca88bdf5765fe092a04eea0bcd48d968c28744e0 Mon Sep 17 00:00:00 2001 From: John Stachurski Date: Mon, 29 Jun 2020 17:23:24 +1000 Subject: [PATCH 12/16] minor corrections --- .../tools_and_techniques/seir_model_sde.rst | 76 +++++++++---------- 1 file changed, 36 insertions(+), 40 deletions(-) diff --git a/source/rst/tools_and_techniques/seir_model_sde.rst b/source/rst/tools_and_techniques/seir_model_sde.rst index 44380d7a..e7767034 100644 --- a/source/rst/tools_and_techniques/seir_model_sde.rst +++ b/source/rst/tools_and_techniques/seir_model_sde.rst @@ -16,12 +16,10 @@ Coauthored with Chris Rackauckas This is a Julia version of code for analyzing the COVID-19 pandemic. -The purpose of these notes is to introduce economists to quantitative modeling -of infectious disease dynamics, and to modeling with ordinary and stochastic differential +The purpose of these notes is to introduce economists to quantitative modeling of infectious disease dynamics, and to modeling with ordinary and stochastic differential equations. -The main objective is to study the impact of suppression through social -distancing on the spread of the infection. +The main objective is to study the impact of suppression through social distancing on the spread of the infection. The focus is on US outcomes but the parameters can be adjusted to study other countries. @@ -41,7 +39,7 @@ See further variations on the classic SIR model in Julia `here 0` leads to a strictly positive flow into the exposed state. -We will discuss this topic (and related topics of irreducibility) in the lecture on continuous-time -markov-chains, as well as the limitations of these approximations when the discretness becomes essential (e.g. continuum approximations are incapable of modeling extinguishing of an outbreak). +We will discuss this topic (and related topics of irreducibility) in the lecture on continuous-time Markov chains, as well as the limitations of these approximations when the discretness becomes essential (e.g. continuum approximations are incapable of modeling extinguishing of an outbreak). Given this system, we choose an initial condition and a timespan, and create a ``ODEProblem`` encapsulating the system. @@ -180,14 +177,14 @@ With this, we can choose an ODE algorithm (e.g. a good default for non-stiff ODE plot(sol, labels = ["s" "e" "i" "r"], title = "SEIR Dynamics", lw = 2) -We did not provide either a set of timesteps or a ``dt`` time stepsize to the ``solve``. Most accurate and high-performance ODE solvers appropriate for this problem use adaptive time-stepping, changing the stepsize based the degree of curvature in the derivatives. +We did not provide either a set of time steps or a ``dt`` time step size to the ``solve``. Most accurate and high-performance ODE solvers appropriate for this problem use adaptive time-stepping, changing the step size based the degree of curvature in the derivatives. Or, as an alternative visualization, the proportions in each state over time .. code-block:: julia - areaplot(sol.t, sol', labels = ["s" "e" "i" "r"], title = "SIER Proportions") + areaplot(sol.t, sol', labels = ["s" "e" "i" "r"], title = "SEIR Proportions") While implementing the system of ODEs in :math:`(s, e, i)`, we will extend the basic model to enable some policy experiments and calculations of aggregate values. @@ -195,9 +192,9 @@ While implementing the system of ODEs in :math:`(s, e, i)`, we will extend the b Extending the Model ----------------------- -First, we can consider some additional calculations such as the cumulative caseload (i.e., all those who have or have had the infection) as :math:`c = i + r`. Differentiating that expression and substituing from the time-derivatives of :math:`i(t), r(t)` yields :math:`\frac{d c}{d t} = \sigma e` +First, we can consider some additional calculations such as the cumulative caseload (i.e., all those who have or have had the infection) as :math:`c = i + r`. Differentiating that expression and substituting from the time-derivatives of :math:`i(t), r(t)` yields :math:`\frac{d c}{d t} = \sigma e` -We will assume that the transmission rate follows a process with a reversion to a value :math:`B(t)` which could conceivably be influenced by policy. The intuition is that even if the targetted :math:`B(t)` was changed through social distancing/etc., lags in behavior and implementation would smooth out the transition, where :math:`\eta` governs the speed of :math:`R(t)` moves towards :math:`B(t)`. +We will assume that the transmission rate follows a process with a reversion to a value :math:`B(t)` which could conceivably be influenced by policy. The intuition is that even if the targeted :math:`B(t)` was changed through social distancing/etc., lags in behavior and implementation would smooth out the transition, where :math:`\eta` governs the speed of :math:`R(t)` moves towards :math:`B(t)`. .. math:: \begin{aligned} @@ -215,7 +212,7 @@ Finally, let :math:`m(t)` be the mortality rate, which we will leave constant fo \end{aligned} :label: Mode -While we could conveivably integate the total deaths given the solution to the model, it is more convenient to use the integrator built into the ODE solver. That is, we add :math:`d M(t)/dt` rather than calculating :math:`M(t) = \int_0^t \gamma m(\tau) i(\tau) d \tau` ex-post. +While we could conceivably integrate the total deaths given the solution to the model, it is more convenient to use the integrator built into the ODE solver. That is, we add :math:`d M(t)/dt` rather than calculating :math:`M(t) = \int_0^t \gamma m(\tau) i(\tau) d \tau` ex-post. This is a common trick when solving systems of ODEs. While equivalent in principle to using the appropriate quadrature scheme, this becomes especially important and convenient when adaptive time-stepping algorithms are used to solve the ODEs (i.e. there is no fixed time grid). @@ -313,7 +310,7 @@ Setting initial conditions, we will assume a fixed :math:`i, e`, :math:`r=0`, :m prob = ODEProblem(F, x_0, (0.0, p.T), p) -The ``tspan`` of ``(0.0, p.T)`` determines that the :math:`t` used by the sovler. The time scale needs to be consistent with the arrival +The ``tspan`` of ``(0.0, p.T)`` determines that the :math:`t` used by the solver. The time scale needs to be consistent with the arrival rate of the transition probabilities (i.e. the :math:`\gamma, \sigma` were chosen based on daily data). The time period we investigate will be 550 days, or around 18 months: @@ -327,8 +324,7 @@ Let's run some experiments using this code. sol = solve(prob, Tsit5()) @show length(sol.t); -We see that the adaptive time-stepping used approximately 50 time-steps to solve this problem to the desires accuracy. Evaluating the solver at points outside of those time-steps uses the an interpolator consistent with the -solution to the ODE. +We see that the adaptive time-stepping used approximately 50 time-steps to solve this problem to the desires accuracy. Evaluating the solver at points outside of those time-steps uses an interpolator consistent with the solution to the ODE. See `here `__ for details on analyzing the solution, and `here `__ for plotting tools. The built-in plots for the solutions provide all of the `attributes `__ in `Plots.jl `__. @@ -477,7 +473,7 @@ and 75,000 agents already exposed to the virus and thus soon to be contagious. prob_late = ODEProblem(F, x_0, tspan, p_late) -Unlike the previous examples, the :math:`B(t)` functions have discontinuties which might occur. We can tell the adaptive time-stepping methods to ensure they include those points using ``tstops`` +Unlike the previous examples, the :math:`B(t)` functions have discontinuities which might occur. We can tell the adaptive time-stepping methods to ensure they include those points using ``tstops`` Let's calculate the paths: @@ -506,7 +502,7 @@ if a vaccine is found. Despite its richness, the model above is fully deterministic. The policy :math:`B(t)` could change over time, but only in predictable ways. -One source of randomness which would enter the model is considering the discretness of individuals. This topic, the connection to between SDEs and the Langevin equations typically used in the approximation of chemical reactions in well-mixed media are explored in our lecture on continuous time markov chains. +One source of randomness which would enter the model is considering the discretness of individuals. This topic, the connection to between SDEs and the Langevin equations typically used in the approximation of chemical reactions in well-mixed media are explored in our lecture on continuous time Markov chains. But rather than examining how granularity leads to aggregate fluctuations, we will concentrate on randomness that comes from aggregate changes in behavior or policy. @@ -515,18 +511,18 @@ Introduction to SDEs: Aggregate Shocks to Transmission Rates We will start by extending our model to include randomness in :math:`R(t)`, which makes it a system of Stochastic Differential Equations (SDEs). -Continous-Time Stochastic Processes +Continuous-Time Stochastic Processes ----------------------------------- -In continuous-time, there is an important distinction between randomness that leads to continuous paths vs. those which may have jumps (which are almost surely right-continous). The most tractable of these includes the theory of `Levy Processes `_. +In continuous-time, there is an important distinction between randomness that leads to continuous paths vs. those which may have jumps (which are almost surely right-continuous). The most tractable of these includes the theory of `Levy Processes `_. -**TBD:** Add definition of levy processes and the intuitive connection between statoionary increments and independence of increments. Also, is there intuition why additive processes in continuous time lead to additive SDEs/etc.? Might be useful. +**TBD:** Add definition of levy processes and the intuitive connection between stationary increments and independence of increments. Also, is there intuition why additive processes in continuous time lead to additive SDEs/etc.? Might be useful. Among the appealing features of Levy Processes is that they fit well into the sorts of Markov modeling techniques that economists tend to use in discrete time... Unlike in discrete-time, where a modeller has license to be creative, the rules of continuous-time stochastic processes are much more strict. In practice, there are only two types of Levy Processes that can be used without careful measure theory. -#. `Weiner Processes `__ (as known as Brownian Motion) which leads to a diffusion equations, and is the only continous-time Levy process with continuous paths +#. `Weiner Processes `__ (as known as Brownian Motion) which leads to a diffusion equations, and is the only continuous-time Levy process with continuous paths #. `Poisson Processes `__ with an arrival rate of jumps in the variable. Every other Levy Process you will typically work with is composed of those building blocks (e.g. a `Diffusion Process `__ such as Geometric Brownian Motion is a transformation of a Weiner process, and a `jump diffusion `__ is a diffusion process with a Poisson arrival of jumps). @@ -539,7 +535,7 @@ Shocks to Transmission Rates Consider that the effective transmission rate :math:`R(t)` could depend on degrees of randomness in behavior and implementation. For example, -* Misinformation on facebook spreading non-uniformly +* Misinformation on Facebook spreading non-uniformly * Large political rallies, elections, or protests * Deviations in the implementation and timing of lockdown policy within demographics, locations, or businesses within the system. * Aggregate shocks in opening/closing industries @@ -577,7 +573,7 @@ As the shock only effects :math:`dR`, which is the 5th equation, define the matr Since these are additive shocks, we will not need to modify the :math:`F` from our equation. -First create a new settings generator, and and then define a `SDEProblem `__ with Diagonal Noise. +First create a new settings generator, and then define a `SDEProblem `__ with Diagonal Noise. We solve the problem with the `SRA `__ algorithm (Adaptive strong order 1.5 methods for additive Ito and Stratonovich SDEs) @@ -599,9 +595,9 @@ We solve the problem with the `SRA \bar{i}`, then the available medical resources are exhuasted, leading to quadratically increased death rate. Second, the only social objective measure we can explore with the current framework is minimizing the total deaths. That ignores the possible policy tradeoffs between minimizing deaths and the impact on the general economy. -While a particular planner may decide that the only relevant welfare criteria is aggregate mortality, that leads to implausibly extreme policy (e.g. set :math:`B(t) = 0` forever). Hence, we will add a proxy for economic impact of Covid and the shotdown policy, summarized by :math:`u(t)` for excess covid-related "unemployment". A particular planner can then decide the weighting of the tradeoffs. +While a particular planner may decide that the only relevant welfare criteria is aggregate mortality, that leads to implausibly extreme policy (e.g. set :math:`B(t) = 0` forever). Hence, we will add a proxy for economic impact of COVID and the shutdown policy, summarized by :math:`u(t)` for excess COVID-related "unemployment". A particular planner can then decide the weighting of the tradeoffs. -The policy choice :math:`B(t)` is then made a markov-process rather than current exogeous and deterministic one. +The policy choice :math:`B(t)` is then made a Markov process rather than current exogenous and deterministic one. The inherent discretness of medical innovations and policy changes provides us an opportunity to explore the use of Poisson and jump diffusion. @@ -753,9 +749,9 @@ As vaccines are not instantaneously delivered to all, we can let :math:`\nu` be Unemployment and Lockdown Policy ---------------------------------- -As a simple summary of the economic and social distortions due to the covid, consider that an aggregate state :math:`u(t)` (indicating "excess" unemployment due to covid) increases as the :math:`B(t)`policy deviates from the natural :math:`\bar{R}` level, but can also increase due to flow of deaths from covid. +As a simple summary of the economic and social distortions due to the COVID, consider that an aggregate state :math:`u(t)` (indicating "excess" unemployment due to COVID) increases as the :math:`B(t)`policy deviates from the natural :math:`\bar{R}` level, but can also increase due to flow of deaths from COVID. -The second part of this is an important consideration: if the death rate is extremely high, then opening the economy may not help much as individuals are reluctant to return to normal economic activities. We will assume that weights on these are :math:`\mu_B` and :math:`mu_M` respetivly. +The second part of this is an important consideration: if the death rate is extremely high, then opening the economy may not help much as individuals are reluctant to return to normal economic activities. We will assume that weights on these are :math:`\mu_B` and :math:`mu_M` respectively. To represent the slow process of coming back to normal, we assume that the stock value :math:`u(t)` depreciates at rate :math:`\delta`. Put together gives the ODE, @@ -765,7 +761,7 @@ To represent the slow process of coming back to normal, we assume that the stock \end{aligned} :label: du -The policy choice :math:`B(t)` is markov, and will not consider implementation of forward looking behavior in this lecture. For a simple example, consider a policy driven by myopic political incentives and driven entirely by the death rates. If :math:`\pi_t > \bar{\pi}` then set :math:`B(t) = R_0` and leave it at :math:`B(t) = \bar{R}` otherwise. +The policy choice :math:`B(t)` is Markov, and will not consider implementation of forward looking behavior in this lecture. For a simple example, consider a policy driven by myopic political incentives and driven entirely by the death rates. If :math:`\pi_t > \bar{\pi}` then set :math:`B(t) = R_0` and leave it at :math:`B(t) = \bar{R}` otherwise. Without getting into the measure theory, this sort of bang-bang policy doesn't work as the process ceases to be a Levy Process (and much of the intuitive of considering measurability of stochastic processes and expected present discounted value fail). @@ -786,7 +782,7 @@ To our model, we have added in three new variables (:math:`u, V,` and :math:`B`) Stacking the :math:`u, V, B` at the end of the existing :math:`x` vector, we add or modify using :eq:`seir_system_vacc`, :eq:`Mode_nl`, and :eq:`du`. -Here, we will move from an "out of place" to an inplace differential equation. +Here, we will move from an "out of place" to an in-place differential equation. .. code-block:: julia @@ -812,7 +808,7 @@ Here, we will move from an "out of place" to an inplace differential equation. Note that the ``V, B`` terms do not have a drift as they are pure jump processes. -Next, we need to consider how the variance term of the diffusion changes. With the exception of the new brownian +Next, we need to consider how the variance term of the diffusion changes. With the exception of the new Brownian motion associated with the jump diffusion in :eq:`dmt`, everything else remains unchanged or zero. .. code-block:: julia @@ -909,7 +905,7 @@ Solving and simulating, plot(sol, vars = [8, 10]) TODO: Chris, lets make sure this is right, pick some parameters. -TODO: Lets show the same SIER decomposition we had before as a proportion, but with a vaccination arrival. +TODO: Lets show the same SEIR decomposition we had before as a proportion, but with a vaccination arrival. TODO: Maybe show a slice at time T = 550 or whatever of the distribution M(t) and u(t) as a histogram with 5, 50, and 95 percent confidence intervals? TODO: Then show how those two change as the myopic B(t) policy changes? TODO: Show how M(t) and u(t) change as the eta changes for a given policy? @@ -924,8 +920,8 @@ focus on the machine learning/etc. in the following lecture and not implement an Some thoughts on the SciML lecture which would be separate: 1. Bayesian estimation and model uncertainty. Have uncertainty on the parameter such as the :math:`m_0` parameter. Show uncertainty that comes after the estimation of the parameter. -2. Explain how bayesian priors are a form of regularization. With rare events like this, you will always have +2. Explain how Bayesian priors are a form of regularization. With rare events like this, you will always have more parameters than data. -3. Solve the time-0 optimal policy problem using Flux of choosing the :math:`B` policy in some form. For the objective, consider the uncertainty and the assymetry of being wrong. -4. If we wanted to add in a model element, I think that putting in an asymptomatic state would be very helpful for explaining the role of uncertainty in evaluating the counterfactuals. -5. Solve for the optimal markovian B(t) policy through simulation and using Flux. +3. Solve the time-0 optimal policy problem using Flux of choosing the :math:`B` policy in some form. For the objective, consider the uncertainty and the asymmetry of being wrong. +4. If we wanted to add in a model element, I think that putting in an asymptomatic state would be very helpful for explaining the role of uncertainty in evaluating the counter-factuals. +5. Solve for the optimal Markovian B(t) policy through simulation and using Flux. From fd35f58bfca0ec65ff5c5371c50ab0cc1d179d92 Mon Sep 17 00:00:00 2001 From: Jesse Perla Date: Tue, 30 Jun 2020 12:36:01 -0700 Subject: [PATCH 13/16] Moved to continous_time section --- conf.py | 3 +- source/rst/continuous_time/index.rst | 21 ++++++++++++++ .../seir_model_sde.rst | 28 +++++-------------- source/rst/index_postgrad.rst | 1 + source/rst/index_toc.rst | 1 + source/rst/tools_and_techniques/index.rst | 1 - 6 files changed, 32 insertions(+), 23 deletions(-) create mode 100644 source/rst/continuous_time/index.rst rename source/rst/{tools_and_techniques => continuous_time}/seir_model_sde.rst (97%) diff --git a/conf.py b/conf.py index b379c196..eb357f82 100644 --- a/conf.py +++ b/conf.py @@ -463,7 +463,8 @@ 'more_julia': [], 'multi_agent_models': [], 'time_series_models': [], - 'tools_and_techniques': [] + 'tools_and_techniques': [], + 'continuous_time': [] } jupyter_download_nb_execute=True diff --git a/source/rst/continuous_time/index.rst b/source/rst/continuous_time/index.rst new file mode 100644 index 00000000..cba0b06d --- /dev/null +++ b/source/rst/continuous_time/index.rst @@ -0,0 +1,21 @@ +.. _continuous_time: + +.. include:: /_static/includes/header.raw + +*************************************** +Modeling in Continuous Time +*************************************** + +This section of the course contains foundational mathematical and computational tools for working with continuous time deterministic and stochastic models. + +.. only:: html + + Lectures + ******** + + + +.. toctree:: + :maxdepth: 2 + + seir_model_sde diff --git a/source/rst/tools_and_techniques/seir_model_sde.rst b/source/rst/continuous_time/seir_model_sde.rst similarity index 97% rename from source/rst/tools_and_techniques/seir_model_sde.rst rename to source/rst/continuous_time/seir_model_sde.rst index aa4436ed..4e5439d8 100644 --- a/source/rst/tools_and_techniques/seir_model_sde.rst +++ b/source/rst/continuous_time/seir_model_sde.rst @@ -1,10 +1,12 @@ +.. _seir_model_sde: + .. include:: /_static/includes/header.raw .. highlight:: julia -****************************************************************** +******************************************************************* :index:`Modeling COVID 19 with (Stochastic) Differential Equations` -****************************************************************** +******************************************************************* .. contents:: :depth: 2 @@ -66,7 +68,7 @@ others covered in previous lectures .. code-block:: julia using OrdinaryDiffEq, StochasticDiffEq, DiffEqJump - using Parameters, StaticArrays, Plots + using Parameters, StaticArrays, Plots The SEIR Model @@ -170,7 +172,7 @@ Given this system, we choose an initial condition and a timespan, and create a ` x_0 = [s_0, e_0, i_0, r_0] # initial condition tspan = (0.0, 350.0) # ≈ 350 days - prob = ODEProblem(F_simple, x_0, tspan) + prob = ODEProblem(F_simple, x_0, tspan, nothing) With this, we can choose an ODE algorithm (e.g. a good default for non-stiff ODEs of this sort might be ``Tsit5()``, which is the Tsitouras 5/4 Runge-Kutta method). @@ -718,7 +720,7 @@ In previous versions of model, :eq:`Mode` had deaths as the integration of the d .. math:: \begin{aligned} \pi(m, i) &:= m + \psi \max(0, i - \bar{i})\\ - d M_t &= \pi(m_t, i_t)\gamma i_t dt + d M_t &= \pi(m_t, i_t)\gamma i_t dt\\ \end{aligned} :label: Mode_nl @@ -911,19 +913,3 @@ Solving and simulating, TODO: Chris, lets make sure this is right, pick some parameters, and solve it? TODO: What sorts of experiments afte rthis? - - -ROUGH NOTES ON THE SCIML LECTURE -===================================== - -The hope is that we have all or almost all of the model elements developed here, so that we can -focus on the machine learning/etc. in the following lecture and not implement any new model features (i.e. just shut things down as required). - - -Some thoughts on the SciML lecture which would be separate: - -1. Bayesian estimation and model uncertainty. Have uncertainty on the parameter such as the :math:`m_0` parameter. Show uncertainty that comes after the estimation of the parameter. -2. Explain how bayesian priors are a form of regularization. With rare events like this, you will always have -more parameters than data. -3. Solve the time-0 optimal policy problem using Flux of choosing the :math:`B` policy in some form. For the objective, consider the uncertainty and the assymetry of being wrong. -4. If we wanted to add in a model element, I think that putting in an asymptomatic state would be very helpful for explaining the role of uncertainty in evaluating the counterfactuals. \ No newline at end of file diff --git a/source/rst/index_postgrad.rst b/source/rst/index_postgrad.rst index 487e5007..ec8ade6b 100644 --- a/source/rst/index_postgrad.rst +++ b/source/rst/index_postgrad.rst @@ -26,6 +26,7 @@ more_julia/index tools_and_techniques/index_grad dynamic_programming/index_grad + continuous_time/index multi_agent_models/index_grad time_series_models/index_grad zreferences diff --git a/source/rst/index_toc.rst b/source/rst/index_toc.rst index 09e33d17..96b42393 100644 --- a/source/rst/index_toc.rst +++ b/source/rst/index_toc.rst @@ -27,6 +27,7 @@ more_julia/index tools_and_techniques/index dynamic_programming/index + continuous_time/index multi_agent_models/index time_series_models/index dynamic_programming_squared/index diff --git a/source/rst/tools_and_techniques/index.rst b/source/rst/tools_and_techniques/index.rst index 7019bd82..cefec116 100644 --- a/source/rst/tools_and_techniques/index.rst +++ b/source/rst/tools_and_techniques/index.rst @@ -20,7 +20,6 @@ tools and techniques. :maxdepth: 2 linear_algebra - seir_model_sde orth_proj lln_clt linear_models From 446326a95e399537df5c6c34a39f5822fea07f18 Mon Sep 17 00:00:00 2001 From: Jesse Perla Date: Tue, 30 Jun 2020 12:57:02 -0700 Subject: [PATCH 14/16] Placeholder for next lecture --- source/rst/continuous_time/index.rst | 1 + source/rst/continuous_time/seir_sciml.rst | 67 +++++++++++++++++++++++ 2 files changed, 68 insertions(+) create mode 100644 source/rst/continuous_time/seir_sciml.rst diff --git a/source/rst/continuous_time/index.rst b/source/rst/continuous_time/index.rst index cba0b06d..4ffdd5de 100644 --- a/source/rst/continuous_time/index.rst +++ b/source/rst/continuous_time/index.rst @@ -19,3 +19,4 @@ This section of the course contains foundational mathematical and computational :maxdepth: 2 seir_model_sde + seir_sciml \ No newline at end of file diff --git a/source/rst/continuous_time/seir_sciml.rst b/source/rst/continuous_time/seir_sciml.rst new file mode 100644 index 00000000..e47ac68c --- /dev/null +++ b/source/rst/continuous_time/seir_sciml.rst @@ -0,0 +1,67 @@ +.. _seir_sciml: + +.. include:: /_static/includes/header.raw + +.. highlight:: julia + +******************************************************************* +:index:`Parameter and Model Uncertainty in the COVID 19 Models` +******************************************************************* + +.. contents:: :depth: 2 + +Overview +============= + +Coauthored with Chris Rackauckas + + +Using the model developed in the :doc:`Modeling COVID 19 with (Stochastic) Differential Equations <../continuous_time/seir_model_sde>`, we explore further topics such as: + +* Finding the derivatives of solutions to ODEs/SDEs +* Parameter uncertainty +* Bayesian estimation of the differential equations +* Minimizing loss functions with solutions to ODE/SDEs + + +Setup +------------------ + +.. literalinclude:: /_static/includes/deps_generic.jl + :class: hide-output + +.. code-block:: julia + + using LinearAlgebra, Statistics, Random, SparseArrays + +.. code-block:: julia + :class: Test + + using Test # Put this before any code in the lecture. + +In addition, we will be exploring packages within the `SciML ecosystem `__ and +others covered in previous lectures + +.. code-block:: julia + + using OrdinaryDiffEq, StochasticDiffEq, DiffEqJump, DiffEqFlux, DiffEqBayes + using Parameters, StaticArrays, Plots + + +ROUGH NOTES +============== + +SOME OLDER NOTES: + +The hope is that we have all or almost all of the model elements developed here, so that we can +focus on the machine learning/etc. in the following lecture and not implement any new model features (i.e. just shut things down as required). + + +Some thoughts on the SciML lecture which would be separate: + +1. Bayesian estimation and model uncertainty. Have uncertainty on the parameter such as the :math:`m_0` parameter. Show uncertainty that comes after the estimation of the parameter. +2. Explain how Bayesian priors are a form of regularization. With rare events like this, you will always have +more parameters than data. +3. Solve the time-0 optimal policy problem using Flux of choosing the :math:`B` policy in some form. For the objective, consider the uncertainty and the asymmetry of being wrong. +4. If we wanted to add in a model element, I think that putting in an asymptomatic state would be very helpful for explaining the role of uncertainty in evaluating the counter-factuals. +5. Solve for the optimal Markovian B(t) policy through simulation and using Flux. \ No newline at end of file From 990b19cb4a4041cb43934f4211ce53f0bc8cde8d Mon Sep 17 00:00:00 2001 From: Jesse Perla Date: Tue, 30 Jun 2020 13:16:35 -0700 Subject: [PATCH 15/16] Placeholder for CTMC --- .../continuous_time_markov_chains.rst | 292 ++++++++++++++++++ source/rst/continuous_time/index.rst | 3 +- .../numerical_linear_algebra.rst | 229 -------------- 3 files changed, 294 insertions(+), 230 deletions(-) create mode 100644 source/rst/continuous_time/continuous_time_markov_chains.rst diff --git a/source/rst/continuous_time/continuous_time_markov_chains.rst b/source/rst/continuous_time/continuous_time_markov_chains.rst new file mode 100644 index 00000000..0851d4d7 --- /dev/null +++ b/source/rst/continuous_time/continuous_time_markov_chains.rst @@ -0,0 +1,292 @@ +.. _seir_sciml: + +.. include:: /_static/includes/header.raw + +.. highlight:: julia + +******************************************************************* +:index:`Continuous Time Markov Chains with COVID 19 Applications` +******************************************************************* + +.. contents:: :depth: 2 + +Overview +============= + +Coauthored with Chris Rackauckas + +This lecture develops the theory of using continuous time markov chains, in the same spirit as :doc:`Discrete Time Markov Chains <../tools_and_techniques/finite_markov>` and with some of the same applications as that of :doc:`Discrete State Dynamic Programming <../dynamic_programming/discrete_dp>` + +Here, we will also explore where the ODE approximation of :doc:`Modeling COVID 19 with (Stochastic) Differential Equations <../continuous_time/seir_model_sde>`, comes from + + +Setup +------------------ + +.. literalinclude:: /_static/includes/deps_generic.jl + :class: hide-output + +.. code-block:: julia + + using LinearAlgebra, Statistics, Random, SparseArrays + +.. code-block:: julia + :class: Test + + using Test # Put this before any code in the lecture. + +In addition, we will be exploring packages within the `SciML ecosystem `__ and +others covered in previous lectures + +.. code-block:: julia + + using OrdinaryDiffEq, StochasticDiffEq, DiffEqJump, DiffEqFlux, DiffEqBayes + using Parameters, StaticArrays, Plots + + +ROUGH NOTES +============== + +Notes from John:- +- matrix exponentials +- Q matrices and their relationship to Markov matrices (via matrix exponential, Kolmogorov fwd and backward equations) +- exponential clock interpretation of Q matrices +- ergodicity +- applications (inventory dynamics?) + +Lets of ideas from Chris in https://github.com/QuantEcon/lecture-source-jl/issues/859#issuecomment-632998856 + +The goal for that might me to building up from a Poisson process to show how the SDE approximation works Gillespie's SSA, and then show how in the limit of a +large number of particles we get the ODE approximation that economists are familiar with. + +I also think that getting more familiar with the SDE approximation of a discrete number of agents could have other applciations in economics (e.g. search) + + +Continuous-Time Markov Chains (CTMCs) +===================================== + +In the lecture on :doc:`discrete-time Markov chains `, we saw that the transition probability +between state :math:`x` and state :math:`y` was summarized by the matrix :math:`P(x, y) := \mathbb P \{ X_{t+1} = y \,|\, X_t = x \}`. + +As a brief introduction to continuous time processes, consider the same state space as in the discrete +case: :math:`S` is a finite set with :math:`n` elements :math:`\{x_1, \ldots, x_n\}`. + +A **Markov chain** :math:`\{X_t\}` on :math:`S` is a sequence of random variables on :math:`S` that have the **Markov property**. + +In continuous time, the `Markov Property `_ is more complicated, but intuitively is +the same as the discrete-time case. + +That is, knowing the current state is enough to know probabilities for future states. Or, for realizations :math:`x(\tau)\in S, \tau \leq t`, + +.. math:: + + \mathbb P \{ X(t+s) = y \,|\, X(t) = x, X(\tau) = x(\tau) \text{ for } 0 \leq \tau \leq t \} = \mathbb P \{ X(t+s) = y \,|\, X(t) = x\} + + +Heuristically, consider a time period :math:`t` and a small step forward, :math:`\Delta`. Then the probability to transition from state :math:`i` to +state :math:`j` is + +.. math:: + + \mathbb P \{ X(t + \Delta) = j \,|\, X(t) \} = \begin{cases} q_{ij} \Delta + o(\Delta) & i \neq j\\ + 1 + q_{ii} \Delta + o(\Delta) & i = j \end{cases} + +where the :math:`q_{ij}` are "intensity" parameters governing the transition rate, and :math:`o(\Delta)` is `little-o notation `_. That is, :math:`\lim_{\Delta\to 0} o(\Delta)/\Delta = 0`. + +Just as in the discrete case, we can summarize these parameters by an :math:`N \times N` matrix, :math:`Q \in R^{N\times N}`. + +Recall that in the discrete case every element is weakly positive and every row must sum to one. With continuous time, however, the rows of :math:`Q` sum to zero, where the diagonal contains the negative value of jumping out of the current state. That is, + +- :math:`q_{ij} \geq 0` for :math:`i \neq j` +- :math:`q_{ii} \leq 0` +- :math:`\sum_{j} q_{ij} = 0` + +The :math:`Q` matrix is called the intensity matrix, or the infinitesimal generator of the Markov chain. For example, + +.. math:: + + Q = \begin{bmatrix} -0.1 & 0.1 & 0 & 0 & 0 & 0\\ + 0.1 &-0.2 & 0.1 & 0 & 0 & 0\\ + 0 & 0.1 & -0.2 & 0.1 & 0 & 0\\ + 0 & 0 & 0.1 & -0.2 & 0.1 & 0\\ + 0 & 0 & 0 & 0.1 & -0.2 & 0.1\\ + 0 & 0 & 0 & 0 & 0.1 & -0.1\\ + \end{bmatrix} + +In the above example, transitions occur only between adjacent states with the same intensity (except for a ``bouncing back'' of the bottom and top states). + +Implementing the :math:`Q` using its tridiagonal structure + +.. code-block:: julia + + using LinearAlgebra + α = 0.1 + N = 6 + Q = Tridiagonal(fill(α, N-1), [-α; fill(-2α, N-2); -α], fill(α, N-1)) + +Here we can use ``Tridiagonal`` to exploit the structure of the problem. + +Consider a simple payoff vector :math:`r` associated with each state, and a discount rate :math:`ρ`. Then we can solve for +the expected present discounted value in a way similar to the discrete-time case. + +.. math:: + + \rho v = r + Q v + +or rearranging slightly, solving the linear system + +.. math:: + + (\rho I - Q) v = r + +For our example, exploiting the tridiagonal structure, + +.. code-block:: julia + + r = range(0.0, 10.0, length=N) + ρ = 0.05 + + A = ρ * I - Q + +Note that this :math:`A` matrix is maintaining the tridiagonal structure of the problem, which leads to an efficient solution to the +linear problem. + +.. code-block:: julia + + v = A \ r + +The :math:`Q` is also used to calculate the evolution of the Markov chain, in direct analogy to the :math:`ψ_{t+k} = ψ_t P^k` evolution with the transition matrix :math:`P` of the discrete case. + +In the continuous case, this becomes the system of linear differential equations + +.. math:: + + \dot{ψ}(t) = Q(t)^T ψ(t) + +given the initial condition :math:`\psi(0)` and where the :math:`Q(t)` intensity matrix is allowed to vary with time. In the simplest case of a constant :math:`Q` matrix, this is a simple constant-coefficient system of linear ODEs with coefficients :math:`Q^T`. + +If a stationary equilibrium exists, note that :math:`\dot{ψ}(t) = 0`, and the stationary solution :math:`ψ^{*}` needs to satisfy + +.. math:: + + 0 = Q^T ψ^{*} + + +Notice that this is of the form :math:`0 ψ^{*} = Q^T ψ^{*}` and hence is equivalent to finding the eigenvector associated with the :math:`\lambda = 0` eigenvalue of :math:`Q^T`. + +With our example, we can calculate all of the eigenvalues and eigenvectors + +.. code-block:: julia + + λ, vecs = eigen(Array(Q')) + +Indeed, there is a :math:`\lambda = 0` eigenvalue, which is associated with the last column in the eigenvector. To turn that into a probability, +we need to normalize it. + +.. code-block:: julia + + vecs[:,N] ./ sum(vecs[:,N]) + +Multiple Dimensions +-------------------- + +A frequent case in discretized models is dealing with Markov chains with multiple "spatial" dimensions (e.g., wealth and income). + +After discretizing a process to create a Markov chain, you can always take the Cartesian product of the set of states in order to +enumerate it as a single state variable. + +To see this, consider states :math:`i` and :math:`j` governed by infinitesimal generators :math:`Q` and :math:`A`. + +.. code-block:: julia + + function markov_chain_product(Q, A) + M = size(Q, 1) + N = size(A, 1) + Q = sparse(Q) + Qs = blockdiag(fill(Q, N)...) # create diagonal blocks of every operator + As = kron(A, sparse(I(M))) + return As + Qs + end + + α = 0.1 + N = 4 + Q = Tridiagonal(fill(α, N-1), [-α; fill(-2α, N-2); -α], fill(α, N-1)) + A = sparse([-0.1 0.1 + 0.2 -0.2]) + M = size(A,1) + L = markov_chain_product(Q, A) + L |> Matrix # display as a dense matrix + +This provides the combined Markov chain for the :math:`(i,j)` process. To see the sparsity pattern, + +.. code-block:: julia + + using Plots + gr(fmt = :png); + spy(L, markersize = 10) + +To calculate a simple dynamic valuation, consider whether the payoff of being in state :math:`(i,j)` is :math:`r_{ij} = i + 2j` + +.. code-block:: julia + + r = [i + 2.0j for i in 1:N, j in 1:M] + r = vec(r) # vectorize it since stacked in same order + +Solving the equation :math:`\rho v = r + L v` + +.. code-block:: julia + + ρ = 0.05 + v = (ρ * I - L) \ r + reshape(v, N, M) + +The ``reshape`` helps to rearrange it back to being two-dimensional. + + +To find the stationary distribution, we calculate the eigenvalue and choose the eigenvector associated with :math:`\lambda=0` . In this +case, we can verify that it is the last one. + +.. code-block:: julia + + L_eig = eigen(Matrix(L')) + @assert norm(L_eig.values[end]) < 1E-10 + + ψ = L_eig.vectors[:,end] + ψ = ψ / sum(ψ) + + +Reshape this to be two-dimensional if it is helpful for visualization. + +.. code-block:: julia + + reshape(ψ, N, size(A,1)) + +Irreducibility +-------------- + +As with the discrete-time Markov chains, a key question is whether CTMCs are reducible, i.e., whether states communicate. The problem +is isomorphic to determining whether the directed graph of the Markov chain is `strongly connected `_. + +.. code-block:: julia + + using LightGraphs + α = 0.1 + N = 6 + Q = Tridiagonal(fill(α, N-1), [-α; fill(-2α, N-2); -α], fill(α, N-1)) + +We can verify that it is possible to move between every pair of states in a finite number of steps with + +.. code-block:: julia + + Q_graph = DiGraph(Q) + @show is_strongly_connected(Q_graph); # i.e., can follow directional edges to get to every state + +Alternatively, as an example of a reducible Markov chain where states :math:`1` and :math:`2` cannot jump to state :math:`3`. + +.. code-block:: julia + + Q = [-0.2 0.2 0 + 0.2 -0.2 0 + 0.2 0.6 -0.8] + Q_graph = DiGraph(Q) + @show is_strongly_connected(Q_graph); diff --git a/source/rst/continuous_time/index.rst b/source/rst/continuous_time/index.rst index 4ffdd5de..3fa66085 100644 --- a/source/rst/continuous_time/index.rst +++ b/source/rst/continuous_time/index.rst @@ -19,4 +19,5 @@ This section of the course contains foundational mathematical and computational :maxdepth: 2 seir_model_sde - seir_sciml \ No newline at end of file + seir_sciml + continuous_time_markov_chains \ No newline at end of file diff --git a/source/rst/tools_and_techniques/numerical_linear_algebra.rst b/source/rst/tools_and_techniques/numerical_linear_algebra.rst index 9a2ab59c..382b2a91 100644 --- a/source/rst/tools_and_techniques/numerical_linear_algebra.rst +++ b/source/rst/tools_and_techniques/numerical_linear_algebra.rst @@ -580,235 +580,6 @@ To see this, Keep in mind that a real matrix may have complex eigenvalues and eigenvectors, so if you attempt to check ``Q * Λ * inv(Q) - A`` - even for a positive-definite matrix - it may not be a real number due to numerical inaccuracy. -Continuous-Time Markov Chains (CTMCs) -==================================== - -In the previous lecture on :doc:`discrete-time Markov chains `, we saw that the transition probability -between state :math:`x` and state :math:`y` was summarized by the matrix :math:`P(x, y) := \mathbb P \{ X_{t+1} = y \,|\, X_t = x \}`. - -As a brief introduction to continuous time processes, consider the same state space as in the discrete -case: :math:`S` is a finite set with :math:`n` elements :math:`\{x_1, \ldots, x_n\}`. - -A **Markov chain** :math:`\{X_t\}` on :math:`S` is a sequence of random variables on :math:`S` that have the **Markov property**. - -In continuous time, the `Markov Property `_ is more complicated, but intuitively is -the same as the discrete-time case. - -That is, knowing the current state is enough to know probabilities for future states. Or, for realizations :math:`x(\tau)\in S, \tau \leq t`, - -.. math:: - - \mathbb P \{ X(t+s) = y \,|\, X(t) = x, X(\tau) = x(\tau) \text{ for } 0 \leq \tau \leq t \} = \mathbb P \{ X(t+s) = y \,|\, X(t) = x\} - - -Heuristically, consider a time period :math:`t` and a small step forward, :math:`\Delta`. Then the probability to transition from state :math:`i` to -state :math:`j` is - -.. math:: - - \mathbb P \{ X(t + \Delta) = j \,|\, X(t) \} = \begin{cases} q_{ij} \Delta + o(\Delta) & i \neq j\\ - 1 + q_{ii} \Delta + o(\Delta) & i = j \end{cases} - -where the :math:`q_{ij}` are "intensity" parameters governing the transition rate, and :math:`o(\Delta)` is `little-o notation `_. That is, :math:`\lim_{\Delta\to 0} o(\Delta)/\Delta = 0`. - -Just as in the discrete case, we can summarize these parameters by an :math:`N \times N` matrix, :math:`Q \in R^{N\times N}`. - -Recall that in the discrete case every element is weakly positive and every row must sum to one. With continuous time, however, the rows of :math:`Q` sum to zero, where the diagonal contains the negative value of jumping out of the current state. That is, - -- :math:`q_{ij} \geq 0` for :math:`i \neq j` -- :math:`q_{ii} \leq 0` -- :math:`\sum_{j} q_{ij} = 0` - -The :math:`Q` matrix is called the intensity matrix, or the infinitesimal generator of the Markov chain. For example, - -.. math:: - - Q = \begin{bmatrix} -0.1 & 0.1 & 0 & 0 & 0 & 0\\ - 0.1 &-0.2 & 0.1 & 0 & 0 & 0\\ - 0 & 0.1 & -0.2 & 0.1 & 0 & 0\\ - 0 & 0 & 0.1 & -0.2 & 0.1 & 0\\ - 0 & 0 & 0 & 0.1 & -0.2 & 0.1\\ - 0 & 0 & 0 & 0 & 0.1 & -0.1\\ - \end{bmatrix} - -In the above example, transitions occur only between adjacent states with the same intensity (except for a ``bouncing back'' of the bottom and top states). - -Implementing the :math:`Q` using its tridiagonal structure - -.. code-block:: julia - - using LinearAlgebra - α = 0.1 - N = 6 - Q = Tridiagonal(fill(α, N-1), [-α; fill(-2α, N-2); -α], fill(α, N-1)) - -Here we can use ``Tridiagonal`` to exploit the structure of the problem. - -Consider a simple payoff vector :math:`r` associated with each state, and a discount rate :math:`ρ`. Then we can solve for -the expected present discounted value in a way similar to the discrete-time case. - -.. math:: - - \rho v = r + Q v - -or rearranging slightly, solving the linear system - -.. math:: - - (\rho I - Q) v = r - -For our example, exploiting the tridiagonal structure, - -.. code-block:: julia - - r = range(0.0, 10.0, length=N) - ρ = 0.05 - - A = ρ * I - Q - -Note that this :math:`A` matrix is maintaining the tridiagonal structure of the problem, which leads to an efficient solution to the -linear problem. - -.. code-block:: julia - - v = A \ r - -The :math:`Q` is also used to calculate the evolution of the Markov chain, in direct analogy to the :math:`ψ_{t+k} = ψ_t P^k` evolution with the transition matrix :math:`P` of the discrete case. - -In the continuous case, this becomes the system of linear differential equations - -.. math:: - - \dot{ψ}(t) = Q(t)^T ψ(t) - -given the initial condition :math:`\psi(0)` and where the :math:`Q(t)` intensity matrix is allowed to vary with time. In the simplest case of a constant :math:`Q` matrix, this is a simple constant-coefficient system of linear ODEs with coefficients :math:`Q^T`. - -If a stationary equilibrium exists, note that :math:`\dot{ψ}(t) = 0`, and the stationary solution :math:`ψ^{*}` needs to satisfy - -.. math:: - - 0 = Q^T ψ^{*} - - -Notice that this is of the form :math:`0 ψ^{*} = Q^T ψ^{*}` and hence is equivalent to finding the eigenvector associated with the :math:`\lambda = 0` eigenvalue of :math:`Q^T`. - -With our example, we can calculate all of the eigenvalues and eigenvectors - -.. code-block:: julia - - λ, vecs = eigen(Array(Q')) - -Indeed, there is a :math:`\lambda = 0` eigenvalue, which is associated with the last column in the eigenvector. To turn that into a probability, -we need to normalize it. - -.. code-block:: julia - - vecs[:,N] ./ sum(vecs[:,N]) - -Multiple Dimensions --------------------- - -A frequent case in discretized models is dealing with Markov chains with multiple "spatial" dimensions (e.g., wealth and income). - -After discretizing a process to create a Markov chain, you can always take the Cartesian product of the set of states in order to -enumerate it as a single state variable. - -To see this, consider states :math:`i` and :math:`j` governed by infinitesimal generators :math:`Q` and :math:`A`. - -.. code-block:: julia - - function markov_chain_product(Q, A) - M = size(Q, 1) - N = size(A, 1) - Q = sparse(Q) - Qs = blockdiag(fill(Q, N)...) # create diagonal blocks of every operator - As = kron(A, sparse(I(M))) - return As + Qs - end - - α = 0.1 - N = 4 - Q = Tridiagonal(fill(α, N-1), [-α; fill(-2α, N-2); -α], fill(α, N-1)) - A = sparse([-0.1 0.1 - 0.2 -0.2]) - M = size(A,1) - L = markov_chain_product(Q, A) - L |> Matrix # display as a dense matrix - -This provides the combined Markov chain for the :math:`(i,j)` process. To see the sparsity pattern, - -.. code-block:: julia - - using Plots - gr(fmt = :png); - spy(L, markersize = 10) - -To calculate a simple dynamic valuation, consider whether the payoff of being in state :math:`(i,j)` is :math:`r_{ij} = i + 2j` - -.. code-block:: julia - - r = [i + 2.0j for i in 1:N, j in 1:M] - r = vec(r) # vectorize it since stacked in same order - -Solving the equation :math:`\rho v = r + L v` - -.. code-block:: julia - - ρ = 0.05 - v = (ρ * I - L) \ r - reshape(v, N, M) - -The ``reshape`` helps to rearrange it back to being two-dimensional. - - -To find the stationary distribution, we calculate the eigenvalue and choose the eigenvector associated with :math:`\lambda=0` . In this -case, we can verify that it is the last one. - -.. code-block:: julia - - L_eig = eigen(Matrix(L')) - @assert norm(L_eig.values[end]) < 1E-10 - - ψ = L_eig.vectors[:,end] - ψ = ψ / sum(ψ) - - -Reshape this to be two-dimensional if it is helpful for visualization. - -.. code-block:: julia - - reshape(ψ, N, size(A,1)) - -Irreducibility --------------- - -As with the discrete-time Markov chains, a key question is whether CTMCs are reducible, i.e., whether states communicate. The problem -is isomorphic to determining whether the directed graph of the Markov chain is `strongly connected `_. - -.. code-block:: julia - - using LightGraphs - α = 0.1 - N = 6 - Q = Tridiagonal(fill(α, N-1), [-α; fill(-2α, N-2); -α], fill(α, N-1)) - -We can verify that it is possible to move between every pair of states in a finite number of steps with - -.. code-block:: julia - - Q_graph = DiGraph(Q) - @show is_strongly_connected(Q_graph); # i.e., can follow directional edges to get to every state - -Alternatively, as an example of a reducible Markov chain where states :math:`1` and :math:`2` cannot jump to state :math:`3`. - -.. code-block:: julia - - Q = [-0.2 0.2 0 - 0.2 -0.2 0 - 0.2 0.6 -0.8] - Q_graph = DiGraph(Q) - @show is_strongly_connected(Q_graph); - Banded Matrices =============== From dd992714572f44ed45b52df44739b25ba7918e96 Mon Sep 17 00:00:00 2001 From: Jesse Perla Date: Tue, 30 Jun 2020 13:17:03 -0700 Subject: [PATCH 16/16] renamed reference --- source/rst/continuous_time/continuous_time_markov_chains.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/source/rst/continuous_time/continuous_time_markov_chains.rst b/source/rst/continuous_time/continuous_time_markov_chains.rst index 0851d4d7..68e66d6d 100644 --- a/source/rst/continuous_time/continuous_time_markov_chains.rst +++ b/source/rst/continuous_time/continuous_time_markov_chains.rst @@ -1,4 +1,4 @@ -.. _seir_sciml: +.. _ctmc: .. include:: /_static/includes/header.raw