Oil Shock Model
Update: I added an appendix at the bottom that provides the algorithm for the oil shock model. Per Bubba's suggestion, I am thinking of how to redisplay the URR plot to better guide the eye.
Following up on my critique of Stuart Staniford's estimate of global cumulative endpoint production of oil (aka the URR estimate), I found enough additional misinformation or plain mistakes in SS's model to refine my macro peak oil model (1, 2) further.
The mistake that Staniford made in applying the logistics model to BP's data came from the wrong initial (pre-1965) cumulative production. Look at the following plot from TOD and you see a telling "cusp" at the top of the curve:
That cusp disappears when one applies the right initial cumulative data; Staniford added about 100 Billion Barrels (BBls) too much to this data. This becomes obvious when you notice that his cumulative as of 2004 sits at about 1050 BBls -- 10% over the generally accepted value of 952 BBls from the USGS.
This time I downloaded the data from BP's World Energy 2005 page and tried to make a more decent try at a model fit. Based on comments at The Oil Drum, two issues that I consider very important, (a) using historical data and (b) using realistic yet simple models, apparently tends to rub some people the wrong way. I think I have good arguments covering these concerns.
For some reason, the use of historical data causes certain oil watchers to scoff at us "backward" thinkers; they don't realize that we can learn a lot, both in a Bayesian sense and avoiding repeating the same mistakes over and over again (how shocking). Otherwise, I can only buy the argument up to a point; for example, until we hit the peak, we can't with 100% certainty predict the peak. Having a real inflection point or an obvious concavity does wonders for this class of analytical model.
I also find the simplicity of the model important -- something that Staniford agrees with. Unfortunately, he simplifies to such an extent, by using the Logistics curve, that he only has a single parameter to adjust. I find that much like having to make a prediction given only the mean of some data set. I would at least like to have a couple of parameters, say the mean and variance, to improve my chances. Here, Staniford throws away valuable information that has a good physical basis to achieve a simple empirical fit. I side with a good physical model that will give good intuition and insight to the problem at hand.
With that in mind I used my convolution-based Macro Peak Oil model to fit the BP data along with the USGS cumulative value of 952 BBls. The following uses a symmetric Triangular Discovery1 curve starting in 1944 with a width of 87 years and normalized to an URR value of 2400 BBls.
The early years (pre 1970) of my model feature an oil extraction rate of 6% of total content per year. For an exponential drop-off, this gives a 1/e point of approximately 16 years.
I add three shocks in the years 1973-1974, 1979-1980, and 1991-1992, whereby the oil extraction rate changed dramatically. Because of the stationary properties of the Markov model, I can simply adjust the rates of the normalized exponential term in the middle of the convolution without violating the stochastic nature of oil extraction over time. The shocks tend to momentarily suppress and flatten the production/consumption rate. (As a technical aside, the suppressive dip in the output comes about mathematically from the exponential normalization)
The three shocks correspond to the OPEC embargo, Iranian crises coupled with a deep recession, and the first Gulf war.
No doubt we will get more shocks in the future. The crucial finding in my mind: the shocks serve to delay significantly the onset of peak oil. Before the 70's, we used oil as if it came out of the tap; we have since made significant corrections in the extraction rate and our more conservative use of oil. We will likely get a shock fairly soon -- this will not invalidate the model results and it will not disprove the peak oil hypothesis. It simply will serve to delay the peak a bit more. I really believe (and the model shows) that we would have hit a peak several years ago without these shocks in place.
The final URR fit looks like the following curve. Of course the asymptote hits 2400 BBls because of the triangular distribution I started with. I mainly want to demonstrate how the annoying cusp disappears when plotted correctly.
1 I chose the triangular curve because of the proliferation of the following type of data.
Note that I use the term "discovery" a bit loosely; in my mind it actually indicates when pumping activity actually started on the previously discovered reserves.
This code compiles with the GCC
gnatmakecompiler provided with most recent vintages of Linux (or available for Windows here). To build, save as file "shock_model.adb", then
gnatmake shock_modelto make the executable. Spreadsheets like the CSV-format: "year", "BBls/day", "cumulative BBls"
procedure Shock_Model is
subtype Years is Natural range 1944 .. Natural'Last;
type Shock is record
Year : Years; -- Year of shock
Rate : Float; -- Prod. or depletion rate (fraction / year)
type Shocks is array (Years range <>) of Shock;
Width : constant := 87; -- Triangle window
Volume : constant := 2400.0; -- In BBls
S : constant Shocks :=
((Years'First, 0.06), -- Starting Year
(1974, 0.048), -- OPEC
(1980, 0.035), -- Iran +
(1983, 0.029), -- Recession
(1991, 0.026), -- Gulf I
(2150, 0.0)); -- Ending Year
--##### You don't have to look beyond this point unless you
--##### want to mess with the algorithm or output.
-- General purpose array of floats indexed by year
type Flt is array (Natural range <>) of Float;
-- Safe array get function, truncates values out of range
function Get (Arr : Flt; I : Integer) return Float is
if I in Arr'Range then
return Arr (I);
-- Takes 2 arrays and returns the convolution array
function Convolve (A, B : in Flt) return Flt is
Total : constant Natural := A'Length + B'Length;
C : Flt (0 .. Total);
V : Float;
for J in 0 .. Total loop
V := 0.0;
for I in 0 .. J loop
V := V + Get (A, I) * Get (B, J - I);
C (J) := V;
-- Creates a symmetric window, truncated left by Start
-- The prior accumulated data is provided by Impulse
(Width : in Natural;
Start : in Natural := 0;
Impulse : Float := 0.0)
R : Flt (0 .. Width - 1 - Start);
Half : constant Natural := Width / 2;
Offset : Natural;
for I in R'Range loop
Offset := I + Start;
if Offset < Half then
R (I) := Float (Offset) / Float (Half);
R (I) := Float (Width - Offset) / Float (Half);
R (0) := R (0) + Impulse;
-- An exponential array populated to length
function Exponential (Length : in Natural;
Alpha : in Float) return Flt is
R : Flt (0 .. Length - 1);
for I in 0 .. Length - 1 loop
R (I) := Alpha * Exp (-Alpha * Float (I));
-- The summed difference between two arrays up to Width
-- This is used to renormalize when an Alpha rate changes
function Leftover (A, B : in Flt;
Width : in Natural) return Float is
V1 : Float := 0.0;
V2 : Float := 0.0;
for J in 0 .. Width - 1 loop
V1 := V1 + Get (A, J);
V2 := V2 + Get (B, J);
V1 := V1 + 0.5 * Get (A, Width);
V2 := V2 + 0.5 * Get (B, Width);
return (V1 - V2);
-- Recursive formula to re-convolute on successive rate shocks
-- Returns production curve in array
(S : in Shocks; -- Input stimulus
Left : in Float; -- Pre-history reserve level
Curve : in Flt; -- The starting curve
TS : in Natural) -- TimeSpan of results
if S'Length = 1 then -- returns on completion
W1 : constant Natural := S (S'First).Year - Years'First;
W2 : constant Natural := S (S'First + 1).Year - Years'First;
T : constant Flt := Triangle_Window (Width, W1, Left);
R : constant Flt :=
Convolve (T, Exponential (TS, S (S'First).Rate));
(S => S (S'First + 1 .. S'Last),
Left => Leftover (T, R, W2 - W1),
Curve => Curve & R (1 .. W2 - W1),
TS => TS);
Time_Span : Natural := S (S'Last).Year - S (S'First).Year;
R : constant Flt := Gen (S, 0.0, (0 => 0.0), Time_Span);
Total : Float := 0.0;
K : constant Float := Volume / Float (Width) * 2.0;
for I in R'Range loop
Total := Total + K * R (I);
(Integer (Years'First + I)'Img & "," & -- Year
Float'Image (K * R (I) / 365.0) & "," & -- Per Day
As a check, the red curve displays the output of the above code. The yellow curve shows what would have happened without the oil shocks. Through the roof, baby.