Cournot optimization problem.
NoutHakkesteegt opened this issue · 3 comments
Hi,
My group member Viktorija emailed you last week regarding the model we are trying to recreate. We have since run into a problem while recreating the model from the paper 'International R & D Rivalry and Industrial Strategy'
In the paper they first simply solve the cournot competition to find optimal q_1 and q_2. Then find the level of R&D the paper create a 'gains' function that uses the profit function from the first step to maximize 'gains' with respect to R &D, see function 10 from the paper. we are trying to recreate this function using reaction_3 and reaction_4 in the code below. currently they return a a value but we don't know of anyway to verify these values and were wondering if you had any suggestions?
###last stage: choice of output
##variables
q_1 = 1
q_2 = 1
#quantities:
def Q(q_1, q_2):
return q_1 + q_2
#demand (price):
max_wtp= 10
b = 1 #elasticity of substitution between 1 and 2, paper assumes substitutes (1 for perfect substitutes)
def P(q_1, q_2, max_wtp):
return max_wtp - q_1 - b*q_2
#costs:
#fixed_c_1= not sure if we need these
#fixed_c_2=
mc_1 = 0.5
mc_2 = 0.5
#R&D levels & costs:
q_rd_1 = 0
q_rd_2 = 0
c_rd_1 = 0.5
c_rd_2 = 0.5
#revenue for firm 1 and 2:
def revenue_1(max_wtp, q_1, q_2):
return P(q_1, q_2, max_wtp)*q_1
def revenue_2(max_wtp, q_1, q_2):
return P(q_1, q_2, max_wtp)*q_2
#variable costs for firm 1 and 2: (NOTE: R&D level q_rd_1 and q_rd_2 should be added here to show the effect of cost-reducting R&D and they should depend on quantities)
def variable_c_1(q_rd_1, q_1, x):
if q_rd_1 == 0:
variable_c_1 = 0.5(q_1)*2
else:
variable_c_1 = 0.5(q_1/q_rd_1)**2
return variable_c_1
def variable_c_2(q_rd_2, q_2, x):
if q_2 == 0:
variable_c_2 = 0.5*(q_2)*2
else:
variable_c_2 = 0.5(q_2/q_rd_2)**2
return variable_c_2
#profits for firm 1 and 2: (q_rd_1 and q_rd_2 are squared to denote diminishing returns to scale and so that the code doesnt run forever)
def profit_1(q_1, mc_1, q_rd_1, c_rd_1):
return revenue_1(max_wtp, q_1, q_2) - variable_c_1(mc_1, q_1, x) - (q_rd_1**2)*c_rd_1
def profit_2(q_2, mc_2, q_rd_2, c_rd_2):
return revenue_2(max_wtp, q_1, q_2) - variable_c_2(mc_2, q_2, x) - (q_rd_2**2)*c_rd_2
#reaction functions for firm 1 and 2: (I assume that these are the FOC's and R&D will be taken into account)
def reaction_1(mc_1, q_rd_1, c_rd_1):
q_1 = optimize.brute(lambda q_1: -profit_1(q_1, mc_1, q_rd_1, c_rd_1), ((0,1,),)) # brute minimizes -profit_1 the function
return q_1[0]
def reaction_2(mc_2, q_rd_2, c_rd_2):
q_2 = optimize.brute(lambda q_2: -profit_2(q_2, mc_2, q_rd_2, c_rd_2), ((0,1,),)) # brute minimizes -profit_2 the function
return q_2[0]
#Nash Equilibrium (not sure if its correct, should depend on R&D variables)
parameters_1 = (mc_1, q_rd_1, c_rd_1)
parameters_2 = (mc_2, q_rd_2, c_rd_2)
parameters_total = (parameters_1, parameters_2)
def vector_reaction(q, parameters_1, parameters_2):
return np.array(q)-np.array([reaction_1(q[1], q_rd_1, c_rd_1),reaction_2(q[0], q_rd_2, c_rd_2)])
initial_guess = [0.3, 0.3]
cournot_equilibrium = optimize.fsolve(vector_reaction, initial_guess, args = (parameters_total))
#now for rd
def reaction_3(q_1, mc_1, q_rd_1, c_rd_1):
q_rd_1 = optimize.brute(lambda q_1: -profit_1(q_1, mc_1, q_rd_1, c_rd_1), ((0,1,),)) # brute minimizes -profit_1 the function
return q_rd_1[0]
def reaction_4(q_2, mc_2, q_rd_2, c_rd_2):
q_rd_2 = optimize.brute(lambda q_2: -profit_2(q_2, mc_2, q_rd_2, c_rd_2), ((0,1,),)) # brute minimizes -profit_1 the function
return q_rd_2[0]
q_rd_1_test = reaction_3(cournot_equilibrium[0], mc_1, q_rd_1, c_rd_1)
q_rd_2_test = reaction_4(cournot_equilibrium[1], mc_2, q_rd_2, c_rd_2)
print(q_rd_1_test)
print(q_rd_2_test)
print(cournot_equilibrium)
print(profit_1(cournot_equilibrium[0], mc_1, q_rd_1, c_rd_1))
print(profit_2(cournot_equilibrium[1], mc_2, q_rd_2, c_rd_2))
Some comments on your code to simplify things:
I am not sure why in the last stage you set =q_1, q_2= equal to 1. Is this as initial values in an optimization routine?
you define a function =Q= that you do not seem to use?
the function =P= that you define is correct only for firm 1 (in case b < 1)
If you want to consider b < 1: the definition of =revenue_2= is actually not correct; should be =P(q_2,q_1,max_wtp)=
Indeed, I would not use fixed costs; the R&D costs already generate fixed costs for the firm
as above: why are you setting =q_rd_1, q_rd_2=?
your definition of =variable_c_1= is counter-intuitive. If a firm invests 0 in R&D, its costs are lower than when it invests, say, 0.5 in R&D
I understand that you do not want to divide by 0; but you can fix this by working with =(q_1/(1+q_rd_1))**2=
it is, actually, easier to work with 0.5/(1+q_rd_1)*q_1**2
Basically, the final (Cournot) stage should give you profits for both firms as a function of =q_rd_1, q_rd_2= (solved numerically as the intersection of the reaction functions and then the resulting outputs are substituted into the profit functions)
With these functions you can specify the first stage problem where, say, firm 1 maximizes:
profits(q_rd_1,q_rd_2) - cost(q_rd_1)
The structure of this problem is similar to a standard 2nd stage Cournot model, except that the profits function is a bit more complicated (as it is the solution to the Cournot output problem).
Is this helpful?
Jan.
Thank you for your suggestions, It cleaned up the code quite a bit. Unfortunately we are still running into a bit of a problem when it comes to optimizing profits (q_rd_1,q_rd_2) - cost(q_rd_1). as you can see below we are using reaction3 and 4 for this, but our results for q_rd_1 and q_rd_2 are only for the q_1 and q_2 found in the cournot equilibrium where there is no R&D yet. Because we need both quantity and R&D to be optimized at the same time but in the current code only R&D is. we feel like we are either missing a step or overseeing a very obvious mistake. One of the solutions we came up with is to run a while loop to test profit for increasing levels of R&D until profits no longer increase, but that's very slow and gimmicky with the exit break conditions. Hopefully you can help us hank you in advance.
##variables, initial values
q_1 = 1
q_2 = 1
#demand (price):
max_wtp= 10
b = 1 #elasticity of substitution between 1 and 2, paper assumes substitutes (1 for perfect substitutes)
def P1(q_1, q_2, max_wtp):
return max_wtp - q_1 - bq_2
def P2(q_2, q_1, max_wtp):
return max_wtp - q_2 - bq_1
#R&D levels & costs:
q_rd_1 = 0
q_rd_2 = 0
c_rd_1 = 0.5
c_rd_2 = 0.5
#revenue for firm 1 and 2:
def revenue_1(max_wtp, q_1, q_2):
return P1(q_1, q_2, max_wtp)*q_1
def revenue_2(max_wtp, q_2, q_1):
return P2(q_2, q_1, max_wtp)*q_2
#variable costs for firm 1 and 2:
def variable_c_1(q_rd_1, q_1):
return ((q_1/(1+q_rd_1))**2)
def variable_c_2(q_rd_2, q_2):
return ((q_2/(1+q_rd_2))**2)
#profits for firm 1 and 2:
def profit_1(q_1, q_2, q_rd_1, c_rd_1):
return revenue_1(max_wtp, q_1, q_2) - variable_c_1(q_rd_1, q_1) - (q_rd_1)*c_rd_1
def profit_2(q_2, q_1, q_rd_2, c_rd_2):
return revenue_2(max_wtp, q_2, q_1) - variable_c_2(q_rd_2, q_2) - (q_rd_2)*c_rd_2
#reaction functions for firm 1 and 2: (I assume that these are the FOC's and R&D will be taken into account)
def reaction_1(q_1, q_2, q_rd_1, c_rd_1):
q_1 = optimize.brute(lambda q_1: -profit_1(q_1, q_2, q_rd_1, c_rd_1), ((0,1,),)) # brute minimizes -profit_1 the function
return q_1[0]
def reaction_2(q_2, q_1, q_rd_2, c_rd_2):
q_2 = optimize.brute(lambda q_2: -profit_2(q_2, q_1, q_rd_2, c_rd_2), ((0,1,),)) # brute minimizes -profit_2 the function
return q_2[0]
#Nash Equilibrium (not sure if its correct, should depend on R&D variables)
parameters_1 = [q_1, q_2, q_rd_1, c_rd_1]
parameters_2 = [q_2, q_1, q_rd_2, c_rd_2]
#parameters_total = [parameters_1, parameters_2]
def vector_reaction(q, parameters_1, parameters_2):
return np.array(q)-np.array([reaction_1(q[1], parameters_1[1], parameters_1[2], parameters_1[3]),reaction_2(q[0], parameters_2[1], parameters_2[2], parameters_2[3])])
initial_guess = [0.1, 0.1]
cournot_equilibrium = optimize.fsolve(vector_reaction, initial_guess, args = (parameters_1, parameters_2))
print(cournot_equilibrium)
#now for rd
def reaction_3(q_1, q_2, q_rd_1, c_rd_1):
q_rd_1 = optimize.brute(lambda q_rd_1: -profit_1(q_1, q_2, q_rd_1, c_rd_1), ((0,1,),)) # brute minimizes -profit_1 the function
return q_rd_1[0]
def reaction_4(q_2, q_1, q_rd_2, c_rd_2):
q_rd_2 = optimize.brute(lambda q_rd_2: -profit_2(q_2, q_1, q_rd_2, c_rd_2), ((0,1,),)) # brute minimizes -profit_1 the function
return q_rd_2[0]
def vector_reaction_rd(q_rd, parameters_1, parameters_2):
return np.array(q_rd)-np.array([reaction_3(parameters_1[0], parameters_1[1], q_rd[1], parameters_1[3]),reaction_4(parameters_2[0], parameters_2[1], q_rd[0], parameters_2[3])])
rd_equilibrium = optimize.fsolve(vector_reaction_rd, initial_guess, args = (parameters_3, parameters_4))
print(rd_equilibrium)
You are moving in the right direction; well done!
You have defined the relevant functions to calculate the Cournot equilibrium in the last stage as a function of c_1, c_2 (I simplify notation here to make this more readable).
Hence by gathering these functions together, can you define a function profits(c_1,c_2)? This function calculates the Cournot equilibrium profits for firm 1 if it has costs c_1 while its competitor has costs c_2. The equilibrium profits for 2 are then given by profits(c_2,c_1)
Once you have these functions, the stage 1 profit function for firm 1 becomes, say, Pi(r_1,r_2) = profits(c/r_1,c/r_2) - 0.5*r_1**2 with a similar expression for firm 2.
Now you can define reaction function r_1 = f(r_2) anf r_2 = f(r_1) where f(r_2) maximizes Pi(r_1,r_2) over r_1 for a given value of r_2. For this maximization you can use optimize.brute as you do.
As in the last Cournot stage, a symmetric equilibrium in R&D investments is a fixed point such that r = f(r).
Does this help?
Jan.