# ECON 514 Julia Tutorial

Couple of things before we start. If this kind of coding is new to you (i.e. not straight up using one line from packages like Stata or some of R), in my experience I found several things quite useful.

- Try to change my code and see what works or what breaks. Does it do what you intend it to?
- If there is a block of code that does not run or is not giving what you expected, maybe use @show to check out what specific steps are doing. You can also do this if there are lines for which you are not sure what they are doing.
- A significant part of knowing how to code is knowing how to google.
- There is a ton of Julia documentation available online. They are mostly readable, so give it a go if you are stuck.


This tutorial is not an exhaustive guide to Julia. You are not meant to be able to right away solve problems. This tutorial is meant to give you some idea of

- what a function is
- how they fit together to do stuff
- solve a system of equations in SymPy
- pickup some syntax(?)

We are going to do this by solving a game from Running Man, a Korean gameshow.

## Problem (Copied from Mike's message on Discord)

There are 8 rooms and 8 contestants.  Each person must choose to go to one of the rooms (assume the choice is once for all and each person decides simultaneously).  Once everyone is in their room, a prize in each room is split among all the players in that room.  The prize is room 1 is 40K (won I assume), the prize in rooms 2,3 and 4 is 50K won, the prize in rooms 5,6 and 7 is 70K won, while the prize in room 8 is 100K won. So if all 8 people go to room 8, each of them will get 100/8 K won. We want to

### i) Find a pure strategy equilibrium

### ii) Find the unique symmetric mixed equilibrium.


We start by defining some parameters.


```julia
valuations = [40,50,50,50,70,70,70,100] 
options = 1:1:8

# These are arrays/vetors. Note that they have square brackets. If we use normal brackets, 
# they become a list of objects ("tuples"). For example,

tup = (1,"Wednesday 8:30 class...",[1 3 2;1 2 3])

# On one hand they do have something in common. To access the second element of say valuations, we can do 

@show valuations[2]
tup[2]
```

    valuations[2] = 50





    "Wednesday 8:30 class..."



At other times they have to be handled quite differently. For arrays of numbers, we can do calculations. Note that the default operation is linear algebraic. In R, when you have an array and you do something like valuations^2, you will square each element. Here, it will not make sense. Instead, we have to put a dot to say that we want to apply a function element-by-element. 

This will NOT run

       valuations^2


```julia
# Element by Element
@show valuations.^2

# v' means v transpose. Therefore this is the dot product with itself
valuations' * valuations
```

    valuations .^ 2 = [1600, 2500, 2500, 2500, 4900, 4900, 4900, 10000]





    33800



In general, stuff like multiplication and addition with tuples will not work. However, functions may work depending on how they are written.


```julia
a = (1,2,3)
b = (4,5,6)

@show typeof(a)

a'*b
```




We solve for pure NE using the following algorithm. 

(i) Start off by placing all of them in one room. 

(ii) In each round, we let players 1 to 8 move in turn into rooms that maximize their own payoff given others' current position. 

(iii) Repeat (ii) until nobody wants to move.


To implement step 2, we need to calculate how strategies map to payoffs, and figure out the best responses.


```julia
function room_population(s)

    """ 
        Input: s = vector of dimension 8x1. Each component is an integer from 1 to 8 representating a player's room choice.

        Output: pop = vector of dimension 8x1 representing the number of people in each room 
    """
    
    pop = zeros(8)

    for i =1:8
        pop[i] = sum(s.==i)
    end

    return pop
end
```




    room_population (generic function with 1 method)




```julia
@show 3 == 2
@show [1,2,3] == 2

s = [1,1,2,3,4,6,7,8]

@show s.==1

room_population(s)
```

    3 == 2 = false
    [1, 2, 3] == 2 = false
    s .== 1 = Bool[1, 1, 0, 0, 0, 0, 0, 0]





    8-element Vector{Float64}:
     2.0
     1.0
     1.0
     1.0
     0.0
     1.0
     1.0
     1.0




