Solution sketch¶

Solutions to EMEA Ch. 6

In [ ]:
import numpy as np
import matplotlib.pyplot as plt
plt.style.use('seaborn-poster')

/var/folders/q_/twndxknx06n5v1nllkjtz0km0000gn/T/ipykernel_5195/2492506434.py:3: MatplotlibDeprecationWarning: The seaborn styles shipped by Matplotlib are deprecated since 3.6, as they no longer correspond to the styles shipped by seaborn. However, they will remain available as 'seaborn-v0_8-<style>'. Alternatively, directly use the seaborn API instead.
  plt.style.use('seaborn-poster')
In [ ]:
# step size
h = 0.1

# define grid
Q = np.arange(.01, 20, h) # why not zero?

pi = np.exp(Q**2/(Q+Q**2))*Q - Q*np.log(Q)
In [ ]:
pi
Out[ ]:
array([ 0.0561512 ,  0.36425957,  0.57753619,  0.75583253,  0.91391813,
        1.05831744,  1.19251765,  1.31859144,  1.43786254,  1.55122399,
        1.65930612,  1.76257168,  1.86137264,  1.95598537,  2.04663334,
        2.13350201,  2.21674897,  2.29651089,  2.37290838,  2.44604955,
        2.51603252,  2.58294729,  2.64687716,  2.70789979,  2.76608802,
        2.82151052,  2.87423227,  2.92431503,  2.97181759,  3.01679611,
        3.05930433,  3.09939374,  3.13711379,  3.172512  ,  3.20563411,
        3.23652416,  3.26522462,  3.29177649,  3.31621932,  3.33859137,
        3.3589296 ,  3.37726978,  3.39364653,  3.40809337,  3.4206428 ,
        3.43132631,  3.44017442,  3.44721676,  3.45248209,  3.45599833,
        3.45779259,  3.45789121,  3.4563198 ,  3.45310326,  3.44826581,
        3.441831  ,  3.43382178,  3.42426045,  3.41316877,  3.4005679 ,
        3.38647849,  3.37092064,  3.35391397,  3.3354776 ,  3.31563017,
        3.2943899 ,  3.27177454,  3.24780143,  3.2224875 ,  3.1958493 ,
        3.16790297,  3.13866431,  3.10814874,  3.07637134,  3.04334687,
        3.00908975,  2.97361409,  2.93693369,  2.89906207,  2.86001244,
        2.81979777,  2.77843072,  2.73592372,  2.69228892,  2.64753825,
        2.60168339,  2.55473578,  2.50670664,  2.45760698,  2.40744758,
        2.35623904,  2.30399172,  2.25071582,  2.19642133,  2.14111806,
        2.08481564,  2.02752353,  1.969251  ,  1.91000717,  1.849801  ,
        1.78864128,  1.72653666,  1.66349562,  1.59952651,  1.53463754,
        1.46883676,  1.4021321 ,  1.33453137,  1.26604222,  1.19667219,
        1.1264287 ,  1.05531905,  0.98335042,  0.91052986,  0.83686434,
        0.76236069,  0.68702565,  0.61086585,  0.53388782,  0.45609799,
        0.37750269,  0.29810816,  0.21792054,  0.13694588,  0.05519014,
       -0.0273408 , -0.11064116, -0.19470522, -0.27952737, -0.36510206,
       -0.45142383, -0.53848729, -0.62628713, -0.7148181 , -0.80407506,
       -0.89405289, -0.98474657, -1.07615116, -1.16826175, -1.26107352,
       -1.35458171, -1.44878162, -1.54366862, -1.63923813, -1.73548563,
       -1.83240666, -1.92999683, -2.02825179, -2.12716724, -2.22673896,
       -2.32696275, -2.42783449, -2.5293501 , -2.63150554, -2.73429684,
       -2.83772006, -2.94177132, -3.04644678, -3.15174264, -3.25765517,
       -3.36418065, -3.47131542, -3.57905588, -3.68739844, -3.79633957,
       -3.90587577, -4.0160036 , -4.12671964, -4.23802051, -4.34990288,
       -4.46236344, -4.57539893, -4.68900613, -4.80318184, -4.9179229 ,
       -5.03322619, -5.14908862, -5.26550715, -5.38247873, -5.50000038,
       -5.61806915, -5.73668211, -5.85583635, -5.97552901, -6.09575726,
       -6.21651828, -6.33780929, -6.45962756, -6.58197035, -6.70483496,
       -6.82821874, -6.95211904, -7.07653324, -7.20145877, -7.32689304,
       -7.45283353, -7.57927773, -7.70622314, -7.8336673 , -7.96160778])
In [ ]:
fig, ax = plt.subplots()
ax.plot(Q,pi)
Out[ ]:
[<matplotlib.lines.Line2D at 0x7fa779e82fd0>]
In [ ]:
dpi = np.diff(pi)/h 
Q_diff = Q[:-1]
In [ ]:
fig, ax = plt.subplots()
ax.plot(Q_diff,dpi)
ax.axhline(0, color="r")
Out[ ]:
<matplotlib.lines.Line2D at 0x7fa768ebea00>
In [ ]:
absMin = np.min(np.abs(dpi))
absMin
Out[ ]:
0.0009861971553171145
In [ ]:
Q_diff[dpi==absMin]
Out[ ]:
array([5.01])

Solution, Problem 4¶

In [ ]:
rng = np.random.default_rng()
x = np.linspace(10,100, 100)
#x = rng.normal(50,20, 100)  #let x be the price of oil 
beta = .02 #(since this is a small value) we can interpret
#roughly as that if the oil price increases by 1, t
# then the probability of drilling will increase by about 2%

thetas = np.exp(beta*x)
drills = rng.poisson(thetas)
drills
Out[ ]:
array([ 2,  1,  0,  3,  0,  3,  1,  0,  2,  3,  0,  2,  0,  2,  4,  1,  0,
        1,  1,  1,  1,  1,  0,  4,  3,  1,  1,  2,  5,  3,  2,  5,  1,  0,
        4,  2,  2,  3,  3,  3,  0,  2,  4,  3,  1,  2,  2,  0,  1,  2,  3,
        6,  3,  7,  6,  5,  6,  7,  2,  3,  3,  3,  4,  5,  4,  5,  7,  6,
        6,  3,  6,  3,  3,  5, 10,  2,  3,  7,  7,  5,  3, 11,  3,  3,  6,
        5,  3,  7,  9,  7,  4,  9,  4,  5,  2,  8,  6,  7,  9,  7])
In [ ]:
In [ ]:
# Define a poisson log-likelihood
from scipy.stats import poisson

def PoissonloglikFunction(Y, X, b):
    theta_hats = np.exp(b*X)
    loglik = np.sum(np.log(poisson.pmf(Y, theta_hats)))
    return loglik
In [ ]:
import scipy.optimize as opt
result = opt.minimize_scalar(lambda b: -PoissonloglikFunction(drills, x, b), bounds=(0, .05), method='bounded')
print(result)

     fun: 194.21353099501928
 message: 'Solution found.'
    nfev: 10
     nit: 10
  status: 0
 success: True
       x: 0.020299127628022103
In [ ]: