# OK to post homework
# RDB, 2022-02-06, HW5

# 1. By the way, this is not Bertrand Russel, but some other guy named
# Bertrand. Unfortunate, but true.

# Anyway, now the firms choose a *price*, and the consumers determine the
# demand.

# If Firm 1 chooses price p1 and Firm 2 chooses price p2, then the demands are

#       Qi(p1, p2) = a - pi + b * p(1 - i),

# for some parameters a, b > 0. The marginal cost of production is a constant
# c. So, the profit here is

#       Pi(p1, p2) = pi * Qi(p1, p2) - c Qi(p1, p2)
#                  = (pi - c) Qi(p1, p2).

# This is, just like like in the Cournot model, some quadratic.

Q1 := (p1, p2) -> a - p1 + b * p2:
Q2 := (p1, p2) -> a - p2 + b * p1:
P1 := (p1 - c) * Q1(p1, p2):
P2 := (p2 - c) * Q2(p1, p2):

bestResponseOne := solve(diff(P1, p1), p1):
bestResponseTwo := solve(diff(P2, p2), p2):

solve({p1 = bestResponseOne, p2 = bestResponseTwo}, {p1, p2});

# Solutions:

#       p1 = p2 = (a + c) / (2 - b).

# We can confirm this graphically by plotting the lines p1 = bestResponseOne
# and p2 = bestResponseTwo in the Cartesian plane.

with(plots):
bestResponseOneN := subs([b=3/2, a=10, c=5], bestResponseOne):
bestResponseTwoN := subs([b=3/2, a=10, c=5], bestResponseTwo):

evalf(solve({p1 = bestResponseOneN, p2 = bestResponseTwoN}, {p1, p2}));

implicitplot({p1 = bestResponseOneN, p2 = bestResponseTwoN}, p1=0..50, p2=0..50);

(*

 50+                                                               HH
   +                                                            HHHH
   +                                                          HHH        HHHH
   +                                                       HHH      HHHHHH
 40+                                                    HHHH   HHHHHH
   +                                                  HHH  HHHHH
   +                                               HHHHHHHH
   +                                            HHHHHHH
 30+                                        HHHHHH
   +                                    HHHHHH
   +                               HHHHHHHHH
   +                          HHHHHH HHHH
 20+                     HHHHHH    HHH
   +                 HHHHH      HHHH
   +            HHHHH        HHHH
   +       HHHHHH          HHH
 10+  HHHHHH            HHH
   *HHH              HHHH
   +              HHHH
   +            HHH
  -+--+--+--+-***-+--+--+--+-+--+--+--+--+--+--+--+-+--+--+--+--+--+--+--+--+-
                 10            20             30            40             50

*)

# 5.

u1 := q1 * (1 - q1 - q2)^d1:
u2 := q2 * (1 - q1 - q2)^d2:

bestResponse1 := solve(diff(u1, q1), q1):
bestResponse2 := solve(diff(u2, q2), q2):

NE := solve({q1 = bestResponse1, q2 = bestResponse2}, {q1, q2}):

nashQ1 := rhs(NE[1]);
nashQ2 := rhs(NE[2]);

print(`Check:`);
subs({d1=1, d2=1}, [nashQ1, nashQ2]);

# 6.

# What quantities maximize the total payoff?

total := u1 + u2;
gradient := simplify([diff(total, q1), diff(total, q2)]);
soln := solve(gradient, {q1, q2});

# There are three "extremal" regions for the total payoff:

#   1.
#       a. If d1 != d2, then the point

#               (q1, q2) = (1 / (d1 - d2), -1 / (d1 - d2)).

#       b. If d1 = d2, then the line q1 + q2 = 1 / (d1 + 1).

#   2. The line q1 = 1 - q2.
#   3. The boundaries q1 = 0 and q2 = 0.

# In reverse:

#   3. q1 = 0 gives total = q2(1 - q2)^d2, and the maximum here occurs at q2 =
#   1 / (d2 + 1), with a value of

#       d2^d2 / (d2 + 1)^(d2 + 1).

#   Similarly, q2 = 0 will give a maximum value of

#       d1^d1 / (d1 + 1)^(d1 + 1).

#   2. The line q1 = 1 - q2 gives a total payoff of 0.
#   1. a. The point (q1, q2) = (1 / (d1 - d2), -1 / (d1 - d2)) gives a total
#         payoff of 0.
#      b. The total payoff on this line is

#           d1^d1 / (d1 + 1)^(d1 + 1).

# So, in general, the maximum total payoff is

#       max(di^di / (di + 1)^(di + 1)).

# We have to compare these to the total pay-off at the Nash equilibrium.

maximumPayoff := (d1, d2) -> max(d1^d1 / (d1 + 1)^(d1 + 1), d2^d2 / (d2 + 1)^(d2 + 1)):
nashPayoffExpr := simplify(subs({q1=nashQ1, q2=nashQ2}, total));
nashPayoff := (d1N, d2N) -> simplify(subs({d1=d1N, d2=d2N}, nashPayoffExpr)):

# Now we want to know when nashPayoff(d1, d2) < maximumPayoff(d1, d2).

# We could probably figure this out rigorously, but that seems like an awful
# lot of work before we know what the answer is. So I made a 3D plot. (See
# attached plot.)

# It seems like the answer is this:

#       The Nash equilibrium is inefficient unless d1 = 0 or d2 = 0.

# It's pretty easy to check by hand that the Nash equilibrium is efficient if
# d1 = 0 or d2 = 0.

# I leave the remaining case as an exercise to the reader.