# OK to post homework
# Vikrant, Feb 06 2022, Assignment 5

# ================================================================================
# 0. Code that has been given.
# ================================================================================

read("C5.txt"):

# ================================================================================
# 1. Read and understand 1.2B of Gibbons.
#    Write code for payoffs.
#    Find best response and Nash Equilibrium, and confirm graphically.
# ================================================================================

(*

BertrandU1:= (p1,p2) -> (a-p1+b*p2)*(p1-c);
BertrandU2:= (p1,p2) -> (a-p2+b*p1)*(p2-c);
# This is a best response as it's feasible (>= 0), and U1 is a parabola with negative leading coefficient.
BertrandBR__1:= solve({diff(U1(p1,p2),p1)=0},{p1});
# This is a best response as it's feasible (>= 0), and U2 is a parabola with negative leading coefficient.
BertrandBR__2:= solve({diff(U2(p1,p2),p2)=0},{p2});
BertrandNE:= solve({diff(U1(p1,p2),p1)=0,diff(U2(p1,p2),p2)=0},{p1,p2});

with(plots,implicitplot,pointplot):
a:= 1: b:= 1.5: c:= 0.5:
curve:= implicitplot([diff(U1(p1,p2),p1)=0,diff(U2(p1,p2),p2)=0], p1=0..10, p2=0..10):
inter:= pointplot([rhs(BertrandNE[1])],[rhs(BertrandNE[2])],symbol=cross,symbolsize=20):
plots[display]({curve,inter});
a:= 'a': b:= 'b': c:= 'c':

*)

# ================================================================================
# 2. Read and understand C5.txt. 
#    Find Nash Equilibria, welfare maximizing pair, and BetterForBoth for three games.
# ================================================================================

(*

randomize(1900788349):

for i from 1 to 3 do
    G:= RG(30,30,10000):
    print(cat("Game ", i));
    print(cat("NEs: ", NE(G), "."));
    print(cat("Best welfare: ", BestTot(G)[1], "."));
    for ne in NE(G) do
        print(cat(ne, " is worse for both than ", BetterForBoth(G,ne[1],ne[2]), "."));
    od:
    print();
od:

*)

# ================================================================================
# 3. BeatNE.
# ================================================================================

BeatNE:=
    proc(a,b)
        local G:
        {seq(`if`(nops(NE(G)) = 1 and BetterForBoth(G,NE(G)[1][1],NE(G)[1][2]) <> {},G,NULL),G in AllGames(a,b))}:
    end:

(*

BeatNE(2,2);
# {[[[1, 1], [3, 2]], [[4, 3], [2, 4]]], [[[1, 1], [3, 4]], [[2, 3], [4, 2]]], [[[1, 2], [3, 4]], [[2, 3], [4, 1]]], [[[1, 3], [3, 4]], [[2, 2], [4, 1]]], [[[1, 4], [3, 3]], [[2, 2], [4, 1]]], [[[1, 4], [4, 3]], [[2, 2], [3, 1]]], [[[1, 4], [4, 3]], [[3, 2], [2, 1]]], [[[2, 1], [3, 2]], [[4, 3], [1, 4]]], [[[2, 2], [3, 1]], [[1, 4], [4, 3]]], [[[2, 2], [4, 1]], [[1, 3], [3, 4]]], [[[2, 2], [4, 1]], [[1, 4], [3, 3]]], [[[2, 3], [4, 1]], [[1, 2], [3, 4]]], [[[2, 3], [4, 2]], [[1, 1], [3, 4]]], [[[2, 4], [4, 3]], [[3, 2], [1, 1]]], [[[3, 1], [2, 2]], [[4, 3], [1, 4]]], [[[3, 2], [1, 1]], [[2, 4], [4, 3]]], [[[3, 2], [2, 1]], [[1, 4], [4, 3]]], [[[3, 3], [1, 4]], [[4, 1], [2, 2]]], [[[3, 4], [1, 1]], [[4, 2], [2, 3]]], [[[3, 4], [1, 2]], [[4, 1], [2, 3]]], [[[3, 4], [1, 3]], [[4, 1], [2, 2]]], [[[4, 1], [2, 2]], [[3, 3], [1, 4]]], [[[4, 1], [2, 2]], [[3, 4], [1, 3]]], [[[4, 1], [2, 3]], [[3, 4], [1, 2]]], [[[4, 2], [2, 3]], [[3, 4], [1, 1]]], [[[4, 3], [1, 4]], [[2, 1], [3, 2]]], [[[4, 3], [1, 4]], [[3, 1], [2, 2]]], [[[4, 3], [2, 4]], [[1, 1], [3, 2]]]}

BeatNE(2,3);
# [Length of output exceeds limit of 1000000]
# Not going to adjust precision. Instead: 
nops(BeatNE(2,3));
# 32400

*)

