Monte Carlo Analysis Python Oil and Gas Volumetrics

Question

I am trying to teach myself python and I wanted to start with learning how to do a monte carlo analysis (I am a scientist by trade who uses MCA alot). I am trying to write a program in that will perform a montecarlo simulation of 7 variables to calculate the range of possible outcomes from a given formula.

I am EXTREMELY new at python. I have experience in VBA, but am still learning python.

All of the examples I find online relate to finance, and unfortunately I am having a hard time canablazing these codes as they are not very relatable to what I am trying to do.

I feel like this should be a very simple thing as the equation is quite simple. However, I cant seem to make any headway. Could someone please take a look at the code I have right now and point me in the right direction? If you have a non finance example of a montecarlo analysis that is easy for a beginner to understand, please point me that way. I want to learn python, but I have found the best way I learn is to look at examples from other people.

What I am trying to accomplish is to have a list of 7 variables each with normal distributions. I would like python to perform 10,000 iterations the formula in the code below, pulling a different set of the 7 different variables each time it calculates. I would eventually like to print the P90, P50, and P10 values of the calculation. At this time I dont really care about looking at any graphs (that will come later).

Below is the code I have come up with so far. Again, I am looking for direction on how to accomplish this. I know my syntax is probably messed up but I am really struggling with this. Any help would greatly be appreciated.

from scipy.stats import *
import numpy as np


n = 10000

for i in range(n):
    Area = norm(200,50)
    Thickness = norm(100,25)
    NTG = norm(.85,.2)
    POR = norm(.32,.02)
    GS = norm(.80,.2)
    BG= norm(.0024,.0001)
    Feather = 1
    return ((((Area*Thickness*NTG*POR*GS)/BG)*43560)*Feather)/1000000000


Result = ((((Area*Thickness*NTG*POR*GS)/BG)*43560)*Feather)/1000000000


print ('Result is ', Result, 'ft')

Answer 1

There are several problems with your code, at the least:

• You cannot place return outside of a function

• Your objects are random number generators, not random numbers; you cannot multiply them

Presumably, you want something like this:

n = 10000
Area = norm(200,50).rvs(n)
Thickness = norm(100,25).rvs(n)
NTG = norm(.85,.2).rvs(n)
POR = norm(.32,.02).rvs(n)
GS = norm(.80,.2).rvs(n)
BG= norm(.0024,.0001).rvs(n)
Feather = 1
Results = Area*Thickness*NTG*POR*GS/BG*43560*Feather/1000000000
print(Results)

Output:

array([  43.88063752,   88.94160248,   46.89368561, ...,   87.32369654,
    210.16899452,   32.21191486])

The Results array could then be used to perform statistical analysis, such as P10, P50, and P90 calculations.

Notice the absence of loops, which you should try to avoid in numerical python. This already created n results.

I must say I think you'd be more productive reading a Python tutorial first.

Answer 2

I don't know what you are trying to accomplish, but I can offer some suggestions:

• It looks like there is no covariance between realizations from your simulation, so to speed up the computation and make the code a little clearer, you can draw the data from a joint multivariate normal distribution.

• I hate for loops, so vectorization is better whenever it can be done.

Here is some code

import numpy as np
import matplotlib.pyplot as plt

N = 10000
mean_vector = np.array([200, #Area
                        100, #Thickness
                        0.85, #NTG  
                        0.32, #POR
                        0.80, #GS
                        0.0024 #BG
                        ])

covariance = np.diag([50,
                      25,
                      0.2,
                      0.02,
                      0.2,
                      0.0001])

simulated = np.random.multivariate_normal(mean_vector, covariance, size = N)

area, thickness, ntg, por, gs, bg = simulated.T

result= (area*thickness*ntg*por*gs*43560)/bg

The result array contains the results of your simulation. Now you can plot a histogram or compute means if needed.