gt03_correlated_equilibria.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <center>Correlated equilibria</center>\n",
    "### <center>Alfred Galichon (NYU & Sciences Po) and Antoine Jacquet (Sciences Po)</center>\n",
    "## <center>'math+econ+code' masterclass series</center>\n",
    "#### <center>With python code examples</center>\n",
    "© 2018–2025 by Alfred Galichon. \n",
    "Past and present support from NSF grant DMS-1716489, ERC grant CoG-866274 are acknowledged, as well as inputs from contributors listed [here](http://www.math-econ-code.org/team).\n",
    "\n",
    "**If you reuse material from this masterclass, please cite as:**<br>\n",
    "Alfred Galichon, 'math+econ+code' masterclass series. https://www.math-econ-code.org/\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### References\n",
    "\n",
    "* Aumann (1974). \"Subjectivity and correlation in randomized strategies\". *Journal of Mathematical Economics*.\n",
    "* Forges, F. (1990). \"Correlated Equilibrium in Two-Person Zero-Sum Games\". *Econometrica*.\n",
    "\n",
    "### Learning objectives\n",
    "\n",
    "* Correlated equilibrium and linear programming formulation\n",
    "* Computation using Gurobi and simplex algorithm\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "#!pip install mec --upgrade\n",
    "import numpy as np\n",
    "import gurobipy as grb\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Game of chicken?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Correlated equilibrium in a two-player bimatrix game\n",
    "\n",
    "We consider the same set-up as in the lecture on Nash equilibria: a two-player game where player 1 plays $i \\in \\mathcal I$, player 2 plays $j \\in \\mathcal J$, and their payoffs are respectively $A_{ij}$ and $B_{ij}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A *correlated equilibrium* is a vector of probabilities $(\\pi_{ij})_{i \\in \\mathcal I, j \\in \\mathcal J}$ such that $\\sum_{ij} \\pi_{ij} = 1$ and\n",
    "\n",
    "\\begin{align}\n",
    "\\textstyle\\sum_j \\pi_{ij} A_{ij} &\\geq \\textstyle\\sum_j \\pi_{ij} A_{i' j} \\quad \\text{for all $i, i'$},  \\\\\n",
    "\\textstyle\\sum_i \\pi_{ij} B_{ij} &\\geq \\textstyle\\sum_i \\pi_{ij} B_{i j'} \\quad \\text{for all $j, j'$}.\n",
    "\\end{align}\n",
    "\n",
    "This can be interpreted as follows: with probability $\\pi_{ij}$, a mediator sends to player 1 a private signal instructing him to play $i$, and to player 2 a private signal instructing her to play $j$.\n",
    "Then $(\\pi_{ij})_{i \\in \\mathcal I, j \\in \\mathcal J}$ is a correlated equilibrium if no player has an ex-ante incentive to deviate from the instruction they received, assuming that the other player is going to follow it.\n",
    "\n",
    "*Remark.* Considering $\\pi$ as an $\\mathcal I \\times \\mathcal J$ matrix, the equilibrium conditions can be rewritten:\n",
    "\n",
    "$(A \\pi^\\top)_{ii} \\geq (A \\pi^\\top)_{i' i}$ for all $i, i'$  \n",
    "$(B^\\top \\pi)_{jj} \\geq (B^\\top \\pi)_{j' j}$ for all $j, j'$.\n",
    "\n",
    "*Exercise.* Show that if $(p,q)$ is a Nash equilibrium, then the vector $(\\pi_{ij})$ defined by $\\pi_{ij} = p_i q_j$ is a correlated equilibrium."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We import the `Bimatrix_game` class that we built to study Nash equilibria.\n",
    "(We use again the data from Palacios–Huerta on penalty kicks.)\n",
    "\n",
    "We then build a method `is_Aumann_eq` to check whether a vector $(\\pi_{ij})$ is a correlated equilibrium, and we test it on a Nash equilibrium."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A_i_j =\n",
      " [[ 53.21  71.35  93.8 ]\n",
      " [ 90.26  42.81  86.12]\n",
      " [ 96.88 100.    75.43]]\n",
      "B_i_j =\n",
      " [[ 96.79  28.65   6.2 ]\n",
      " [  9.74 107.19  13.88]\n",
      " [  3.12   0.    74.57]]\n"
     ]
    }
   ],
   "source": [
    "from mec.gt import Bimatrix_game\n",
    "\n",
    "penalty_nonzero_sum = Bimatrix_game.create_penalty_game(zero_sum=False)\n",
    "print('A_i_j =\\n', penalty_nonzero_sum.A_i_j)\n",
    "print('B_i_j =\\n', penalty_nonzero_sum.B_i_j)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Bimatrix_game_is_AumannEq(self, π_i_j, tol=1e-5):\n",
    "    if abs(π_i_j.sum() - 1) > tol :\n",
    "        print('π_i_j does not sum to 1.')\n",
    "        return False\n",
    "    for i in range(self.nbi) :\n",
    "        for k in range(self.nbi) :\n",
    "            if self.A_i_j[i,:] @ π_i_j[i,:] < self.A_i_j[k,:] @ π_i_j[i,:] - tol :\n",
    "                print('Player 1 can deviate from i = ' + str(i) + ' with k = ' + str(k) + '.')\n",
    "                return False\n",
    "    for j in range(self.nbj) :\n",
    "        for l in range(self.nbj) :\n",
    "            if self.B_i_j[:,j] @ π_i_j[:,j] < self.B_i_j[:,l] @ π_i_j[:,j] - tol :\n",
    "                print('Player 2 can deviate from j = ' + str(j) + ' with l = ' + str(l) + '.')\n",
    "                return False\n",
    "    return True\n",
    "\n",
    "Bimatrix_game.is_AumannEq = Bimatrix_game_is_AumannEq"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Nash_sol = penalty_nonzero_sum.lemke_howson_solve()\n",
    "p_i, q_j = Nash_sol['p_i'], Nash_sol['q_j']\n",
    "\n",
    "penalty_nonzero_sum.is_AumannEq(np.outer(p_i, q_j))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Linear programming formulation\n",
    "\n",
    "Let's write in matrix form the two inequalities which define the correlated equilibrium.\n",
    "\n",
    "The matrix corresponding to the first inequality is\n",
    "\n",
    "\\begin{equation}\n",
    "    \\left[ \\begin{smallmatrix}\n",
    "    A_{11}-A_{11} & A_{12}-A_{12} & \\dots & A_{1J}-A_{1J} & & & & & & & & & \\\\\n",
    "    A_{11}-A_{21} & A_{12}-A_{22} & \\dots & A_{1J}-A_{2J} & & & & & & & & & \\\\\n",
    "    \\vdots        & \\vdots        &       & \\vdots        & & & & & & & & & \\\\\n",
    "    A_{11}-A_{I1} & A_{12}-A_{I2} & \\dots & A_{1J}-A_{IJ} & & & & & & & & & \\\\\n",
    "    & & & & A_{21}-A_{11} & A_{22}-A_{12} & \\dots & A_{2J}-A_{1J} & & & & & \\\\\n",
    "    & & & & A_{21}-A_{21} & A_{22}-A_{22} & \\dots & A_{2J}-A_{2J} & & & & & \\\\\n",
    "    & & & & \\vdots        & \\vdots        &       & \\vdots        & & & & & \\\\\n",
    "    & & & & A_{21}-A_{I1} & A_{22}-A_{I2} & \\dots & A_{2J}-A_{IJ} & & & & & \\\\\n",
    "    & & & &               &               &       &               & \\ddots & & & & \\\\\n",
    "    & & & & & & & & & A_{I1}-A_{11} & A_{I2}-A_{12} & \\dots & A_{IJ}-A_{1J} \\\\\n",
    "    & & & & & & & & & A_{I1}-A_{21} & A_{I2}-A_{22} & \\dots & A_{IJ}-A_{2J} \\\\\n",
    "    & & & & & & & & & \\vdots        & \\vdots        &       & \\vdots        \\\\\n",
    "    & & & & & & & & & A_{I1}-A_{I1} & A_{I2}-A_{I2} & \\dots & A_{IJ}-A_{IJ} \\\\\n",
    "    \\end{smallmatrix} \\right].\n",
    "\\end{equation}\n",
    "\n",
    "Using Kronecker products, we are able to rewrite it as\n",
    "\n",
    "\\begin{equation}\n",
    "M_A = \\big[ (\\text I_{\\mathcal I} \\otimes 1_{\\mathcal J}^\\top) \\, \\text{diag}(A) \\big] \\otimes 1_{\\mathcal I} - \\text I_{\\mathcal I} \\otimes A\n",
    "\\end{equation}\n",
    "\n",
    "where $\\text{diag}(A)$ is the diagonal matrix of size $IJ \\times IJ$ obtained by lining the entries of the matrix $A$ along the diagonal (in row-major order).\n",
    "\n",
    "- Alternative with Hadamard product:\n",
    "\\begin{equation}\n",
    "M_A = \\big[ (\\text I_{\\mathcal I} \\otimes 1_{\\mathcal J}^\\top) \\odot ( 1_{\\mathcal I}^\\top \\otimes A) \\big] \\otimes 1_{\\mathcal I} - \\text I_{\\mathcal I} \\otimes A\n",
    "\\end{equation}\n",
    "\n",
    "- Alternative wih vectorization:  \n",
    "Let $M_1 = I_{\\mathcal I} \\otimes 1_{\\mathcal I} \\otimes I_{\\mathcal J}$,  \n",
    "$M_2 = 1_{\\mathcal I} \\otimes I_{\\mathcal I}  \\otimes I_{\\mathcal J}$,  \n",
    "so that $\\sum_{i, i', j} \\pi_{ij} (A_{ij} - A_{i' j}) = (M_1 \\text{vec}\\, \\pi)^\\top (M_1 - M_2)\\text{vec}\\, A$.  \n",
    "We want the sum to be only on $j$, for all $i, i'$ fixed:\n",
    "\\begin{equation}\n",
    "\\sum_j \\pi_{ij} (A_{ij} - A_{i' j}) = (M_1 \\text{vec}\\, \\pi)^\\top (S_i \\otimes S_{i'} \\otimes I_{\\mathcal J}) (M_1 - M_2)\\text{vec}\\, A\n",
    "\\end{equation}\n",
    "where $S_i$ is the matrix with 0 everywhere except at index $ii$, where it is 1.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The matrix corresponding to the second inequality is\n",
    "\n",
    "\\begin{equation}\n",
    "    \\left[ \\begin{smallmatrix}\n",
    "    B_{11}-B_{11} &               &        &               &\n",
    "    B_{21}-B_{21} &               &        &               &       &\n",
    "    B_{I1}-B_{I1} &               &        &                \\\\\n",
    "    B_{11}-B_{12} &               &        &               &\n",
    "    B_{21}-B_{22} &               &        &               &       &\n",
    "    B_{I1}-B_{I2} &               &        &                \\\\\n",
    "    \\vdots        &               &        &               &\n",
    "    \\vdots        &               &        &               & \\dots &\n",
    "    \\vdots        &               &        &                \\\\\n",
    "    B_{11}-B_{1J} &               &        &               &\n",
    "    B_{21}-B_{2J} &               &        &               &       &\n",
    "    B_{I1}-B_{IJ} &               &        &                \\\\\n",
    "                  & B_{12}-B_{11} &        &               &\n",
    "                  & B_{22}-B_{21} &        &               &       &\n",
    "                  & B_{I2}-B_{I1} &        &                \\\\\n",
    "                  & B_{12}-B_{12} &        &               &\n",
    "                  & B_{22}-B_{22} &        &               &       &\n",
    "                  & B_{I2}-B_{I2} &        &                \\\\\n",
    "                  & \\vdots        &        &               &\n",
    "                  & \\vdots        &        &               & \\dots &\n",
    "                  & \\vdots        &        &                \\\\\n",
    "                  & B_{12}-B_{1J} &        &               &\n",
    "                  & B_{22}-B_{2J} &        &               &       &\n",
    "                  & B_{I2}-B_{IJ} &        &                \\\\\n",
    "                  &               & \\ddots &               &\n",
    "                  &               & \\ddots &               &       &\n",
    "                  &               & \\ddots &                \\\\\n",
    "                  &               &        & B_{1J}-B_{11} &\n",
    "                  &               &        & B_{2J}-B_{21} &       &\n",
    "                  &               &        & B_{IJ}-B_{I1}  \\\\\n",
    "                  &               &        & B_{1J}-B_{12} & \n",
    "                  &               &        & B_{2J}-B_{22} &       &\n",
    "                  &               &        & B_{IJ}-B_{I2}  \\\\\n",
    "                  &               &        & \\vdots        & \n",
    "                  &               &        & \\vdots        & \\dots &\n",
    "                  &               &        & \\vdots         \\\\\n",
    "                  &               &        & B_{1J}-B_{1J} & \n",
    "                  &               &        & B_{2J}-B_{2J} &       &\n",
    "                  &               &        & B_{IJ}-B_{IJ}\n",
    "    \\end{smallmatrix} \\right].\n",
    "\\end{equation}\n",
    "\n",
    "As it is, there is no simple way to rewrite this using Kronecker products.\n",
    "However, if we reorder the lines as $jj' = 11, 21, \\dots, J1, 12, \\dots$, then we are able to rewrite this matrix as\n",
    "\\begin{equation}\n",
    "M_B = 1_{\\mathcal J} \\otimes \\big[ (1_{\\mathcal I}^\\top \\otimes \\text I_{\\mathcal J}) \\, \\text{diag}(B) \\big] -  B^\\top \\otimes \\text I_{\\mathcal J}\n",
    "\\end{equation}\n",
    "\n",
    "where $\\text{diag}(B)$ is defined similarly as $\\text{diag}(A)$ (also using the row-major order)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finding a correlated equilibrium is a linear feasibility problem, which we can always reformulate as a linear programming problem:\n",
    "\n",
    "\\begin{align}\n",
    "\\max_{\\pi_{ij} \\geq 0} ~ & 1_{\\mathcal I \\mathcal J}^\\top \\pi \\\\\n",
    "\\text{s.t.} ~ & 1_{\\mathcal I \\mathcal J}^\\top \\pi \\leq 1 \\\\\n",
    "& M_A \\, \\pi \\geq 0_{\\mathcal I^2} \\\\\n",
    "& M_B \\, \\pi \\geq 0_{\\mathcal J^2}\n",
    "\\end{align}\n",
    "\n",
    "where $\\pi$ is understood as a vector.\n",
    "We know that the value of this program is 1 (since it is attained for $\\pi$ constructed from a Nash equilibrium, for instance), and therefore the constraint $1_{\\mathcal I \\mathcal J}^\\top \\pi \\leq 1$ must bind."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Zero-sum games\n",
    "\n",
    "In zero-sum games, the margins of correlated equilibria constitute Nash equilibrium strategies, or optimal strategies.\n",
    "\n",
    "**Proposition.** If $(\\pi_{ij})$ is a correlated equilibrium, then $p_i = \\sum_j \\pi_{ij}$ is an optimal strategy for player 1, and $q_j = \\sum_i \\pi_{ij}$ is an optimal strategy for player 2.\n",
    "\n",
    "*Proof.*\n",
    "Recall that in zero-sum games, optimal strategies are solutions to the max-min problem\n",
    "\n",
    "\\begin{equation*}\n",
    "\\max_{p \\in \\Delta_I} \\, \\min_{q \\in \\Delta_J} ~ p^\\top A q,\n",
    "\\end{equation*}\n",
    "\n",
    "and any couple of solutions $(p,q)$ forms a Nash equilibrium.\n",
    "\n",
    "Here, the correlated equilibrium conditions rewrite as:  \n",
    "$\\sum_j \\pi_{ij} A_{ij} \\geq \\sum_j \\pi_{ij} A_{i' j}$ for all $i$, $i'$  \n",
    "$\\sum_i \\pi_{ij} A_{ij} \\leq \\sum_i \\pi_{ij} A_{i j'}$ for all $j$, $j'$.\n",
    "\n",
    "By summing these conditions on $i$ and $j$ respectively, and defining $p_i = \\sum_j \\pi_{ij}$ and $q_j = \\sum_i \\pi_{ij}$, we obtain that\n",
    "\n",
    "\\begin{equation*}\n",
    "\\sum_j q_j A_{i' j} \\leq \\sum_{ij} \\pi_{ij} A_{ij} \\leq \\sum_i p_i A_{i j'} \\quad \\text{for all $i'$ and $j'$},\n",
    "\\end{equation*}\n",
    "and as a consequence \n",
    "\\begin{equation*}\n",
    "\\max_{\\tilde p \\in \\Delta_I} \\sum_{ij} \\tilde p_i A_{ij} q_j\n",
    "\\leq \\sum_{ij} \\pi_{ij} A_{ij}\n",
    "\\leq \\min_{\\tilde q \\in \\Delta_J} \\sum_{ij} p_i A_{ij} \\tilde q_j,\n",
    "\\end{equation*}\n",
    "\n",
    "with equality when $\\tilde p = p$ and $\\tilde q = q$.\n",
    "Thus $p$ and $q$ are best responses to each other, so they are optimal strategies. *Q.E.D.*\n",
    "\n",
    "An immediate consequence is that in a zero-sum game, the payoff of the players at a correlated equilibrium is equal to the value of the game:\n",
    "\n",
    "\\begin{equation*}\n",
    "\\sum_{ij} \\pi_{ij} A_{ij}\n",
    "= \\sum_{ij} p_i A_{ij} q_j.\n",
    "\\end{equation*}\n",
    "\n",
    "We also have the following result:\n",
    "\n",
    "**Proposition (Forges 1990).**\n",
    "Let $\\pi$ be a correlated equilibrium of a zero-sum game, and let $j$ such that $q_j := \\sum_i \\pi_{ij} > 0$.\n",
    "Then the strategy $(\\pi_{ij}/q_j)_{i \\in \\mathcal I}$ is a Nash equilibrium strategy for player 1.\n",
    "\n",
    "*Proof.*\n",
    "Let $j$ such that $q_j > 0$ and denote $V$ the value of the game.\n",
    "For player 1, the expected payoff of playing $(\\pi_{ij}/q_j)_{i \\in \\mathcal I}$ against $q$ is\n",
    "\n",
    "\\begin{equation*}\n",
    "\\sum_{il} \\frac{\\pi_{ij}}{q_j} q_l A_{il}\n",
    "= \\sum_i \\frac{\\pi_{ij}}{q_j} \\underbrace{\\sum_l q_l A_{il}}_{= V \\text{ if $p_i > 0$}}\n",
    "= \\sum_i \\frac{\\pi_{ij}}{q_j} V = V\n",
    "\\qquad \\text{since $\\pi_{ij} > 0$ implies $p_i > 0$.}\n",
    "\\end{equation*}\n",
    "\n",
    "(Note that the bracketed equality is true because all pure strategies played with positive probability in a Nash equilibrium have the same expected payoff.)  \n",
    "Hence $(\\pi_{ij}/q_j)_{i \\in \\mathcal I}$ is a best response to the optimal strategy $q$, and therefore it is itself an optimal strategy. *Q.E.D.*\n",
    "\n",
    "An immediate consequence of the proposition above is: \n",
    "\n",
    "**Corollary.** If a zero-sum game has a unique Nash equilibrium, then it also has a unique correlated equilibrium (which must then be given by $\\pi_{ij} = p_i q_j$).\n",
    "\n",
    "*Proof.*\n",
    "If $q_j > 0$ we must have $\\pi_{ij}/q_j = p_i$ by uniqueness of player 1's optimal strategy, hence $\\pi_{ij} = p_i q_j$.  \n",
    "If $q_j = 0$ then necessarily $\\pi_{ij} = 0$ (since $\\sum_i \\pi_{ij} = q_j$) so we also have $\\pi_{ij} = p_i q_j$.\n",
    "*Q.E.D.*\n",
    "\n",
    "-------"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, the correlated equilibrium conditions rewrite as:  \n",
    "$\\sum_j \\pi_{ij} A_{ij} \\geq \\sum_j \\pi_{ij} A_{i' j}$ for all $i$, $i'$  \n",
    "$\\sum_i \\pi_{ij} A_{ij} \\leq \\sum_i \\pi_{ij} A_{i j'}$ for all $j$, $j'$.\n",
    "\n",
    "$\\max_{\\pi_{ij} \\geq 0} \\min_{\\alpha_{ii'}, \\beta_{jj'} \\geq 0}\n",
    "~\\alpha_{ii'} (\\sum_j \\pi_{ij} A_{ij} - \\sum_j \\pi_{ij} A_{i' j})\n",
    "- \\beta_{jj'} (\\sum_i \\pi_{ij} A_{ij} - \\sum_i \\pi_{ij} A_{i j'})$\n",
    "\n",
    "$\\sum_j \\pi_{ij} > 0 \\implies \\sum_j \\pi_{ij} A_{ij} = V_1$\n",
    "\n",
    "Recall the linear program associated with finding player 1's Nash strategy in a zero-sum game:\n",
    "\n",
    "\\begin{equation}\n",
    "V_1 = \\max_{ p_i \\geq 0, U} ~ U \n",
    "\\qquad \\text{s.t.} ~ U \\leq \\textstyle\\sum_i A_{ij} p_i \\quad (\\forall j),\n",
    "\\quad \\textstyle\\sum_i p_i = 1.\n",
    "\\end{equation}\n",
    "\n",
    "We can introduce new variables $\\pi_{ij} \\geq 0$ verifying $\\sum_j \\pi_{ij} = p_i$:\n",
    "\n",
    "\\begin{align}\n",
    "V_1 &= \\max_{ p_i \\geq 0, \\pi_{ij} \\geq 0, U} ~ U \n",
    "\\qquad \\text{s.t.} ~ U \\leq \\textstyle\\sum_i A_{ij} p_i \\quad (\\forall j),\n",
    "\\quad \\textstyle\\sum_i p_i = 1,\n",
    "\\quad \\textstyle\\sum_j \\pi_{ij} = p_i \\\\\n",
    "&= \\max_{ p_i \\geq 0, \\pi_{ij} \\geq 0, U} ~ U \n",
    "\\qquad \\text{s.t.} ~ U \\leq \\textstyle\\sum_i A_{ij'} (\\textstyle\\sum_j \\pi_{ij}) \\quad (\\forall j'),\n",
    "\\quad \\textstyle\\sum_{ij} \\pi_{ij} = 1,\n",
    "\\quad \\textstyle\\sum_j \\pi_{ij} = p_i \\\\\n",
    "&= \\max_{ \\pi_{ij} \\geq 0, U} ~ U \n",
    "\\qquad \\text{s.t.} ~ U \\leq \\textstyle\\sum_j \\big(\\textstyle\\sum_i A_{ij'} \\pi_{ij}\\big) \\quad (\\forall j'),\n",
    "\\quad \\textstyle\\sum_{ij} \\pi_{ij} = 1.\n",
    "\\end{align}\n",
    "\n",
    "We can consider another program, which is more constrained:\n",
    "\n",
    "\\begin{align}\n",
    "\\max_{ \\pi_{ij} \\geq 0, U} ~ U \n",
    "\\qquad \\text{s.t.} ~ \\big(\\textstyle\\sum_i \\pi_{ij'}\\big) U \\leq \\textstyle\\sum_i A_{ij'} \\pi_{ij} \\quad (\\forall j, j'),\n",
    "\\quad \\textstyle\\sum_{ij} \\pi_{ij} = 1.\n",
    "\\end{align}\n",
    "\n",
    "$\\min_{\\pi_{ij}, U_i, V_j} ~ \\sum_i (\\sum_j \\pi_{ij}) U_i - \\sum_j (\\sum_i \\pi_{ij}) V_j$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{align}\n",
    "\\max_{ \\pi_{ij}, U_i, V_j } &~ \\sum_i (\\sum_j \\pi_{ij}) U_i - \\sum_j (\\sum_i \\pi_{ij}) V_j \\\\\n",
    "\\qquad \\text{s.t.} &~ \\big(\\textstyle\\sum_j \\pi_{ij}\\big) U_i \\geq \\textstyle\\sum_j \\pi_{ij} A_{ij'} \\quad (\\forall i, i'),\n",
    "\\qquad \\text{s.t.} &~ \\big(\\textstyle\\sum_i \\pi_{ij}\\big) V_j \\leq \\textstyle\\sum_i \\pi_{ij} A_{ij'} \\quad (\\forall j, j'),\n",
    "\\quad \\textstyle\\sum_{ij} \\pi_{ij} = 1.\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "--------\n",
    "\n",
    "In zero-sum games, the margins of correlated equilibria coincide with Nash equilibria, in the sense that if $(\\pi_{ij})$ is a correlated equilibrium, then $p_i = \\sum_j \\pi_{ij}$ and $q_j = \\sum_i \\pi_{ij}$ constitute a Nash equilibrium.\n",
    "An immediate consequence is that in a zero-sum game, the payoff of the players at a correlated equilibrium is equal to the value of the game.\n",
    "\n",
    "Recall the conditions for $(p,q)$ to be a Nash equilibrium:\n",
    "\n",
    "$\\sum_{ij} p_i A_{ij} q_j \\geq \\sum_j A_{i'j} q_j$ for all $i'$\n",
    "\n",
    "$\\sum_{ij} p_i B_{ij} q_j \\geq \\sum_j B_{ij'} q_j$ for all $j'$.\n",
    "\n",
    "Let's show that the first condition holds.\n",
    "Fix $i'$, then for all $j$ we have\n",
    "\n",
    "\\begin{align}\n",
    "\\sum_i p_i A_{ij}\n",
    "&= \\sum_{il} A_{ij} \\pi_{il} \\\\\n",
    "&\\geq \\sum_{il} A_{il} \\pi_{il}\n",
    "\\quad \\text{since } \\sum_i B_{il} \\pi_{il} \\geq \\sum_i B_{ij} \\pi_{il} \\text{ and } B_{ij} = -A_{ij} \\\\\n",
    "&\\geq \\sum_{il} A_{i' l} \\pi_{il}\n",
    "\\quad \\text{since } \\sum_l A_{il} \\pi_{il} \\geq \\sum_l A_{i' l} \\pi_{il} \\\\\n",
    "&= \\sum_l A_{i' l} q_l\n",
    "\\end{align}\n",
    "hence\n",
    "\\begin{equation}\n",
    "\\sum_{ij} p_i A_{ij} q_j\n",
    "= \\sum_j \\left( \\sum_i p_i A_{ij} \\right) q_j\n",
    "\\geq \\left( \\sum_l A_{i' l} q_l \\right) \\sum_j q_j\n",
    "= \\sum_l A_{i' l} q_l.\n",
    "\\end{equation}\n",
    "\n",
    "The second condition holds as well by symmetry.\n",
    "\n",
    "**Question.** If we have a Nash $(p,q)$, does any $\\pi$ such that its margins are $p$ and $q$ constitute a correlated equilibrium?\n",
    "\n",
    "**Proposition (Forges 1990).**\n",
    "Let $\\pi$ be a correlated equilibrium of a zero-sum game, and let $j$ such that $q_j := \\sum_i \\pi_{ij} > 0$.\n",
    "Then the strategy $(\\pi_{ij}/q_j)_{i \\in \\mathcal I}$ is a Nash equilibrium strategy for player 1.\n",
    "\n",
    "An immediate consequence of the proposition above is that, if a zero-sum game has a unique Nash equilibrium, then it also has a unique correlated equilibrium (which must then be $\\pi_{ij} = p_i q_j$).\n",
    "Indeed, if $q_j > 0$ we must have $\\pi_{ij}/q_j = p_i$ by uniqueness of player 1's Nash equilibrium strategy, hence $\\pi_{ij} = p_i q_j$.\n",
    "If $q_j = 0$ then $\\pi_{ij} = 0$ since $\\sum_i \\pi_{ij} = q_j$ so we also have $\\pi_{ij} = p_i q_j$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Computing correlated equilibria"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Computation using Gurobi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Bimatrix_game_solve_AumannEq(self, verbose=0):\n",
    "    model = grb.Model()\n",
    "    model.Params.OutputFlag = 0\n",
    "    π = model.addMVar(shape=(self.nbi,self.nbj))\n",
    "    model.setObjective(π.sum(), grb.GRB.MAXIMIZE)\n",
    "    model.addConstr(π.sum() <= 1)\n",
    "    for i in range(self.nbi) :\n",
    "        for k in range(self.nbi) :\n",
    "            if k != i : model.addConstr(self.A_i_j[i,:] @ π[i,:] >= self.A_i_j[k,:] @ π[i,:])\n",
    "    for j in range(self.nbj) :\n",
    "        for l in range(self.nbj) :\n",
    "            if l != j : model.addConstr(self.B_i_j[:,j] @ π[:,j] >= self.B_i_j[:,l] @ π[:,j])\n",
    "    model.optimize() \n",
    "    π_i_j = np.array(model.getAttr('x')).reshape(self.nbi,self.nbj)\n",
    "    if verbose > 0: print('π_i_j =\\n', π_i_j)\n",
    "    return {'π_i_j': π_i_j,\n",
    "            'val1': self.A_i_j.flatten() @ π_i_j.flatten(),\n",
    "            'val2': self.B_i_j.flatten() @ π_i_j.flatten()}\n",
    "\n",
    "Bimatrix_game.solve_AumannEq = Bimatrix_game_solve_AumannEq"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Set parameter Username\n",
      "Academic license - for non-commercial use only - expires 2025-12-03\n",
      "π_i_j =\n",
      " [[0.0739268  0.03428835 0.22921855]\n",
      " [0.0545915  0.02532035 0.16926722]\n",
      " [0.09056711 0.04200638 0.28081374]]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "{'π_i_j': array([[0.0739268 , 0.03428835, 0.22921855],\n",
       "        [0.0545915 , 0.02532035, 0.16926722],\n",
       "        [0.09056711, 0.04200638, 0.28081374]]),\n",
       " 'val1': 82.62606457928054,\n",
       " 'val2': 36.37698012002112}"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sol = penalty_nonzero_sum.solve_AumannEq(verbose=1)\n",
    "π_i_j, val1, val2 = sol['π_i_j'], sol['val1'], sol['val2']\n",
    "\n",
    "sol"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "penalty_nonzero_sum.is_AumannEq(π_i_j)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "penalty_nonzero_sum.is_NashEq(π_i_j.sum(axis=1), π_i_j.sum(axis=0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.0739268 , 0.03428835, 0.22921855],\n",
       "       [0.0545915 , 0.02532035, 0.16926722],\n",
       "       [0.09056711, 0.04200638, 0.28081374]])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.outer(π_i_j.sum(axis=1), π_i_j.sum(axis=0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Computation using the simplex\n",
    "\n",
    "Recall the three inequalities we need to implement:\n",
    "\n",
    "$1_{\\mathcal{IJ}}^\\top \\, \\pi \\leq 1$\n",
    "\n",
    "$\\big( \\big[ (\\text I_{\\mathcal I} \\otimes 1_{\\mathcal J}^\\top) \\, \\text{diag}(A) \\big] \\otimes 1_{\\mathcal I} - \\text I_{\\mathcal I} \\otimes A\\big) \\, \\pi \\geq 0$\n",
    "\n",
    "$\\big( 1_{\\mathcal J} \\otimes \\big[ (1_{\\mathcal I}^\\top \\otimes \\text I_{\\mathcal J}) \\, \\text{diag}(B) \\big] -  B^\\top \\otimes \\text I_{\\mathcal J} \\big) \\, \\pi \\geq 0$.\n",
    "\n",
    "The associated linear program is degenerate (because some inequalities have a constant 0 on the right-hand side) so we solve an approximation with the simplex."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "from mec.lp import Tableau\n",
    "\n",
    "def Bimatrix_game_simplex_solve_AumannEq(self, verbose=0, eps=1e-5):\n",
    "    M_1 = np.ones(self.nbi*self.nbj)\n",
    "    M_A = np.kron( np.kron(np.eye(self.nbi), np.ones((1, self.nbj))) @ np.diag(self.A_i_j.flatten()),\n",
    "                   np.ones((self.nbi, 1))\n",
    "                 ) - np.kron(np.eye(self.nbi), self.A_i_j)\n",
    "    M_B = np.kron( np.ones((self.nbj, 1)),\n",
    "                   np.kron(np.ones((1, self.nbi)), np.eye(self.nbj)) @ np.diag(self.B_i_j.flatten())\n",
    "                 ) - np.kron(self.B_i_j.T, np.eye(self.nbj))\n",
    "    tab = Tableau(A_i_j = np.vstack([M_1, -M_A, -M_B]),\n",
    "                  b_i = np.eye(1+self.nbi**2+self.nbj**2)[0] + eps*(1-np.eye(1+self.nbi**2+self.nbj**2)[0]),\n",
    "                  c_j = np.ones(self.nbi*self.nbj),\n",
    "                  decision_var_names_j = ['π_'+str(i)+'_'+str(j) for j in range(self.nbj) for i in range(self.nbi)])\n",
    "    πstar, _, _ = tab.simplex_solve()\n",
    "    π_i_j = np.array(πstar.reshape(self.nbi,self.nbj))\n",
    "    if verbose > 0: print('π_i_j = \\n', π_i_j)\n",
    "    return {'π_i_j': π_i_j,\n",
    "            'val1': self.A_i_j.flatten() @ π_i_j.flatten(),\n",
    "            'val2': self.B_i_j.flatten() @ π_i_j.flatten()}\n",
    "\n",
    "Bimatrix_game.simplex_solve_AumannEq = Bimatrix_game_simplex_solve_AumannEq"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "π_i_j = \n",
      " [[0.07392655 0.03428756 0.22921896]\n",
      " [0.05459142 0.02532056 0.16926736]\n",
      " [0.09056694 0.04200653 0.28081412]]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "{'π_i_j': array([[0.07392655, 0.03428756, 0.22921896],\n",
       "        [0.05459142, 0.02532056, 0.16926736],\n",
       "        [0.09056694, 0.04200653, 0.28081412]]),\n",
       " 'val1': 82.62607459750402,\n",
       " 'val2': 36.37698676485681}"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sol = penalty_nonzero_sum.simplex_solve_AumannEq(verbose=1, eps=1e-5)\n",
    "π_i_j, val1, val2 = sol['π_i_j'], sol['val1'], sol['val2']\n",
    "\n",
    "sol"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Player 1 can deviate from i = 0 with k = 1.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "penalty_nonzero_sum.is_AumannEq(π_i_j)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Pure strategy i = 0 beats p_i.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "penalty_nonzero_sum.is_NashEq(π_i_j.sum(axis=1), π_i_j.sum(axis=0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Benchmarking"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "I, J = 10, 7\n",
    "random_game = Bimatrix_game(A_i_j = 100*np.random.rand(I,J),\n",
    "                            B_i_j = 100*np.random.rand(I,J))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Simplex"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "sol = random_game.simplex_solve_AumannEq(eps=1e-8)\n",
    "π_i_j, val1, val2 = sol['π_i_j'], sol['val1'], sol['val2']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "random_game.is_AumannEq(π_i_j)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Pure strategy i = 4 beats p_i.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "random_game.is_NashEq(π_i_j.sum(axis=1), π_i_j.sum(axis=0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Gurobi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "sol = random_game.solve_AumannEq()\n",
    "π_i_j, val1, val2 = sol['π_i_j'], sol['val1'], sol['val2']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "random_game.is_AumannEq(π_i_j)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Pure strategy i = 2 beats p_i.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "random_game.is_NashEq(π_i_j.sum(axis=1), π_i_j.sum(axis=0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}