# Compute evolution of a distribution over time

We have a population of people with different age $a$, time is indexed with $t$. There is a rate at which people die, $d(a, t)$. For simplicity, ignore births. I want to compute the evolution of the distribution of ages over time.

Denote the mass of people at or below age $a$ by $F(a,t)$

$$F(a,t) = \int_0^{a} m(\tilde a,t) d\tilde a$$

Ultimately, I am after some Kolmogorov forward equation, that is, the solution for

$$\partial_t F(a,t)$$

My approach Let $f(a, t)$ denote the density of people at age $a$ and point in time $t$. I will start with a discrete time approximation and let $\Delta$ go to zero. At each discrete point in time,

$$f(a+\Delta, t+\Delta) = (1-P(a, t))f(a, t)$$

where $P(a, t)$ is the discrete time analogue of $d(a,t)$. As I'm going to let $\Delta\to 0$, I can approximate $1-P$ with $1-\Delta d)$:

$$f(a+\Delta, t+\Delta) = (1- \Delta d(a,t))f(a, t)\\ \frac{f(a+\Delta, t+\Delta) -f(a,t)}{\Delta} = -d(a,t))f(a, t)\\ (\partial_t + \partial_a)f(a,t) = \lim_{\Delta\to 0}\frac{f(a+\Delta, t+\Delta) -f(a,t)}{\Delta} = -d(a,t))f(a, t)\\$$

I can integrate both sides w.r.t. a and get

$$\partial_t F(t, a) = - f(t, a) - \int q(t, a) f(t, a) da \\ \partial_t F(t, a) = - \partial_a F(t, a) - \int q(t, a) \partial_a F(t, a) da$$

I know that $\partial_a q(t, a) = q(t, a) (1-q(t, a))$. However, that doesn't really help me with solving the integral. Is there perhaps another angle to attack this problem? Or did I miss something?

• What is $\partial_t F(a,t)$? Is that notation for $\frac{\partial F}{\partial t} (a, t)$? – jmbejara May 11 '15 at 16:28
• it is notation for $\frac{\partial F(a, t)}{\partial t}$; the left hand side is the total derivative. – FooBar May 11 '15 at 16:29
• If you know $m$, why can't you just compute $\frac{\partial F(a, t)}{\partial t} = \int_0^a m_t(\tilde a,t) d \tilde a$ with Leibniz's rule? – jmbejara May 11 '15 at 16:32
• You can look at en.wikipedia.org/wiki/Feynman%E2%80%93Kac_formula – VicAche May 13 '15 at 15:53
• quant.stackexchange.com/questions/10359/… maybe you should ask your question there if this is not enough, actually. They seem far more confident that I am about the magic :) – VicAche May 14 '15 at 14:41

Here's my best guess. I haven't checked to thoroughly if this is right, but maybe it will help.

# Evolution of population density

I understand the model as follows. $f(a,t)$ is the density of people of age $a$ at time $t$. Suppose at time $t=0$, the density of the population is $f_0(a)$. To model the aging process as well as the mortality rate, the density $f$ must evolve over time so that it satisfies the condition that you derived. Thus, $m$ must satisfy the following partial differential equation $$\frac{\partial f}{\partial t} + \frac{\partial f}{\partial a} + d(a, t) f(a,t) = 0$$ and the initial condition $$f(a,0) = f_0(a).$$ I believe that PDEs of this form can have relatively straightforward solutions depending on the functional form chosen for $d$. Here are some cases.

Case 1: Constant mortality rate, $d(a,t) = d_0$.

Suppose a constant mortality rate, $d(a,t) = d_0$. Then (using Mathematica),

DSolve[
{D[f[a, t], t] + D[f[a, t], a] + d0 f[a, t] == 0,
f[a, 0] ==  f0[a]},
f[a, t], {a, t}]


results in

{{f[a, t] -> E^(-a d0 + d0 (a - t)) f0[a - t]}}


So, as we can see, this gives a simple solution $$f(a,t) = \exp\{-a d_0 + d_0 (a-t)\} \cdot f_0(a-t)$$ Case 2: Log Mortality, $d(a,t) = \log(a+1)$.

Here we have

DSolve[
{D[f[a, t], t] + D[f[a, t], a] + Log[a+1] f[a, t] == 0,
f[a, 0] ==  f0[a]},
f[a, t], {a, t}]


which results in

{{f[a, t] -> (1 + a)^(-1 - a) E^t (1 + a - t)^(1 + a - t) f0[a - t]}}


This also gives a simple solution $$f(a,t) = (1+a)^{-a-1} e^t (1+a-t)^{1+a-t} \cdot f_0(a-t)$$

Case 3: The general case

For the case where we do not yet specify $d(a,t)$,

DSolve[
{D[f[a, t], t] + D[f[a, t], a] + d[a, t] f[a, t] == 0,
f[a, 0] ==  f0[a]},
f[a, t], {a, t}]


gives us (in TeXForm for clarity)

\begin{align*} \left\{\left\{f(a,t)\to \\ \text{f0}(a-t) \cdot \\ \exp \left(\int_1^a -d(K[1],K[1]-a+t) \, dK[1]-\int_1^{a-t} -d(K[1],K[1]-a+t) \, dK[1]\right)\right\}\right\} \end{align*}

• As stated in the comments of the question above, they are (Monte-Carlo) approaches that are pretty efficient on the generic case. Why would you assume log or constant MR? – VicAche May 18 '15 at 8:54
• I used those two cases as simple examples. I've put the general case in now. I don't quite see how you would use MC methods. I'll have to look into it more. I understand that Feynman-Kac gives an equivalence between certain PDEs and stochastic processes and that this allows us to solve these PDEs by approximated a conditional expectation with MC methods. However, as @FooBar mentioned, I'm having trouble seeing how to fit this problem into that framework. – jmbejara May 18 '15 at 18:07
• Also, as for efficiency, I believe that the strength of Feynman-Kac is that you can avoid Monte Carlo methods. As I understand it, solution by MC methods is reliable but not particularly efficient. Shreve, in Stochastic Calculus for Finance II (p.268), says "the Euler methods described (an MC method)... for determining this function converges slowly and gives the function value for only one pair (t,x). Numerical algorithms for solving [the PDE] converge quickly in the case of one dimensional $x$ being considered here... and given the function g(t,x) for all values simultaneously." (abridged) – jmbejara May 18 '15 at 18:13
• Your answer is much more convincing in its current form. You get the upvote :) – VicAche May 18 '15 at 18:19
• tho you don't quite answer the question I believe. You're giving the integral we need to compute, not a proposed method (except "using mathematica") to compute it. – VicAche May 18 '15 at 18:20