vintage_cost_calculation.rb

require 'integration'
require 'bigdecimal/math'

# Calculates the weighted costs of ER credits from a list of given vintages and their contained ERs.
# Early vintages have higher prices allowing later vintages to be less effected by initial costs.
class VintageCostCalculation
  # Calculating raises an error if the sum of vintage ERs multiplied with the ER costs in this
  # vintage differs from the given project costs by more than ACCEPTABLE_PROJECT_COST_DIFFERENCE.
  ACCEPTABLE_PROJECT_COST_DIFFERENCE = 16

  # FULL_AREA_POSITION_FOR_NORMAL defines the cut-off point for the normal distribution after which
  # the area is ignored. The area from 0 to this point is considered to be 1 to avoid working with
  # infinity.
  # For standard normal distribution 5 is already at 0.9999994256677729.
  # To be really correct the integrate gem supports the GSL (GNU Scientific Library). See:
  # https://github.com/SciRuby/rb-gsl.
  # Using rb-gsl would mean depending on a new gem and would require GSL to be installed on the
  # system. If all the prerequisites are met full_are_position_for_normal could be set to :infinity
  # when the integrate method is set to method: :qag, instead of e.g. method: :simpson.
  FULL_AREA_POSITION_FOR_NORMAL = 100.0

  BIGMATH_PRECISION = 8

  require 'bigdecimal'

  Struct.new('VintageData', :year, :contained_er, :start_tonnes_position, :end_tonnes_position)
  attr_reader :project_costs, :vintages_data, :method_params

  def initialize(project_costs, vintages_data, method_params = { method: :n_of_cost_on_first_m_of_vintages })
    @vintages_data = vintages_data.each_with_object([]) do |data, memo|
      memo << Struct::VintageData.new(
        data.first,
        data.last,
        memo.sum(&:contained_er),
        memo.sum(&:contained_er) + data.last
      )
    end
    @project_costs = project_costs
    @method_params = method_params
  end

  def cost_per_tonne_by_vintage
    result = vintages_data.each_with_object({}) do |vintage_data, memo|
      memo[vintage_data.year] = cost_per_t_for_vintage(vintage_data)
    end
    if valid_output?(result)
      raise "Calculation not precise enough! Given project_costs were: #{project_costs}. Calculated sum of vintage costs is #{calculated_project_costs(result)}."
    end
    result
  end

  private

  def valid_output?(result)
    (project_costs - calculated_project_costs(result)).abs > ACCEPTABLE_PROJECT_COST_DIFFERENCE
  end

  def calculated_project_costs(result)
    @calculated_project_costs ||= result.sum do |year, costs_per_er|
      costs_per_er * vintages_data.find { |vintage| vintage.year == year }.contained_er
    end.round
  end

  def all_vintage_er
    @all_vintage_er ||= vintages_data.sum(&:contained_er)
  end

  def cost_per_t_for_vintage(vintage_data)
    fraction_of_costs(vintage_data) * project_costs / vintage_data.contained_er
  end

  def fraction_of_costs(vintage_data)
    if method_params.fetch(:method) == :n_of_cost_on_first_m_of_vintages
      integrate(
        vintage_data.start_tonnes_position,
        vintage_data.end_tonnes_position,
        lambda do |position|
          n_of_cost_on_first_m_of_vintages(
            position, method_params[:cost_fraction], method_params[:first_vintages_fraction]
          )
        end
      )
    elsif method_params.fetch(:method) == :normal_distribution
      integrate(
        (vintage_data.start_tonnes_position * FULL_AREA_POSITION_FOR_NORMAL / all_vintage_er),
        vintage_end_position_for_normal(vintage_data),
        ->(position) { normal(position, method_params[:sigma]) }
      ) * 2
      # * 2 because we're only looking at the positive side of the normal distribution.
    end
  end

  def integrate(start_pos, end_pos, func)
    Integration.integrate(start_pos, end_pos, tolerance: 1e-4, method: :simpson) do |position|
      func.call(position)
    end
  end

  def vintage_end_position_for_normal(vintage_data)
    if last_vintage?(vintage_data)
      FULL_AREA_POSITION_FOR_NORMAL
    else
      vintage_data.end_tonnes_position * FULL_AREA_POSITION_FOR_NORMAL / all_vintage_er
    end
  end

  def last_vintage?(vintage_data)
    vintage_data.end_tonnes_position == all_vintage_er
  end

  def n_of_cost_on_first_m_of_vintages(position, cost_fraction = nil, first_vintages_fraction = nil)
    first_vintages_fraction ||= BigDecimal(1) / BigDecimal(3)
    cost_fraction ||= BigDecimal(1) / BigDecimal(2)

    if position <= first_vintages_fraction * all_vintage_er
      cost_fraction / (first_vintages_fraction * all_vintage_er)
    elsif position <= all_vintage_er
      (1 - cost_fraction) / ((1 - first_vintages_fraction) * all_vintage_er)
    else
      0
    end
  end

  def normal(position, sigma = nil)
    # The larger sigma the closer the bell-curve resembles an even distribution.
    # The smaller sigma the higher the weighting of the first vintages.

    sigma ||= 1.0

    1 / (sigma * BigMath.sqrt(2 * BigMath.PI(BIGMATH_PRECISION), BIGMATH_PRECISION)) *
      BigMath.exp(-0.5 * (position / sigma)**2, BIGMATH_PRECISION)
  end
end

vintage_cost_calculation_spec.rb

require_relative '../vintage_cost_calculation.rb'

describe VintageCostCalculation do
  it 'returns an empty list if called with no vintages' do
    expect(described_class.new(0, []).cost_per_tonne_by_vintage).to(eq({}))
  end

  context 'if called with a single vintage' do
    let(:vintages_data) { { '2016' => 10_000 } }

    it 'returns the average cost of an er in the vintage' do
      expect(described_class.new(8000, vintages_data).cost_per_tonne_by_vintage)
        .to include('2016' => a_value_within(0.01).of(0.8))
    end
  end

  context 'with three 10_000t vintages' do
    let(:vintages_data) { { '2016' => 10_000, '2017' => 10_000, '2018' => 10_000 } }

    it 'returns ekos calculations' do
      expect(described_class.new(8000, vintages_data).cost_per_tonne_by_vintage)
        .to include(
          '2016' => a_value_within(0.01).of(0.4),
          '2017' => a_value_within(0.01).of(0.2),
          '2018' => a_value_within(0.01).of(0.2)
        )
    end

    it 'makes sense with custom cost_fraction and first_vintages_fraction' do
      expect(
        described_class
        .new(
          8000,
          vintages_data,
          method: :n_of_cost_on_first_m_of_vintages,
          cost_fraction: 1 / 2.0,
          first_vintages_fraction: 1 / 5.0
        )
        .cost_per_tonne_by_vintage
      ).to include(
        '2016' => a_value_within(0.01).of(0.466),
        '2017' => a_value_within(0.01).of(0.167),
        '2018' => a_value_within(0.01).of(0.167)
      )
    end
  end

  context 'with a contrived example fitting the n-th of cost' do
    let(:vintages_data) do
      { '2015' => 10_000, '2016' => 10_000, '2017' => 10_000, '2018' => 10_000, '2019' => 10_000 }
    end

    it 'makes sense' do
      expect(
        described_class
        .new(
          10_000,
          vintages_data,
          method: :n_of_cost_on_first_m_of_vintages,
          cost_fraction: 1 / 2.0,
          first_vintages_fraction: 1 / 5.0
        )
        .cost_per_tonne_by_vintage
      ).to include(
        '2015' => a_value_within(0.01).of(5000.0 / 10_000.0),
        '2016' => a_value_within(0.01).of(5000.0 / 40_000.0),
        '2017' => a_value_within(0.01).of(5000.0 / 40_000.0),
        '2018' => a_value_within(0.01).of(5000.0 / 40_000.0),
        '2019' => a_value_within(0.01).of(5000.0 / 40_000.0)
      )
    end
  end

  context 'with a cut-off in the first vintage' do
    let(:vintages_data) { { '2016' => 20_000, '2017' => 10_000, '2018' => 10_000 } }

    it 'makes sense with one third cut-off being in one vintage' do
      expect(described_class.new(8000, vintages_data).cost_per_tonne_by_vintage)
        .to include(
          '2016' => a_value_within(0.01).of(0.25),
          '2017' => a_value_within(0.01).of(0.15),
          '2018' => a_value_within(0.01).of(0.15)
        )
    end
  end

  context 'with a cut-off in the second vintage' do
    let(:vintages_data) { { '2016' => 10_000, '2017' => 20_000, '2018' => 10_000 } }

    it 'makes sense with one third cut-off being in second vintage' do
      expect(described_class.new(8000, vintages_data).cost_per_tonne_by_vintage)
        .to include(
          '2016' => a_value_within(0.1).of(0.30),
          '2017' => a_value_within(0.1).of(0.175),
          '2018' => a_value_within(0.1).of(0.15)
        )
    end

    it 'uses normal distribution if called with that option' do
      expect(described_class.new(8000, vintages_data, method: :normal_distribution).cost_per_tonne_by_vintage)
        .to include(
          '2016' => a_value_within(0.1).of(0.79),
          '2017' => a_value_within(0.1).of(0.004),
          '2018' => a_value_within(0.1).of(0.0)
        )
    end
  end

  context 'tembea example' do
    let(:vintages_data) do
      # vintages_data = CarbonProject.find(id).each_with_object({}) do |memo|
      #  vintage.name => vintage.er_contracted
      # end
      {
        '20121' => 10_467,
        '20122' => 9_553,
        '2013' => 26_740,
        '2014' => 35_586,
        '2015' => 58_089,
        '2016' => 85_606,
        '2017' => 107_505,
        '20181' => 100_000,
        '20182' => 29_620,
        '20191' => 100_000,
        '20192' => 38_395,
        '20201' => 100_000,
        '20202' => 46_050,
        '20211' => 100_000,
        '20212' => 53_593,
        '20221' => 100_000,
        '20222' => 57_255,
        '20231' => 100_000,
        '20232' => 43_672,
        '20241' => 145_785,
        '20242' => 45_785
      }
    end
    let(:project_cost) { 250_000 }

    it 'returns tfis calculated values' do
      # p described_class.new(project_cost, vintages_data).cost_per_tonne_by_vintage
      described_class.new(project_cost, vintages_data, method: :normal_distribution, sigma: 25.0).cost_per_tonne_by_vintage.each do |year, costs|
        p "Costs in Vintage #{year} are #{costs.to_f.round(3)}"
      end
    end
  end
end