```julia
function purepayoff(s) 

    """ 
        Input: s = vector of dimension 8x1. Each component is an integer from 1 to 8 representating a player's room choice.

        Output: payoff_from_rooms = 8x8 matrix. ith column represents players payoffs from room i
                payoff = vector of dimension 8x1 each representating the payoff of players

    """

    pop = room_population(s)   
    payoff_from_rooms = zeros(8,8)

    for i =1:8
        # payoff_from_rooms[:,i] = valuations[i]/pop[i] * (s==i)  # This is wrong
        payoff_from_rooms[:,i] = valuations[i]/pop[i] * (s.==i)
    end

    payoff = (sum(payoff_from_rooms,dims = 2))[:,1]
    
    return payoff_from_rooms , payoff
end

function purepayoff(si,s_noti) 
    s = vcat(si,s_noti)
    return purepayoff(s)[2]
end

payoff_from_room,payoff_vector = purepayoff(s)

@show typeof(payoff_from_room)
@show payoff_vector

payoff_from_room

```

    typeof(payoff_from_room) = Matrix{Float64}
    payoff_vector = [20.0, 20.0, 50.0, 50.0, 50.0, 70.0, 70.0, 100.0]





    8×8 Matrix{Float64}:
     20.0   0.0   0.0   0.0  0.0   0.0   0.0    0.0
     20.0   0.0   0.0   0.0  0.0   0.0   0.0    0.0
      0.0  50.0   0.0   0.0  0.0   0.0   0.0    0.0
      0.0   0.0  50.0   0.0  0.0   0.0   0.0    0.0
      0.0   0.0   0.0  50.0  0.0   0.0   0.0    0.0
      0.0   0.0   0.0   0.0  0.0  70.0   0.0    0.0
      0.0   0.0   0.0   0.0  0.0   0.0  70.0    0.0
      0.0   0.0   0.0   0.0  0.0   0.0   0.0  100.0



Once we have the payoff function, we can see what the best response is by checking the payoff of each action, and pick the best one.


```julia
function BR(s_noti) 
    f(x) = purepayoff(x,s_noti)[1]
    rooms = 1:1:8
    payoffs = f.(rooms)
    br_payoff,br = findmax(payoffs)
    return br_payoff,br
end

A = [1,3,5,5,2,5,1]
findmax(A)
```




    (5, 3)




```julia
output = BR([1,2,3,4,5,6,7])
@show typeof(output)

@show output[1]
output[2]
```

    typeof(output) = Tuple{Float64, Int64}
    output[1] = 100.0





    8



Since now that given other player's strategy, we know what is the best response payoff, we can check if a player's action is a best response.


```julia
function is_BR(si,s_noti)
    """
    Input: si = strategy of the player we want to check
           s_noti = strategy of other players. Vector!

    Output: check = Boolean whether player is playing BR
            br_payoff = player's BR payoff
    """
    
    br_payoff,br = BR(s_noti)
    check = (purepayoff(si,s_noti)[1] == br_payoff)
    return check, br_payoff
end

is_BR(2,[1,3,4,3,5,6,7])
```




    (false, 100.0)



We now have what we need to do step 2. 


```julia
function player_move(s,k)
    """
    Input: s = current strategy of players
           k = which player is currently allowed to move
    
    Output: s = updated strategy of players after player k's turn
            move = Boolean whether player k has moved
    """

    si = s[k]
    s_noti = s[1:end .!=k]

    check = is_BR(si,s_noti)[1]

    if check == true
        move = false
    else
        s[k] = BR(s_noti)[2]
        move = true
    end

    return s,move
end

player_move([2,3,4,5,6,2,4,7],2)       
```




    ([2, 8, 4, 5, 6, 2, 4, 7], true)




```julia
function pureNE_search(s0)
    stable = false
    count = 0
    while stable == false
        moves = ones(8)

        for i = 1:8
            s,moves[i] = player_move(s0,i)
            s0 = s
        end

        stable = (maximum(moves)==0)
        count = count+1
    end
    return s0,count
end

#s0 = [8,8,8,8,8,8,8,8]
s0 = [1,1,1,1,1,1,1,1]
        
@show pure_ne = pureNE_search(s0)[1]
@show room_population(pure_ne)
purepayoff(pure_ne)[2]
```

    pure_ne = (pureNE_search(s0))[1] = [8, 5, 6, 7, 2, 3, 4, 8]
    room_population(pure_ne) = [0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 2.0]





    8-element Vector{Float64}:
     50.0
     70.0
     70.0
     70.0
     50.0
     50.0
     50.0
     50.0



To fixed symmetric mixed equilibrium, suppose the probabilites are denoted 

$$
p = (p_1,p_2,p_3,p_4,p_5,p_6,p_7,p_8)
$$

In the set of i such that $p_i>0$, it has to be that any two actions in that set give the same expected payoff when others are playing mixed strategy p. Now, the expected payoff from choosing room i is

$$
   \overset{7}{\underset{j=0}{\sum}} \frac{v_i}{1+j}\times C^7_j (1-p_i)^{7-j}p_i^j
$$

Therefore, if the equilibrium has m non-zero probabilities, then we have m-1 equalities. Our final equation comes from the fact that the probabilites have to sum to 1. Hence there exists a unique solution. 

(Technically, for this to be a mixed equilibrium, we also require that the expected payoff from playing mixed actions are better than actions that are not played.)


```julia
using SymPy

v,p,p1,p2,p3,p4,p5,p6,p7,p8 = symbols("v,p,p1 p2 p3 p4 p5 p6 p7 p8")

EU = sum( v/(1+j)* binomial(7,j)* (1-p)^(7-j)* p^j for j in 0:7);

EU
```




$\begin{equation*}\frac{p^{7} v}{8} + p^{6} v \left(1 - p\right) + \frac{7 p^{5} v \left(1 - p\right)^{2}}{2} + 7 p^{4} v \left(1 - p\right)^{3} + \frac{35 p^{3} v \left(1 - p\right)^{4}}{4} + 7 p^{2} v \left(1 - p\right)^{5} + \frac{7 p v \left(1 - p\right)^{6}}{2} + v \left(1 - p\right)^{7}\end{equation*}$




Notice how in SymPy the variable names are tied with functions. 

We can evaluate the function at different values/variables using =>. We can also define a new function this way.


```julia
EU1 = EU(p=>p1,v=>valuations[1])
EU1
```




$\begin{equation*}5 p_{1}^{7} + 40 p_{1}^{6} \left(1 - p_{1}\right) + 140 p_{1}^{5} \left(1 - p_{1}\right)^{2} + 280 p_{1}^{4} \left(1 - p_{1}\right)^{3} + 350 p_{1}^{3} \left(1 - p_{1}\right)^{4} + 280 p_{1}^{2} \left(1 - p_{1}\right)^{5} + 140 p_{1} \left(1 - p_{1}\right)^{6} + 40 \left(1 - p_{1}\right)^{7}\end{equation*}$





```julia
EU2 = EU(p=>p2,v=>valuations[2])
EU3 = EU(p=>p3,v=>valuations[3])
EU4 = EU(p=>p4,v=>valuations[4])
EU5 = EU(p=>p5,v=>valuations[5])
EU6 = EU(p=>p6,v=>valuations[6])
EU7 = EU(p=>p7,v=>valuations[7])
EU8 = EU(p=>p8,v=>valuations[8])

system_of_eq = [EU2-EU1, EU2-EU3, EU2-EU4, EU2-EU5, EU2-EU6, EU2-EU7, EU2-EU8, 1-p1-p2-p3-p4-p5-p6-p7-p8]

system_of_eq
```




