Multi-Well DCF.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import scipy.optimize as optimize\n",
    "from scipy.interpolate import interp1d\n",
    "\n",
    "def arps_decline(UID, phase, Qi, Dei, Def, b, t, prior_cum, prior_t):\n",
    "    # UID is a unique identifier for the well such as API, must be a number\n",
    "    # phase is 1 = oil, 2 = gas, or 3 = water\n",
    "    # Qi is the initial production rate typically in bbl/day or Mcf/day\n",
    "    # Dei is the initial effective annual decline rate\n",
    "    # Def is the final effective annual decline rate at which point the decline becomes exponential\n",
    "    # b is the b-factor used in hyperbolic or harmonic decline equations\n",
    "    # t is the time as a month integer\n",
    "    # prior_cum is the cumulative amount produced before the start of the decline calcuations\n",
    "    # prior_t is an integer representing the final month from a previous decline segment\n",
    "    \n",
    "    # Calculations to determine decline type\n",
    "    if Dei == Def:\n",
    "        Type = 'exp'\n",
    "    elif Dei > Def and b == 1:\n",
    "        Type = 'har'\n",
    "        Dn = Dei / (1 - Dei)\n",
    "        Qlim = Qi * ((-np.log(1 - Def)) / Dn)\n",
    "        tlim = (((Qi / Qlim) - 1) / Dn) * 12 # output in months\n",
    "    else:\n",
    "        Type = 'hyp'\n",
    "        Dn = (1 / b) * (((1 - Dei) ** -b) - 1)\n",
    "        Qlim = Qi * ((-np.log(1 - Def)) / Dn) ** (1 / b)\n",
    "        tlim = ((((Qi / Qlim) ** b) - 1) / ( b * Dn)) * 12 # output in months\n",
    "    \n",
    "    # Generate volumes\n",
    "    if Type == 'hyp':\n",
    "        Dn_t = Dn / (1 + b * Dn * (t / 12))\n",
    "        De_t = 1 - (1 / ((Dn_t * b) + 1)) ** (1 / b)\n",
    "        if De_t > Def:\n",
    "            q = Qi * (1 + b * Dn * (t / 12)) ** (-1/b)\n",
    "            Np = ((Qi ** b) / (Dn * (1 - b))) * ((Qi ** (1 - b)) - (q ** (1 - b))) * 365\n",
    "        else:\n",
    "            q = Qlim * np.exp(-(-np.log(1 - Def)) * ((t - tlim) / 12))\n",
    "            Np = ((Qlim - q) / (-np.log(1 - Def)) * 365) + (((Qi ** b) / \n",
    "                    (Dn * (1 - b))) * ((Qi ** (1 - b)) - (Qlim ** (1 - b))) * 365)\n",
    "            De_t = Def\n",
    "    elif Type == 'har':\n",
    "        Dn_t = Dn / (1 + Dn * (t / 12))\n",
    "        De_t = 1 - (1 / (Dn_t + 1))\n",
    "        if De_t > Def:\n",
    "            q = Qi / (1 + b * Dn * (t / 12))\n",
    "            Np = (Qi / Dn) * np.log(Qi / q) * 365\n",
    "        else:\n",
    "            q = Qlim * np.exp(-(-np.log(1 - Def)) * ((t - tlim) / 12))\n",
    "            Np = ((Qlim - q) / (-np.log(1 - Def)) * 365) + ((Qi / Dn) * np.log(Qi / Qlim) * 365)\n",
    "            De_t = Def\n",
    "    else:\n",
    "        q = Qi * np.exp(-(-np.log(1 - Dei)) * (t / 12))\n",
    "        Np = (Qi - q) / (-np.log(1 - Dei)) * 365\n",
    "        De_t = Dei\n",
    "    \n",
    "    return UID, phase, t + prior_t, q, De_t, Np + prior_cum\n",
    "\n",
    "# Vectorize the arps_decline function to allow it to work with numpy arrays\n",
    "varps_decline = np.vectorize(arps_decline)"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "source": [
    "# Calculate Dei from Qi and Qf based on exponential decline equation\r\n",
    "\r\n",
    "def exp_Dei(Qi, Qf, duration):\r\n",
    "    # Qi is the initial production rate typically in bbl/day or Mcf/day\r\n",
    "    # Qf is the final production rate typically in bbl/day or Mcf/day\r\n",
    "    # duration is the time interval in months over which you are trying to calculate the exponential decline rate\r\n",
    "    Dei = 1 - np.exp(-np.log(Qi / Qf) / (duration / 12))\r\n",
    "    return Dei"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "source": [
    "# Function to manage multiple segments\r\n",
    "# Segment 1 is the initial incline period and uses Arps exponential equation\r\n",
    "# Segment 2 is the period between the incline and decline periods and uses Arps exponential equation\r\n",
    "# Segment 3 is the decline period\r\n",
    "\r\n",
    "def arps_segments(UID, phase, Q1, Q2, Q3, Dei, Def, b, t1, t2, duration):\r\n",
    "    # Determine valid segment count\r\n",
    "    if t1 > 0 and t2 > 0:\r\n",
    "        segment_ct = 3\r\n",
    "        if Q2 == Q3:\r\n",
    "            Q2 = Q2 * 1.0001\r\n",
    "    elif t1 > 0:\r\n",
    "        segment_ct = 2\r\n",
    "    elif Q3 == 0 or np.isnan(Q3):\r\n",
    "        segment_ct = 0\r\n",
    "    else:\r\n",
    "        segment_ct = 1\r\n",
    "    \r\n",
    "    # 3 segment logic\r\n",
    "    if segment_ct == 3:\r\n",
    "        t_seg1 = np.arange(0, t1 + 1, 1)\r\n",
    "        t_seg2 = np.arange(1, t2 + 1, 1)\r\n",
    "        t_seg3 = np.arange(1, duration - t1 - t2 + 1, 1)\r\n",
    "        Dei1 = exp_Dei(Q1, Q2, t1)\r\n",
    "        Dei2 = exp_Dei(Q2, Q3, t2)\r\n",
    "        seg1 = varps_decline(UID, phase, Q1, Dei1, Dei1, 1.0, t_seg1, 0, 0)\r\n",
    "        prior_cum1 = np.max(seg1[5])\r\n",
    "        seg2 = varps_decline(UID, phase, Q2, Dei2, Dei2, 1.0, t_seg2, prior_cum1, t1)\r\n",
    "        prior_cum2= np.max(seg2[5])\r\n",
    "        seg3 = varps_decline(UID, phase, Q3, Dei, Def, b, t_seg3, prior_cum2, t1 + t2)\r\n",
    "        out_nparr = np.column_stack((seg1, seg2, seg3))\r\n",
    "    elif segment_ct == 2:\r\n",
    "        t_seg1 = np.arange(0, t1 + 1, 1)\r\n",
    "        t_seg3 = np.arange(1, duration - t1 + 1, 1)\r\n",
    "        Dei1 = exp_Dei(Q1, Q3, t1)\r\n",
    "        seg1 = varps_decline(UID, phase, Q1, Dei1, Dei1, 1.0, t_seg1, 0, 0)\r\n",
    "        prior_cum1 = np.max(seg1[5])\r\n",
    "        seg3 = varps_decline(UID, phase, Q3, Dei, Def, b, t_seg3, prior_cum1, t1)\r\n",
    "        out_nparr = np.column_stack((seg1, seg3))\r\n",
    "    elif segment_ct == 1:\r\n",
    "        t_seg3 = np.arange(0, duration + 1, 1)\r\n",
    "        out_nparr = varps_decline(UID, phase, Q3, Dei, Def, b, t_seg3, 0, 0)\r\n",
    "    else:\r\n",
    "        t_nan = np.arange(0, duration + 1, 1)\r\n",
    "        UID_nan = np.full((1, duration + 1), UID)\r\n",
    "        phase_nan = np.full((1, duration + 1), phase)\r\n",
    "        val_nan = np.full((3, duration + 1), 0)\r\n",
    "        out_nparr = np.vstack((UID_nan, phase_nan, t_nan, val_nan))\r\n",
    "    \r\n",
    "    # Add monthly volumes to array\r\n",
    "    \r\n",
    "    Cum_i = out_nparr[5][:-1]\r\n",
    "    Cum_f = out_nparr[5][1:]\r\n",
    "    cum = Cum_f - Cum_i\r\n",
    "    cum[cum < 0] = 0\r\n",
    "    cum = np.insert(cum, 0, 0)\r\n",
    "    out_nparr = np.vstack((out_nparr, cum))[:,1:]\r\n",
    "    \r\n",
    "    return out_nparr"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "source": [
    "# Create MonthDiff function to calculate time difference in months\r\n",
    "\r\n",
    "def MonthDiff(BaseDate, StartDate):\r\n",
    "    MonthDiff = int(((np.datetime64(StartDate) - (np.datetime64(BaseDate) - np.timedelta64(1, 'M'))) / np.timedelta64(1, 'M')) - 1)\r\n",
    "    return MonthDiff"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "source": [
    "# function to calculate monthly cash flow output, no economic truncation applied to output\r\n",
    "\r\n",
    "def monthly_cf(index, uid, wi, nri, weight, prod_wt, inv_wt, shrink, btu, ngl_yield, pri_oil, pri_gas, paj_oil, \r\n",
    "               paj_gas, paj_ngl, stx_oil, stx_gas, stx_ngl, adval, opc_fix, opc_oil, opc_gas, capex):\r\n",
    "    # Calculate cash flow components\r\n",
    "    mos = np.extract(vol_np[index][0] == uid, vol_np[index][1])\r\n",
    "    gr_oil = np.extract(vol_np[index][0] == uid, vol_np[index][2])\r\n",
    "    gr_gas = np.extract(vol_np[index][0] == uid, vol_np[index][3])\r\n",
    "    oil_sales = gr_oil * nri * weight * prod_wt\r\n",
    "    gas_sales = gr_gas * nri * shrink * weight * prod_wt\r\n",
    "    ngl_sales = ngl_yield * gr_gas * nri * weight * prod_wt\r\n",
    "    oil_rev = np.round(oil_sales * (pri_oil + paj_oil), 2)\r\n",
    "    gas_rev = np.round(gas_sales * btu * (pri_gas + paj_gas), 2)\r\n",
    "    ngl_rev = np.round(ngl_sales * (pri_oil * paj_ngl), 2)\r\n",
    "    total_rev = oil_rev + gas_rev + ngl_rev\r\n",
    "    oil_tax = oil_rev * stx_oil\r\n",
    "    gas_tax = gas_rev * stx_gas\r\n",
    "    ngl_tax = ngl_rev * stx_ngl\r\n",
    "    adval_tax = total_rev * adval\r\n",
    "    total_tax = oil_tax + gas_tax + ngl_tax + adval_tax\r\n",
    "    FirstProd = np.min(np.nonzero(gr_oil + gr_gas))\r\n",
    "    opex = np.where(mos >= FirstProd, (wi * weight * inv_wt) * ((gr_oil * prod_wt * opc_oil) + (gr_gas * prod_wt * shrink * opc_gas) + opc_fix), 0)\r\n",
    "    op_cf = total_rev - total_tax - opex\r\n",
    "    net_cf = np.where(mos == FirstProd, op_cf - (capex * weight * inv_wt * wi), op_cf)\r\n",
    "    cum_opcf = np.cumsum(op_cf)\r\n",
    "    cum_ncf = np.cumsum(net_cf)\r\n",
    "\r\n",
    "    # create output array\r\n",
    "    out_full = np.vstack((vol_np[index][:, vol_np[index][0] == uid], oil_sales, gas_sales, ngl_sales, oil_rev, gas_rev, \r\n",
    "                          ngl_rev, total_rev, total_tax, opex, op_cf, cum_opcf, net_cf, cum_ncf))\r\n",
    "    return out_full"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "source": [
    "# function to calculate monthly cash flow output with economic truncation\r\n",
    "\r\n",
    "def econ_cf(index, uid, wi, nri, roy, eloss, weight, prod_wt, inv_wt, shrink, btu, ngl_yield, pri_oil, pri_gas, paj_oil, paj_gas,\r\n",
    "            paj_ngl, stx_oil, stx_gas, stx_ngl, adval, opc_fix, opc_oil, opc_gas, capex, aban):\r\n",
    "    if wi == 0 or nri / wi > (1 - roy):\r\n",
    "        life_cf = monthly_cf(index, uid, 1, 1 - roy, weight, prod_wt, inv_wt, shrink, btu, ngl_yield, pri_oil, \r\n",
    "                             pri_gas, paj_oil, paj_gas, paj_ngl, stx_oil, stx_gas, stx_ngl, adval, opc_fix, opc_oil, \r\n",
    "                             opc_gas, capex)\r\n",
    "        try:\r\n",
    "            life = np.where(life_cf[14] == np.max(life_cf[14]))[0][0]\r\n",
    "        except:\r\n",
    "            life = 1\r\n",
    "    else:\r\n",
    "        life_cf = monthly_cf(index, uid, wi, nri, weight, prod_wt, inv_wt, shrink, btu, ngl_yield, pri_oil, \r\n",
    "                             pri_gas, paj_oil, paj_gas, paj_ngl, stx_oil, stx_gas, stx_ngl, adval, opc_fix, opc_oil, \r\n",
    "                             opc_gas, capex)\r\n",
    "        try:\r\n",
    "            life = np.where(life_cf[14] == np.max(life_cf[14]))[0][0]\r\n",
    "        except:\r\n",
    "            life = 1\r\n",
    "    \r\n",
    "    result = monthly_cf(index, uid, wi, nri, weight, prod_wt, inv_wt, shrink, btu, ngl_yield, pri_oil, pri_gas, paj_oil, \r\n",
    "                        paj_gas, paj_ngl, stx_oil, stx_gas, stx_ngl, adval, opc_fix, opc_oil, opc_gas, capex)\r\n",
    "    \r\n",
    "    # add abandonment costs to end of life\r\n",
    "    net_cf = np.where(result[1] == life + eloss, result[15] - (aban * weight * inv_wt * wi), result[15])\r\n",
    "    cum_ncf = np.cumsum(net_cf)\r\n",
    "    \r\n",
    "    # replace net cf array\r\n",
    "    result = np.vstack((result[0:15], net_cf, cum_ncf))\r\n",
    "    \r\n",
    "    # zero out values past life\r\n",
    "    monthly_econ = np.where(result[1] <= life + eloss, result[:], np.vstack((result[0:2], result[2:] * 0)))\r\n",
    "    return monthly_econ"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "source": [
    "# function to calculate npv, assumes monthly cash flow and annualized discount rate\r\n",
    "\r\n",
    "def npv(disc_rate, ncf, n):\r\n",
    "    # disc_rate is an annualized discount rate\r\n",
    "    # ncf is the monthly cashflow array\r\n",
    "    # n is the array of month integers\r\n",
    "    rate = disc_rate / 12 # Corrects effective discount rate from annual to monthly for discounting\r\n",
    "    factor = 1 / ((1 + rate) ** -0.5) # Adjusts discounting to mid-period similar to Aries output\r\n",
    "    pv = np.sum(ncf / ((1 + rate) ** n)) * factor\r\n",
    "    return pv"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "source": [
    "# Import data for decline parameters\n",
    "prop_list = pd.read_csv('https://raw.githubusercontent.com/jshumway0475/Petroleum/main/property_table.csv')\n",
    "prop_list = prop_list.query('INCLUDE == 1')\n",
    "prop_list = prop_list.fillna(0)\n",
    "prop_list.reset_index(drop = True, inplace = True)\n",
    "C = prop_list.keys() # Columns\n",
    "R = prop_list.index # Rows"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "source": [
    "# Loop through DataFrame and output monthly oil volumes\r\n",
    "\r\n",
    "duration = 600\r\n",
    "\r\n",
    "oil = lambda w: arps_segments(prop_list.loc[w, 'PROPNUM'], 1, prop_list.loc[w, 'OIL_IP'], prop_list.loc[w, 'OIL_IP2'], \r\n",
    "                              prop_list.loc[w, 'OIL_IP3'], prop_list.loc[w, 'OIL_DI'] / 100, prop_list.loc[w, 'Oil_Def'] / 100, \r\n",
    "                              prop_list.loc[w, 'OIL_B'], prop_list.loc[w, 'OIL_SEG1_TIME'], prop_list.loc[w, 'OIL_SEG2_TIME'], \r\n",
    "                              duration)\r\n",
    "v_oil = np.vectorize(oil, otypes = [object])\r\n",
    "oil_nparr = v_oil(R)\r\n",
    "\r\n",
    "# Loop through DataFrame and output monthly gas volumes\r\n",
    "\r\n",
    "gas = lambda w: arps_segments(prop_list.loc[w, 'PROPNUM'], 1, prop_list.loc[w, 'GAS_IP'], prop_list.loc[w, 'GAS_IP2'], \r\n",
    "                              prop_list.loc[w, 'GAS_IP3'], prop_list.loc[w, 'GAS_DI'] / 100, prop_list.loc[w, 'GAS_Def'] / 100, \r\n",
    "                              prop_list.loc[w, 'GAS_B'], prop_list.loc[w, 'GAS_SEG1_TIME'], prop_list.loc[w, 'GAS_SEG2_TIME'], \r\n",
    "                              duration)\r\n",
    "v_gas = np.vectorize(gas, otypes = [object])\r\n",
    "gas_nparr = v_gas(R)"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "source": [
    "# Import price files\n",
    "\n",
    "BaseDate = '2021-07'\n",
    "MaxDate = pd.to_datetime(prop_list['Start_Date']).max()\n",
    "MaxDate_np = MaxDate.to_datetime64()\n",
    "add_months = int(round((MaxDate - np.datetime64(BaseDate)) / np.timedelta64(1, 'M'), 0)) + 1\n",
    "str_periods = duration + add_months\n",
    "dates = pd.DataFrame(pd.date_range(BaseDate, periods = str_periods, freq = 'MS'), columns = ['Date'])\n",
    "strip_price = pd.read_csv('https://raw.githubusercontent.com/jshumway0475/Petroleum/main/STR-071521.csv')\n",
    "strip_price['Date'] = pd.to_datetime(strip_price['Date'])\n",
    "\n",
    "# Create diff array\n",
    "pepl = {BaseDate: -0.137, '01/2022': -0.17, '01/2023': -0.249, '01/2024': -0.26, '01/2025': -0.214, '01/2026': -0.212}\n",
    "pepl_pd = pd.DataFrame(pepl.items(), columns = ['Date', 'paj_gas'])\n",
    "pepl_pd['Date'] = pepl_pd['Date'].astype('datetime64[ns]')\n",
    "\n",
    "# Create numpy arrays for cash flow calcs\n",
    "strip_price = pd.merge(dates, strip_price, how = 'left', on = 'Date')\n",
    "strip_price = pd.merge(strip_price, pepl_pd, how = 'left', on = 'Date')\n",
    "strip_price.fillna(method = 'ffill', inplace = True)\n",
    "oil_price_full = np.transpose(strip_price[['OilPrice']].to_numpy())[0]\n",
    "gas_price_full = np.transpose(strip_price[['GasPrice']].to_numpy())[0]\n",
    "pepl_diff_full = np.transpose(strip_price[['paj_gas']].to_numpy())[0]"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "source": [
    "# Generate price and volume arrays with proper indexing\r\n",
    "\r\n",
    "def oil_price(x):\r\n",
    "    # Create shift integer for arrays\r\n",
    "    StartDate = prop_list.loc[x, 'Start_Date']\r\n",
    "    DateDiff = MonthDiff(BaseDate, StartDate)\r\n",
    "    \r\n",
    "    # Shift price arrays\r\n",
    "    uid = np.full((duration + DateDiff,), oil_nparr[x][0][0])\r\n",
    "    oil_pri = oil_price_full[:duration + DateDiff]\r\n",
    "    oil_price = np.vstack((uid, oil_pri))\r\n",
    "    return oil_price\r\n",
    "\r\n",
    "voil_price = np.vectorize(oil_price, otypes = [object])                         \r\n",
    "oilprice = voil_price(R)\r\n",
    "\r\n",
    "def gas_price(x):\r\n",
    "    # Create shift integer for arrays\r\n",
    "    StartDate = prop_list.loc[x, 'Start_Date']\r\n",
    "    DateDiff = MonthDiff(BaseDate, StartDate)\r\n",
    "    \r\n",
    "    # Shift price arrays\r\n",
    "    uid = np.full((duration + DateDiff,), gas_nparr[x][0][0])\r\n",
    "    gas_pri = gas_price_full[:duration + DateDiff]\r\n",
    "    gas_price = np.vstack((uid, gas_pri))\r\n",
    "    return gas_price\r\n",
    "\r\n",
    "vgas_price = np.vectorize(gas_price, otypes = [object])                         \r\n",
    "gasprice = vgas_price(R)\r\n",
    "\r\n",
    "def gas_diff(x):\r\n",
    "    # Create shift integer for arrays\r\n",
    "    StartDate = prop_list.loc[x, 'Start_Date']\r\n",
    "    DateDiff = MonthDiff(BaseDate, StartDate)\r\n",
    "    \r\n",
    "    # Shift price arrays\r\n",
    "    uid = np.full((duration + DateDiff,), gas_nparr[x][0][0])\r\n",
    "    pepl = pepl_diff_full[:duration + DateDiff]\r\n",
    "    gas_diff = np.vstack((uid, pepl))\r\n",
    "    return gas_diff\r\n",
    "\r\n",
    "vgas_diff = np.vectorize(gas_diff, otypes = [object])                         \r\n",
    "gasdiff = vgas_diff(R)\r\n",
    "\r\n",
    "def volarray(x):\r\n",
    "    uid = oil_nparr[x][0]\r\n",
    "    StartDate = prop_list.loc[x, 'Start_Date']\r\n",
    "    DateDiff = MonthDiff(BaseDate, StartDate)\r\n",
    "    \r\n",
    "    # create well specific volume arrays\r\n",
    "    month = np.round(oil_nparr[x][2] + DateDiff - 1, 0)\r\n",
    "    oil_vol = np.round(oil_nparr[x][6], 4)\r\n",
    "    gas_vol = np.round(gas_nparr[x][6], 4)\r\n",
    "    nvol_np = np.vstack((uid, month, oil_vol, gas_vol))\r\n",
    "    prop = np.full((DateDiff,), uid[0])\r\n",
    "    delay = np.arange(DateDiff)\r\n",
    "    oil_zeros = np.zeros((DateDiff),)\r\n",
    "    gas_zeros = oil_zeros\r\n",
    "    shift = np.vstack((prop, delay, oil_zeros, gas_zeros))\r\n",
    "    vol_arr = np.column_stack((shift, nvol_np))\r\n",
    "    \r\n",
    "    return vol_arr\r\n",
    "\r\n",
    "vvolarray = np.vectorize(volarray, otypes = [object])\r\n",
    "vol_np = vvolarray(R)"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "source": [
    "# Generate full monthly cash flow arrays\r\n",
    "    \r\n",
    "# define constant input parameters\r\n",
    "eloss = 0\r\n",
    "weight = 1.0\r\n",
    "prod_wt = 1.0\r\n",
    "inv_wt = 1.0\r\n",
    "stx_oil = 0.0795\r\n",
    "stx_gas = 0.0795\r\n",
    "stx_ngl = 0.0795\r\n",
    "adval = 0\r\n",
    "aban = 150000\r\n",
    "\r\n",
    "# Create function for slicing the volume array and calculating the monthly cash flow\r\n",
    "def econ_ncf_iter(r):    \r\n",
    "    econ_ncf_iter = econ_cf(index = r, uid = prop_list.loc[r, 'PROPNUM'], wi = prop_list.loc[r, 'WI'], \r\n",
    "                            nri = prop_list.loc[r, 'NRI'], roy = prop_list.loc[r, 'Royalty'], eloss = eloss, \r\n",
    "                            weight = weight, prod_wt = prod_wt, inv_wt = inv_wt, \r\n",
    "                            shrink = np.round(prop_list.loc[r, 'SHRINK'] / 100, 6), \r\n",
    "                            btu = np.round(prop_list.loc[r, 'BTU'] / 1000, 6), \r\n",
    "                            ngl_yield = np.round(prop_list.loc[r, 'NGL/GAS'], 6), \r\n",
    "                            pri_oil = np.extract(oilprice[r][0] == prop_list.loc[r, 'PROPNUM'], oilprice[r][1]),\r\n",
    "                            pri_gas = np.extract(gasprice[r][0] == prop_list.loc[r, 'PROPNUM'], gasprice[r][1]),\r\n",
    "                            paj_oil = prop_list.loc[r, 'PAJ_OIL'], \r\n",
    "                            paj_gas = np.extract(gasdiff[r][0] == prop_list.loc[r, 'PROPNUM'], gasdiff[r][1]), \r\n",
    "                            paj_ngl = prop_list.loc[r, 'PAJ_NGL'], stx_oil = stx_oil, stx_gas = stx_gas, stx_ngl = stx_ngl,\r\n",
    "                            adval = adval, opc_fix = np.round(prop_list.loc[r, 'OPC/T'], 2), \r\n",
    "                            opc_oil = np.round(prop_list.loc[r, 'OIL_OPEX'], 2), \r\n",
    "                            opc_gas = np.round(prop_list.loc[r, 'GAS_OPEX'], 2), \r\n",
    "                            capex = np.round(prop_list.loc[r, 'CAPITAL'] * 1000, 2), aban = aban)\r\n",
    "    return econ_ncf_iter\r\n",
    "\r\n",
    "# generate net cash flow array\r\n",
    "econ_ncf = lambda r: econ_ncf_iter(r)\r\n",
    "vecon_ncf = np.vectorize(econ_ncf_iter, otypes = [object])\r\n",
    "ncf_arr_packed = vecon_ncf(R)"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "source": [
    "ncf_pd_dflist = []\r\n",
    "columns = ['UID', 'Month', 'Grs Oil', 'Grs Gas', 'Net Oil', 'Net Gas', 'Net NGL', 'Oil Revenue', 'Gas Revenue', \r\n",
    "           'NGL Revenue', 'Total Revenue', 'Total Tax', 'OPEX', 'Operating Income', 'Cumulative Op CF', 'Net Cashflow',\r\n",
    "           'Cumulative Net CF']\r\n",
    "\r\n",
    "pbar = tqdm(len(R))\r\n",
    "for r in R:\r\n",
    "    ncf_pd_dflist.append(pd.DataFrame(np.transpose(ncf_arr_packed[r])))\r\n",
    "    pbar.update()\r\n",
    "ncf_pd = pd.concat(ncf_pd_dflist)\r\n",
    "ncf_pd.columns = columns\r\n",
    "pbar.close()"
   ],
   "outputs": [
    {
     "output_type": "stream",
     "name": "stderr",
     "text": [
      "199it [00:00, 13259.20it/s]\n"
     ]
    }
   ],
   "source": [
    "ncf_pd_dflist = []\n",
    "columns = ['UID', 'Month', 'Grs Oil', 'Grs Gas', 'Net Oil', 'Net Gas', 'Net NGL', 'Oil Revenue', 'Gas Revenue', \n",
    "           'NGL Revenue', 'Total Revenue', 'Total Tax', 'OPEX', 'Operating Income', 'Cumulative Op CF', 'Net Cashflow',\n",
    "           'Cumulative Net CF']\n",
    "\n",
    "for r in R:\n",
    "    ncf_pd_dflist.append(pd.DataFrame(np.transpose(ncf_arr_packed[r])))\n",
    "ncf_pd = pd.concat(ncf_pd_dflist)\n",
    "ncf_pd.columns = columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "source": [
    "# Truncate pandas dataframe\r\n",
    "monthly_out = 360 # months of cash flow output to export\r\n",
    "ncf_pd_trunc = ncf_pd.query(f'Month < {monthly_out}')\r\n",
    "ncf_pd_trunc.to_csv('monthly_ncf.csv')"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "source": [
    "# loop though wells in array and create oneline output of economic metrics\r\n",
    "oneline_cat = prop_list.iloc[:, :9]\r\n",
    "disc_rate = [0.05, 0.08, 0.09, 0.1, 0.12, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]\r\n",
    "propnum = []\r\n",
    "ror = []\r\n",
    "life = []\r\n",
    "payout = []\r\n",
    "roi = []\r\n",
    "droi = []\r\n",
    "result_nparr = np.empty([6, 0])\r\n",
    "npv_list = []\r\n",
    "\r\n",
    "pbar = tqdm(len(R))\r\n",
    "for r in R:\r\n",
    "    propID = prop_list.loc[r, 'PROPNUM']\r\n",
    "    wi = prop_list.loc[r, 'WI']\r\n",
    "    capex = prop_list.loc[r, 'CAPITAL']\r\n",
    "    ncf_r = np.extract(ncf_arr_packed[r][0] == propID, ncf_arr_packed[r][15])\r\n",
    "    month_r = np.extract(ncf_arr_packed[r][0] == propID, ncf_arr_packed[r][1])\r\n",
    "    cum_ocf_r = np.extract(ncf_arr_packed[r][0] == propID, ncf_arr_packed[r][14])\r\n",
    "    cum_ncf_r = np.extract(ncf_arr_packed[r][0] == propID, ncf_arr_packed[r][14])\r\n",
    "    \r\n",
    "    # calculate npv at all discount rates\r\n",
    "    pv_calc = lambda i: npv(i, ncf_r, month_r.astype(int))\r\n",
    "    pv = list(map(pv_calc, disc_rate))\r\n",
    "    npv_list.append(pv)\r\n",
    "    \r\n",
    "    # calculate irr\r\n",
    "    if wi == 0 or capex == 0:\r\n",
    "        irr = max(max(disc_rate), 1)\r\n",
    "    else:\r\n",
    "        f = lambda x: npv(x, ncf_r, month_r)\r\n",
    "        r = optimize.root(f, [0])\r\n",
    "        irr = min(r.x[0], 1)\r\n",
    "    ror.append(irr)\r\n",
    "    \r\n",
    "    # calculate life\r\n",
    "    try:\r\n",
    "        life_calc = np.where(cum_ocf_r == np.max(cum_ocf_r))[0][0] + eloss\r\n",
    "    except:\r\n",
    "        life_calc = 1 + eloss\r\n",
    "    life.append(np.round(life_calc / 12, 2))\r\n",
    "        \r\n",
    "    # calculate payout\r\n",
    "    ncf_cum = cum_ncf_r[:life_calc]\r\n",
    "    month_arr = month_r[:life_calc]\r\n",
    "    try:\r\n",
    "        payout_interp = interp1d(ncf_cum, month_arr, kind = 'cubic')\r\n",
    "        payout_calc = payout_interp(0) / 12\r\n",
    "    except:\r\n",
    "        payout_calc = life_calc / 12\r\n",
    "    payout.append(np.round(payout_calc, 2))\r\n",
    "        \r\n",
    "    # calculate ROI and DROI\r\n",
    "    net_capex = capex * weight * inv_wt * wi\r\n",
    "    if net_capex == 0:\r\n",
    "        roi_calc = 0\r\n",
    "        droi_calc = 0\r\n",
    "    else:\r\n",
    "        roi_calc = np.round(((net_capex + np.sum(ncf_r)) / net_capex), 2)\r\n",
    "        droi_calc = np.round(((net_capex + pv[3]) / net_capex), 2)\r\n",
    "    roi.append(roi_calc)\r\n",
    "    droi.append(droi_calc)\r\n",
    "       \r\n",
    "    propnum.append(propID)\r\n",
    "    oneline_nparr = np.array((propnum, ror, life, payout, roi, droi))\r\n",
    "    pbar.update()\r\n",
    "\r\n",
    "result_nparr = np.column_stack((result_nparr, oneline_nparr))\r\n",
    "npv_pd = pd.DataFrame(npv_list, columns = disc_rate)\r\n",
    "pbar.close()"
   ],
   "outputs": [
    {
     "output_type": "stream",
     "name": "stderr",
     "text": [
      "199it [00:00, 1816.34it/s]\n"
     ]
    }
   ],
   "source": [
    "# loop though wells in array and create oneline output of economic metrics\n",
    "oneline_cat = prop_list.iloc[:, :9]\n",
    "disc_rate = [0.05, 0.08, 0.09, 0.1, 0.12, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]\n",
    "propnum = []\n",
    "ror = []\n",
    "life = []\n",
    "payout = []\n",
    "roi = []\n",
    "droi = []\n",
    "result_nparr = np.empty([6, 0])\n",
    "npv_list = []\n",
    "\n",
    "for r in R:\n",
    "    propID = prop_list.loc[r, 'PROPNUM']\n",
    "    wi = prop_list.loc[r, 'WI']\n",
    "    capex = prop_list.loc[r, 'CAPITAL']\n",
    "    ncf_r = np.extract(ncf_arr_packed[r][0] == propID, ncf_arr_packed[r][15])\n",
    "    month_r = np.extract(ncf_arr_packed[r][0] == propID, ncf_arr_packed[r][1])\n",
    "    cum_ocf_r = np.extract(ncf_arr_packed[r][0] == propID, ncf_arr_packed[r][14])\n",
    "    cum_ncf_r = np.extract(ncf_arr_packed[r][0] == propID, ncf_arr_packed[r][14])\n",
    "    \n",
    "    # calculate npv at all discount rates\n",
    "    pv_calc = lambda i: npv(i, ncf_r, month_r.astype(int))\n",
    "    pv = list(map(pv_calc, disc_rate))\n",
    "    npv_list.append(pv)\n",
    "    \n",
    "    # calculate irr\n",
    "    if wi == 0 or capex == 0:\n",
    "        irr = max(max(disc_rate), 1)\n",
    "    else:\n",
    "        f = lambda x: npv(x, ncf_r, month_r)\n",
    "        r = optimize.root(f, [0])\n",
    "        irr = min(r.x[0], 1)\n",
    "    ror.append(irr)\n",
    "    \n",
    "    # calculate life\n",
    "    try:\n",
    "        life_calc = np.where(cum_ocf_r == np.max(cum_ocf_r))[0][0] + eloss\n",
    "    except:\n",
    "        life_calc = 1 + eloss\n",
    "    life.append(np.round(life_calc / 12, 2))\n",
    "        \n",
    "    # calculate payout\n",
    "    ncf_cum = cum_ncf_r[:life_calc]\n",
    "    month_arr = month_r[:life_calc]\n",
    "    try:\n",
    "        payout_interp = interp1d(ncf_cum, month_arr, kind = 'cubic')\n",
    "        payout_calc = payout_interp(0) / 12\n",
    "    except:\n",
    "        payout_calc = life_calc / 12\n",
    "    payout.append(np.round(payout_calc, 2))\n",
    "        \n",
    "    # calculate ROI and DROI\n",
    "    net_capex = capex * weight * inv_wt * wi\n",
    "    if net_capex == 0:\n",
    "        roi_calc = 0\n",
    "        droi_calc = 0\n",
    "    else:\n",
    "        roi_calc = np.round(((net_capex + np.sum(ncf_r)) / net_capex), 2)\n",
    "        droi_calc = np.round(((net_capex + pv[3]) / net_capex), 2)\n",
    "    roi.append(roi_calc)\n",
    "    droi.append(droi_calc)\n",
    "       \n",
    "    propnum.append(propID)\n",
    "    oneline_nparr = np.array((propnum, ror, life, payout, roi, droi))\n",
    "\n",
    "result_nparr = np.column_stack((result_nparr, oneline_nparr))\n",
    "npv_pd = pd.DataFrame(npv_list, columns = disc_rate)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "source": [
    "# Create oneline output export to csv\n",
    "result_pd = pd.DataFrame(np.transpose(result_nparr), columns = ['UID', 'IRR', 'Life', 'Payout', 'ROI', 'DROI'])\n",
    "result_pd = pd.merge(result_pd, npv_pd, how = 'inner', left_index = True, right_index = True)\n",
    "result_pd.to_csv('multiwell_dcf_oneline.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "source": [
    "# Calculate breakevens and add to pandas output\n",
    "equiv_ratio = 20\n",
    "disc1 = 0.1\n",
    "disc2 = 0.2\n",
    "pri = 50 # breakeven price guess\n",
    "pajgas = -0.25 # gas differential\n",
    "run_breakevens = 1\n",
    "\n",
    "if run_breakevens == 0:\n",
    "    pass\n",
    "\n",
    "else:\n",
    "    # Create function for slicing the volume array and calculating the monthly cash flow\n",
    "    def econ_be(pri, paj_gas, disc, equiv_ratio, r):    \n",
    "        econ_be = econ_cf(index = r, uid = prop_list.loc[r, 'PROPNUM'], wi = prop_list.loc[r, 'WI'], \n",
    "                           nri = prop_list.loc[r, 'NRI'], roy = prop_list.loc[r, 'Royalty'], eloss = eloss, \n",
    "                           weight = weight, prod_wt = prod_wt, inv_wt = inv_wt, \n",
    "                           shrink = np.round(prop_list.loc[r, 'SHRINK'] / 100, 6), \n",
    "                           btu = np.round(prop_list.loc[r, 'BTU'] / 1000, 6), \n",
    "                           ngl_yield = np.round(prop_list.loc[r, 'NGL/GAS'], 6), \n",
    "                           pri_oil = pri, \n",
    "                           pri_gas = pri / equiv_ratio, \n",
    "                           paj_oil = prop_list.loc[r, 'PAJ_OIL'], \n",
    "                           paj_gas = pajgas, \n",
    "                           paj_ngl = prop_list.loc[r, 'PAJ_NGL'], stx_oil = stx_oil, stx_gas = stx_gas, stx_ngl = stx_ngl,\n",
    "                           adval = adval, opc_fix = np.round(prop_list.loc[r, 'OPC/T'], 2), \n",
    "                           opc_oil = np.round(prop_list.loc[r, 'OIL_OPEX'], 2), \n",
    "                           opc_gas = np.round(prop_list.loc[r, 'GAS_OPEX'], 2), \n",
    "                           capex = np.round(prop_list.loc[r, 'CAPITAL'] * 1000, 2), aban = aban)\n",
    "        pv_beo = npv(disc, econ_be[15], econ_be[1])\n",
    "        return pv_beo\n",
    "\n",
    "    # calculate breakeven prices\n",
    "    def econ_be_iter1(r):\n",
    "        pv_beo1 = lambda p: econ_be(p, pajgas, disc1, equiv_ratio, r)\n",
    "        r_beo1 = optimize.root(pv_beo1, pri)\n",
    "        return round(r_beo1.x[0], 2)\n",
    "\n",
    "    vpv_beo1 = np.vectorize(econ_be_iter1)\n",
    "    beo1 = vpv_beo1(R)\n",
    "\n",
    "    def econ_be_iter2(r):\n",
    "        pv_beo2 = lambda p: econ_be(p, pajgas, disc2, equiv_ratio, r)\n",
    "        r_beo2 = optimize.root(pv_beo2, pri)\n",
    "        return round(r_beo2.x[0], 2)\n",
    "\n",
    "    vpv_beo2 = np.vectorize(econ_be_iter2)\n",
    "    beo2 = vpv_beo2(R)\n",
    "\n",
    "    # Add results to oneline output\n",
    "    result_pd['Oil_BE1'] = beo1\n",
    "    result_pd['Oil_BE2'] = beo2\n",
    "    result_pd.to_csv('multiwell_dcf_oneline.csv')"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [],
   "outputs": [],
   "metadata": {}
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "base",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.15"
  },
  "vscode": {
   "interpreter": {
    "hash": "a3ead6ff9e9c1be087d7e40f893bad9c1e8f35d013c1ee3985ae6710c98a7135"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}