Certain Equivalents and Sensitivity Analysis

By: Chengyi (Jeff) Chen

%load_ext autotime
%load_ext nb_black
%matplotlib notebook
%config InlineBackend.figure_format = 'retina'

# Graphing
import matplotlib.pyplot as plt
import seaborn as sns

# Plotting Design Configs
sns.set(
    font="Verdana",
    rc={
        "axes.axisbelow": False,
        "axes.edgecolor": "lightgrey",
        "axes.facecolor": "None",
        "axes.grid": True,
        "grid.color": "lightgrey",
        "axes.labelcolor": "dimgrey",
        "axes.spines.right": False,
        "axes.spines.top": False,
        "figure.facecolor": "white",
        "lines.solid_capstyle": "round",
        "patch.edgecolor": "w",
        "patch.force_edgecolor": True,
        "text.color": "dimgrey",
        "xtick.bottom": False,
        "xtick.color": "dimgrey",
        "xtick.direction": "out",
        "xtick.top": False,
        "ytick.color": "dimgrey",
        "ytick.direction": "out",
        "ytick.left": False,
        "ytick.right": False,
        "figure.dpi": 100,
        "figure.figsize": (8, 6),
    },
)
sns.set_context(
    "notebook", rc={"font.size": 10, "axes.titlesize": 12, "axes.labelsize": 10}
)
sns.color_palette(palette="Spectral")

# General
import numpy as np
import scipy.stats as sp
import pandas as pd
from functools import partial
time: 864 ms (started: 2021-03-26 18:57:59 +08:00)
def ρ(γ):
    """Risk tolerance"""
    return 1 / γ
time: 294 µs (started: 2021-03-26 18:57:59 +08:00)
def r(γ):
    """Risk Odds"""
    return np.log(γ)
time: 280 µs (started: 2021-03-26 18:57:59 +08:00)
def u(x, γ: np.array = np.array([0.0]), a: float = 0.0, b: float = 1.0):
    """Assuming user satisifies the Δ property, 
    calculates the U-values using either a 
    Piecewise Linear u = a + bx if γ = 0 (Risk-Neutral) 
    else Exponential u = a + be^(-xγ) U-curve U(x)
    
    Args:
        x (np.array): Payoffs of prospects matrix, Shape = (Number of deals, Number of prospects in each deal)
        γ (np.array): Risk-aversion coefficients, Shape = (Number of different risk-aversion coefficients for sensitivity analysis), Default = 0.0 (Risk-neutral) [γ > 0: Risk-averse, γ < 0: Risk-seeking]
        a (float): Constant for U-curve
        b (float): Coefficient for payoff variable
        
    Returns:
        np.array:
            U-values, Shape = (Number of Risk Aversion Coefficients γ, Number of deals, Number of prospects in each deal)
    """
    assert (
        x.ndim == 2
    ), "Payoffs require 2 dimensions, first dim is number of deals, second is number of prospects in each deal."
    γ = np.array([γ]) if np.isscalar(γ) else np.array(γ)
    return np.array(
        [
            np.apply_along_axis(
                func1d=lambda x, γ=γ: a + b * x if γ == 0 else a + b * np.exp(-x * γ),
                axis=-1,
                arr=x,
                γ=γ_i,
            )
            for γ_i in γ
        ]
    )
time: 597 µs (started: 2021-03-26 18:57:59 +08:00)
def eu(u, p):
    """Assuming user satisifies the Δ property, calculates the 
    hadamard product of the u values matrix and probabilities 
    matrix (respective probabilities associated with each 
    prospect in the u matrix)
    
    Args:
        u (np.array): U-values matrix, Shape = (Number of Risk Aversion Coefficients γ, Number of deals, Number of prospects in each deal)
        p (np.array): Probabilities of each prospect matrix, Shape = (Number of Risk Aversion Coefficients γ, Number of deals, Number of prospects in each deal)
        
    Returns:
        np.array:
            E-values, Shape = (Number of Risk Aversion Coefficients γ, Number of deals) 
    """
    assert (
        u.shape == p.shape
    ), "U-values matrix must be the same shape as the probabilities matrix."

    assert (
        u.ndim == 3
    ), "Both matrices must have 3 dimensions, Shape = (Number of Risk Aversion Coefficients γ, Number of deals, Number of prospects in each deal)."

    return np.sum(u * p, axis=-1)
time: 473 µs (started: 2021-03-26 18:57:59 +08:00)
def u_inv(eu, γ: np.array = np.array([0.0]), a: float = 0, b: float = 1):
    """Piecewise Inverse of Linear if γ = 0 (Risk-Neutral) 
    else Inverse of Exponential U-curve Certain equivalent 
    function U^{-1}(x)
    
    Args:
        eu (np.array): E[U-Values] / Expectation over u-values, AKA E-values, Shape = (Number of Risk Aversion Coefficients γ, Number of deals)
        γ (np.array): Risk-aversion coefficients, Shape = (Number of different risk-aversion coefficients for sensitivity analysis), Default = 0.0 (Risk-neutral) [γ > 0: Risk-averse, γ < 0: Risk-seeking]
        a (float): Constant for U-curve
        b (float): Coefficient for payoff variable
        
    Returns:
        np.array:
             Certainty Equivalent values, Shape = (Number of Risk Aversion Coefficients γ, Number of deals)
    """
    γ = np.array([γ]) if np.isscalar(γ) else np.array(γ)

    for eu_i, γ_i in zip(eu, γ):
        if not np.isclose(γ_i, 0):
            assert np.alltrue(
                eu_i > 0
            ), "E[U-Values] / Expectation over u-values must be positive for γ > 0 and γ < 0 in order for inverse of Exponential U-curve to work."

    assert (
        eu.shape[0] == γ.shape[0]
    ), "E-values first dimension must be the same as γ first dimension for inverse operation."

    return np.array(
        [
            list(
                map(
                    lambda eu, γ=γ_i: (eu - a) / b if γ == 0 else -(1 / γ) * np.log(eu),
                    eu_i,
                )
            )
            for eu_i, γ_i in zip(eu, γ)
        ]
    )
time: 819 µs (started: 2021-03-26 18:58:00 +08:00)
def certainty_equivalent_values_calculator(
    u,
    u_inv,
    x: np.array = None,
    p: np.array = None,
    N: int = 10,
    payout_lb: float = 0.0,
    payout_ub: float = 100.0,
    γ: np.array = np.array([0.0]),
):
    """Assuming user satisifies the Δ property, calculates the Certainty Equivalent values for a given 
    payoff matrix `x` and probability matrix `p`. If both matrices
    are not provided, random matrices will be sampled to simulate
    calculations.

    Args:
        u (function): U-curve function
        u_inv (function): Inverse U-curve function
        x (np.array): Payoffs of prospects matrix, Shape = (Number of deals, Number of prospects in each deal)
        p (np.array): Probabilities of each prospect matrix, Shape = (Number of deals, Number of prospects in each deal)
        N (int): Number of prospects
        payout_lb (float): Lower Bound of Payout for simulation, Default = 0.0
        payout_ub (float): Upper Bound of Payout for simulation, Default = 100.0
        γ (np.array): Risk-aversion coefficients, Shape = (Number of different risk-aversion coefficients for sensitivity analysis), Default = 0.0 (Risk-neutral) [γ > 0: Risk-averse, γ < 0: Risk-seeking]
        
    Returns:
        np.array:
            Certainty Equivalent values, Shape = (Number of Risk Aversion Coefficients γ, Number of deals)
    """
    γ = np.array([γ]) if np.isscalar(γ) else np.array(γ)

    # Case 1: Simulator - Both x and p arent provided
    if p is None and x is None:

        # Preferential probabilities for each prospect
        p = np.random.dirichlet(np.ones(N), size=1)

        # Payouts for each prospect
        x = np.expand_dims(np.random.randint(payout_lb, payout_ub, size=N), axis=0)

    # Case 2: Simulator - Only payouts provided
    elif p is None:

        # Payouts for each prospect
        x = np.expand_dims(
            np.random.randint(payout_lb, payout_ub, size=p.shape[-1]), axis=0
        )

    # Case 3: Simulator - Only probabilities provided
    elif x is None:

        # Preferential probabilities for each prospect
        p = np.random.dirichlet(np.ones(x.shape[-1]), size=1)

    # Case 4: Calculator - Both are provided
    else:

        pass

    # Check that payoffs and probability assignments are the same shape
    assert (
        x.shape == p.shape
    ), "Payoffs `x` and Probabilities `p` must be the same shape=(Number of deals, Number of prospects in each deal)."

    # Calculates U-values Shape = (Number of Risk Aversion Coefficients γ, Number of deals, Number of prospects in each deal)
    u_values = u(x, γ=γ)

    # Reshape p to match U-values shape
    p_reshaped = np.array([p for _ in range(u_values.shape[0])])

    # Expectation of U values Shape = (Number of Risk Aversion Coefficients γ, Number of deals)
    eu_values = eu(u=u_values, p=p_reshaped)

    # Certainty equivalent value of deal
    ce_values = u_inv(eu_values, γ=γ)

    return ce_values
time: 958 µs (started: 2021-03-26 18:58:00 +08:00)

1) Simple Certain Equivalent Calculator

Build a certain equivalent calculator for a decision maker with an exponential u-curve facing a deal with 5 possible prospects with pay-offs. Build a model that allows a sensitivity analysis to the risk aversion coefficient.

# Number of prospects
N = 5

# Risk-aversion coefficient
γ = np.linspace(-1, 1, num=100)

fig, ax = plt.subplots(1, 1)
ax.plot(
    γ,
    certainty_equivalent_values_calculator(
        u, u_inv, N=5, payout_lb=0, payout_ub=100, γ=γ
    ),
)
ax.set_xlabel("γ")
ax.set_ylabel("CE")
ax.set_title("CE(Decision Maker)")
plt.show()
time: 35.4 ms (started: 2021-03-26 18:58:00 +08:00)

2) 1000-degree Certain Equivalent Calculator

Build a certain equivalent calculator for a decision maker with an exponential u-curve facing a deal represented by a

  • Discretized Uniform distribution on the range 1 ~ 1000

# Number of prospects
N = 100

# Payouts for each prospect
x = np.expand_dims(
    np.array([sp.uniform.rvs(loc=0, scale=1000) for _ in range(N)]), axis=0
)

# Preferential probabilities for each prospect
p = np.array([sp.uniform.pdf(x=rv, loc=0, scale=1000) for rv in x])

# Risk-aversion coefficient
γ = np.linspace(-1, 1, num=100)


fig, ax = plt.subplots(1, 1)
ax.plot(
    γ,
    certainty_equivalent_values_calculator(
        u, u_inv, x, p, N, payout_lb=None, payout_ub=None, γ=γ
    ),
)

ax.set_xlabel("γ")
ax.set_ylabel("CE")
ax.set_title(r"CE(Decision Maker) with $p(X) \sim U(0, 1000)$")
plt.show()
time: 34.4 ms (started: 2021-03-26 18:58:00 +08:00)
<ipython-input-4-9a6997586247>:24: RuntimeWarning: overflow encountered in exp
  func1d=lambda x, γ=γ: a + b * x if γ == 0 else a + b * np.exp(-x * γ),
  • Discretized Beta distribution on the range 1 ~ 1000

# Number of prospects
N = 100

# Preferential probabilities for each prospect
x = np.expand_dims(
    np.array([sp.beta.rvs(a=2, b=2, loc=0, scale=1000) for _ in range(N)]), axis=0
)

# Payouts for each prospect
p = np.array([sp.beta.pdf(x=rv, a=2, b=2, loc=0, scale=1000) for rv in x])

# Risk-aversion coefficient
γ = np.linspace(-1, 1, num=100)


fig, ax = plt.subplots(1, 1)
ax.plot(
    γ,
    certainty_equivalent_values_calculator(
        u, u_inv, x, p, N, payout_lb=None, payout_ub=None, γ=γ
    ),
)

ax.set_xlabel("γ")
ax.set_ylabel("CE")
ax.set_title(r"CE(Decision Maker) with $p(X) \sim Beta(a=2, b=2, 0, 1000)$")
plt.show()
time: 40.4 ms (started: 2021-03-26 18:58:00 +08:00)
<ipython-input-4-9a6997586247>:24: RuntimeWarning: overflow encountered in exp
  func1d=lambda x, γ=γ: a + b * x if γ == 0 else a + b * np.exp(-x * γ),

3) Heat Maps: Probability vs Severity Plots

  1. Assuming a Risk-Neutral decision maker, and starting with a deal that has a 0.9 probability of losing \(\$1\) Million and a 0.1 probability of a zero pay-off, find the value of \(X\) in the following deals that makes you indifferent to this deal.

Let \(u(x, \gamma) = e^{-x\gamma} \) and \(u^{-1}(u, \gamma) = x = -\frac{1}{\gamma}\ln u\) be the \(u\)-value function and it’s inverse.

Let \(N, M\) be the number of prospects in deal \(A\) and \(B\) respectively.

Let \(p^{(A)} \in \mathbb{R}^N: \mathbb{1}^\top p^{(A)} = 1\) and \(p^{(B)} \in \mathbb{R}^M: \mathbb{1}^\top p^{(B)} = 1\) be the vector of probabilities of getting different prospects \(a\), \(b\), for deals \(A\) and \(B\) respectively:

(32)\[\begin{align} u(a, \gamma)^\top p^{(A)} &= u(b, \gamma)^\top p^{(B)} \\ b_j &= u^{-1}\left(\frac{1}{p^{(B)}_j}\left[u(a, \gamma)^\top p^{(A)} - \sum^M_{i, i \not= j} u(b_i, \gamma)p^{(B)}_i\right], \gamma\right) \\ \end{align}\]
b_j = lambda a, p_A, b_excl_j, p_B_excl_j, p_B_j, γ=0: u_inv(
    (np.sum(u(a, γ=γ) * p_A, axis=-1) - np.sum(u(b_excl_j, γ=γ) * p_B_excl_j, axis=-1))
    / p_B_j,
    γ=γ,
)
time: 1.25 ms (started: 2021-03-26 18:58:00 +08:00)

a) A Deal with 0.8 Probability of \(X\) and 0.2 Probability of Zero

b_j(
    a=np.array([[-1e6, 0]]),
    p_A=np.array([[0.9, 0.1]]),
    b_excl_j=np.array([[0]]),
    p_B_excl_j=np.array([[0.2]]),
    p_B_j=np.array([[0.8]]),
)[0][0]
-1125000.0
time: 3.24 ms (started: 2021-03-26 18:58:00 +08:00)

b) A Deal with 0.7 Probability of \(X\) and 0.3 Probability of Zero

b_j(
    a=np.array([[-1e6, 0]]),
    p_A=np.array([[0.9, 0.1]]),
    b_excl_j=np.array([[0]]),
    p_B_excl_j=np.array([[0.3]]),
    p_B_j=np.array([[0.7]]),
)[0][0]
-1285714.2857142857
time: 3.54 ms (started: 2021-03-26 18:58:00 +08:00)

c) A Deal with 0.6 Probability of \(X\) and 0.4 Probability of Zero

b_j(
    a=np.array([[-1e6, 0]]),
    p_A=np.array([[0.9, 0.1]]),
    b_excl_j=np.array([[0]]),
    p_B_excl_j=np.array([[0.4]]),
    p_B_j=np.array([[0.6]]),
)[0][0]
-1500000.0
time: 3.91 ms (started: 2021-03-26 18:58:00 +08:00)

d) A Deal with 0.5 Probability of \(X\) and 0.5 Probability of Zero

b_j(
    a=np.array([[-1e6, 0]]),
    p_A=np.array([[0.9, 0.1]]),
    b_excl_j=np.array([[0]]),
    p_B_excl_j=np.array([[0.5]]),
    p_B_j=np.array([[0.5]]),
)[0][0]
-1800000.0
time: 3.37 ms (started: 2021-03-26 18:58:00 +08:00)

e) A Deal with 0.4 Probability of \(X\) and 0.6 Probability of Zero

b_j(
    a=np.array([[-1e6, 0]]),
    p_A=np.array([[0.9, 0.1]]),
    b_excl_j=np.array([[0]]),
    p_B_excl_j=np.array([[0.6]]),
    p_B_j=np.array([[0.4]]),
)[0][0]
-2250000.0
time: 2.96 ms (started: 2021-03-26 18:58:00 +08:00)

f) A Deal with 0.3 Probability of \(X\) and 0.7 Probability of Zero

