<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Introduction to Computational General Equilibrium (CGE) Models
In this notebook we'll cover some basic ideas that will introduce you to CGE modelling, i.e., of solving Walrasian GE models in a computer. To make things simple we will focus on the easiest CGE model possible
Step1:  Let's compute the utility for "all" values of $x$ in $[0,10]^M$ and plot them.
Step2:  Demands
Since we computed the actual demand functions $x^_i=x^(p, \omega_i)$ for each good $j$ for each individual $i$, in this case we can simply compute the optimal demand for each individual $i$ given her initial endowments $\omega_i$ and prices $p$.
Step3:  This is an equilibrium, but of course in an economy with only 1 individual, $(p^, x^)=(p^,\omega)$ is an equilibrium for any $p^$.
Step4:  Here we have excess supply for good 1 and excess demand for good 2
Step5:  Again we have excess supply for good 1 and excess demand for good 2
Step6:  We have excess supply for goods 1 & 3 and excess demand for goods 2 & 4
Step7:  Again we have excess supply for good 1 and excess demand for good 2
Step8:  Now we have excess supply for goods 1, 7, 8, and 10, and excess demand for the other
How can we find the equilibrium price?
We need excess demand to be equal to zero!
So, we only need to use our computer to find the zero of an equation. Luckily, this is a very well studied problem (see this Wikipedia entry). In python we can use the scipy package.
Step9:  But before we use this to try to find an equilibrium, we need to remember a few things.
Equlibrium prices only determine relative prices, not absolute ones, so we need to decide on a numeraire. Various choices, e.g. chose $p_j=1$ for some good $j$. Another common option is to normalize the price vector so it has length 1.
Walras Law
Step10:  Simple $2\times2$ Economy
Step11:  Notice that the excess demand with normalized prices generates the same excess demand as with the original price vector
Step12:  Similar $2\times 4 $ Economy
Step13:  $6\times 2$ Economy
Let's try with a larger economy
Step14:  And even larger economy $10\times5$
<ASSISTANT_TASK:>
Python Code:
import numpy as np
from scipy.misc import derivative
import itertools
%matplotlib widget
import matplotlib.pyplot as plt
def CobbDouglas(x, alpha, h=1e-10, deriv=False):
    '''
    Compute the utility of an individual with Cobb-Douglas preferences
    Additionally it returns the exact and numerical derivatives
    x: (1 x N) vector
    alpha:(1 x N) vector
    '''
    lx = len(x)
    try:
        u = (x ** alpha).prod(axis=1)
    except:
        u = (x ** alpha).prod()
    if deriv:
        uprime = alpha * u / x
        nuprime = (((x + (np.zeros([lx, lx]) + h * np.eye(lx)))**alpha).prod(axis=1) - ((x + (np.zeros([lx, lx]) - h * np.eye(lx)))**alpha).prod(axis=1)) / (2 * h)
        du = [uprime, nuprime]
        return [u, du]
    else:
        return u
Explanation: Introduction to Computational General Equilibrium (CGE) Models
In this notebook we'll cover some basic ideas that will introduce you to CGE modelling, i.e., of solving Walrasian GE models in a computer. To make things simple we will focus on the easiest CGE model possible: a pure endowment economy with $M$ individuals and $N$ commodities. When $M=2$ and $N=2$, this is an economy that can be represented in an Edgeworth box.
Preferences
We will focus on an economy in which all individuals have Cobb-Douglas preferences, i.e., we'll assume that each individual $i$ has utility
$$
u_i(x_{i1}, x_{i2}, \dots, x_{iN})= x_{i1}^{\alpha_{i1}} x_{i2}^{\alpha_{i2}} \cdots x_{iN}^{\alpha_{iN}}
$$
where $x_{ij}$ denoted the amount of good $j=1,\ldots,M$ she consumes. We assume that $\alpha_{ij}\ge0$ for all $j=1,2,\ldots,M$, $\sum_{j=1}^M\alpha_{ij}\le1$ to ensure the utility function is strictly quasi-concave. Let's denote the consumption plan of individual $i$ by $x_i\equiv(x_{i1}, x_{i2}, \ldots, x_{iM})$. 
Given a price vector $p=(p_1, p_2, \ldots,p_M)$, the value of her total expenditures is
$$
p\cdot x_i = \sum_{j=1}^M p_i\cdot x_{ij}.
$$
Endowments
Each individual in the economy is endowed with a quantity of each good $\omega_{ij}\ge0$. We will denote by $\omega_i=(\omega_{i1}, \omega_{i2}, \ldots, \omega_{iM})$ the vector of endowments of an individual $i$.
Given her endowments vector and a price vector $p=(p_1, p_2, \ldots,p_M)$, the value of her resources is
$$
y_i \equiv p\cdot \omega_i = \sum_{j=1}^M p_i\cdot\omega_{ij}.
$$
Individual's Optimal Consumption/Demand
Given her initial endowments $\omega_i$ and a price vector $p$, each individual $i$ choses her consumption basket $x^*_i$, which solves the following optimization problem
$$
\max_{x_i} u_i(x_{i1}, x_{i2}, \dots, x_{iN})= x_{i1}^{\alpha_{i1}} x_{i2}^{\alpha_{i2}} \cdots x_{iN}^{\alpha_{iN}}
$$
subject to
$$
p\cdot x_i\le y_i.
$$
Since the utility function is strictly quasiconcave and the restriction is convex, we can apply the Kuhn-Tucker or Arrow-Enthoven Theorems to find the optimal choices of the individual. But, since her utility is strictly increasing in all goods, we know the restriction will hold with equality and we can just use Langrange's method (if you need to review this material, check out La Fuente Ch.7 or Simon and Blume Ch.18). The Langrangian is 
$$
\mathcal{L}(x_i, \lambda) = u_i(x_i)-\lambda(p\cdot x_i-y_i)
$$
so, her optimal consumption choice $x^*_i$ must satisfy
$$
\frac{\partial u_i(x^*i)}{\partial x{ij}}=\lambda p_j, \quad j=1,\ldots,M
$$
and the budget constraint. Given our assumption about the functional form of $u_i(x_i)$, this implies it must satisfy
$$
\frac{\alpha_{ij}}{x^_{ij}}u_i(x^_i)=\lambda p_j, \quad j=1,\ldots,M,
$$
which implies that for any two goods $j,j'=1,\ldots,M$
$$
\frac{\alpha_{ij}}{p_j \cdot x^{ij}}=\frac{\alpha{ij'}}{p_{j'} \cdot x^_{ij'}}.
$$
Replacing in the budget contraint we have
$$
p\cdot x_i = \sum_{j=1}^M p_j \cdot x^{ij} = \sum{j=1}^M \left(\frac{\alpha_{ij}}{\alpha_{i1}}p_1 x^{i1}\right) = p_1 x^*{i1}\frac{\sum_{j=1}^M \alpha_{ij}}{\alpha_{i1}}=y_i,
$$
so that the optimal demand for good $j$ satisfies
$$
p_j x^{ij}=\frac{\alpha{ij}}{\sum_{j=1}^M \alpha_{ij}}y_i \iff x^{ij}=\frac{1}{p_j}\frac{\alpha{ij}}{\sum_{j=1}^M \alpha_{ij}}y_i.
$$
Notice that since our utility function is homothetic, we knew (right?) that the share of income spent on good $j$ had to be equal to
$$
\frac{\alpha_{ij}}{\sum_{j=1}^M \alpha_{ij}}.
$$
Equilibrium
A Walrasian equilibrium for this endowment economy is a set of prices $p^$ and an allocation $x^=(x^_1, \ldots, x^_N)$ such that
$x^_i$ maximizes the utility of individual $i$ given the value of her endowmwents (at those prices $p^$),
$x^*$ is feasible, i.e., 
$$
\sum_{i=1}^N x^*{ij}\le \sum{i=1}^N \omega_{ij}, \quad \text{for all } j=1,\ldots,M.
$$
all markets clear, i.e., 
$$
\sum_{i=1}^N x^*{ij}=\sum{i=1}^N \omega_{ij}, \quad \text{for all } j=1,\ldots,M.
$$
Moreover, it must be the case that $p_j\ge0$ and
$$
p_j\left(\sum_{i=1}^N x^*{ij}-\sum{i=1}^N \omega_{ij}\right)=0, \quad  \text{for all } j=1,\ldots,M.
$$
Computation
Given our assumptions on the utility functions of individuals, we know there exists a unique Walrasian equilibrium for our economy. This is important, since if we want to compute an equilibrium, we need to know it exists (otherwise we would be trying the impossible!) and, additionally, that what we find is indeed the only solution (if there are multiple isolated equilibria we could try to find them, but this is only possible if it is a finite set).
Let us create some functions that will allow us to solve this problem.
End of explanation
N = 1 # Individuals
M = 2 # Goods
alpha = np.array([0.5, 0.5])
x = np.linspace(0, 10, num=100)
X = np.array([i for i in itertools.product(x, repeat=M)])
u = CobbDouglas(X, alpha)
u
fig, ax = plt.subplots(figsize=(15,10))
plt.contourf(x, x, u.reshape(len(x), len(x)), 30)
plt.xlabel('Good 1', size=18)
plt.ylabel('Good 2', size=18)
plt.title('Indiference Map', size=22)
from mpl_toolkits.mplot3d import Axes3D
xx, yy = np.meshgrid(x, x)
fig = plt.figure(figsize=(15, 10))
ax = fig.add_subplot(111, projection='3d')
surf = ax.plot_surface(xx, yy, u.reshape(len(x), len(x)), rstride=1, cstride=1, cmap='viridis')
ax.set_xlabel('Good 1', size=18)
ax.set_ylabel('Good 2', size=18)
ax.set_zlabel('Utility', size=18)
plt.show()
Explanation: Let's compute the utility for "all" values of $x$ in $[0,10]^M$ and plot them.
End of explanation
def CobbDouglasDemand(p, w, alpha):
    '''
    Compute the demand vector of an individual with Cobb-Douglas Utility with parameter vector alpha,
    given her initial endowments and prices for all goods
    '''
    # Total income
    y = w.dot(p)
    ly = len(y)
    x = ((alpha / alpha.sum(axis=1, keepdims=True) * y).T / p).T
    return x
def ExcessDemand(p, w, alpha):
    '''
    Compute excess demand function for each good
    '''
    z = CobbDouglasDemand(p, w, alpha).sum(axis=0) - w.sum(axis=0)
    return z
# Let's test to make sure function is well written
alpha = np.array([[0.5, 0.5]])
w = np.array([[1, 2]])
p = np.array([[2, 1]]).T
CobbDouglasDemand(p, w, alpha)
ExcessDemand(p, w, alpha)
Explanation: Demands
Since we computed the actual demand functions $x^_i=x^(p, \omega_i)$ for each good $j$ for each individual $i$, in this case we can simply compute the optimal demand for each individual $i$ given her initial endowments $\omega_i$ and prices $p$.
End of explanation
alpha = np.array([[0.5, 0.5], [0.25, 0.75]])
w = np.array([[1, 2], [2, 2]])
p = np.array([[2, 1]]).T
CobbDouglasDemand(p, w, alpha)
ExcessDemand(p, w, alpha)
Explanation: This is an equilibrium, but of course in an economy with only 1 individual, $(p^, x^)=(p^,\omega)$ is an equilibrium for any $p^$.
End of explanation
alpha = np.array([[0.5, 0.5], [0.25, 0.75], [0.5, 0.5]])
w = np.array([[1, 2], [2, 2], [1, 2]])
p = np.array([[2, 1]]).T
CobbDouglasDemand(p, w, alpha)
ExcessDemand(p, w, alpha)
Explanation: Here we have excess supply for good 1 and excess demand for good 2
End of explanation
# Economy with 2 individuals of each type and 4 goods 2 groups of identical goods
alpha = np.array([[0.5, 0.5] * 2, [0.25, 0.75] * 2])
w = np.array([[1, 2] * 2, [2, 2] * 2] )
p = np.array([[2, 1] * 2]).T
CobbDouglasDemand(p, w, alpha)
ExcessDemand(p, w, alpha)
Explanation: Again we have excess supply for good 1 and excess demand for good 2
End of explanation
# Same economy with 6 individuals of each type
alpha = np.array([[0.5, 0.5], [0.25, 0.75]] * 6)
w = np.array([[1, 2], [2, 2]] * 6)
p = np.array([[2, 1]]).T
CobbDouglasDemand(p, w, alpha)
ExcessDemand(p, w, alpha)
Explanation: We have excess supply for goods 1 & 3 and excess demand for goods 2 & 4
End of explanation
# Generate a random economy
np.random.seed(123456)
N = 5 # Individuals
M = 10 # Goods
alpha = np.random.uniform(0, 1, size=(N, M))
w = np.random.uniform(0, 1, size=(N, M))
p = np.random.uniform(0, 1, size=(1, M)).T
CobbDouglasDemand(p, w, alpha)
ExcessDemand(p, w, alpha)
Explanation: Again we have excess supply for good 1 and excess demand for good 2
End of explanation
from scipy.optimize import root
Explanation: Now we have excess supply for goods 1, 7, 8, and 10, and excess demand for the other
How can we find the equilibrium price?
We need excess demand to be equal to zero!
So, we only need to use our computer to find the zero of an equation. Luckily, this is a very well studied problem (see this Wikipedia entry). In python we can use the scipy package.
End of explanation
def ExcessDemand2(p, w, alpha, normalize=False):
    '''
    Compute excess demand function for all but the last good
    Price of good 1 = 1
    Prices are normalized to length of vector is equal to 1
    '''
    # Ensure p has the right shape
    p = p.reshape((np.max(p.shape), 1))
    p = np.vstack([[1], p])
    if normalize:
        p = (p / np.linalg.norm(p))
    z = CobbDouglasDemand(p, w, alpha).sum(axis=0) - w.sum(axis=0)
    return z[1:]
Explanation: But before we use this to try to find an equilibrium, we need to remember a few things.
Equlibrium prices only determine relative prices, not absolute ones, so we need to decide on a numeraire. Various choices, e.g. chose $p_j=1$ for some good $j$. Another common option is to normalize the price vector so it has length 1.
Walras Law: If $M-1$ markets clear, then the $M^{th}$ market also clears.
So let's define a new excess demand function that takes these issues into account.
End of explanation
alpha = np.array([[0.5, 0.5], [0.25, 0.75]])
w = np.array([[1, 2], [2, 2]])
p = np.array([[1, 0.5]]).T
Explanation: Simple $2\times2$ Economy
End of explanation
CobbDouglasDemand(p, w, alpha)
ExcessDemand(p, w, alpha)
ExcessDemand2(p[1:], w, alpha)
sol = root(ExcessDemand2, np.array([1]), args=(w, alpha))
sol
pstar = np.vstack([[1], sol.x])
ExcessDemand(pstar, w, alpha)
Explanation: Notice that the excess demand with normalized prices generates the same excess demand as with the original price vector
End of explanation
# Economy with 2 individuals of each type and 4 goods 2 groups of identical goods
alpha = np.array([[0.5, 0.5] * 2, [0.25, 0.75] * 2])
w = np.array([[1, 2] * 2, [2, 2] * 2] )
p = np.array([[2, 1] * 2]).T
p = p / p[0]
ExcessDemand(p, w, alpha)
ExcessDemand2(p[1:], w, alpha)
sol = root(ExcessDemand2, p[1:], args=(w, alpha))
sol
pstar = sol.x
pstar = pstar.reshape((np.max(pstar.shape), 1))
pstar = np.vstack([[1], pstar])
ExcessDemand(pstar, w, alpha)
Explanation: Similar $2\times 4 $ Economy
End of explanation
# Same economy with 6 individuals of each type
alpha = np.array([[0.5, 0.5], [0.25, 0.75]] * 6)
w = np.array([[1, 2], [2, 2]] * 6)
p = np.array([[2, 1]]).T
sol = root(ExcessDemand2, p[1:], args=(w, alpha))
sol
pstar = sol.x
pstar = pstar.reshape((np.max(pstar.shape), 1))
pstar = np.vstack([[1], pstar])
ExcessDemand(pstar, w, alpha)
Explanation: $6\times 2$ Economy
Let's try with a larger economy
End of explanation
# Generate a random economy
np.random.seed(123456)
N = 5 # Individuals
M = 10 # Goods
alpha = np.random.uniform(0, 1/M, size=(N, M))
w = np.random.uniform(0, 1/M, size=(N, M))
p = np.random.uniform(0, 1/M, size=(1, M)).T
alpha.sum(axis=1)
ExcessDemand2(np.ones((1, M-1)).T, w, alpha)
sol = root(ExcessDemand2, np.ones((1, M-1)).T, args=(w, alpha))
sol
pstar = sol.x
pstar = pstar.reshape((np.max(pstar.shape), 1))
pstar = np.vstack([[1], pstar])
ExcessDemand(pstar, w, alpha)
Explanation: And even larger economy $10\times5$
End of explanation
<END_TASK>