Oil Shock Model - Continuous

I bit the bullet and updated the oil shock model to instead do a grind-it-out numerical integration. As it stood, the original code got a little too intricate for my tastes; yet it served it's purpose in providing me some valuable insight. Overall, the results duplicate that of the previous algorithm, plus or minus some numerical error.

The code at the bottom of the post solves the following differential equation:

R(t+dt) = R(t) + (T(t) - R(t) * E(t)) * dt

P(t) = E(t) * R(t)

where

R(t) = Current reserves
T(t) = Triangular discovery curve
E(t) = Extraction rate (yearly or daily)
P(t) = Yearly (or daily) Production

The equation basically states that the reserves accumulate by discovery but deplete by extraction proportional to the amount available for extraction. This latter proportionality allowed me to do an exponential convolution, but the book-keeping of a variable extraction rate caused by temporal oil shocks made the algorithm a bit unwieldy. With a numerical integration, I can apply the proverbial mathematical sledge hammer to get at the results.

One can also now see how this differs from the logistic formulation favored by other analysts:

Q(t+dt)= Q(t) + k*Q(t)*(1-Q(t))*dt

(where Q is the cumulative production as a proportion of the ultimately recovered resource, and k is a constant that sets the width of the peak) Yet I still don't understand what this means in terms of a physical model.

I kept the same temporal oil shocks in the model but I made their duration match the span of the political crisis in real terms.


It helps to know the dates at which the oil shocks occurred. The oil embargo started in late 1973 and lasted until the middle of 1974. The Iranian hostage crisis started in late 1979 and the early 80's recession officially ended by 1984. The gulf war started in 1990 and it's associated recession ended in 1992. The change in extraction rate was interpolated over each of these intervals to make the suppressions in the curve a bit smoother.

with Text_IO;
procedure Shock_Model_Continuous is

  type Shock is record
    Year : Float; -- Year of shock
    Rate : Float; -- Rate of extraction
  end record;
  type Shocks is array (Natural range <>) of Shock;

  Triangle_Width : constant := 87.0;    -- Triangle window
  Volume         : constant := 2400.0;  -- In BBls

  S : constant Shocks :=
   ((1944.0, 0.060), -- Start of data
    (1974.0, 0.060), -- Start of oil embargo
    (1974.5, 0.051), -- End of oil embargo
    (1980.0, 0.051), -- Start of Iran hostage crisis
    (1984.0, 0.034), -- End of recession
    (1990.0, 0.034), -- Start of Gulf War
    (1992.0, 0.030), -- End of recesssion
    (2100.0, 0.030));

  -- Make the shock a continuous function
  -- by interpolating between points in the list
  function Interpolate_Shocks (Year : Float) return Float is
    K : Integer := S'Last - 1;
  begin
    for J in  S'First + 1 .. S'Last loop
      if S (J).Year > Year then
        K := J - 1;
        exit;
      end if;
    end loop;
    return S (K).Rate +
           (S (K + 1).Rate - S (K).Rate) * (Year - S (K).Year) /
           (S (K + 1).Year - S (K).Year);
  end Interpolate_Shocks;

  -- Creates a symmetric window for dicovered reserves
  function Triangle_Window
   (Year   : Float;
    Start  : Float;
    Length : Float)
    return   Float
  is
    Y    : constant Float := Year - Start;
    Half : constant Float := Length / 2.0;
  begin
    if Y < 0.0 or Y > Length then
      return 0.0;
    end if;
    if Y < Half then
      return Y / Half;
    else
      return (Length - Y) / Half;
    end if;
  end Triangle_Window;

  C      : Float          := 0.0;  -- Reserve
  R      : Float;         -- Extraction Rate
  P      : Float;         -- Production Rate
  DT     : constant Float := 0.01; -- 1/100th of year
  Start  : constant Float := S (S'First).Year;
  Finish : constant Float := S (S'Last).Year;
  Time   : Float          := Start;
  K      : constant Float := Volume / Triangle_Width * 2.0;
  V      : Float          := 0.0; -- Volume of oil extracted
begin
  loop
    R := Interpolate_Shocks (Time);
    C := C +
         K * Triangle_Window (Time, Start, Triangle_Width) * DT -
         C * R * DT;
    P := R * C / 365.0; -- BBls/day
    V := V + R * C * DT;
    if Time > Float (Integer (Time)) - DT / 2.0 and
       Time < Float (Integer (Time)) + DT / 2.0
    then  -- Print output only once per year
      Text_IO.Put_Line
       (Integer'Image (Integer (Time)) & "," & P'Img & "," & V'Img);
    end if;
    Time := Time + DT;
    exit when Time > Finish;
  end loop;
end Shock_Model_Continuous;

The output plotted against the BP data.