b_j(
    a=np.array([[-1e6, 0]]),
    p_A=np.array([[0.9, 0.1]]),
    b_excl_j=np.array([[0]]),
    p_B_excl_j=np.array([[0.7]]),
    p_B_j=np.array([[0.3]]),
)[0][0]
-3000000.0
time: 3.16 ms (started: 2021-03-26 18:58:00 +08:00)

g) A Deal with 0.2 Probability of \(X\) and 0.8 Probability of Zero

b_j(
    a=np.array([[-1e6, 0]]),
    p_A=np.array([[0.9, 0.1]]),
    b_excl_j=np.array([[0]]),
    p_B_excl_j=np.array([[0.8]]),
    p_B_j=np.array([[0.2]]),
)[0][0]
-4500000.0
time: 3.17 ms (started: 2021-03-26 18:58:00 +08:00)

h) A Deal with 0.1 Probability of \(X\) and 0.9 Probability of Zero

b_j(
    a=np.array([[-1e6, 0]]),
    p_A=np.array([[0.9, 0.1]]),
    b_excl_j=np.array([[0]]),
    p_B_excl_j=np.array([[0.9]]),
    p_B_j=np.array([[0.1]]),
)[0][0]
-9000000.0
time: 2.98 ms (started: 2021-03-26 18:58:00 +08:00)

Plot the (Probability of \(X\)) vs. \(X\).

fig, ax = plt.subplots(1, 1)
ax.plot(
    np.array(
        [
            b_j(
                a=np.array([[-1e6, 0]]),
                p_A=np.array([[0.9, 0.1]]),
                b_excl_j=np.array([[0]]),
                p_B_excl_j=np.array([[1 - p_B_j]]),
                p_B_j=np.array([[p_B_j]]),
            )[0][0]
            for p_B_j in np.linspace(0.01, 0.99, num=100)
        ]
    ),
    np.linspace(0.01, 0.99, num=100),
    label="γ = 0",
)
ax.set_xlabel(r"$X$")
ax.set_ylabel(r"$p(X)$")
ax.set_title(r"$p(X)$ Vs. $X$ Heat Map for Risk-Neutral (γ=0) Decision Maker")
plt.show()
time: 55.4 ms (started: 2021-03-26 18:58:00 +08:00)
  1. Repeat for a decision maker with a Risk Aversion Coefficient of \(10^{-6}\).

fig, ax = plt.subplots(1, 1)
ax.plot(
    np.array(
        [
            b_j(
                a=np.array([[-1e6, 0]]),
                p_A=np.array([[0.9, 0.1]]),
                b_excl_j=np.array([[0]]),
                p_B_excl_j=np.array([[1 - p_B_j]]),
                p_B_j=np.array([[p_B_j]]),
                γ=0,
            )[0][0]
            for p_B_j in np.linspace(0.01, 0.99, num=100)
        ]
    ),
    np.linspace(0.01, 0.99, num=100),
    label="γ = 0",
)

ax.plot(
    np.array(
        [
            b_j(
                a=np.array([[-1e6, 0]]),
                p_A=np.array([[0.9, 0.1]]),
                b_excl_j=np.array([[0]]),
                p_B_excl_j=np.array([[1 - p_B_j]]),
                p_B_j=np.array([[p_B_j]]),
                γ=1e-06,
            )[0]
            for p_B_j in np.linspace(0.01, 0.99, num=100)
        ]
    ),
    np.linspace(0.01, 0.99, num=100),
    label="γ = 10^-6",
)

ax.set_xlabel(r"$X$")
ax.set_ylabel(r"$p(X)$")
ax.set_title(
    r"$p(X)$ Vs. $X$ Heat Map for Risk-Neutral (γ = 0) and Risk-Averse (γ = 10^-6) Decision Maker"
)
plt.legend()
plt.show()
time: 86.5 ms (started: 2021-03-26 18:58:00 +08:00)

4) Risk-Sensitivity Profile

Plot the Risk Sensitivity profile for the following two deals.

The Risk-Sensitivity profile is a plot of the certain equivalent of a deal vs, the risk aversion coefficient, \(\gamma\). Draw the plot from \(\gamma = −1\) to \(\gamma = +1\)

../../_images/lab_assignment_q4.png
# Risk-aversion coefficient
γ = np.linspace(-1, 1, num=100)


fig, ax = plt.subplots(1, 1)
ax.plot(
    γ,
    certainty_equivalent_values_calculator(
        u,
        u_inv,
        np.array([[80, 70, 0]]),
        np.array([[0.3, 0.2, 0.5]]),
        N,
        payout_lb=None,
        payout_ub=None,
        γ=γ,
    ),
    label=r"CE(Deal A)",
)

ax.plot(
    γ,
    certainty_equivalent_values_calculator(
        u,
        u_inv,
        np.array([[90, 50, 20]]),
        np.array([[0.1, 0.3, 0.6]]),
        N,
        payout_lb=None,
        payout_ub=None,
        γ=γ,
    ),
    label=r"CE(Deal B)",
)
ax.set_xlabel("γ")
ax.set_ylabel("CE")
ax.set_title("Risk Sensitivity Profile for Deal A and B")
plt.legend()
plt.show()
time: 43 ms (started: 2021-03-26 18:58:00 +08:00)

5) Large Monetary Prospects Using an Exponential U-Curve

In the deal below, assume \(p = 0.5\) and a decision maker with an exponential u-curve and a risk-aversion coefficient \(\gamma = 0.001\). I.e. the equation for the u-curve is

(33)\[\begin{align} u(x) &= a + b{e}^{-\gamma x} ,b \not= 0. \\ \end{align}\]

Plot the certain equivalent of the deal vs. \(X\) for the range where \(X=\$0\) to \(X=\$10,000\) in increments of \({\$500}\).

Explain and interpret the results.

../../_images/lab_assignment_q5.png
# Preferential probabilities for each prospect
x = np.array(
    list(zip(np.arange(0, 1e4 + 1, 500), np.zeros(len(np.arange(0, 1e4 + 1, 500)))))
)

# Payouts for each prospect
p = np.array([[0.5, 0.5] for _ in x])

# Risk-aversion coefficient
γ = np.array([0, 0.001, -0.001])

fig, ax = plt.subplots(1, 1)

ax.plot(
    x[:, 0],
    certainty_equivalent_values_calculator(
        u, u_inv, x, p, N=None, payout_lb=None, payout_ub=None, γ=γ
    )[0],
    label="γ = 0 (Risk-Neutral Individual)",
)

ax.plot(
    x[:, 0],
    certainty_equivalent_values_calculator(
        u, u_inv, x, p, N=None, payout_lb=None, payout_ub=None, γ=γ
    )[1],
    label="γ = 0.001 (Risk-Averse Individual)",
)

ax.plot(
    x[:, 0],
    certainty_equivalent_values_calculator(
        u, u_inv, x, p, N=None, payout_lb=None, payout_ub=None, γ=γ
    )[2],
    label="γ = -0.001 (Risk-Seeking Individual)",
)

ax.set_xlabel(r"$X$")
ax.set_ylabel("CE")
ax.set_title(
    r"CE(Decision Maker) $U(X) = a + be^{-\gamma x}$, $a = 0, b = 1$ Vs. $X$(\$)"
)
plt.legend()
plt.show()
time: 33.3 ms (started: 2021-03-26 18:58:00 +08:00)

We notice that the risk-averse individual (\(\gamma = 0.001 > 0\))’s certainty equivalent value for all gambles are lower than the risk-neutral individual (\(\gamma = 0\)), and much lower than the risk-seeking individual (\(\gamma = -0.001 < 0\)). For example, when presented with a deal to win \(\$4,000\) or nothing with equal probability, the risk-neutral individual is indifferent to playing the deal or simply getting the traditional expected value of \(\mathbb{E}\left[X\right] = \$2,000.00\). The risk-averse individual will be indifferent to that deal or getting \(\$675.00\) while the risk-seeking individual will be indifferent to that deal or getting \(\$3,325.00\).


6) EXTRA CREDIT: Sensitivity to Value of Clairvoyance

In the party problem, the following table represents the probabilities, prospects, and \(u\)-values for Kim, Mary, and Jane.

../../_images/lab_assignment_q6.png

a) Kim’s risk Odd for \(\$50\) are \(r_{50} = 2 :1\). What is her risk tolerance?

Assuming that Kim satisfies the \(\Delta\) property,

