Skip to content
This repository has been archived by the owner on Sep 21, 2021. It is now read-only.

WIP: Covid sde #867

Merged
merged 22 commits into from
Jul 22, 2020
Merged

WIP: Covid sde #867

merged 22 commits into from
Jul 22, 2020

Conversation

jlperla
Copy link
Member

@jlperla jlperla commented Jun 27, 2020

This is a work in progress, but it is time to start examining it. There has been virtually no copy editing/etc.

@jstac @ChrisRackauckas @thomassargent30

@jstac
Copy link
Contributor

jstac commented Jun 28, 2020

Implements some of the discussion in #859

@jstac
Copy link
Contributor

jstac commented Jun 29, 2020

Very interesting. I added only minor corrections. The stochastic time series map better to the data than the deterministic ones. So there must be significant aggregate shocks,

covid

I'll review more carefully when I have time but everything is clear and nicely explained.

@ChrisRackauckas
Copy link
Member

I plan to go through this in detail, but let's say this weekend. I had a paper resubmission come up.

@jlperla
Copy link
Member Author

jlperla commented Jun 30, 2020

@jstac The only real chance on this PR from the one you modified is that I move this into its own section on continuous-time - where we can build up the other lecture slowly, and I moved the raw notes at the bottom intended for the 2nd lecture into #868

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`
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Not sure what v is supposed to be referring to.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sorry, not on computer, but it might be old notation for the death probability in the python ? Maybe I changed it to m(t) for mortality?

* 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 the :math:`\sqrt{R}` ensures that the process can never go negative since the variance converges to zero as R goes to zero.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Not sure it actually ends up as a provably positive process. At least the argument and intuition aren't necessarily correct. u' = -sqrt(abs(u)) is an example where the rate goes to zero but this process hits exactly zero in finite time, at which point its solution is non-unique but it's either 0 or it goes negative. I think the fact that the diffusion length scales by another sqrt might make this work out, but it's not exactly trivial.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I added in the scaling by the sqrt when we discussed this previously, and you had suggested it was a good way to keep things positive.

But if I have learned anythign over the years, it is that continuous time stochastic processes are is often counterintutive, and that the details matter.

For this one, I think the finance crew calls it the CIR process: https://en.wikipedia.org/wiki/Cox%E2%80%93Ingersoll%E2%80%93Ross_model

It says there that the scaling avoids the possibility of negative interest rates for all a and b, and an absorbing state of zero is precluded as long as the drift is strong enough relative to the variance. I think this comes out of the ergodic distribution being Gamma for the right combination of parameters: https://en.wikipedia.org/wiki/Cox%E2%80%93Ingersoll%E2%80%93Ross_model#Distribution

With that said, I don't think it means that the SDE simulators would necessarily always keep things positive? Which might be what I ran into later in the code (which you might see in https://github.com/QuantEcon/lecture-source-jl/blob/459e806e3acaa3371013112ed10eae7c5f38ce4a/source/rst/continuous_time/seir_model_sde.rst#implementing-the-complete-system where I fudge the variance to avoid a sqrt of a negative. Probably better ways to do this.

I think a little adhoc on the SDE side for this is OK as long as we are clea rwhat is going on.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It says there that the scaling avoids the possibility of negative interest rates for all a and b, and an absorbing state of zero is precluded as long as the drift is strong enough relative to the variance.

Alright, cool. I think we should mention that the explanation given is really just a heuristic, and to know it's going to be positive requires more details.

With that said, I don't think it means that the SDE simulators would necessarily always keep things positive?

Yeah, being that close to zero can cause numerical difficulties

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Comment added


First create a new settings generator, and then define a `SDEProblem <https://docs.sciml.ai/stable/tutorials/sde_example/#Example-2:-Systems-of-SDEs-with-Diagonal-Noise-1>`__ with Diagonal Noise.

We solve the problem with the `SRA <https://docs.sciml.ai/stable/solvers/sde_solve/#Full-List-of-Methods-1>`__ algorithm (Adaptive strong order 1.5 methods for additive Ito and Stratonovich SDEs)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this is not an additive noise SDE

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I see. This is because of the scaling of the variance? If so, maybe putting in a note on what the general goto algorithms are for additive vs. other noise couldn't hurt. Economists (and definetely me!) are utterly ignorant of this stuff. Everything is homegrown Euler-Maruyama at best.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, additive noise is what dW is only multiplied by constants. Here it's just diagonal noise.


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 (and related topics of irreducibility) in the lecture on continuous-time Markov chains, as well as the limitations of these approximations when the discreteness becomes essential (e.g. continuum approximations are incapable of modeling extinguishing of an outbreak).
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Perhaps this is not true once aggregate shocks are added, if, say, the paths admit jumps. Might be worth qualifying?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I spent a bunch of time looking into this and clearing this up. This turns out to be a very subtle point, but what I had written before was correct. It would only be true if R_0 > 0 for this sort of model, I think.

The latest updates though have changed a lot so this won't be an issue.

* :math:`\gamma` is called the *recovery rate* (the rate at which infected people recover or die).


In addition, the transmission rate can be re-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.
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@jstac I took this from the python lecture, but I am not sure it is correct in epidemilogy lingo. This might be the R_0 instead (i.e. the "Basic reproduction number") instead. See https://en.wikipedia.org/wiki/Basic_reproduction_number and other stuff.

But understanding from this is that the base reproduction number is that where we make everyone susceptible, while the effective reproduction number takes into account the current equilibrium?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I believe in your implementation you may talk about R0 instead? https://python.quantecon.org/sir_model.html#Implementation

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Finally, for the "effective reproduction rate", which might actually be useful to plot as well, I think it might be simply R_e(t) = R_0 * s(t) i.e. the derivative of the R_0 * I/N * S flow wrt I. Or the marginal number of new infections for an increase in the number of infected?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@jstac I implemented these sorts of changes in the latest version, and there are a significant number of changes along these topics.

@jlperla
Copy link
Member Author

jlperla commented Jul 22, 2020

@jstac @ChrisRackauckas @thomassargent30

I just pushed a major set of changes

The summary of changes is:

  • I moved the jump diffusion stuff to a do later. There was enough content for the SDE stuff that it would be overwhelming otherwise.
  • Also, I split the deterministic and SDE versions into separate lectures for the same reason. It will also mean we don't have to fully nest the determinstic and SDE models.
  • One benefit of this is that without the jump-diffusions, it will be much easier to push these lectures out since they don't require deploying new packages.
  • The SDE model now uses an SIRD model. This lets us teach people about a different type of covid approach (and it is referenced in a paper by jesus and chad), and simplifies those equations a bit.
  • There is a lot of new material on the R_0 and it is now "correct" (or at least consistent with much of the literature).
  • There are a bunch of new experiments and results in teh SDE lecture. At this point, I think the lectures are complete.

I am going to merge this PR because it is the only way to initiate a preview build for quantecon. We can then review the preview and do any copy editing. Or we do it in waves (i.e. get it live, then copy edit over time).

@jlperla jlperla merged commit 7038f71 into master Jul 22, 2020
@jlperla jlperla deleted the covid-sde branch July 22, 2020 23:48
Sign up for free to subscribe to this conversation on GitHub. Already have an account? Sign in.
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants