Solution sketch, Integration

MET 430

NTNU Business School

Johannes Mauritzen

Solution to EMEA ch.10

import numpy as np
import matplotlib.pyplot as plt

4. Computational Integration: consumer and Producer surplus

Say that the total inverse demand function for an economy for a certain product can be written:

$$P(Q) = \frac{150}{(ln(Q))}$$

And that the total inverse supply function for an economy can be written:

$$P(Q) = \frac{50}{(e^{1/Q})}$$

a.) Discretize and plot these two functions between Q=2 and Q=50, with n=1000. (We start at 2 here for convenience since $log(0)=\infty$ and $log(1)=0$)

b.) Find the approximate equilibrium price and quantity

c.) Write the consumer and producer surpluses as integrals. Calculate approximations of the consumer and producer surpluses (from Q=2), using both the trapezoid rule and simpsons rule.

a = 2
b = 50
n = 1000
h = (b - a) / (n - 1)
Q = np.linspace(a, b, n)
P_d = 150/np.log(Q)
P_s = 50/np.exp(1/Q)

fig, ax = plt.subplots()
ax.plot(Q, P_s, label="Supply")
ax.plot(Q, P_d, label="Demand")
ax.set_ylim(0,100)
ax.legend()

<matplotlib.legend.Legend at 0x7fa898edc970>

#to find the (approximate) equilibrium price and quantity, we can take the difference 
#between the two functions and find the absolute minimum

diff = P_d - P_s
P_star = P_s[diff==np.min(np.abs(diff))]
Q_star = Q[diff==np.min(np.abs(diff))]
print("Equilibrium Q", Q_star)
print("Equilibrium P", P_star)

Equilibrium Q [22.94894895]
Equilibrium P [47.86803841]

dCS = P_d - P_star
dPS = P_star - P_s

dCS = dCS[Q<=Q_star]
dPS = dPS[Q<=Q_star]
Q_s = Q[Q<=Q_star]

from scipy.integrate import trapz
from scipy.integrate import simpson

print("Consumer Surplus (trapezoid)", trapz(dCS,Q_s))
print("Producer Surplus (trapezoid)", trapz(dPS,Q_s))

Consumer Surplus (trapezoid) 470.55918074370743
Producer Surplus (trapezoid) 66.88781638170128

print("Consumer Surplus (simpson)", simpson(dCS,Q_s))
print("Producer Surplus (simpson)", simpson(dPS,Q_s))

Consumer Surplus (simpson) 470.5293103605592
Producer Surplus (simpson) 66.88637550623774

5.) Integrating a logistic diffusion model

T = 5000
ts = np.linspace(0,T, T)
L = 10000
k = .002
t0 = ts.mean()
R = L/(1 + np.exp(-k*(ts-t0)))

fig, ax = plt.subplots()
ax.plot(ts, R)

[<matplotlib.lines.Line2D at 0x7fa898f8d9a0>]

Suggested solution

We'll start by putting our logistic function in a function

def logisticFunc(ts, L, k, t0):
    return L/(1 + np.exp(-k*(ts-t0)))

T = 5000
ts = np.linspace(0,T, T)
L = 10000
k = .002
t0 = ts.mean()
carSales = logisticFunc(ts, L, k, t0)

Let's use the cumtrapz methodology for numerical integration since we want to get the cumulative integration.

from scipy.integrate import cumtrapz


carFleet = cumtrapz(logisticFunc(ts, L, k, t0), ts)
fig, ax = plt.subplots()
ax.plot(ts[0:-1], carFleet)

[<matplotlib.lines.Line2D at 0x7fa88b0abfa0>]

One possibility for improving this model is to also think about how a car exits the car park. That is, we should also consider that older electric cars will eventually get scrapped. We could introduce this into the model to get a better model of the total number of electric cars on the market.

6.) calculating the expected value of future revenue

# import normal pdf
from scipy.stats import norm
norm.pdf(3, loc=8, scale=3)

0.03315904626424956

def demandFunc(P):
    return 1000/np.log(P)

def revenueFunction(P):
    return demandFunc(P)*P

def priceProb(P):
    return norm.pdf(P, loc=10, scale=3)

Below we'll do a calculation that integrates over the function

$\int_{P=1.01}^{19} R(P)*prob(P)$

where $R(P)$ is the revenue function and $prob(P)$ is the probability distribution of P. The limits are chosen to be approximate +/- 3 standard deviations of the normal distribution, which represents about 99.8% of the probability.

Notice the use of the lambda function below. This is just a quick way of writing a function

from scipy.integrate import quad 

expected_revenue, est_err = quad(lambda P: revenueFunction(P)*priceProb(P), 1.01, 19)

expected_revenue

4332.740917852739

1000/np.log(10)*10

4342.944819032517