(34)\[\begin{align} r_{50} &= \frac{2}{1} \\ r &= {\left(\frac{2}{1}\right)}^{\frac{1}{50}} \\ \text{Risk Tolerance } \rho &= \frac{1}{\ln (r)} \\ &= \frac{50}{\ln (2)} \\ \end{align}\]

Also, we can find her risk tolerance by comparing the forms of the exponential \(U\)-curves:

(35)\[\begin{align} a + br^{-x} &= a + be^{-\frac{x}{\rho}} \\ \frac{4}{3} - \frac{4}{3} {\left(\frac{1}{2}\right)}^{\frac{x}{50}} &= \frac{4}{3} - \frac{4}{3}e^{-\frac{x}{\rho}} \\ {2}^{-\frac{x}{50}} &= e^{-\frac{x}{\rho}} \\ \frac{x}{50}\ln (2) &= \frac{x}{\rho} \\ \rho &= \frac{50}{\ln (2)} \\ \end{align}\]

b) Build a spreadsheet to calculate the Value of Clairvoyance for each person.

Assuming a person satisfies the \(\Delta\) property,

(36)\[\begin{align} \text{Value of Clairvoyance} &= \text{Value of Free Clairvoyance} - \text{Value of No Clairvoyance} \\ \end{align}\]
def value_of_clairvoyance(u, u_inv, x, p):
    """Calculates the Value of Clairvoyance
    = Value of Free Clairvoyance - Value of No Clairvoyance
    AKA How much more value we get when we have access
    to an oracle that can tell us the actual outcome of a
    random event.
    
    Args:
        u (function): U-curve function
        u_inv (function): Inverse U-curve function
        x (np.array): Payoffs of prospects matrix, Shape = (Number of deals, Number of prospects in each deal)
        p (np.array): Probabilities of each prospect matrix, Shape = (Number of deals, Number of prospects in each deal)
        
    Returns:
        float:
            Value of the Clairvoyance
    """

    # Check that payoffs and probability assignments are the same shape
    assert (
        x.shape == p.shape
    ), "Payoffs `x` and Probabilities `p` must be the same shape=(Number of deals, Number of prospects in each deal)."

    # Value of Free Clairvoyance
    vfc = certainty_equivalent_values_calculator(
        u=u,
        u_inv=u_inv,
        x=np.array([x.max(axis=0)]),
        p=np.array([p[np.argmax(x, axis=0), np.array(range(x.shape[1]))]]),
        γ=np.array([u.keywords["γ"]]),
    )[0][0]

    # Value of No Clairvoyance
    vnc = np.max(
        certainty_equivalent_values_calculator(
            u=u, u_inv=u_inv, x=x, p=p, γ=np.array([u.keywords["γ"]]),
        )
    )

    # Value of Clairvoyance
    return np.round(vfc - vnc, 2)
time: 2 ms (started: 2021-03-26 18:58:00 +08:00)

Kim’s Value of Clairvoyance

# Probabilities of Sun and Rain
p = np.array([[0.4, 0.6]])

# Prospect Valuation
x_kim = pd.DataFrame(
    np.array([[100, 90, 40], [0, 20, 50]]),
    index=["Sun", "Rain"],
    columns=["Outdoors", "Porch", "Indoors"],
)

value_of_clairvoyance(
    u=partial(u, γ=np.log(2) / 50, a=4 / 3, b=-4 / 3),
    u_inv=partial(u_inv, γ=np.log(2) / 50, a=4 / 3, b=-4 / 3),
    x=x_kim.T.values,
    p=np.repeat(p, repeats=x_kim.T.shape[0], axis=0),
)
-50.0
time: 6.39 ms (started: 2021-03-26 18:58:00 +08:00)

Jane’s Value of Clairvoyance

# Probabilities of Sun and Rain
p = np.array([[0.4, 0.6]])

# Prospect Valuation
x_jane = pd.DataFrame(
    np.array([[100, 90, 40], [0, 20, 50]]),
    index=["Sun", "Rain"],
    columns=["Outdoors", "Porch", "Indoors"],
)

value_of_clairvoyance(
    u=partial(u, γ=0, a=0, b=1),
    u_inv=partial(u_inv, γ=0, a=0, b=1),
    x=x_jane.T.values,
    p=np.repeat(p, repeats=x_jane.T.shape[0], axis=0),
)
22.0
time: 7.24 ms (started: 2021-03-26 18:58:00 +08:00)

