Cost and Benefit Analysis - Basic simulation with SymPy - Econometrics with Python Programming
" ~~~~~~ Initial setup ~~~~~~ "
# importing liberaries
import math
import numpy as np
from scipy.optimize import newton
import matplotlib.pyplot as plt
# SymPy for polynomial intepretation
import sympy as smp
# Showing the plots bigger
#plt.rcParams["figure.figsize"] = (18,9) # Bigger
plt.rcParams["figure.figsize"] = (9,5)
" ~~~~~~ Profit and Loss (Revenue and Cost) ~~~~~~ "
# Qauntity demanded
Q=np.linspace(1,1000,1000) # Numerical scale
q =smp.symbols('q',real=True, positive=True) # Defining an algebraic symbol
# Demand = Average Revenue in Polynomial Symbolic
ar = 2000*(q**0) -2*(q**1)
display("Demand = Average Revenue =",ar)
AR=[ar.subs([(q,Q[k])]) for k in range(len(Q))]
# Supply ≈ Average Cost
ac = 0.01*(q-300)**2
AC=[ac.subs([(q,Q[k])]) for k in range(len(Q))]
C=200; #min(AC)
ac = ac+C
ac = smp.expand(ac)
display("Supply ≈ Average Cost =",ac)
AC=[ac.subs([(q,Q[k])]) for k in range(len(Q))]
# Plot both
plt.figure;
plt.plot(AR, label="Demand = Average Revenue (AR)")
plt.plot(AC, label="Supply ≈ Average Cost (AC)")
plt.legend(loc="best")
plt.xlabel('Passengers')
plt.ylabel('Price')
plt.show(block=False)
q_nonProfit=smp.solve(ar-ac)[-1]
p_nonProfit=ar.subs(q,q_nonProfit)
display(f"Break-even (neither profit nor loss) when the price is {round(p_nonProfit)} and the quantity is {round(q_nonProfit)}.")
# Total Revenue
tr = smp.expand(ar * q)
display("Total Revenue = AR * q =",tr)
TR=[tr.subs([(q,Q[k])]) for k in range(len(Q))]
# Total Cost
tc = smp.expand(ac * q)
display("Total Cost = AC * q =",tc)
TC=[tc.subs([(q,Q[k])]) for k in range(len(Q))]
# Plot both
plt.figure;
plt.plot(TR, label="Total Revenue")
plt.plot(TC, label="Total Cost")
plt.legend(loc="best")
plt.xlabel('Passengers')
plt.ylabel('Price')
plt.show(block=False)
# Marginal Revenue
mr = smp.diff(tr,q)
display("Marginal Revenue = dTR/dq =",mr)
MR=[mr.subs([(q,Q[k])]) for k in range(len(Q))]
# Marginal Cost
mc = smp.diff(tc,q)
display("Marginal Cost = dTC/dq =",mc)
MC=[mc.subs([(q,Q[k])]) for k in range(len(Q))]
def DispProfit(y, x):
display(f"The price maximising the profit is {round(y)} with the quantity {round(x)} .")
# Finding the quantity maximising the profit
# where the marginal revenue is equal to the marginal cost
mr_mc=mr-mc
display("MR -MC =",mr_mc)
MR_MC_0=smp.solve(mr-mc)
q_opt=MR_MC_0[-1]
p_opt=ar.subs(q,q_opt)
print("Find with MR=MC")
DispProfit(p_opt, q_opt)
GraphCut=600
plt.figure;
plt.plot(Q[:GraphCut],AR[:GraphCut], label="Average Revenue")
plt.plot(Q[:GraphCut],AC[:GraphCut], label="Average Cost")
plt.plot(Q[:GraphCut],MR[:GraphCut], label="Marginal Revenue")
plt.plot(Q[:GraphCut],MC[:GraphCut], label="Marginal Cost")
plt.axvline(x=q_opt, color='Grey')
plt.legend(loc="best")
plt.xlabel('Passenger')
plt.ylabel('Price')
plt.show(block=False)
# Finding the quantity maximising the profit
# with the function max() and index()
profit=tr-tc
display("Profit Function =",profit)
Profit=[profit.subs([(q,Q[k])]) for k in range(len(Q))]
q_opt_=round(Profit.index(max(Profit))+min(Q))
print("Find with the function max() and index() ")
DispProfit(p_opt, q_opt)
" ~~~~~~ Capital and Investment ~~~~~~ "
# Cash flow
apl=ar-ac
P=apl.subs([(q,q_opt)])*q_opt
print("Profit per day is",P)
# Yearly nominal interest rate to lend
r=0.005324 # Interest rate #0.000324
d=0.005 # Discount rate
Delta=365
# Initial investment: Bond issued
I=1e9 # Investment
print("The initial investment (Debt) is", I)
# Calculating the repayment plan
def B_j(r,d,Delta,I):
# Borrowing after the initial construction
Y_I=2 # Years spent for the construction
B_0=I*math.exp(Y_I*(r-d));
B_j=[B_0]
for j in range(10*Delta):
B_j.append((B_j[j]-P)*math.exp((j*(r-d))/Delta)-P)
if B_j[-1] <= 0:
#print(B_j[-1])
break
else:
continue
if B_j[-1]>0:
print("Unable to repay!!")
return j,B_j,B_j[-1]
B_j=B_j(r,d,Delta,I)
t=round(B_j[0]/Delta,2)
print(f'It needs {math.floor(t)} years and {math.ceil(t%12)} months to repay')
plt.plot(B_j[1])