$\left[ \begin{array}{r}- 5 p_{1}^{7} - 40 p_{1}^{6} \left(1 - p_{1}\right) - 140 p_{1}^{5} \left(1 - p_{1}\right)^{2} - 280 p_{1}^{4} \left(1 - p_{1}\right)^{3} - 350 p_{1}^{3} \left(1 - p_{1}\right)^{4} - 280 p_{1}^{2} \left(1 - p_{1}\right)^{5} - 140 p_{1} \left(1 - p_{1}\right)^{6} + \frac{25 p_{2}^{7}}{4} + 50 p_{2}^{6} \left(1 - p_{2}\right) + 175 p_{2}^{5} \left(1 - p_{2}\right)^{2} + 350 p_{2}^{4} \left(1 - p_{2}\right)^{3} + \frac{875 p_{2}^{3} \left(1 - p_{2}\right)^{4}}{2} + 350 p_{2}^{2} \left(1 - p_{2}\right)^{5} + 175 p_{2} \left(1 - p_{2}\right)^{6} - 40 \left(1 - p_{1}\right)^{7} + 50 \left(1 - p_{2}\right)^{7}\\\frac{25 p_{2}^{7}}{4} + 50 p_{2}^{6} \left(1 - p_{2}\right) + 175 p_{2}^{5} \left(1 - p_{2}\right)^{2} + 350 p_{2}^{4} \left(1 - p_{2}\right)^{3} + \frac{875 p_{2}^{3} \left(1 - p_{2}\right)^{4}}{2} + 350 p_{2}^{2} \left(1 - p_{2}\right)^{5} + 175 p_{2} \left(1 - p_{2}\right)^{6} - \frac{25 p_{3}^{7}}{4} - 50 p_{3}^{6} \left(1 - p_{3}\right) - 175 p_{3}^{5} \left(1 - p_{3}\right)^{2} - 350 p_{3}^{4} \left(1 - p_{3}\right)^{3} - \frac{875 p_{3}^{3} \left(1 - p_{3}\right)^{4}}{2} - 350 p_{3}^{2} \left(1 - p_{3}\right)^{5} - 175 p_{3} \left(1 - p_{3}\right)^{6} + 50 \left(1 - p_{2}\right)^{7} - 50 \left(1 - p_{3}\right)^{7}\\\frac{25 p_{2}^{7}}{4} + 50 p_{2}^{6} \left(1 - p_{2}\right) + 175 p_{2}^{5} \left(1 - p_{2}\right)^{2} + 350 p_{2}^{4} \left(1 - p_{2}\right)^{3} + \frac{875 p_{2}^{3} \left(1 - p_{2}\right)^{4}}{2} + 350 p_{2}^{2} \left(1 - p_{2}\right)^{5} + 175 p_{2} \left(1 - p_{2}\right)^{6} - \frac{25 p_{4}^{7}}{4} - 50 p_{4}^{6} \left(1 - p_{4}\right) - 175 p_{4}^{5} \left(1 - p_{4}\right)^{2} - 350 p_{4}^{4} \left(1 - p_{4}\right)^{3} - \frac{875 p_{4}^{3} \left(1 - p_{4}\right)^{4}}{2} - 350 p_{4}^{2} \left(1 - p_{4}\right)^{5} - 175 p_{4} \left(1 - p_{4}\right)^{6} + 50 \left(1 - p_{2}\right)^{7} - 50 \left(1 - p_{4}\right)^{7}\\\frac{25 p_{2}^{7}}{4} + 50 p_{2}^{6} \left(1 - p_{2}\right) + 175 p_{2}^{5} \left(1 - p_{2}\right)^{2} + 350 p_{2}^{4} \left(1 - p_{2}\right)^{3} + \frac{875 p_{2}^{3} \left(1 - p_{2}\right)^{4}}{2} + 350 p_{2}^{2} \left(1 - p_{2}\right)^{5} + 175 p_{2} \left(1 - p_{2}\right)^{6} - \frac{35 p_{5}^{7}}{4} - 70 p_{5}^{6} \left(1 - p_{5}\right) - 245 p_{5}^{5} \left(1 - p_{5}\right)^{2} - 490 p_{5}^{4} \left(1 - p_{5}\right)^{3} - \frac{1225 p_{5}^{3} \left(1 - p_{5}\right)^{4}}{2} - 490 p_{5}^{2} \left(1 - p_{5}\right)^{5} - 245 p_{5} \left(1 - p_{5}\right)^{6} + 50 \left(1 - p_{2}\right)^{7} - 70 \left(1 - p_{5}\right)^{7}\\\frac{25 p_{2}^{7}}{4} + 50 p_{2}^{6} \left(1 - p_{2}\right) + 175 p_{2}^{5} \left(1 - p_{2}\right)^{2} + 350 p_{2}^{4} \left(1 - p_{2}\right)^{3} + \frac{875 p_{2}^{3} \left(1 - p_{2}\right)^{4}}{2} + 350 p_{2}^{2} \left(1 - p_{2}\right)^{5} + 175 p_{2} \left(1 - p_{2}\right)^{6} - \frac{35 p_{6}^{7}}{4} - 70 p_{6}^{6} \left(1 - p_{6}\right) - 245 p_{6}^{5} \left(1 - p_{6}\right)^{2} - 490 p_{6}^{4} \left(1 - p_{6}\right)^{3} - \frac{1225 p_{6}^{3} \left(1 - p_{6}\right)^{4}}{2} - 490 p_{6}^{2} \left(1 - p_{6}\right)^{5} - 245 p_{6} \left(1 - p_{6}\right)^{6} + 50 \left(1 - p_{2}\right)^{7} - 70 \left(1 - p_{6}\right)^{7}\\\frac{25 p_{2}^{7}}{4} + 50 p_{2}^{6} \left(1 - p_{2}\right) + 175 p_{2}^{5} \left(1 - p_{2}\right)^{2} + 350 p_{2}^{4} \left(1 - p_{2}\right)^{3} + \frac{875 p_{2}^{3} \left(1 - p_{2}\right)^{4}}{2} + 350 p_{2}^{2} \left(1 - p_{2}\right)^{5} + 175 p_{2} \left(1 - p_{2}\right)^{6} - \frac{35 p_{7}^{7}}{4} - 70 p_{7}^{6} \left(1 - p_{7}\right) - 245 p_{7}^{5} \left(1 - p_{7}\right)^{2} - 490 p_{7}^{4} \left(1 - p_{7}\right)^{3} - \frac{1225 p_{7}^{3} \left(1 - p_{7}\right)^{4}}{2} - 490 p_{7}^{2} \left(1 - p_{7}\right)^{5} - 245 p_{7} \left(1 - p_{7}\right)^{6} + 50 \left(1 - p_{2}\right)^{7} - 70 \left(1 - p_{7}\right)^{7}\\\frac{25 p_{2}^{7}}{4} + 50 p_{2}^{6} \left(1 - p_{2}\right) + 175 p_{2}^{5} \left(1 - p_{2}\right)^{2} + 350 p_{2}^{4} \left(1 - p_{2}\right)^{3} + \frac{875 p_{2}^{3} \left(1 - p_{2}\right)^{4}}{2} + 350 p_{2}^{2} \left(1 - p_{2}\right)^{5} + 175 p_{2} \left(1 - p_{2}\right)^{6} - \frac{25 p_{8}^{7}}{2} - 100 p_{8}^{6} \left(1 - p_{8}\right) - 350 p_{8}^{5} \left(1 - p_{8}\right)^{2} - 700 p_{8}^{4} \left(1 - p_{8}\right)^{3} - 875 p_{8}^{3} \left(1 - p_{8}\right)^{4} - 700 p_{8}^{2} \left(1 - p_{8}\right)^{5} - 350 p_{8} \left(1 - p_{8}\right)^{6} + 50 \left(1 - p_{2}\right)^{7} - 100 \left(1 - p_{8}\right)^{7}\\- p_{1} - p_{2} - p_{3} - p_{4} - p_{5} - p_{6} - p_{7} - p_{8} + 1\end{array} \right]$





```julia
mixed_ne = sympy.nsolve( system_of_eq , (p1,p2,p3,p4,p5,p6,p7,p8) , (1,1,1,1,1,1,1,1) )
```




$\left[ \begin{array}{r}0.0005900395709166406840137419570967888870576042685833651064880671073501502121675\\0.06553637288318371000730577629947820555476889335433512302866178831665458158983\\0.06553637288318371000730577629947820555476889335433512302866178831665458158983\\0.06553637288318371000730577629947820555476889335433512302866178831665458158983\\0.1696779680897165722324230717817383914913708924507470226023318700539585865871\\0.1696779680897165722324230717817383914913708924507470226023318700539585865871\\0.1696779680897165722324230717817383914913708924507470226023318700539585865871\\0.293766937510382512596799713799253419974523038316170198000530957780810345257\end{array}\right]$




Let's look up the syntax for the function above.


```julia
?sympy.nsolve
```




    
        Solve a nonlinear equation system numerically: ``nsolve(f, [args,] x0,
        modules=['mpmath'], **kwargs)``.
    
        Explanation
        ===========
    
        ``f`` is a vector function of symbolic expressions representing the system.
        *args* are the variables. If there is only one variable, this argument can
        be omitted. ``x0`` is a starting vector close to a solution.
    
        Use the modules keyword to specify which modules should be used to
        evaluate the function and the Jacobian matrix. Make sure to use a module
        that supports matrices. For more information on the syntax, please see the
        docstring of ``lambdify``.
    
        If the keyword arguments contain ``dict=True`` (default is False) ``nsolve``
        will return a list (perhaps empty) of solution mappings. This might be
        especially useful if you want to use ``nsolve`` as a fallback to solve since
        using the dict argument for both methods produces return values of
        consistent type structure. Please note: to keep this consistent with
        ``solve``, the solution will be returned in a list even though ``nsolve``
        (currently at least) only finds one solution at a time.
    
        Overdetermined systems are supported.
    
        Examples
        ========
    
        >>> from sympy import Symbol, nsolve
        >>> import mpmath
        >>> mpmath.mp.dps = 15
        >>> x1 = Symbol('x1')
        >>> x2 = Symbol('x2')
        >>> f1 = 3 * x1**2 - 2 * x2**2 - 1
        >>> f2 = x1**2 - 2 * x1 + x2**2 + 2 * x2 - 8
        >>> print(nsolve((f1, f2), (x1, x2), (-1, 1)))
        Matrix([[-1.19287309935246], [1.27844411169911]])
    
        For one-dimensional functions the syntax is simplified:
    
        >>> from sympy import sin, nsolve
        >>> from sympy.abc import x
        >>> nsolve(sin(x), x, 2)
        3.14159265358979
        >>> nsolve(sin(x), 2)
        3.14159265358979
    
        To solve with higher precision than the default, use the prec argument:
    
        >>> from sympy import cos
        >>> nsolve(cos(x) - x, 1)
        0.739085133215161
        >>> nsolve(cos(x) - x, 1, prec=50)
        0.73908513321516064165531208767387340401341175890076
        >>> cos(_)
        0.73908513321516064165531208767387340401341175890076
    
        To solve for complex roots of real functions, a nonreal initial point
        must be specified:
    
        >>> from sympy import I
        >>> nsolve(x**2 + 2, I)
        1.4142135623731*I
    
        ``mpmath.findroot`` is used and you can find their more extensive
        documentation, especially concerning keyword parameters and
        available solvers. Note, however, that functions which are very
        steep near the root, the verification of the solution may fail. In
        this case you should use the flag ``verify=False`` and
        independently verify the solution.
    
        >>> from sympy import cos, cosh
        >>> f = cos(x)*cosh(x) - 1
        >>> nsolve(f, 3.14*100)
        Traceback (most recent call last):
        ...
        ValueError: Could not find root within given tolerance. (1.39267e+230 > 2.1684e-19)
        >>> ans = nsolve(f, 3.14*100, verify=False); ans
        312.588469032184
        >>> f.subs(x, ans).n(2)
        2.1e+121
        >>> (f/f.diff(x)).subs(x, ans).n(2)
        7.4e-15
    
        One might safely skip the verification if bounds of the root are known
        and a bisection method is used:
    
        >>> bounds = lambda i: (3.14*i, 3.14*(i + 1))
        >>> nsolve(f, bounds(100), solver='bisect', verify=False)
        315.730061685774
    
        Alternatively, a function may be better behaved when the
        denominator is ignored. Since this is not always the case, however,
        the decision of what function to use is left to the discretion of
        the user.
    
        >>> eq = x**2/(1 - x)/(1 - 2*x)**2 - 100
        >>> nsolve(eq, 0.46)
        Traceback (most recent call last):
        ...
        ValueError: Could not find root within given tolerance. (10000 > 2.1684e-19)
        Try another starting point or tweak arguments.
        >>> nsolve(eq.as_numer_denom()[0], 0.46)
        0.46792545969349058
    
        



Suppose we want to see if there is an equilibrium where players only mix between room 2-8. 


```julia
solution = sympy.nsolve((EU2-EU3, EU2-EU4, EU2-EU5, EU2-EU6,EU2-EU7,EU2-EU8,1-p2-p3-p4-p5-p6-p7-p8),(p2,p3,p4,p5,p6,p7,p8),(1,1,1,1,1,1,1))

@show EU2(solution[1])
@show EU3(solution[2])
@show EU4(solution[3])
@show EU5(solution[4])
@show EU6(solution[5])
@show EU7(solution[6])
@show EU8(solution[7])

EU1(0)
```

    EU2(solution[1]) = 39.90698417785192965555710132508902709270987682969983827549323312327389406646
    EU3(solution[2]) = 39.90698417785192965555710132508902709270987682969983827549323312327389406646
    EU4(solution[3]) = 39.90698417785192965555710132508902709270987682969983827549323312327389406646
    EU5(solution[4]) = 39.90698417785192965555710132508902709270987682969983827549323312327389406646
    EU6(solution[5]) = 39.90698417785192965555710132508902709270987682969983827549323312327389406646
    EU7(solution[6]) = 39.90698417785192965555710132508902709270987682969983827549323312327389406646
    EU8(solution[7]) = 39.90698417785192965555710132508902709270987682969983827549323312327389406646





$\begin{equation*}40\end{equation*}$




The expected payoff by deviating to room 1 is 40, which is higher than any actions we are mixing over. So this is not an equilibrium.

### Another way to do part 1

We first write a function to check if a vector of pure strategies is an equilibrium. Then, suppose we have a set of potential candidates, we go through the set to find the first one that works.


```julia
function is_eqm(s)
    check = zeros(8)
    for i = 1:8
        s_i = s[i]
        s_noti = s[1:end .!=i]
        check[i] = is_BR(s_i,s_noti)[1]
    end
    return (minimum(check)==1)
end

@show is_eqm([1,2,3,4,5,6,7,8])
is_eqm(pure_ne)
```

    is_eqm([1, 2, 3, 4, 5, 6, 7, 8]) = false





    true




```julia
function brute_search(S)
    found = false
    for i in eachindex(S)  
                cs = S[i]       
                check = is_eqm(cs)
                
                if check == true
                    return cs
                end

    end  
    println("Did not find any eqm.")
end
```




    brute_search (generic function with 1 method)



We can ues this to search over alll possible combinations of strategies. We can set

        full_S = [[x1,x2,x3,x4,x5,x6,x7,x8] for x1 in 1:8 for x2 in 1:8 for x3 in 1:8 for x4 in 1:8 for x5 in 1:8 for x6 in 1:8 for x7 in 1:8 for x8 in 1:8]

Then run

        search(full_S)

But this will take a very long time. One of the reaseon is it takes up a lot of memory. Also, we know the equilibrium does not contain any 1's. So that means going through the first 8^7 candidates is totally useless. Instead, what if we just generate a bunch of 8x1 vectors with values from 1 to 8, and search through that set? 


```julia
function gen_candidates(n)
    candidates = [rand(1:8,8)]
    count = 1
    for i=2:n
        candidates = vcat(candidates,[rand(1:8,8)])
    end
    return candidates
end

n = 1000

candidates = gen_candidates(n);

brute_search(candidates)
```

    Did not find any eqm.


The probability of successfully finding a NE with n


```julia
p = (factorial(8)/2)/(8^8)
(1-p)^n
```




    0.3004865831455566



231142578d8c308c6403828562f4c650f9fdf523
ipynb
Jupyter Notebook
Julia Intro, Running Man Room Game.ipynb
cytthomas/ECON514_2022_Tutorials
ce9877aad10d1d0b98c07eb21f56631fa12ae321
MIT
2022-01-24T17:49:01.000Z
2022-01-24T17:49:01.000Z
Julia Intro, Running Man Room Game.ipynb
ACrawley15/ECON514_2022_Tutorials
ce9877aad10d1d0b98c07eb21f56631fa12ae321
MIT
Julia Intro, Running Man Room Game.ipynb
ACrawley15/ECON514_2022_Tutorials
ce9877aad10d1d0b98c07eb21f56631fa12ae321
MIT
2022-01-19T16:32:58.000Z
2022-01-19T17:29:03.000Z
Qwen/Qwen-72B
1. YES
2. YES