F(t) = Discovery stimulus
R1(t) = Reserve emerging from fallow state, Rate = a
R2(t) = Reserve emerging from construction state, Rate = b
R3(t) = Reserve emerging from maturation state, Rate = c
R(t) = Reserve emerging from production state, Rate = d
P(t) = Production curve
The stochastic differential equations look like:
dR1/dt = F(t) - a*R1(t)
dR2/dt = R1(t) - b*R2(t)
dR3/dt = R2(t) - c*R3(t)
dR/dt = R3(t) - d*R(t)
P(t) = d*R(t)
This forms a set of linear differential equations. If we take the Laplace transform of this set and do the transitive substitution, we can get the production curve in s-space.
r1(s) = f(s)/(s+a)If we assume a single delta for discoveries, then
r2(s) = a*r1(s)/(s+b)
r3(s) = b*r2(s)/(s+c)
r(s) = c*r3(s)/(s+d)
p(s) = f(s)*a*b*c*d/(s+a)/(s+b)/(s+c)/(s+d)
f(s)=1
. The inverse Laplace transform gives the following (unscaled) time-domain expressionFor values of rates very near 2.0, the production curve looks like this:
Remember to make sure that no two rates identically equate or else the solution becomes degenerate as the multiple poles form singularities. I mention this because the formulation as shown should prove useful in an optimization setting. By scanning through the ranges of the set of (a,b,c,d) one can quickly zero in on a first-order fit for a known discovery date and corresponding production data.
Up to now, I have used a numerical integration scheme to solve these equations, but the straight derivation provides a bit of insight into how the phased time constants arithmetically combine the exponentials into forming the asymmetric production profile. However, the simplistic assumption of a delta discovery and constant rates prevent me from recommending the close-form solution for complex, highly-featured real-world production curves.