Mary’s Value of Clairvoyance

# Probabilities of Sun and Rain
p = np.array([[0.4, 0.6]])

# Prospect Valuation
x_mary = pd.DataFrame(
    np.array([[100, 95, 57], [0, 32, 67]]),
    index=["Sun", "Rain"],
    columns=["Outdoors", "Porch", "Indoors"],
)

value_of_clairvoyance(
    u=partial(u, γ=0, a=0, b=1),
    u_inv=partial(u_inv, γ=0, a=0, b=1),
    x=x_mary.T.values,
    p=np.repeat(p, repeats=x_mary.T.shape[0], axis=0),
)
17.2
time: 8.09 ms (started: 2021-03-26 18:58:00 +08:00)

c) Conduct a sensitivity analysis of the Value of Clairvoyance to the probability of Sun for each person and plot the three-curves on one graph.

# Range of P(Sun) values to perform sensitivity analysis of VOC over
p_sun_values = np.linspace(0.001, 0.999, num=100)

# VOC values for each P(Sun) value
kim_voc_values, jane_voc_values, mary_voc_values = [
    [
        value_of_clairvoyance(
            u=u,
            u_inv=u_inv,
            x=x,
            p=np.repeat(np.array([[p_sun, 1 - p_sun]]), repeats=x.shape[0], axis=0),
        )
        for p_sun in p_sun_values
    ]
    for u, u_inv, x in [
        [
            partial(u, γ=np.log(2) / 50, a=4 / 3, b=-4 / 3),
            partial(u_inv, γ=np.log(2) / 50, a=4 / 3, b=-4 / 3),
            x_kim.T.values,
        ],
        [partial(u, γ=0, a=0, b=1), partial(u_inv, γ=0, a=0, b=1), x_jane.T.values],
        [partial(u, γ=0, a=0, b=1), partial(u_inv, γ=0, a=0, b=1), x_mary.T.values],
    ]
]

fig, ax = plt.subplots(1, 1)

ax.plot(
    p_sun_values, kim_voc_values, label="Kim",
)

ax.plot(
    p_sun_values, jane_voc_values, label="Jane",
)

ax.plot(
    p_sun_values, mary_voc_values, label="Mary",
)

ax.set_xlabel(r"$p(Sun)$")
ax.set_ylabel("Value of Clairvoyance")
ax.set_title(r"Value of Clairvoyance Vs. $p(Sun)$")
plt.legend()
plt.show()
time: 188 ms (started: 2021-03-26 18:58:00 +08:00)

d) Conduct a sensitivity analysis of VOC to both probability of Sun and the risk aversion coefficient for Kim.

from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm, colors
from matplotlib.ticker import LinearLocator

fig, ax = plt.subplots(subplot_kw={"projection": "3d"})

# Because EU values need to be positive for
# the inverse of an exponential to be defined,
# Kim's U-curve demands that γ values should
# be chosen such that e^{-xγ} <= 1. Hence, we wil
# only perfrom sensitivity analysis over γ > 0
γ = np.linspace(0.0001, 0.0099, num=100)
X, Y = np.meshgrid(p_sun_values, γ)
Z = np.array(
    [
        [
            value_of_clairvoyance(
                u=partial(u, γ=γ_i, a=4 / 3, b=-4 / 3),
                u_inv=partial(u_inv, γ=γ_i, a=4 / 3, b=-4 / 3),
                x=x_kim.T.values,
                p=np.repeat(
                    np.array([[p_sun, 1 - p_sun]]), repeats=x_kim.T.shape[0], axis=0
                ),
            )
            for p_sun in p_sun_values
        ]
        for γ_i in γ
    ]
)

# Plot the surface.
surf = ax.plot_surface(X, Y, Z, cmap=cm.coolwarm, linewidth=0, antialiased=False,)
ax.set_xlabel(r"$p(Sun)$")
ax.set_ylabel(r"Risk-aversion Coefficient $\gamma$")
ax.set_zlabel("VOC")
ax.set_title("Sensitivity Analysis of Value of Clairvoyance (VOC)")

# Add a color bar which maps values to colors.
fig.colorbar(surf, shrink=0.5, aspect=5)

plt.show()
time: 8.56 s (started: 2021-03-26 18:58:01 +08:00)