Probability density function that evolves with time according to a delay differential Consider a real valued variable $X(t)$ that evolves with time according to the delay differential
$\frac{dX(t)}{dt} = \alpha X(t-t_0) \int_{t_0}^\infty f(y) h(t-t_0,y) dy - \beta X(t) \int_{t_0}^\infty g(y) h(t,y) dy$ 
where $h(t,z)$ is a time-evolving probability density over $z$ that is determined by $X$ 
$h(t,z) \propto \exp{[-\beta \int_{t_0}^z g(y) dy]} \alpha X(t-z)\int_{t_0}^\infty f(y) h(t-z,y) dy$.  
Assume that functions $f()$ and $g()$ behave nicely.
My question is: does this make sense? It seems rather circular. In particular I'm concerned that in the limit $z\rightarrow 0$ the expression for $h(.,.)$ approaches a function of itself at an identical time. Perhaps I need to make it $h(t+dt,z)$ on the left hand side? If so, is there an intuitive way to understand this that will help me solve these equations numerically?
A curiosity (sign I'm doing something wrong?) of this formulation is that $z$ seems to sometimes behave like a random variable and sometimes not.   
 A: So I was on the right track, but my PDE made no sense.
Let's simplify by letting $t_0=0$ and removing the delay. Then
$$ \frac{\partial X}{\partial t}(t,z) + \frac{\partial X}{\partial z}(t,z) = -\beta g(z) X(t,z) $$
fully specifies the decay along so-called characteristics (loci where $t=z + c$) and 
$$ X(t,0) = \alpha \int_0^\infty f(y)X(t,y) dy $$
is a time varying, non-local (because of the integral) boundary condition which governs growth / renewal. 
Apparently this style of model was originally conceived by Scottish physician and amateur mathematician A. G. McKendrick in 1926, by all accounts a remarkable fellow. 
The PDE reduces to an ODE along characteristics, greatly simplifying its solution.
Thanks to my friend George Leigh for help with this one.
A: Ok, here's my own stab at this.
First thing to realise is that you're dealing with a partial differential equation. To see this, consider the initial condition required: the distribution $h$ over $z$ must be completely specified, not just at a single point. Of course you also have delays making this a delay partial differential: the full initial condition is actually a function of both $z$ (from $t_0$ to $\infty$) and $t$ (from $t-t_0$ to $t$). 
The second thing to realise is that it makes more sense to consider $X$ alone as the unknown function of $t$ and $z$, and just use integrals of $X$ to achieve the work done by $h$ when needed. This simplifies the situation considerably, as well as removing the random - non-random variable weirdness:
$\frac{\partial X}{\partial t}(t,z) = \alpha \int_{t_0}^\infty X(t-z,y)f(y)dy$
$\frac{\partial X}{\partial z}(t,z) = -\beta \int_{t_0}^z X(t,y)g(y)dy$
The intuitive picture here is of a rectangular two dimensional sheet moving smoothly forward along one dimension (the time dimension) and smoothly changing shape as it does so with the "front corner" increasing as a function of the "back edge" and and the front edge decaying from that increase as a function of the front edge.  
As far as numerical solutions go, this is an integro delay partial differential, so it might be hard, but note that the rates of change appear to be independent of each other so there are probably some simplifying tricks.
