WIP: Covid sde #867
WIP: Covid sde #867
Conversation
|
Implements some of the discussion in #859 |
|
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,
I'll review more carefully when I have time but everything is clear and nicely explained. |
|
I plan to go through this in detail, but let's say this weekend. I had a paper resubmission come up. |
| 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` |
There was a problem hiding this comment.
Not sure what v is supposed to be referring to.
There was a problem hiding this comment.
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. |
There was a problem hiding this comment.
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.
There was a problem hiding this comment.
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.
There was a problem hiding this comment.
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
|
|
||
| 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) |
There was a problem hiding this comment.
this is not an additive noise SDE
There was a problem hiding this comment.
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.
There was a problem hiding this comment.
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). |
There was a problem hiding this comment.
Perhaps this is not true once aggregate shocks are added, if, say, the paths admit jumps. Might be worth qualifying?
There was a problem hiding this comment.
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. |
There was a problem hiding this comment.
@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?
There was a problem hiding this comment.
I believe in your implementation you may talk about R0 instead? https://python.quantecon.org/sir_model.html#Implementation
There was a problem hiding this comment.
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?
There was a problem hiding this comment.
@jstac I implemented these sorts of changes in the latest version, and there are a significant number of changes along these topics.
|
@jstac @ChrisRackauckas @thomassargent30 I just pushed a major set of changes The summary of changes is:
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). |
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