我试图在具有给定需求函数的市场中寻找均衡，一家公司最大化其利润，政府设定消费税以最大化税收收入。我为公司和政府用来最大化利润/税收的最大化问题编写了代码。现在我不确定下一步要找到平衡是什么。谁能帮我下一步？
我想我必须定义一个不动点，并用它来求解平衡。在我大学上课时，我们这样做只是为了与多家公司的库尔诺均衡。我不确定在这种情况下如何通过公司和政府设置消费税来做到这一点。稍后我想向模型添加第二个产品(替代或补充)，但我需要弄清楚如何首先基于最佳响应在这个更简单的模型中求解均衡价格和消费税。
这是我编写的代码，我正在寻找下一步，根据他们的最佳反应解决均衡价格和消费税。
from scipy import stats, optimize
alpha1 = 10
beta1 = -0.5
c1 = 0.5
def Q1(P1,t1):
return alpha1+beta1*P1*(1+t1)
def profits1(P1,t1,c1):
return Q1(P1,t1)*(P1-c1)
def taxes(t1,P1):
return Q1(P1,t1)*P1*t1
def reactionfirm1(t1,c1):
P1 = optimize.fminbound(lambda x: -profits1(x,t1,c1),0,100)
return P1
def reactiongovernment1(P1):
t1 = optimize.fminbound(lambda x: -taxes(x,P1),0,100)
return t1
用笔和纸解决这个问题得到t1=3.73和P1=2.36(都是四舍五入的)。
要用纸笔解决，我首先要找到企业和政府的最佳回应。对于公司来说，最好的反应是通过取利润函数相对于 P1 的一阶导数，将其等于 0 并将 P1 写为 t1 的函数来找到的。这给了：
P1 = (10.25+0.25t1)/(1+t)
同样，为了获得政府的最佳响应，可以通过取税收收入相对于 t1 的一阶导数，将其等于 0 并将 t1 写为 P1 的函数，得到：
t1 = 10/P1 -0.5
要找到 P1 和 t1，您可以将 t1 插入 P1或将 P1 插入 t1 并求解它。
所以我将 P1 插入 t1 以获得：
t1 = (10+10吨)/(10.25+0.25吨) -0.5
这给出 t1 = 3.73
通过在公司的最佳响应中插入 t1 = 3.73，我得到 P1 = 2.36
我希望用Python做的是计算和解决公司和政府找到t1和P1的最佳响应。因为我想扩展模型，这使得用笔和纸解决起来更加复杂。此外，通过在python中使用它，可以更轻松地更改参数alpha1，beta1和c1，以查看最佳P1和t1如何变化。
我在论文中使用了这个python代码来寻找库尔诺均衡(最大化利润的数量)。我还把我用来检查凹陷的代码留给你。我为 8 家公司这样做，但您可以将其折叠为两家。此外，成本更复杂，您可以取出字母并只保留可变成本，因为这些字母代表成本的线性降低：
import sympy as sym
import numpy as np
import matplotlib.pyplot as plt
# Define symbols
q1, q2, q3, q4, q5, q6, q7, q8, a1, b2, c3, d4, e1, e5, e6, e7, e8, f2, f5, f6, f7, f8, g3, g5, g6, g7, g8, h4, h5, h6, h7, h8, n = sym.symbols('q1 q2 q3 q4 q5 q6 q7 q8 a1 b2 c3 d4 e1 e5 e6 e7 e8 f2 f5 f6 f7 f8 g3 g5 g6 g7 g8 h4 h5 h6 h7 h8 n')
# Define parameters
#alpha = sym.Symbol('alpha')
#cost = sym.Symbol('cost')
#s = sym.Symbol('s')
#each time you change network also change prices, costs, utility
alpha=10
cost=5
s=0.8
# Define profit functions
P1 = alpha - q1 - s*q2 - s*q8
P2 = alpha - q2 - s*q3 - s*q1
P3 = alpha - q3 - s*q4 - s*q2
P4 = alpha - q4 - s*q5 - s*q3
P5 = alpha - q5 - s*q6 - s*q4
P6 = alpha - q6 - s*q7 - s*q5
P7 = alpha - q7 - s*q8 - s*q6
P8 = alpha - q8 - s*q1 - s*q7
C1 = cost - a1 - n*e1
C2 = cost - b2 - n*f2
C3 = cost - c3 - n*g3
C4 = cost - d4 - n*h4
C5 = cost - e5 - e1 - n*a1 - e6 - e7 - e8 - n*f5 - n*f6 - n*f7 - n*f8 - n*g5 - n*g6 - n*g7 - n*g8 - n*h5 - n*h6 - n*h7 - n*h8
C6 = cost - f6 - f2 - n*b2 - f5 - f7 - f8 - n*e5 - n*e6 - n*e7 - n*e8 - n*g5 - n*g6 - n*g7 - n*g8 - n*h5 - n*h6 - n*h7 - n*h8
C7 = cost - g7 - g3 - n*c3 - g5 - g6 - g8 - n*f5 - n*f6 - n*f7 - n*f8 - n*e5 - n*e6 - n*e7 - n*e8 - n*h5 - n*h6 - n*h7 - n*h8
C8 = cost - h8 - h4 - n*d4 - h5 - h6 - h7 - n*f5 - n*f6 - n*f7 - n*f8 - n*g5 - n*g6 - n*g7 - n*g8 - n*e5 - n*e6 - n*h7 - n*e8
prof1 = (P1-C1)*q1 - 2*a1**2
prof2 = (P2-C2)*q2 - 2*b2**2
prof3 = (P3-C3)*q3 - 2*c3**2
prof4 = (P4-C4)*q4 - 2*d4**2
prof5 = (P5-C5)*q5 - 10*e1**2 - 10*e5**2 - 10*e6**2 - 10*e7**2 - 10*e8**2
prof6 = (P6-C6)*q6 - 10*f2**2 - 10*f5**2 - 10*f6**2 - 10*f7**2 - 10*f8**2
prof7 = (P7-C7)*q7 - 10*g3**2 - 10*g5**2 - 10*g6**2 - 10*g7**2 - 10*g8**2
prof8 = (P8-C8)*q8 - 10*h4**2 - 10*h5**2 - 10*h6**2 - 10*h7**2 - 10*h8**2
# Take derivatives
dq1 = sym.diff(prof1, q1)
dq2 = sym.diff(prof2, q2)
dq3 = sym.diff(prof3, q3)
dq4 = sym.diff(prof4, q4)
dq5 = sym.diff(prof5, q5)
dq6 = sym.diff(prof6, q6)
dq7 = sym.diff(prof7, q7)
dq8 = sym.diff(prof8, q8)
# Set derivatives to 0 and solve for q1, q2, q3
soln = sym.solve((dq1, dq2, dq3, dq4, dq5, dq6, dq7, dq8), (q1, q2, q3, q4, q5, q6, q7, q8))
# Print optimal values of q1, q2, q3
print("Optimal value of q1:", soln[q1])
print("Optimal value of q2:", soln[q2])
print("Optimal value of q3:", soln[q3])
print("Optimal value of q4:", soln[q4])
print("Optimal value of q5:", soln[q5])
print("Optimal value of q6:", soln[q6])
print("Optimal value of q7:", soln[q7])
print("Optimal value of q8:", soln[q8])
p1 = P1.subs({q1: soln[q1], q2: soln[q2], q3: soln[q3], q4: soln[q4], q5: soln[q5], q6: soln[q6], q7: soln[q7], q8: soln[q8]})
p2 = P2.subs({q1: soln[q1], q2: soln[q2], q3: soln[q3], q4: soln[q4], q5: soln[q5], q6: soln[q6], q7: soln[q7], q8: soln[q8]})
p3 = P3.subs({q1: soln[q1], q2: soln[q2], q3: soln[q3], q4: soln[q4], q5: soln[q5], q6: soln[q6], q7: soln[q7], q8: soln[q8]})
p4 = P4.subs({q1: soln[q1], q2: soln[q2], q3: soln[q3], q4: soln[q4], q5: soln[q5], q6: soln[q6], q7: soln[q7], q8: soln[q8]})
p5 = P5.subs({q1: soln[q1], q2: soln[q2], q3: soln[q3], q4: soln[q4], q5: soln[q5], q6: soln[q6], q7: soln[q7], q8: soln[q8]})
p6 = P6.subs({q1: soln[q1], q2: soln[q2], q3: soln[q3], q4: soln[q4], q5: soln[q5], q6: soln[q6], q7: soln[q7], q8: soln[q8]})
p7 = P7.subs({q1: soln[q1], q2: soln[q2], q3: soln[q3], q4: soln[q4], q5: soln[q5], q6: soln[q6], q7: soln[q7], q8: soln[q8]})
p8 = P8.subs({q1: soln[q1], q2: soln[q2], q3: soln[q3], q4: soln[q4], q5: soln[q5], q6: soln[q6], q7: soln[q7], q8: soln[q8]})
co1 = C1.subs({q1: soln[q1], q2: soln[q2], q3: soln[q3], q4: soln[q4], q5: soln[q5], q6: soln[q6], q7: soln[q7], q8: soln[q8]})
co2 = C2.subs({q1: soln[q1], q2: soln[q2], q3: soln[q3], q4: soln[q4], q5: soln[q5], q6: soln[q6], q7: soln[q7], q8: soln[q8]})
co3 = C3.subs({q1: soln[q1], q2: soln[q2], q3: soln[q3], q4: soln[q4], q5: soln[q5], q6: soln[q6], q7: soln[q7], q8: soln[q8]})
co4 = C4.subs({q1: soln[q1], q2: soln[q2], q3: soln[q3], q4: soln[q4], q5: soln[q5], q6: soln[q6], q7: soln[q7], q8: soln[q8]})
co5 = C5.subs({q1: soln[q1], q2: soln[q2], q3: soln[q3], q4: soln[q4], q5: soln[q5], q6: soln[q6], q7: soln[q7], q8: soln[q8]})
co6 = C6.subs({q1: soln[q1], q2: soln[q2], q3: soln[q3], q4: soln[q4], q5: soln[q5], q6: soln[q6], q7: soln[q7], q8: soln[q8]})
co7 = C7.subs({q1: soln[q1], q2: soln[q2], q3: soln[q3], q4: soln[q4], q5: soln[q5], q6: soln[q6], q7: soln[q7], q8: soln[q8]})
co8 = C8.subs({q1: soln[q1], q2: soln[q2], q3: soln[q3], q4: soln[q4], q5: soln[q5], q6: soln[q6], q7: soln[q7], q8: soln[q8]})
Prof1 = (p1 - co1)*soln[q1] - 2*a1**2
Prof2 = (p2 - co2)*soln[q2] - 2*b2**2
Prof3 = (p3 - co3)*soln[q3] - 2*c3**2
Prof4 = (p4 - co4)*soln[q4] - 2*d4**2
Prof5 = (p5 - co5)*soln[q5] - 10*e1**2 - 10*e5**2 - 10*e6**2 - 10*e7**2 - 10*e8**2
Prof6 = (p6 - co6)*soln[q6] - 10*f2**2 - 10*f5**2 - 10*f6**2 - 10*f7**2 - 10*f8**2
Prof7 = (p7 - co7)*soln[q7] - 10*g3**2 - 10*g5**2 - 10*g6**2 - 10*g7**2 - 10*g8**2
Prof8 = (p8 - co8)*soln[q8] - 10*h4**2 - 10*h5**2 - 10*h6**2 - 10*h7**2 - 10*h8**2
#########################################################
print("Profits with optimal quantities for firm 1:", Prof1)
print("Profits with optimal quantities for firm 2:", Prof2)
print("Profits with optimal quantities for firm 3:", Prof3)
print("Profits with optimal quantities for firm 4:", Prof4)
print("Profits with optimal quantities for firm 5:", Prof5)
print("Profits with optimal quantities for firm 6:", Prof6)
print("Profits with optimal quantities for firm 7:", Prof7)
print("Profits with optimal quantities for firm 8:", Prof8)
# Compute the second derivative of the function with respect to a1, a2
f_double_prime_Prof1_a1 = Prof1.diff(a1, 2)
f_double_prime_Prof2_b2 = Prof2.diff(b2, 2)
f_double_prime_Prof3_c3 = Prof3.diff(c3, 2)
f_double_prime_Prof4_d4 = Prof4.diff(d4, 2)
f_double_prime_Prof5_e1 = Prof5.diff(e1, 2)
f_double_prime_Prof5_e5 = Prof5.diff(e5, 2)
f_double_prime_Prof5_e6 = Prof5.diff(e6, 2)
f_double_prime_Prof5_e7 = Prof5.diff(e7, 2)
f_double_prime_Prof5_e8 = Prof5.diff(e8, 2)
f_double_prime_Prof6_f2 = Prof5.diff(f2, 2)
f_double_prime_Prof6_f5 = Prof6.diff(f5, 2)
f_double_prime_Prof6_f6 = Prof6.diff(f6, 2)
f_double_prime_Prof6_f7 = Prof6.diff(f7, 2)
f_double_prime_Prof6_f8 = Prof6.diff(f8, 2)
f_double_prime_Prof7_g3 = Prof7.diff(g3, 2)
f_double_prime_Prof7_g5 = Prof7.diff(g5, 2)
f_double_prime_Prof7_g6 = Prof7.diff(g6, 2)
f_double_prime_Prof7_g7 = Prof7.diff(g7, 2)
f_double_prime_Prof7_g8 = Prof7.diff(g8, 2)
f_double_prime_Prof8_h4 = Prof8.diff(h4, 2)
f_double_prime_Prof8_h5 = Prof8.diff(h5, 2)
f_double_prime_Prof8_h6 = Prof8.diff(h6, 2)
f_double_prime_Prof8_h7 = Prof8.diff(h7, 2)
f_double_prime_Prof8_h8 = Prof8.diff(h8, 2)
# Values of n to check for concavity
n_values = np.arange(0, 1.01, 0.01)
# Check concavity for each value of n
for n_val in n_values:
# Substitute n with the current value .subs(n, n_val)
f_double_prime_Prof1_a1_n = f_double_prime_Prof1_a1.subs({n: n_val})
f_double_prime_Prof2_b2_n = f_double_prime_Prof2_b2.subs({n: n_val})
f_double_prime_Prof3_c3_n = f_double_prime_Prof3_c3.subs({n: n_val})
f_double_prime_Prof4_d4_n = f_double_prime_Prof4_d4.subs({n: n_val})
f_double_prime_Prof5_e1_n = f_double_prime_Prof5_e1.subs({n: n_val})
f_double_prime_Prof5_e5_n = f_double_prime_Prof5_e5.subs({n: n_val})
f_double_prime_Prof5_e6_n = f_double_prime_Prof5_e6.subs({n: n_val})
f_double_prime_Prof5_e7_n = f_double_prime_Prof5_e7.subs({n: n_val})
f_double_prime_Prof5_e8_n = f_double_prime_Prof5_e8.subs({n: n_val})
f_double_prime_Prof6_f2_n = f_double_prime_Prof6_f2.subs({n: n_val})
f_double_prime_Prof6_f5_n = f_double_prime_Prof6_f5.subs({n: n_val})
f_double_prime_Prof6_f6_n = f_double_prime_Prof6_f6.subs({n: n_val})
f_double_prime_Prof6_f7_n = f_double_prime_Prof6_f7.subs({n: n_val})
f_double_prime_Prof6_f8_n = f_double_prime_Prof6_f8.subs({n: n_val})
f_double_prime_Prof7_g3_n = f_double_prime_Prof7_g3.subs({n: n_val})
f_double_prime_Prof7_g5_n = f_double_prime_Prof7_g5.subs({n: n_val})
f_double_prime_Prof7_g6_n = f_double_prime_Prof7_g6.subs({n: n_val})
f_double_prime_Prof7_g7_n = f_double_prime_Prof7_g7.subs({n: n_val})
f_double_prime_Prof7_g8_n = f_double_prime_Prof7_g8.subs({n: n_val})
f_double_prime_Prof8_h4_n = f_double_prime_Prof8_h4.subs({n: n_val})
f_double_prime_Prof8_h5_n = f_double_prime_Prof8_h5.subs({n: n_val})
f_double_prime_Prof8_h6_n = f_double_prime_Prof8_h6.subs({n: n_val})
f_double_prime_Prof8_h7_n = f_double_prime_Prof8_h7.subs({n: n_val})
f_double_prime_Prof8_h8_n = f_double_prime_Prof8_h8.subs({n: n_val})
# Check if the second derivatives are non-positive
if (f_double_prime_Prof1_a1_n.is_nonpositive and f_double_prime_Prof2_b2_n.is_nonpositive
and f_double_prime_Prof3_c3_n.is_nonpositive and f_double_prime_Prof4_d4_n.is_nonpositive
and f_double_prime_Prof5_e1_n.is_nonpositive and f_double_prime_Prof5_e5_n.is_nonpositive
and f_double_prime_Prof5_e6_n.is_nonpositive and f_double_prime_Prof5_e7_n.is_nonpositive
and f_double_prime_Prof5_e8_n.is_nonpositive and f_double_prime_Prof6_f5_n.is_nonpositive
and f_double_prime_Prof6_f2_n.is_nonpositive and f_double_prime_Prof7_g3_n.is_nonpositive
and f_double_prime_Prof6_f6_n.is_nonpositive and f_double_prime_Prof6_f7_n.is_nonpositive
and f_double_prime_Prof6_f8_n.is_nonpositive and f_double_prime_Prof7_g5_n.is_nonpositive
and f_double_prime_Prof7_g6_n.is_nonpositive and f_double_prime_Prof7_g7_n.is_nonpositive
and f_double_prime_Prof7_g8_n.is_nonpositive and f_double_prime_Prof8_h4_n.is_nonpositive
and f_double_prime_Prof8_h5_n.is_nonpositive and f_double_prime_Prof8_h6_n.is_nonpositive
and f_double_prime_Prof8_h7_n.is_nonpositive and f_double_prime_Prof8_h8_n.is_nonpositive):
print(f"All functions are concave for n = {n_val}")
else:
print(f"Some functions are not concave for n = {n_val}")