<jupyter_start><jupyter_text>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:<jupyter_code>rho = 2.0
beta = 0.96
gamma = 0.1
w = 2.0
b = 1.0
delta = 0.1<jupyter_output><empty_output><jupyter_text>The **relevant levels of human capital** are:<jupyter_code>h_vec = np.linspace(0.1,1.5,100)<jupyter_output><empty_output><jupyter_text>**Question 1:** Solve the model in period 2 and illustrate the solution (including labor supply as a function of human capital).<jupyter_code>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")
#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.')
#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()<jupyter_output><empty_output><jupyter_text>**Question 2:** Solve the model in period 1 and illustrate the solution (including labor supply as a function of human capital).<jupyter_code>#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}')
#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()<jupyter_output><empty_output><jupyter_text>**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.<jupyter_code>#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()<jupyter_output><empty_output><jupyter_text>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:<jupyter_code>par = {}

par['alpha'] = 5.76
par['h'] = 0.5
par['b'] = 0.5
par['phi'] = 0
par['gamma'] = 0.075<jupyter_output><empty_output><jupyter_text>**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}$.<jupyter_code>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

AD = sm.Eq(pt,(1/(h*alpha)*(vt-(1+b*alpha)*yt)))
AD
#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)
#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
#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<jupyter_output><empty_output><jupyter_text>Solving with parameter values:<jupyter_code>solutiony21_par = sm.solve(SRASandAD(par['h'], par['alpha'], par['b'], par['phi'], par['gamma']), yt)
solutiony21_par
solutionpi21_par = (1/(par['h']*par['alpha'])*(vt-(1+par['b']*par['alpha'])*solutiony21_par[0]))
solutionpi21_par<jupyter_output><empty_output><jupyter_text>**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$.<jupyter_code>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
#solution for pi_t
solutionpi22_par = (1/(par['h']*par['alpha'])*(vt-(1+par['b']*par['alpha'])*solutiony22_par[0]))
solutionpi22_par
#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
#Solution for PI
solutionpi22_vt = (1/(par['h']*par['alpha'])*(vt-(1+par['b']*par['alpha'])*solutiony22_vt[0]))
solutionpi22_vt
#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()<jupyter_output><empty_output><jupyter_text>**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:<jupyter_code>par['delta'] = 0.80
par['omega'] = 0.15<jupyter_output><empty_output><jupyter_text>**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$?<jupyter_code>#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()<jupyter_output><empty_output><jupyter_text>**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:<jupyter_code>par['sigma_x'] = 3.492
par['sigma_c'] = 0.2<jupyter_output><empty_output><jupyter_text>**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})$<jupyter_code>#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() 
    

print(f'The variance is (y, pi) = {np.var(yminus), np.var(piminus)}')
print(f'The correlation coefficient between y and pi is: {np.corrcoef(yminus, piminus)}')
print(f'The autocorrelation between yt and y_(t-1) is: {np.corrcoef(np.array([yminus[:-1], yminus[1:]]))}')
print(f'The autocorrelation between pi_t and pi_(t-1) is: {np.corrcoef(np.array([piminus[:-1], piminus[1:]]))}')<jupyter_output>The autocorrelation between pi_t and pi_(t-1) is: [[1.         0.97795762]
 [0.97795762 1.        ]]<jupyter_text>**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.<jupyter_code>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


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}')<jupyter_output>The value of phi, such that the correlation coefficient is closest to 0.31 is: 0.9562298823595026<jupyter_text>**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$<jupyter_code>#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)<jupyter_output><empty_output><jupyter_text>3. Exchange economy Consider an **exchange economy** with1. 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:<jupyter_code># 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<jupyter_output><empty_output><jupyter_text>**Question 1:** Plot the histograms of the budget shares for each good across agents.<jupyter_code>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()<jupyter_output><empty_output><jupyter_text>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.<jupyter_code>#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.<jupyter_output><empty_output><jupyter_text>**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.<jupyter_code>#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)<jupyter_output>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<jupyter_text>**Question 4:** Plot the distribution of utility in the Walras-equilibrium and calculate its mean and variance.<jupyter_code>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()}')
print(f'The variance of the utility is: {utility_wal.var()}')
#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()<jupyter_output><empty_output><jupyter_text>**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?<jupyter_code>#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)

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()
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()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.19544767344196204, 'fraction_numerical': 0.04598521753737611, 'mean_word_length': 4.65360873694207, 'pattern_counts': {'":': 2, '<': 86, '<?xml version=': 0, '>': 91, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '21', 'dir': 'jupyter-structured-clean-dedup', 'id': '449145', 'license': 'MIT', 'n_tokens_mistral': 10221, 'n_tokens_neox': 9667, 'n_words': 2472, 'path': 'examproject/exam_2019.ipynb', 'repo_name': 'NumEconCopenhagen/projects-2019-aloise3', 'size': '23812'}