# ================================================================================
# 4. EstBeatNE.
# ================================================================================

EstBeatNE:=
    proc(a,b,K,M)
        local i,G:
        local Gs:= {seq(`if`(nops(NE(G)) = 1,G,NULL),G in {seq(RG(a,b,K),i=1..M)})}:
        evalf(
            [
                nops(Gs)/M,
                nops([seq(`if`(BestTot(G)[2] > convert(G[op(NE(G)[1])],`+`),1,NULL),G in Gs)])/M,
                nops([seq(`if`(BetterForBoth(G,NE(G)[1][1],NE(G)[1][2]) <> {},1,NULL),G in Gs)])/M
            ]
        ):
    end:

(*

randomize(2146731577):

for i from 1 to 5 do
    print(EstBeatNE(10,10,1000,1000));
od:

# Experiment 1
# [0.4060000000, 0.3390000000, 0.1060000000]
# Experiment 2
# [0.4130000000, 0.3350000000, 0.1060000000]
# Experiment 3
# [0.4260000000, 0.3630000000, 0.1180000000]
# Experiment 4
# [0.3920000000, 0.3240000000, 0.1030000000]
# Experiment 5
# [0.4270000000, 0.3530000000, 0.1070000000]

*)

# ================================================================================
# 5. Nash Equilibrium for {q1*(1-q1-q2)^(d1), q2*(1-q1-q2)^(d2)}.
# ================================================================================

(*

solve({diff(q1*(1-q1-q2)^(d1),q1)=0,diff(q2*(1-q1-q2)^(d2),q2)=0},{q1,q2});

*)

# ================================================================================
# 6. Conditions for when Nash Equilibria for {q1*(A-q1-q2)^(d1), q2*(A-q1-q2)^(d2)} are (not) welfare maximizing.
# ================================================================================

(*
    # We claim that Nash Equilibria in this game are never welfare maximizing.

    # Proof:
        # We prove this by casework on d1 =?= d2.

        # Case 1: d1 <> d2
            # For every Nash Equilibrium, we will show a solution with higher welfare.
            # Suppose (q1,q2) form a Nash Equilibrium (note that 0 < q1,q2; q1+q2 < A).
            # Let q = q1+q2.
            # Then the welfare is q1 * (A-q)^(d1) + q2 * (A-q)^(d2).
            # Assuming wlog that d1 > d2, we can see that:
                # If A-q < 1,
                    # Then (0,q) is a solution with higher welfare than (q1,q2).
                # If A-q = 1,
                    # Then (q,0) and (0,q) are solutions which achieve the same welfare as (q1,q2).
                    # From hw4 we know that (A/(d1+1),0) and (0,A/(d2+1)) are the unique welfare maximizing solutions to monopolies for firms 1 and 2 respectively.
                    # Using these facts we can deduce:
                        # q is different from at least one of A/(d1+1) or A/(d2+1).
                        # Therefore, the welfare for at least one of (A/(d1+1),0) or (0,A/(d2+1)) is strictly higher than that for (q1,q2).
                # If A-q > 1,
                    # Then (q,0) is a solution with higher welfare than (q1,q2).

        # Case 2: d1 = d2 (= d)
            # The welfare is (q1+q2) * (A-(q1+q2))^d,
            # and so from hw4 we know that the Nash Equilibrium welfare is 2*A/(d+2) * (d*A/(d+2))^d,
            # and the maximum welfare, by reducing to monopoly, is A/(d+1) * (d*A/(d+1))^d.
            # These cannot be equated when d > 0:
                #     2*A/(d+2) * (d*A/(d+2))^d = A/(d+1) * (d*A/(d+1))^d
                # iff 2/(d+2)^(d+1) = 1/(d+1)^(d+1)
                # iff 2 = (1 + 1/(d+1))^(d+1)
                # which, by Bernoulli's inequality, does not have a solution.
            # We are using Bernoulli's inequality in the following form: (1+x)^r > 1+rx for all r > 1 and x > 0.
            # Letting x = 1/(d+1) and r = d+1, we can then observe the insolvability in d of 2 = (1 + 1/(d+1))^(d+1).
*)

# ================================================================================
# oo. Helpers from past homeworks.
# ================================================================================

with(combinat, permute):

AllGames:=
    proc(a,b)
        local S,pi1,pi2:
        S:= permute(a*b):
        {seq(seq(MakeGame(pi1,pi2,a,b),pi2 in S),pi1 in S)}: 
    end:

MakeGame:=
    proc(pi1,pi2,a,b)
        local i,j:
        [seq([seq([pi1[b*i+j],pi2[b*i+j]],j=1..b)],i=0..a-1)]:
    end: