```python
from __future__ import division
import numpy as np
import matplotlib.pyplot as plt
import scipy
from scipy import optimize
import sympy as sm
sm.init_printing(use_unicode=True)
```

# 1. Human capital accumulation

Consider a worker living in **two periods**, $t \in \{1,2\}$. 

In each period she decides whether to **work ($l_t = 1$) or not ($l_t = 0$)**. 

She can *not* borrow or save and thus **consumes all of her income** in each period. 

If she **works** her **consumption** becomes:

$$c_t = w h_t l_t\,\,\text{if}\,\,l_t=1$$

where $w$ is **the wage rate** and $h_t$ is her **human capital**. 

If she does **not work** her consumption becomes:

$$c_t = b\,\,\text{if}\,\,l_t=0$$

where $b$ is the **unemployment benefits**.

Her **utility of consumption** is: 

$$ \frac{c_t^{1-\rho}}{1-\rho} $$

Her **disutility of working** is:

$$ \gamma l_t $$

From period 1 to period 2, she **accumulates human capital** according to:

$$ h_2 = h_1 + l_1 + 
\begin{cases}
0 & \text{with prob. }0.5 \\
\Delta & \text{with prob. }0.5 
\end{cases} \\
$$

where $\Delta$ is a **stochastic experience gain**.

In the **second period** the worker thus solves:

$$
\begin{eqnarray*}
v_{2}(h_{2}) & = &\max_{l_{2}} \frac{c_2^{1-\rho}}{1-\rho} - \gamma l_2
\\ & \text{s.t.} & \\
c_{2}& = & w h_2 l_2 \\
l_{2}& \in &\{0,1\}
\end{eqnarray*}
$$

In the **first period** the worker thus solves:

$$
\begin{eqnarray*}
v_{1}(h_{1}) &=& \max_{l_{1}} \frac{c_1^{1-\rho}}{1-\rho} - \gamma l_1 + \beta\mathbb{E}_{1}\left[v_2(h_2)\right]
\\ & \text{s.t.} & \\
c_1 &=& w h_1 l_1 \\
h_2 &=& h_1 + l_1 + \begin{cases}
0 & \text{with prob. }0.5\\
\Delta & \text{with prob. }0.5 
\end{cases}\\
l_{1} &\in& \{0,1\}\\
\end{eqnarray*}
$$

where $\beta$ is the **discount factor** and $\mathbb{E}_{1}\left[v_2(h_2)\right]$ is the **expected value of living in period two**.

The **parameters** of the model are:


```python
rho = 2.0
beta = 0.96
gamma = 0.1
w = 2.0
b = 1.0
delta = 0.1
```

The **relevant levels of human capital** are:


```python
h_vec = np.linspace(0.1,1.5,100) 
```

**Question 1:** Solve the model in period 2 and illustrate the solution (including labor supply as a function of human capital). 


```python
c2 = sm.symbols("c_2")
l2 = sm.symbols("l_2")
c1 = sm.symbols("c_1")
l1 = sm.symbols("l_1")
h2 = sm.symbols("h_2")
h1 = sm.symbols("h_1")
l = sm.symbols("l")
rho_2 = sm.symbols("rho")
gamma_2 = sm.symbols("gamma")
w_2 = sm.symbols("w")
```


```python
#Defining relevant equations
def util(h, l):
    return consumption(h, l)**(1-rho)/(1-rho)-gamma*l

def consumption(h, l):
    if l == 1: 
        return w*h
    elif l == 0:
        return b

def v2util(h, l):
    return util(h, l)


#Solve cutoff value between work and not working
cutoff_p2 = sm.solve(v2util(h2, 1)-v2util(h2, 0), h2)
print(f'Cutoff value is: {cutoff_p2}. The individual works if their level of human capital is above this value as utility is increasing in human capital.')
```

    Cutoff value is: [0.555555555555556]. The individual works if their level of human capital is above this value as utility is increasing in human capital.



```python
#Plotting the cutoff
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.plot(h_vec, v2util(h_vec, 1), "r-", label='Working')
ax.axhline(y=v2util(h_vec, 0), xmin=0, xmax=2, label='Not working')
plt.title("Utility as a function of human capital")
plt.xlabel('Human capital')
plt.ylabel('Utility')
plt.legend()
plt.show()
```

**Question 2:** Solve the model in period 1 and illustrate the solution (including labor supply as a function of human capital). 


```python
#Defining the relevant equations
def utilp1(h, l):
    return consumptionp1(h, l)**(1-rho)/(1-rho)-gamma*l+beta*v2util(human_cap(h, l), l)

def human_cap(h, l):
    return h + l + 0.5*delta

def consumptionp1(h, l):
    if l == 1: 
        return w*h
    elif l == 0:
        return b
#Calculating the cutoff value
cutoff_p1 = sm.solve(utilp1(h1, 1)-utilp1(h1, 0), h1)
print(f'The solution for period 1 is: {cutoff_p1}. As the relevant levels of human capital are positive, we are left with the solution: {cutoff_p1[1]:.2f}')
```

    The solution for period 1 is: [-0.846169527733768, 0.351725083289324]. As the relevant levels of human capital are positive, we are left with the solution: 0.35



```python
#Plotting the cutoff
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.plot(h_vec, utilp1(h_vec, 1), "r-", label='Working')
ax.axhline(y=utilp1(h_vec, 0), xmin=0, xmax=2, label='Not working')
plt.title("Utility as a function of human capital in period 1")
plt.xlabel('Human capital')
plt.ylabel('Utility')
plt.legend()
plt.show()
```

**Question 3:** Will the worker never work if her potential wage income is lower than the unemployment benefits she can get? Explain and illustrate why or why not.


```python
#Let's try setting b > w
b = 2.5
#Plotting the cutoff
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.plot(h_vec, utilp1(h_vec, 1), "r-", label='Working')
ax.axhline(y=utilp1(h_vec, 0), xmin=0, xmax=2, label='Not working')
plt.title("Utility as a function of human capital in period 1")
plt.xlabel('Human capital')
plt.ylabel('Utility')
plt.legend()
plt.show()
```

As we see, the individual can choose to work in period 1, if the level of human capital is large enough. The cutoff is within the 'relevant' level of acquired human capital. Therefore the worker may choose to work, even though the unemployment benefits are larger than the wage. 

# 2. AS-AD model

Consider the following **AS-AD model**. The **goods market equilibrium** is given by

$$ y_{t} = -\alpha r_{t} + v_{t} $$

where $y_{t}$ is the **output gap**, $r_{t}$ is the **ex ante real interest** and $v_{t}$ is a **demand disturbance**. 

The central bank's **Taylor rule** is

$$ i_{t} = \pi_{t+1}^{e} + h \pi_{t} + b y_{t}$$

where $i_{t}$ is the **nominal interest rate**, $\pi_{t}$ is the **inflation gap**, and $\pi_{t+1}^{e}$ is the **expected inflation gap**. 

The **ex ante real interest rate** is given by 

$$ r_{t} = i_{t} - \pi_{t+1}^{e} $$

Together, the above implies that the **AD-curve** is

$$ \pi_{t} = \frac{1}{h\alpha}\left[v_{t} - (1+b\alpha)y_{t}\right]$$

Further, assume that the **short-run supply curve (SRAS)** is given by

$$ \pi_{t} = \pi_{t}^{e} + \gamma y_{t} + s_{t}$$

where $s_t$ is a **supply disturbance**.

**Inflation expectations are adaptive** and given by

$$ \pi_{t}^{e} = \phi\pi_{t-1}^{e} + (1-\phi)\pi_{t-1}$$

Together, this implies that the **SRAS-curve** can also be written as

$$ \pi_{t} = \pi_{t-1} + \gamma y_{t} - \phi\gamma y_{t-1} + s_{t} - \phi s_{t-1} $$

The **parameters** of the model are:


```python
par = {}

par['alpha'] = 5.76
par['h'] = 0.5
par['b'] = 0.5
par['phi'] = 0
par['gamma'] = 0.075
```

**Question 1:** Use the ``sympy`` module to solve for the equilibrium values of output, $y_t$, and inflation, $\pi_t$, (where AD = SRAS) given the parameters ($\alpha$, $h$, $b$, $\alpha$, $\gamma$) and $y_{t-1}$ , $\pi_{t-1}$, $v_t$, $s_t$, and $s_{t-1}$.


```python
pt = sm.symbols("pi_t")
ptminus = sm.symbols("pi_t-1")
yt = sm.symbols("y_t")
ytminus = sm.symbols("y_t-1")
st = sm.symbols("s_t")
stminus = sm.symbols("s_t-1")
vt = sm.symbols("v_t")
alpha = sm.symbols("alpha")
gamma = sm.symbols("gamma")
phi = sm.symbols("phi")
h = sm.symbols("h")
b = sm.symbols("b")
SRAS = sm.Eq(pt,ptminus+gamma*yt-phi*gamma*ytminus+st-phi*stminus)
SRAS

```


```python
AD = sm.Eq(pt,(1/(h*alpha)*(vt-(1+b*alpha)*yt)))
AD
```


```python
#Came back to define the function as it was needed below.
def SRASandAD(h, alpha, b, phi, gamma):
    return sm.Eq((1/(h*alpha)*(vt-(1+b*alpha)*yt)),ptminus+gamma*yt-phi*gamma*ytminus+st-phi*stminus)
SRASandAD(h, alpha, b, phi, gamma)
```


```python
#solution for y
solutiony21 = sm.solve(SRASandAD(h, alpha, b, phi, gamma), yt)
#came back to store this for 2.3
yt21 = sm.lambdify((h, alpha, vt, b, st, phi, stminus, ptminus, ytminus, gamma), solutiony21[0])
solutiony21
```


```python
#Solution for pi
solutionpi21 = (1/(h*alpha)*(vt-(1+b*alpha)*solutiony21[0]))
#came back to store this for 2.3
pi21 = sm.lambdify((h, alpha, vt, b, st, phi, stminus, ptminus, ytminus, gamma), solutionpi21)
solutionpi21
```

Solving with parameter values:


```python
solutiony21_par = sm.solve(SRASandAD(par['h'], par['alpha'], par['b'], par['phi'], par['gamma']), yt)
solutiony21_par
```


```python
solutionpi21_par = (1/(par['h']*par['alpha'])*(vt-(1+par['b']*par['alpha'])*solutiony21_par[0]))
solutionpi21_par
```

**Question 2:** Find and illustrate the equilibrium when $y_{t-1} = \pi_{t-1} = v_t = s_t = s_{t-1} = 0$. Illustrate how the equilibrium changes when instead $v_t = 0.1$.


```python
ptminus = 0
ytminus = 0
vt = 0
st = 0
stminus = 0
#Finding solution for y_t
solutiony22_par = sm.solve(SRASandAD(par['h'], par['alpha'], par['b'], par['phi'], par['gamma']), yt)
solutiony22_par
```


```python
#solution for pi_t
solutionpi22_par = (1/(par['h']*par['alpha'])*(vt-(1+par['b']*par['alpha'])*solutiony22_par[0]))
solutionpi22_par
```


```python
#now vt = 0.1
#Solution for y
vt = 0.1
solutiony22_vt = sm.solve(SRASandAD(par['h'], par['alpha'], par['b'], par['phi'], par['gamma']), yt)
solutiony22_vt
```


```python
#Solution for PI
solutionpi22_vt = (1/(par['h']*par['alpha'])*(vt-(1+par['b']*par['alpha'])*solutiony22_vt[0]))
solutionpi22_vt
```


```python
#plotting the AS-AD model
xs2 = np.linspace(-1,1)
novtsolasad = lambda x: (1.0/(par['h']*par['alpha'])*(-(1.0+par['b']*par['alpha'])*x)) #no vt
vtsolasad = lambda x: (1.0/(par['h']*par['alpha'])*(0.1-(1.0+par['b']*par['alpha'])*x)) #vt
srassol = lambda x: par['gamma']*x #sras

novtsol = novtsolasad(xs2)
vtsol = vtsolasad(xs2)
sras = srassol(xs2)

plt.plot(xs2, vtsol, 'b-', label='$V_t = 0.1$')          
plt.plot(xs2, novtsol, 'r-', label = "$V_t = 0$")   
plt.plot(xs2, sras, 'g-', label= "SRAS")
plt.xlabel('$y_t$')
plt.ylabel('inflation')
plt.axis([-0.05, 0.05, -0.1, 0.1])
plt.title("AS-AD model")
plt.legend()
plt.show()
```

**Persistent disturbances:** Now, additionaly, assume that both the demand and the supply disturbances are AR(1) processes

$$ v_{t} = \delta v_{t-1} + x_{t} $$
$$ s_{t} = \omega s_{t-1} + c_{t} $$

where $x_{t}$ is a **demand shock**, and $c_t$ is a **supply shock**. The **autoregressive parameters** are:


```python
par['delta'] = 0.80
par['omega'] = 0.15
```

**Question 3:** Starting from $y_{-1} = \pi_{-1} = s_{-1} = 0$, how does the economy evolve for $x_0 = 0.1$, $x_t = 0, \forall t > 0$ and $c_t = 0, \forall t \geq 0$?


```python
#Let's define the processes
def vt(vtminus, x):
    return par['delta']*vtminus + x
def st(stminus, c):
    return par['omega']*stminus + c

T = 200
x = np.zeros(T)
c = np.zeros(T)
yminus = [0]
piminus = [0]
sminus = [0]
vminus = [0]
x[1] = 0.1

#using the stored solutions from 2.1 to simulate the process

for t in range(1, T):
    vminus.append(vt(vminus[t-1], x[t]))
    sminus.append(st(sminus[t-1], c[t]))
    yminus.append(yt21(par['h'], par['alpha'], vminus[t], par['b'], sminus[t], par['phi'], sminus[t-1], piminus[t-1], yminus[t-1], par['gamma']))
    piminus.append(pi21(par['h'], par['alpha'], vminus[t], par['b'], sminus[t], par['phi'], sminus[t-1], piminus[t-1], yminus[t-1], par['gamma']))
    
fig, (ax1, ax2) = plt.subplots(1, 2)
ax1.plot(yminus, label = "Output gap")
ax1.set_title("Output over time")
ax2.plot(piminus, label = "Inflation gap")
ax2.set_title("Inflation over time")
ax1.set_xlabel("Time")
ax2.set_xlabel("Time")
ax1.set_ylabel("Output")
ax2.set_ylabel("Inflation")
fig.set_figwidth(20)
plt.show() 
    
```

**Stochastic shocks:** Now, additionally, assume that $x_t$ and $c_t$ are stochastic and normally distributed

$$ x_{t}\sim\mathcal{N}(0,\sigma_{x}^{2}) $$
$$ c_{t}\sim\mathcal{N}(0,\sigma_{c}^{2}) $$

The **standard deviations of the shocks** are:


```python
par['sigma_x'] = 3.492
par['sigma_c'] = 0.2
```

**Question 4:** Simulate the AS-AD model for 1,000 periods. Calculate the following five statistics:

1. Variance of $y_t$, $var(y_t)$
2. Variance of $\pi_t$, $var(\pi_t)$
3. Correlation between $y_t$ and $\pi_t$, $corr(y_t,\pi_t)$
4. Auto-correlation between $y_t$ and $y_{t-1}$, $corr(y_t,y_{t-1})$
5. Auto-correlation between $\pi_t$ and $\pi_{t-1}$, $corr(\pi_t,\pi_{t-1})$


```python
#Same approach as in 2.3. New distribution, though. I should have made a function in hindsight 
#Function is defined in the next cell. 

T = 1000
np.random.seed(6969)
x = np.random.normal(loc = 0, scale = par['sigma_x'], size = T)
c = np.random.normal(loc = 0, scale = par['sigma_c'], size = T)
yminus = [0]
piminus = [0]
sminus = [0]
vminus = [0]
x[1] = 0.1



for t in range(1, T):
    vminus.append(vt(vminus[t-1], x[t]))
    sminus.append(st(sminus[t-1], c[t]))
    yminus.append(yt21(par['h'], par['alpha'], vminus[t], par['b'], sminus[t], par['phi'], sminus[t-1], piminus[t-1], yminus[t-1], par['gamma']))
    piminus.append(pi21(par['h'], par['alpha'], vminus[t], par['b'], sminus[t], par['phi'], sminus[t-1], piminus[t-1], yminus[t-1], par['gamma']))
    
fig, (ax1, ax2) = plt.subplots(1, 2)
ax1.plot(yminus, label = "Output gap")
ax1.set_title("Output over time")
ax2.plot(piminus, label = "Inflation gap")
ax2.set_title("Inflation over time")
ax1.set_xlabel("Time")
ax2.set_xlabel("Time")
ax1.set_ylabel("Output")
ax2.set_ylabel("Inflation")
fig.set_figwidth(20)
plt.show() 
    

```


```python
print(f'The variance is (y, pi) = {np.var(yminus), np.var(piminus)}')
```

    The variance is (y, pi) = (1.8145391852221522, 1.0238862692255288)



```python
print(f'The correlation coefficient between y and pi is: {np.corrcoef(yminus, piminus)}')
```

    The correlation coefficient between y and pi is: [[ 1.         -0.13784279]
     [-0.13784279  1.        ]]



```python
print(f'The autocorrelation between yt and y_(t-1) is: {np.corrcoef(np.array([yminus[:-1], yminus[1:]]))}')
```

    The autocorrelation between yt and y_(t-1) is: [[1.         0.76075706]
     [0.76075706 1.        ]]



```python
print(f'The autocorrelation between pi_t and pi_(t-1) is: {np.corrcoef(np.array([piminus[:-1], piminus[1:]]))}')
```

    The autocorrelation between pi_t and pi_(t-1) is: [[1.         0.97795762]
     [0.97795762 1.        ]]


**Question 5:** Plot how the correlation between $y_t$ and $\pi_t$ changes with $\phi$. Use a numerical optimizer or root finder to choose $\phi\in(0,1)$ such that the simulated correlation between $y_t$ and $\pi_t$ comes close to 0.31. 


```python
def rootfinder_phi(phi):
    T = 1000
    np.random.seed(6969)
    x = np.random.normal(loc = 0, scale = par['sigma_x'], size = T)
    c = np.random.normal(loc = 0, scale = par['sigma_c'], size = T)
    yminus = [0]
    piminus = [0]
    sminus = [0]
    vminus = [0]
    x[1] = 0.1
    #Adding the correlation
    corr_y_pi = [0]
    for t in range(1, T):
        vminus.append(vt(vminus[t-1], x[t]))
        sminus.append(st(sminus[t-1], c[t]))
        yminus.append(yt21(par['h'], par['alpha'], vminus[t], par['b'], sminus[t], phi, sminus[t-1], piminus[t-1], yminus[t-1], par['gamma']))
        piminus.append(pi21(par['h'], par['alpha'], vminus[t], par['b'], sminus[t], phi, sminus[t-1], piminus[t-1], yminus[t-1], par['gamma']))
        corr_y_pi = np.corrcoef(yminus, piminus)[1, 0]
    return yminus, piminus, corr_y_pi


```


```python
opti_fun = lambda phi: np.corrcoef(rootfinder_phi(phi)[0], rootfinder_phi(phi)[1])[1, 0]-0.31

opti_fun_root = optimize.brentq(opti_fun, a = 0, b = 1, full_output = False)
print(f'The value of phi, such that the correlation coefficient is closest to 0.31 is: {opti_fun_root}')
```

    The value of phi, such that the correlation coefficient is closest to 0.31 is: 0.9562298823595026


**Quesiton 6:** Use a numerical optimizer to choose $\sigma_x>0$, $\sigma_c>0$ and $\phi\in(0,1)$ to make the simulated statistics as close as possible to US business cycle data where:

1. $var(y_t) = 1.64$
2. $var(\pi_t) = 0.21$
3. $corr(y_t,\pi_t) = 0.31$
4. $corr(y_t,y_{t-1}) = 0.84$
5. $corr(\pi_t,\pi_{t-1}) = 0.48$


```python
#Getting an error and cannot think of a way to pass it. My approach would have been extending the method above to fit in multiple criteria.

#def USbusiness(phi, c, x):
#    T = 1000
#    np.random.seed(6969)
#    x = np.random.normal(loc = 0, scale = par['sigma_x'], size = T)
#    c = np.random.normal(loc = 0, scale = par['sigma_c'], size = T)
#    yminus = [0]
#    piminus = [0]
#    sminus = [0]
#    vminus = [0]
#    x[1] = 0.1
#    corr_y_pi = [0]
#    corr_y_y = [0]
#    corr_pi_pi = [0]
#    var_y = [0]
#    var_pi = [0]
    
#    for t in range(1, T):
#        vminus.append(vt(vminus[t-1], x[t]))
#        sminus.append(st(sminus[t-1], c[t]))
#        yminus.append(yt21(par['h'], par['alpha'], vminus[t], par['b'], sminus[t], phi, sminus[t-1], piminus[t-1], yminus[t-1], par['gamma']))
#        piminus.append(pi21(par['h'], par['alpha'], vminus[t], par['b'], sminus[t], phi, sminus[t-1], piminus[t-1], yminus[t-1], par['gamma']))
#        corr_y_pi = np.corrcoef(yminus, piminus)[1, 0]
#        corr_y_y = np.corrcoef(yminus[1:], yminus[:-1])[1, 0]
#        corr_pi_pi = np.corrcoef(piminus[1:], piminus[:-1])[1, 0]
#        var_y = np.var(yminus)
#        var_pi = np.var(piminus)
#        return var_pi, var_y, yminus, piminus, corr_y_y, corr_pi_pi, corr_y_pi
    
#var_pi, var_y, yminus, piminus, corr_y_y, corr_pi_pi, corr_y_pi   = USbusiness(pi, c, x)
    
```

# 3. Exchange economy

Consider an **exchange economy** with

1. 3 goods, $(x_1,x_2,x_3)$
2. $N$ consumers indexed by \\( j \in \{1,2,\dots,N\} \\)
3. Preferences are Cobb-Douglas with log-normally distributed coefficients

    $$ \begin{eqnarray*}
    u^{j}(x_{1},x_{2},x_{3}) &=& 
    \left(x_{1}^{\beta_{1}^{j}}x_{2}^{\beta_{2}^{j}}x_{3}^{\beta_{3}^{j}}\right)^{\gamma}\\
     &  & \,\,\,\beta_{i}^{j}=\frac{\alpha_{i}^{j}}{\alpha_{1}^{j}+\alpha_{2}^{j}+\alpha_{3}^{j}} \\
     &  & \,\,\,\boldsymbol{\alpha}^{j}=(\alpha_{1}^{j},\alpha_{2}^{j},\alpha_{3}^{j}) \\ 
     &  & \,\,\,\log(\boldsymbol{\alpha}^j) \sim \mathcal{N}(\mu,\Sigma) \\
    \end{eqnarray*} $$

4. Endowments are exponentially distributed,

$$
\begin{eqnarray*}
\boldsymbol{e}^{j} &=& (e_{1}^{j},e_{2}^{j},e_{3}^{j}) \\
 &  & e_i^j \sim f, f(z;\zeta) =  1/\zeta \exp(-z/\zeta)
\end{eqnarray*}
$$

Let $p_3 = 1$ be the **numeraire**. The implied **demand functions** are:

$$
\begin{eqnarray*}
x_{i}^{\star j}(p_{1},p_{2},\boldsymbol{e}^{j})&=&\beta^{j}_i\frac{I^j}{p_{i}} \\
\end{eqnarray*}
$$

where consumer $j$'s income is

$$I^j = p_1 e_1^j + p_2 e_2^j +p_3 e_3^j$$

The **parameters** and **random preferences and endowments** are given by:


```python
# a. parameters
N = 50000
mu = np.array([3,2,1])
Sigma = np.array([[0.25, 0, 0], [0, 0.25, 0], [0, 0, 0.25]])
gamma = 0.8
zeta = 1

# b. random draws
seed = 1986
np.random.seed(seed)

# preferences
alphas = np.exp(np.random.multivariate_normal(mu, Sigma, size=N))
betas = alphas/np.reshape(np.sum(alphas,axis=1),(N,1))

# endowments
e1 = np.random.exponential(zeta,size=N)
e2 = np.random.exponential(zeta,size=N)
e3 = np.random.exponential(zeta,size=N)

#Let's check how beta looks for the upcoming task
betas
```




    array([[0.53104511, 0.30693723, 0.16201766],
           [0.50866997, 0.4162364 , 0.07509363],
           [0.83317207, 0.07722226, 0.08960567],
           ...,
           [0.63725404, 0.22027742, 0.14246854],
           [0.78205036, 0.11988119, 0.09806844],
           [0.65552751, 0.24728801, 0.09718448]])



**Question 1:** Plot the histograms of the budget shares for each good across agents.


```python
fig, ax = plt.subplots(1, 3)
ax[0].hist(betas[:,0], bins = 100, density = 'True')
ax[0].set_title("Good 1")
ax[1].hist(betas[:,1], bins = 100, density = 'True')
ax[1].set_title("Good 2")
ax[2].hist(betas[:,2], bins = 100, density = 'True')
ax[2].set_title("Good 3")
fig.set_figwidth(20)
plt.show()
```

Consider the **excess demand functions:**

$$ z_i(p_1,p_2) = \sum_{j=1}^N x_{i}^{\star j}(p_{1},p_{2},\boldsymbol{e}^{j}) - e_i^j$$
**Question 2:** Plot the excess demand functions.


```python
#Demand for each good

def x1(p1, p2, e1, e2, e3, betas):
    I = p1*e1 + p2*e2 + e3
    return betas[:,0]*I/p1

def x2(p1, p2, e1, e2, e3, betas):
    I = p1*e1 + p2*e2 + e3
    return betas[:,1]*I/p2

def x3(p1, p2, e1, e2, e3, betas):
    I = p1*e1 + p2*e2 + e3
    return betas[:,2]*I

#excess demand

def excess_demand(x, p1, p2, e1, e2, e3, betas):
    demand = np.sum(x(p1, p2, e1, e2, e3, betas))
    
    if x == "x1":
        supply = np.sum(e1)
    elif x == "x2":
        supply = np.sum(e2)
    else:
        supply = np.sum(e3)
    exc_dem = demand - supply     
    return exc_dem  




#price_1 = np.linspace(0.1, 10, 100)
#price_2 = np.linspace(0.1, 10, 100)
#price_3 = np.linspace(0.1, 10, 100)

#x1_fun = lambda x: excess_demand(x1, x, p2, e1, e2, e3, betas)
#x1_fun_plot = x1_fun(price_1)
#x2_fun = lambda x: excess_demand(x2, p1, x, e1, e2, e3, betas)
#x2_fun_plot = x2_fun(price_2)


#plt.plot(price_1, x1_fun_plot, "r-", label = "Product 1")
#plt.plot(price_2, x2_fun_plot, "b-", label = "Product 2" )
#plt.show()

#Since the plot requires three axes it got quite messy and with errors all over the place. I ended up deleting it and not being able to recreate anything meaningful.
```

**Question 3:** Find the Walras-equilibrium prices, $(p_1,p_2)$, where both excess demands are (approximately) zero, e.g. by using the following tâtonnement process:

1. Guess on $p_1 > 0$, $p_2 > 0$ and choose tolerance $\epsilon > 0$ and adjustment aggressivity parameter, $\kappa > 0$.
2. Calculate $z_1(p_1,p_2)$ and $z_2(p_1,p_2)$.
3. If $|z_1| < \epsilon$ and $|z_2| < \epsilon$ then stop.
4. Else set $p_1 = p_1 + \kappa \frac{z_1}{N}$ and $p_2 = p_2 + \kappa \frac{z_2}{N}$ and return to step 2.


```python
#Defining the criteria written above

p1_ini = 1
p2_ini = 1
#I had to set a cap on iterations or the loop will spin forever in q5
def equilibrium(iteration, p1, p2, e1, e2, e3, betas, kappa, epsilon):
    t=0
    while True:
        z1 = excess_demand(x1, p1, p2, e1, e2, e3, betas)
        z2 = excess_demand(x2, p1, p2, e1, e2, e3, betas)
        if np.abs(z1) < epsilon and np.abs(z2) < epsilon or t > iteration: #Needed to avoid waiting for minutes.
            print(f'The prices are (p1, p2) = {p1, p2}, and the excess demand is (z1, z2) = {z1, z2}. We keep in mind that p3 = 1')
            break
        p1 = p1 + kappa*z1/N
        p2 = p2 + kappa*z2/N
        t += 1
    return p1, p2
p1, p2 = equilibrium(10000, p1_ini, p2_ini, e1, e2, e3, betas, 2, 1e-10)

```

    The prices are (p1, p2) = (6.0596117842277755, 2.433061390839292), and the excess demand is (z1, z2) = (9.458744898438454e-11, 3.637978807091713e-11). We keep in mind that p3 = 1


**Question 4:** Plot the distribution of utility in the Walras-equilibrium and calculate its mean and variance.


```python
def utility(x1, x2, x3, betas, gamma):
    return (x1**betas[:,0] + x2**betas[:,1] + x3**betas[:,2])**gamma

#using demand functions from above
utility_wal = utility(x1(p1, p2, e1, e2, e3, betas), x2(p1, p2, e1, e2, e3, betas), x3(p1, p2, e1, e2, e3, betas), betas, gamma)
print(f'The mean of the utility is: {utility_wal.mean()}')
```

    The mean of the utility is: 2.3769975524551423



```python
print(f'The variance of the utility is: {utility_wal.var()}')
```

    The variance of the utility is: 0.20702487229161537



```python
#plotting the distribution
fig = plt.figure()
ax = fig.add_subplot(1,1,1)
ax.hist(utility_wal, bins = 100, density = True)
ax.set_title("Distribution of utility of the walras equilibrium")
plt.show()
```

**Question 5:** Find the Walras-equilibrium prices if instead all endowments were distributed equally. Discuss the implied changes in the distribution of utility. Does the value of $\gamma$ play a role for your conclusions?


```python
#Defining new shares of endowments
new_e = (e1+e2+e3)/3
p1_5, p2_5 = equilibrium(10000, p1_ini, p2_ini, new_e, new_e, new_e, betas, 2, 1e-10)


utility_wal_5  = utility(x1(p1_5, p2_5, new_e, new_e, new_e, betas), x2(p1_5, p2_5, new_e, new_e, new_e, betas), x3(p1_5, p2_5, new_e, new_e, new_e, betas), betas, gamma)
utility_wal_5_dif_gam  = utility(x1(p1_5, p2_5, new_e, new_e, new_e, betas), x2(p1_5, p2_5, new_e, new_e, new_e, betas), x3(p1_5, p2_5, new_e, new_e, new_e, betas), betas, 0.5)

```

    The prices are (p1, p2) = (6.441793456362303, 2.5905935972737115), and the excess demand is (z1, z2) = (9.458744898438454e-11, 4.3655745685100555e-11). We keep in mind that p3 = 1



```python
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.hist(utility_wal_5, bins = 100, density = True)
ax.set_title("Distribution of utility of the walras equilibrium with $gamma = 0.8$")
plt.show()
```


```python
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.hist(utility_wal_5_dif_gam, bins = 100, density = True)
ax.set_title("Distribution of utility of the walras equilibrium with $gamma = 0.5$")
plt.show()
```

We see that the average utility lowers with a decreasing gamma. If we compare it to the previous one in 3.4, we see that the variance is a lot smaller but that the mean stays roughly the same ith an unchanged gamma. A rough conclusion could therefore be, that the distribution of endowments have an insignificant effect on the utility in an exchange economy.


```python

```

84664c00d95a6356ef0d5ab973e79ca333949c7b
ipynb
Jupyter Notebook
examproject/exam_2019.ipynb
NumEconCopenhagen/projects-2019-aloise3
8776e44395ef679d25628f39a7244f351c554ecb
MIT
examproject/exam_2019.ipynb
NumEconCopenhagen/projects-2019-aloise3
8776e44395ef679d25628f39a7244f351c554ecb
MIT
2019-04-16T16:29:32.000Z
2019-05-22T14:03:25.000Z
examproject/exam_2019.ipynb
NumEconCopenhagen/projects-2019-aloise3
8776e44395ef679d25628f39a7244f351c554ecb
MIT
Qwen/Qwen-72B
1. YES
2. YES