# 1. Human capital accumulation

## Question 1
First we will define each variable and parameters with Sympy.symbols.


```python
# import packages
import numpy as np
import sympy as sm
sm.init_printing(use_unicode=True) # for pretty printing
```


```python
# Define parameters and variables as symbols
rho = sm.symbols("rho")
beta = sm.symbols("beta")
gamma = sm.symbols("gamma")
w = sm.symbols("w")
b = sm.symbols("b")
Delta = sm.symbols("Delta")
h_2 = sm.symbols("h_2")
l_2 = sm.symbols("l_2")
c_2 = sm.symbols("c_2")
l_1 = sm.symbols("l_1")
h_1 = sm.symbols("h_1")
```

The utility the worker receives from consumption depends on the variable $l_t$. If the value is [1], consumption will be $w*h_2*l_2$, and for the value [0], $(1-l_2)*b$. Hence we define the function as the following:
$c_2 = (w*h_2*l_2)+(1-l_2)*b$


```python
# Define consumption function according to the pres-ented formula
# Includes, wage w, unemployment benefit b, human capital for period 2 h_2 and labor supply for period 2 l_2
# Returns the consumption
def consumption_2(w, h_2, l_2, b):
    return (w*h_2*l_2)+(1-l_2)*b
```

Next we will try to define the utility function.


```python
# Define utility function according to the presented formula
utility = c_2 ** (1-rho)/(1-rho) - gamma * l_2
utility
```

Then we will use sympy to substitute the consumption function into the utility function


```python
# Set consumption equal to the variables and parameters
consumption2 = w*h_2*l_2+(1-l_2)*b
consumption2
```


```python
# Use sympy.subs to substitute consumption_2 into the utility function.
utility_subs = utility.subs(c_2,consumption2)
utility_subs
```

This is the final utility function, which we will define below.


```python
# Define final utility function.
def utility_2(gamma, w, b, rho, h, l):
    return ((1-l)*b+(h*l*w))**(1-rho)/(1-rho)-gamma*l
```

Next we will try to solve the model visually.


```python
# Import the following packages
%matplotlib inline
import matplotlib.pyplot as plt
```

Next we will use  np.argmax to solve the model for period 2 by using the given vector for $h_2$ and an empty vector (100,2), "choices2".


```python
# Vectors
choices2 = np.empty((100,2))
h2_vec = np.linspace(0.1,1.5,100) 

# Parameters
rho = 2
beta = 0.96
gamma = 0.1
w = 2
b = 1
Delta = 0.1

# loop over each value of h for l=0 and l=1 for the utility function to fill up choices2
for index, value in enumerate(h2_vec):
    choices2[index, 0] = utility_2(gamma, w, b, rho, value, 0)
    choices2[index, 1] = utility_2(gamma, w, b, rho, value, 1)
```

Now we can find the solutions for labor supply $l_2$, consumption $c_2$ and value $v_2$.


```python
# Find the argument, which maximise the value of l2 given choices2
l2_sol = np.argmax(choices2, axis=1)
# Insert the above sol into consumption.
c2_sol = consumption_2(w, b ,l2_sol,h2_vec) 
# Max choises2 to find v2
v2_sol = np.max(choices2, axis=1)
```

At last we can show the solution in the below subplots.


```python
# Begin plot
fig = plt.figure(figsize=(18,10))

# Add the three subplots
ax = fig.add_subplot(3,3,1)
ax.plot(h2_vec,l2_sol)
ax.grid()
ax.set_xlabel('$h_2$')
ax.set_ylabel('$l_2$')
ax.set_title('Labor supply function in period 2');

ax = fig.add_subplot(3,3,2)
ax.plot(h2_vec,c2_sol, color='green')
ax.grid()
ax.set_xlabel('$h_2$')
ax.set_ylabel('$c_2$')
ax.set_title('Consumption function in period 2');

ax = fig.add_subplot(3,3,3)
ax.plot(h2_vec,v2_sol, color='orange')
ax.grid()
ax.set_xlabel('$h_2$')
ax.set_ylabel('$v_2$')
ax.set_title('Value function in period 2');
```

We see how there is a kink at a certain value of $h_2$, approx. 0.55. At this value, the worker gets more utility from working than living of unemployment benefits. Hence, the worker's consumtion and value increases.

# Question 2
Now the worker has a chance of accumulating human capital.


```python
# Define consumption and utility function (same as period 2)
def consumption_1(w, h_1, l_1, b):
    return (w*h_1*l_1)+(1-l_1)*b

def utility_1(gamma, w, b, rho, h_1, l_1):
    return ((1-l_1)*b+(h_1*l_1*w))**(1-rho)/(1-rho)-gamma*l_1
```

Next we will create a value function, which inputs consits of the possibilities for the worker to accumulate human capital


```python
def v1(gamma, b, w, rho, Delta, beta, h_1, l_1, v2_interp):
    #  h2 value, if no acc
    h2_no_acc = h_1 + l_1
    
    # b2 value, if no acc
    v2_no_acc = v2_interp([h2_no_acc])[0]
    
    # h2 value, if acc
    h2_acc = h_1 + l_1 + Delta
    
    # v2 value if acc
    v2_acc = v2_interp([h2_acc])[0]
    
    # expected v2 value
    v2 = 0.5 * v2_no_acc + 0.5 * v2_acc
    
    # return the total value
    return utility_1(gamma, rho, b, w, h_1, l_1) + beta * v2
```

Construct interpolar


```python
from scipy import interpolate

# Constructing interpolar
v2_interp = interpolate.RegularGridInterpolator([h2_vec], v2_sol,
                                                bounds_error=False,fill_value=None)
```

Similar to how we did in Question 1, we create a empty array and loop over each value of $h$ 


```python
# Vectors
choices1 = np.empty((100,2))


# loop over each value of h for l=0 and l=1 to fill up choices 1
for index, value in enumerate(h2_vec):
    choices1[index, 0] = v1(gamma, b, w, rho, Delta, beta, value, 0, v2_interp)
    choices1[index, 1] = v1(gamma, b, w, rho, Delta, beta, value, 1, v2_interp)

```

Now we can find the following solutions.


```python
# Find the argument, which maximise the value of l1 given choices1
l1_sol = np.argmax(choices1, axis=1)

# Insert the above sol into consumption.
c1_sol = consumption_1(w,b,l1_sol,h2_vec)

# Max choises1 to find v1
v1_sol = np.max(choices1, axis=1)
```


```python
# Begin plotting
fig = plt.figure(figsize=(18,10))

# Add subplots
ax = fig.add_subplot(3,3,1)
ax.plot(h2_vec,l1_sol)
ax.grid()
ax.set_xlabel('$h_1$')
ax.set_ylabel('$l_1$')
ax.set_title('Labor supply in period 1');

ax = fig.add_subplot(3,3,2)
ax.plot(h2_vec,c1_sol, color ='green')
ax.grid()
ax.set_xlabel('$h_1$')
ax.set_ylabel('$c_1$')
ax.set_title('Consumption in period 1');

ax = fig.add_subplot(3,3,3)
ax.plot(h2_vec,v1_sol, color = 'orange')
ax.grid()
ax.set_xlabel('$h_1$')
ax.set_ylabel('$v_1$')
ax.set_title('Value function in period 1');
```

We see that when the worker has a chance of accumulating human capital, the kink happens at a smaller $h$ value, approx. 0.35.

# Question 3
In this question we will explore if the worker will still prefer working even though $b>w$. In order for us to do so, we will research if the worker chooses $l_t$ = [1] for a certain value of $h_t$. For your convenience this will be done only for period 2.

The method is similar to the previous questions.


```python
# Import packages
import matplotlib.pyplot as plt
import ipywidgets as widgets
from ipywidgets import interactive

# Define interactive figure
def interactive_fig_w(w_t):
    
    # Find solution with the method similar to previous questions.
    choices_w = np.empty((100,2))
    
    for index, value in enumerate(h2_vec):
        choices_w[index, 0] = utility_2(gamma, w_t, b, rho, value, 0)
        choices_w[index, 1] = utility_2(gamma, w_t, b, rho, value, 1)
        
    lw_sol = np.argmax(choices_w, axis=1)
    
    # Plot figure
    
    plt.plot(h2_vec,lw_sol)
    plt.title('Labor supply in period 2')
    plt.xlabel('$h_2$')
    plt.ylabel('$l_2$')
    plt.axis([0,1.6,-0.1,1.1])
    plt.show()


interactive_plot_w = interactive(interactive_fig_w, 
                               w_t=widgets.FloatSlider(description="$w$", 
                                                       min=0.0,
                                                       max=1.0,
                                                       step=0.01,
                                                       value=0.9,
                                                       continuous_update=True) 
                              )

interactive_plot_w
```


    interactive(children=(FloatSlider(value=0.9, description='$w$', max=1.0, step=0.01), Output()), _dom_classes=(…


What the above figure tells us is that when the wage is $w = 0.75$ or above, the worker gets more utility from working than relying on unemployment benefits $b = 1$. So yes, the worker will still work ($l_1 = 1$), when the worker's wage is below the unemployment benefits. 

# 2. AS-AD model


```python
# Create dictionary for parameters and their values
par = {}

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

# Question 1
First we will define the following variables and parameters as symbols.


```python
# Defining variables and parameters as symbols
alpha = sm.symbols("alpha")
b = sm.symbols("b")
h = sm.symbols("h")
phi = sm.symbols("phi")
gamma = sm.symbols("gamma")

y_t = sm.symbols("y_t")
y_t_1 = sm.symbols("y_{t-1}")
s_t = sm.symbols("s_t")
s_t_1 = sm.symbols("s_{t-1}")
v_t = sm.symbols("v_t")
pi_t = sm.symbols("pi_t")
pi_t_1 = sm.symbols("\pi_{t-1}")
```

Then we are able to construct the following equations for the SRAS curve and the AD curve.


```python
# Constructing SRAS and AD curve
SRAS = sm.Eq(pi_t,pi_t_1 + gamma * y_t  - phi * gamma * y_t_1 + s_t - phi * s_t_1)
AD = sm.Eq(pi_t, 1/(h * alpha) * (v_t - (1 + b * alpha) * y_t))
SRAS, AD
```

To solve the model, we will have to set the two curves equal to eachother.


```python
SRAS_AD = sm.Eq(
                pi_t_1 + gamma * y_t  - phi * gamma * y_t_1 + s_t - phi * s_t_1, 
                1/(h * alpha) * (v_t - (1 + b * alpha) * y_t))
SRAS_AD
```

Then we can solve for $y_t$


```python
# Solve the equation for y_t
yt_sol = sm.solve(SRAS_AD,y_t)[0]
print(f'The solution for output y_t is:')
yt_sol
```

Then insert the above solution into the AD-curve to find $\pi_t$


```python
# Sub yt_sol into the AD equation.
pit_sol = 1/(h * alpha) * (v_t - (1 + b * alpha) * yt_sol)
print(f'The solution for inflaton pi_t is:')
pit_sol
```

To find the equilibrium values for  output, $y_t$, and inflation, $\pi_t$,given the parameters ($\alpha$, $h$, $b$, $\alpha$, $\gamma$) and $y_{t-1}$ , $\pi_{t-1}$, $v_t$, $s_t$, and $s_{t-1}$, we will substitute the given parameters into the solution.


```python
yt_val = yt_sol.subs({alpha: par['alpha'],
                      h: par['h'],
                      b: par['b'],
                      phi: par['phi'],
                      gamma: par['gamma']})
print(f'The equilibrium for output y_t given the above parameters is:')
yt_val
```


```python
pit_val = pit_sol.subs({alpha: par['alpha'],
                      h: par['h'],
                      b: par['b'],
                      phi: par['phi'],
                      gamma: par['gamma']})
print(f'The equilibrium for inflation pi_t given the above parameters is:')
pit_val
```

# Question 2
We will first the initial values of $y_{t-1} = \pi_{t-1} = v_t = s_t = s_{t-1} = 0$. Next we will transform the sympy equations to lambda functions to calculate numerical values for $y_t$ and $\pi_t$


```python
# Create dict for var and set initials value to 0

var = {}

var['pi_t_1_val'] = 0
var['v_t_val'] = 0
var['y_t_1_val'] = 0
var['s_t_val'] = 0
var['s_t_1_val'] = 0

# Transform the sympy equations into lambda functions

_yt_sol_func = sm.lambdify((alpha, h, b, phi, gamma, pi_t_1, v_t, y_t_1, s_t, s_t_1),yt_sol)

def yt_sol_func(pi_t_1, v_t, y_t_1, s_t, s_t_1,
                    alpha = par['alpha'],
                    h = par['h'],
                    b= par['b'],
                    phi= par['phi'],
                    gamma= par['gamma']):

    return _yt_sol_func(alpha, h, b, phi, gamma, pi_t_1, v_t, y_t_1, s_t, s_t_1)


# Calculate the numerical values
yt_sol_val = yt_sol_func(
                         pi_t_1 = var['pi_t_1_val'], 
                         v_t = var['v_t_val'], 
                         y_t_1 = var['y_t_1_val'], 
                         s_t = var['s_t_val'], 
                         s_t_1 = var['s_t_1_val']
                        )

f"The initial equilibrium value of y_t given the parameters and variables is  {yt_sol_val:.3f}"
```




    'The initial equilibrium value of y_t given the parameters and variables is  0.000'




```python
# Transform the sympy equations into lambda functions
_pi_sol_func = sm.lambdify((alpha, h, b, phi, gamma, pi_t_1, v_t, y_t_1, s_t, s_t_1, y_t),pit_sol)

def pi_sol_func(pi_t_1, v_t, y_t_1, s_t, s_t_1, y_t,
                    alpha = par['alpha'],
                    h = par['h'],
                    b= par['b'],
                    phi= par['phi'],
                    gamma= par['gamma']):

    return _pi_sol_func(alpha, h, b, phi, gamma, pi_t_1, v_t, y_t_1, s_t, s_t_1, y_t)

pi_sol_val = pi_sol_func(
                         pi_t_1 = var['pi_t_1_val'], 
                         v_t = var['v_t_val'], 
                         y_t_1 = var['y_t_1_val'], 
                         s_t = var['s_t_val'], 
                         s_t_1 = var['s_t_1_val'],
                         y_t = yt_sol_val)
f"The initial equilibrium value of pi_t given the parameters and variables is  {pi_sol_val:.3f}"
```




    'The initial equilibrium value of pi_t given the parameters and variables is  0.000'



Now we will try to illustrate the AD_curve and SRAS_curve in a AS-AD diagram. We will also create a slider to simulate the changing value of $v_t$


```python
# First define the two curves as functions.
def AD_curve(alpha, h, b, v_t, y_t):
    return 1/(h * alpha) * (v_t - (1 + b * alpha) * y_t)

def SRAS_curve(gamma, phi, pi_t_1, y_t_1, s_t, s_t_1, y_t):
    return pi_t_1 + gamma * y_t  - phi * gamma * y_t_1 + s_t - phi * s_t_1
```

Creating the plot


```python
# Define interactive figure
def interactive_fig(v_t_val):

    # Range for output, y_t
    y = np.linspace(-0.1, 0.1, 100)
    
    # Plot the optimal values for y_t and pi_t
    y_opt_val = yt_sol_func(
                         pi_t_1 = var['pi_t_1_val'], 
                         v_t = v_t_val, 
                         y_t_1 = var['y_t_1_val'], 
                         s_t = var['s_t_val'], 
                         s_t_1 = var['s_t_1_val'],
                        )
    
    pi_opt_val = pi_sol_func(
                         pi_t_1 = var['pi_t_1_val'], 
                         v_t = v_t_val, 
                         y_t_1 = var['y_t_1_val'], 
                         s_t = var['s_t_val'], 
                         s_t_1 = var['s_t_1_val'],
                         y_t = y
                        )

    # xlabel, ylabel and title
    plt.xlabel(r'$y_t$', )
    plt.ylabel(r'$\pi_{t}$', fontsize = 15)
    plt.title(r'Equilibrium for AD and SRAS', fontsize = 15)

    # Shorten
    AD_c = AD_curve(
                alpha = par["alpha"],
                h = par['h'],
                b = par['b'], 
                v_t = v_t_val,
                y_t = y
               )
    SRAS_c = SRAS_curve(
                gamma = par["gamma"],
                phi = par['phi'], 
                pi_t_1 = var['pi_t_1_val'], 
                y_t_1 = var['y_t_1_val'], 
                s_t = var['s_t_val'], 
                s_t_1 = var['s_t_1_val'], 
                y_t = y
                )
    # Plot the two curves
    plt.plot(y, AD_c, label = "AD-Curve")
    plt.plot(y, SRAS_c, label = "SRAS_Curve")

    # Create a dashed lines to show optimal values for y_t and pi_t
    plt.axvline(y_opt_val, 
                linestyle = 'dashed', 
                color = "red", 
                label = '$y^*$')
    
    plt.axhline(pi_opt_val, 
                linestyle = 'dashed', 
                color = "green", 
                label = '$\pi^*$')

    plt.legend()
    
    plt.axis([-0.15,0.15,-0.18,0.18])
    
    # Create arrow, which points at the equilibrium
    plt.annotate('Equilibrium', xy=(y_opt_val, pi_opt_val), xytext=(y_opt_val-0.10, pi_opt_val-0.05),
            arrowprops=dict(facecolor='black', shrink=0.01))

    plt.show()

# create slider and plot
interactive_plot = interactive(interactive_fig, 
                               v_t_val=widgets.FloatSlider(description="$v_t$", 
                                                       min=0.0,
                                                       max=0.1,
                                                       step=0.1,
                                                       value=0,
                                                       continuous_update=True) 
                              )

interactive_plot
```


    interactive(children=(FloatSlider(value=0.0, description='$v_t$', max=0.1), Output()), _dom_classes=('widget-i…


We see how the $v_t$ shock moves the AD-curve up and how output increases.

# Question 3
For this question we will begin by setting the initial values to 0. Then create an array for the two shock variables and filling them up with the autoregressive parameters to simulate the economy.


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


```python
# Initial values
y_1 = 0
pi_0 = 0
s_0 = 0
v_0 = 0
x_0 = 0.1
x_t = 0
c_t = 0

# Set periods to 100
period = 100

# Create array and fill in value for initial period
v_t = [ par["delta"] * v_0 + x_0 ]
s_t = [ par["omega"] * s_0 + c_t ]

# Fill up arrays for each period
for i in range(0, period):
    v_t.append(par['delta'] * v_t[i] + x_t)
    s_t.append(par['omega'] * s_t[i] + c_t)
```

Next we want to find all the numeric solutions for $y_t$ and $\pi_t$. We will create an empty array for $y_t$ and $\pi_t$ and then fill it up with the solutions from the simulation.


```python
# Create arrays
y_vals = [0]
pi_vals = [0]

# Loop over each period  and append the solution to the array.
for i in range(1, period):
    
    # y
    y_vals_i = yt_sol_func(
        v_t=v_t[i],
        s_t=s_t[i],
        pi_t_1=pi_vals[i-1],
        y_t_1=y_vals[i-1],
        s_t_1=s_t[i-1]
        )
    y_vals.append(y_vals_i)
    # Pi
    pi_vals_i = pi_sol_func(
        v_t=v_t[i],
        s_t=s_t[i],
        y_t=y_vals[i],
        pi_t_1=pi_vals[i-1],
        y_t_1=y_vals[i-1],
        s_t_1=s_t[i-1]
        )
    pi_vals.append(pi_vals_i)
```

Now we will illustrate this shock to the economy. First we will define functions for the AD curve and the SRAS curve. For the first period, the function will take the initial input.


```python
# Define functions, which returns the initial numeric solution for the first period and otherwise returns 
# the value for the simulated periods to come.
def AD_curves_func(y_t,period,
             h = par['h'],
             b = par['b'],
             alpha = par['alpha'],
             v_t = v_t): 
    if period == 0: 
        return 1 / (h * alpha) * (v_0 - (1 + b * alpha) * y_t)
    else:
        return 1 / (h * alpha) * (v_t[period - 1] - (1 + b * alpha) * y_t)
    
def SRAS_curves_func(y_t,period,
               gamma = par['gamma'],
               phi = par['phi'],
               s_t = s_t): 

    if period == 0: 
        return pi_0 + gamma * y_t - phi * y_1 + s_t[period] - phi * s_0
    else:
        y_values = [0]
        pi_values = [0]
        for i in range(1,period): 
            y_val_i = yt_sol_func(
            v_t=v_t[i],
            s_t=s_t[i],
            pi_t_1=pi_vals[i-1],
            y_t_1=y_vals[i-1],
            s_t_1=s_t[i-1]
        )
            y_values.append(y_val_i)
            pi_val_i = pi_sol_func(
            v_t = v_t[i],
            s_t = s_t[i],
            y_t = y_vals[i],
            pi_t_1 = pi_vals[i-1],
            y_t_1 = y_vals[i-1],
            s_t_1 = s_t[i-1])
            pi_values.append(pi_val_i)
            
        return pi_values[period-1] + gamma * y_t - phi * gamma * y_values[period - 1] + s_t[period] - phi * s_t[period-1] 
    
```


```python
#  Constructing graph

def simulate_economy(vt,st): 
    
    # range for output, y
    
    y = np.linspace(-0.05,0.05,100)
    
    # Create AD-curves
    AD_0 = AD_curves_func(y, 0, v_t = v_t)
    AD_1 = AD_curves_func(y, 1, v_t = v_t)
    AD_2 = AD_curves_func(y, 2, v_t = v_t)
    AD_5 = AD_curves_func(y, 5, v_t = v_t)
    AD_10 = AD_curves_func(y, 10, v_t = v_t)
    AD_100 = AD_curves_func(y, 100, v_t = v_t)
    
    # Create SRAS-curves
    SRAS_0 = SRAS_curves_func(y, 0, s_t = s_t)
    SRAS_1 = SRAS_curves_func(y, 1, s_t = s_t)
    SRAS_2 = SRAS_curves_func(y, 2, s_t = s_t)
    SRAS_5 = SRAS_curves_func(y, 5, s_t = s_t)
    SRAS_10 = SRAS_curves_func(y, 10, s_t = s_t)
    SRAS_100 = SRAS_curves_func(y, 10 , s_t = s_t)
    
    
    # Plot the figure
    plt.figure(figsize=(12,6))
    
    # Plot AD curves
    plt.plot(y,AD_0,label=('$AD_0$'),linewidth=5,color='#1843ce')
    plt.plot(y,AD_1,label='$AD_1$',color='#1863ce')
    plt.plot(y,AD_2,label='$AD_2$',color='#1879ce')
    plt.plot(y,AD_5,label='$AD_5$',color='#189dce')
    plt.plot(y,AD_10,label='$AD_{10}$',color='#18b5ce')
    plt.plot(y,AD_100,label='$AD_{100}$',color='#38c0e2')
    
    # Plot SRAS curves
    plt.plot(y,SRAS_0,label='$SRAS_0$',linewidth=5,color='#961618')
    plt.plot(y,SRAS_1,label='$SRAS_1$',color='#ba1a1c')
    plt.plot(y,SRAS_2,label='$SRAS_2$',color='#d62022')
    plt.plot(y,SRAS_5,label='$SRAS_{5}$',color='#f92527')
    plt.plot(y,SRAS_10,label='$SRAS_{10}$',color='#f94345')
    plt.plot(y,SRAS_100,label='$SRAS_{100}$',color='#fc585a')
    
    # label, title 
    plt.title('AD and SRAS', fontsize = 15)
    plt.xlabel('$y_t$', fontsize = 15)
    plt.ylabel('$\pi_t$', fontsize = 15)

    # d. Create a dashed lines to show optimal values for y_t and pi_t 
    plt.axvline(yt_sol_val, linestyle='dashed', label='$y^*$')
    plt.axhline(pi_sol_val, linestyle='dashed', label='$\pi^*$')
    
    plt.legend()
    
    # Creates an arrow, which point at the way the AD-curve transition back to equilibrium.
    plt.annotate('', xy=(0, 0), xytext=(0.0, 0.035),
            arrowprops=dict(facecolor='black', shrink=0.01))
    
simulate_economy(v_t,s_t)
```

# Question 4
For this question we choose a random seed numer, 189, to generate $x_t$ and $c_t$. 


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


```python
# Choose 1000 periods and set seed number to 189, then generate x_t and c_t 
period = 1000
seed = 189
np.random.seed(seed)
x_t = np.random.normal(0,par['sigma_x'],period)
c_t = np.random.normal(0,par['sigma_c'],period)
```

Next we will create an empty array for the two shock variables and fill them up as done previously.


```python
# Create empty vectors
v_vec = [0]
s_vec = [0]


# Loop over each x_t and c_t
for i in range(1, period):
    v_i = par["delta"] * v_vec[i-1] + x_t[i]
    v_vec.append(v_i)
    s_i = par["omega"] * s_vec[i-1] + c_t[i]
    s_vec.append(s_i)
```

Next we want to find the new numeric solution to arrays for $y_t$ and $\pi_t$


```python
# 
y_vec = [0]
pi_vec = [0]

# Loop over each period and find the numeric solutions.
for i in range(1, period): 
    
    y_val_i = yt_sol_func(
        v_t=v_vec[i],
        s_t=s_vec[i],
        pi_t_1=pi_vec[i-1],
        y_t_1=y_vec[i-1],
        s_t_1=s_vec[i-1]
        )
    
    y_vec.append(y_val_i)
    
    pi_val_i = pi_sol_func(
        pi_t_1=pi_vec[i-1],
        v_t=v_vec[i],
        y_t_1=y_vec[i-1],
        s_t=s_vec[i],
        s_t_1=s_vec[i-1],
        y_t=y_vec[i],
        )
    pi_vec.append(pi_val_i)
```

Now we can find the mean, variance, correlation etc.


```python
# Variance for y
y_variance = np.var(y_vec)
f'1, The variance of y_t is: {y_variance:.3f}'
```




    '1, The variance of y_t is: 1.993'




```python
# Variance for pi
pi_variance = np.var(pi_vec)
f'2. The variance of pi_t is: {pi_variance:.3f}'
```




    '2. The variance of pi_t is: 1.206'




```python
# Correlation between y and pi
y_pi_corr = np.corrcoef(y_vec,pi_vec)
f'3. The correlation between y_t and pi_t is: {y_pi_corr[1,0]:.3f}'
```




    '3. The correlation between y_t and pi_t is: -0.234'




```python
# Auto-correlation between y_t_1 and y_t
autocoef_y = np.corrcoef(np.array([y_vec[:-1],y_vec[1:]]))
f'4. The auto-correlation between y_t and y_t_1 is: {autocoef_y[1,0]:.3f}'
```




    '4. The auto-correlation between y_t and y_t_1 is: 0.791'




```python
# Auto-correlation between pi_t_1 and pi_t
autocoef_pi = np.corrcoef(np.array([pi_vec[:-1],pi_vec[1:]]))
f'5. The auto-correlation between pi_t and pi_t_1 is: {autocoef_pi[1,0]:.3f}'
```




    '5. The auto-correlation between pi_t and pi_t_1 is: 0.980'



# Question 5 
For this question, we will begin by resetting the variables by setting them equal to their initial values. Then we will define a function, which first create two empty arrays and then loop over each period to fill them up with the numeric solutions. The function will return the correlation coeffecient between $y_t$ and $\pi_t$


```python
# Set variables to their initial values
s_0 = 0
v_0 = 0
x_0 = 0.1
x_t = 0
c_t = 0

# Define function
def corr_sim(corr):
    
    # Create arrays
    y_sim_vec = [0]
    pi_sim_vec = [0]
    
    
    # Loop over each period and find numerical solutions.
    for i in range (1, period):
        
        y_sim = yt_sol_func(
        v_t=v_vec[i],
        s_t=s_vec[i],
        pi_t_1=pi_sim_vec[i-1],
        y_t_1=y_sim_vec[i-1],
        s_t_1=s_vec[i-1],
        phi = corr
        )
        y_sim_vec.append(y_sim)
        
        
        pi_sim = pi_sol_func(
        pi_t_1=pi_sim_vec[i-1],
        v_t=v_vec[i],
        y_t_1=y_sim_vec[i-1],
        s_t=s_vec[i],
        s_t_1=s_vec[i-1],
        y_t=y_sim_vec[i],
        phi = corr
        )
        pi_sim_vec.append(pi_sim)
        
    # Return the correlation coeffecient
    return np.corrcoef(y_sim_vec,pi_sim_vec)[1,0]
```

Now we will plot how the correlation increases for increasing values of $\phi$.


```python
# Create vec for phi
phi_vec = np.linspace(0,1,100)

# Create empty vec for correlation and loop over each correlation between y_t and pi_t
correlation_vec = np.empty(100)
for i, val in enumerate(phi_vec):
    correlation_vec[i] = corr_sim(val)

# plot
plt.figure(figsize=(8,5))
plt.plot(phi_vec,correlation_vec)
plt.title('Correlation between $y_t$ and $\pi_t$')
plt.xlabel('$\phi$')
plt.ylabel('Correlation coef between $y_t$ and $\pi_t$')
plt.grid(True)
```

We see from the above plot how the correlation increases for increasing values of $\phi$
Now we will make use of scipy.optimize to find the value of $\phi$, which ensure that the correlation coef between $y_t$ and $pi_t$ is 0.31


```python
# import optimize
from scipy import optimize
# Set x_0 = 0
x_0 = 0

# Create a lambda function
obj = lambda corr: (corr_sim(corr)-0.31)**2

# minimize_scalar with the bounded method
result = optimize.minimize_scalar(obj,x_0,method='bounded',bounds=[0.5,1])
print(f'"The optimal value of phi is {result.x:.3} to ensure the correlation coef between y_t and pi_t is 0.31')
```

    "The optimal value of phi is 0.979 to ensure the correlation coef between y_t and pi_t is 0.31


# Quesiton 6
Unfortunately, we ran into some issues when trying to solve the last question. 
We tried to enforce the similar method as in the previous question, but was unable to do it correctly.
Anyways, we notice that the $var(y_t)$, $corr(y_t,\pi_t$) and $corr(y_{t-1},y_t)$ somewhat matches the US economy, while the others does not.

# 3. Exchange economy


```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)
```

# Question 1
As the utility function is presented as a Cobb-douglas, we know that $beta_j$ is the budget share for each good.


```python
# Create arrays for each agent's preferences of the three goods
beta1=betas[:,0]
beta2=betas[:,1]
beta3=betas[:,2]
```


```python
# Create figure
fig = plt.figure(figsize=(15,10))

# Add a subplot for each good
ax1 = fig.add_subplot(2,3,1)
ax1.hist(beta1, bins=200, color ="#f7d56f")
ax1.set_title('Good #1')
ax1.set_xlabel('$\\beta_1$')
ax1.set_ylabel('Number of agents')
plt.axis([0, 1, 0, 1300])
        
ax2 = fig.add_subplot(2,3,2)
ax2.hist(beta2, bins=200, color ="#6ab2f2")
ax2.set_title('Good #2')
ax2.set_xlabel('$\\beta_2$')
plt.axis([0, 1, 0, 1300])

ax3 = fig.add_subplot(2,3,3)
ax3.hist(beta3, bins=200, color ="#66f47e")
ax3.set_title('Good #3')
ax3.set_xlabel('$\\beta_3$')
plt.axis([0, 1, 0, 1300])


plt.show()    
```

# Question 2
We will begin by defining the demand function for each of the three goods. 


```python
# The function takes the two price vectors p1, p2 and three endowments e1, e2, e3 and the betas.
# The income is defined as price times endowments
# The function returns the demand for each good 

def demand_good_1_func(p1,p2,e1,e2,e3,betas):
    I = p1*e1+p2*e2+e3
    return beta1*I/p1

def demand_good_2_func(p1,p2,e1,e2,e3,betas):
    I = p1 * e1 + p2 * e2 + e3
    return beta2*I/p2

def demand_good_3_func(p1,p2,e1,e2,e3,betas):
    I = p1 * e1 + p2 * e2 + e3
    return beta3*I
```

Then we will define the excess demand functions for each good.


```python
def excess_demand_good_1_func(p1,p2,e1,e2,e3,betas):

    # a. demand - Sum of all demand for good 1
    demand = np.sum(demand_good_1_func(p1,p2,e1,e2,e3,betas))
    
    # b. supply - Sum of all endowments
    supply = np.sum(e1)
    
    # c. excess demand - Demand less the supply
    excess_demand = demand-supply
    
    return excess_demand

def excess_demand_good_2_func(p1,p2,e1,e2,e3,betas):
    
    # a. demand - Sum of all demand for good 2
    demand = np.sum(demand_good_2_func(p1,p2,e1,e2,e3,betas))
    
    # b. upply - Sum of all endowments
    supply = np.sum(e2)
    
    # c. excess demand - Demand less the supply
    excess_demand = demand-supply
    
    return excess_demand


def excess_demand_good_3_func(p1,p2,e1,e2,e3,betas):
    
    # a. demand - Sum of all demand for good 3
    demand = np.sum(demand_good_3_func(p1,p2,e1,e2,e3,betas))
    
    # b. Sum of all endowments
    supply = np.sum(e3)
    
    # c. excess demand - Demand less the supply
    excess_demand = demand-supply
    
    return excess_demand
```

Next we will first define the two price vectors for p1 and p2, which we define as a linear vector from [0.5,10]. Then we will make use of np.meshgrid to return grids from the price vectors.
We will then create three empty grids for each excess grid and then fill them up with each value in the price vectors given the demand excess functions.


```python
# choose a size for the vectors
Num = 100


# Create shape tuple
shape_tuple = (Num,Num)

# Create price vectors
p1_vec = np.linspace(0.5,10,Num)
p2_vec = np.linspace(0.5,10,Num)


# Create meshgrid from price vectors
p1_grid, p2_grid = np.meshgrid(p1_vec, p2_vec) 


# Define three empty excess grids 
excess1 = np.empty((shape_tuple))
excess2 = np.empty((shape_tuple))
excess3 = np.empty((shape_tuple))


# Iterate over each value in the price vector to fill up the excess grids
for i, p1 in enumerate(p1_vec):
    for j, p2 in enumerate(p2_vec):
        excess1[i,j]=excess_demand_good_1_func(p1,p2,e1,e2,e3,betas)
        excess2[i,j]=excess_demand_good_2_func(p1,p2,e1,e2,e3,betas)
        excess3[i,j]=excess_demand_good_3_func(p1,p2,e1,e2,e3,betas)

```


```python
# Import packages
%matplotlib inline
import matplotlib.pyplot as plt # baseline modul
from mpl_toolkits.mplot3d import Axes3D # for 3d figures
from matplotlib import cm # for color
```


```python

fig = plt.figure(figsize=(15,5)) # create the figure

# For each good, plot the price meshes and the excess demand. 

ax1 = fig.add_subplot(1,3,1,projection='3d') # create a 3d axis in the figure 
ax1.plot_surface(p1_grid,p2_grid,excess1,cmap=cm.jet) # create surface plot in the axis
ax1.invert_xaxis()
ax1.set_title('Excess demand of good 1')
ax1.set_xlabel('$p_2$')
ax1.set_ylabel('$p_1$')
ax1.set_zlabel('$z_1$') 
ax1.xaxis.pane.fill = False
ax1.yaxis.pane.fill = False
ax1.zaxis.pane.fill = False

ax2 = fig.add_subplot(1,3,2,projection='3d') # create a 3d axis in the figure 
ax2.plot_surface(p1_grid,p2_grid,excess2,cmap=cm.jet)
ax2.invert_xaxis()
ax2.set_title('Excess demand of good 2')
ax2.set_xlabel('$p_2$')
ax2.set_ylabel('$p_1$')
ax2.set_zlabel('$z_2$'); 
ax2.xaxis.pane.fill = False
ax2.yaxis.pane.fill = False
ax2.zaxis.pane.fill = False

ax3 = fig.add_subplot(1,3,3,projection='3d') # create a 3d axis in the figure 
ax3.plot_surface(p1_grid,p2_grid,excess3,cmap=cm.jet)
ax3.invert_xaxis()
ax3.set_title('Excess demand of good 3')
ax3.set_xlabel('$p_2$')
ax3.set_ylabel('$p_1$')
ax3.set_zlabel('$z_3$')
ax3.xaxis.pane.fill = False
ax3.yaxis.pane.fill = False
ax3.zaxis.pane.fill = False;
```

# Question 3
We will begin this question by defining a function for the tâtonnement process, which takes p1, p2, e1, e2, e3, betas, kappa (Which adjust the aggresivivity), eps (the tolerance), maxiter (Max iterations).


```python
def find_equilibrium(p1,p2,e1,e2,e3,betas,kappa=0.9,eps=1e-8,maxiter=3000):
    
    t = 0
    while True:

        # a. excess demand for good 1 and good 2
        Z1 = excess_demand_good_1_func(p1,p2,e1,e2,e3,betas)
        Z2 = excess_demand_good_2_func(p1,p2,e1,e2,e3,betas)
        
        # b: If the excess demand for good 1 or good 2 is lower (or maxiter has been reached) than the tolerance
        #  then break.
        if  np.abs(Z1) < eps and np.abs(Z2) < eps or t >= maxiter:
            print(f'{t:3d}: p1 = {p1:12.8f} -> excess demand -> {Z1:14.8f} p2 = {p2:12.8f} -> excess demand -> {Z2:14.8f}')
            break    
    
        # c. If not, then set the prices as below
        p1 = p1 + kappa*Z1/betas.size
        p2 = p2 + kappa*Z2/betas.size
        
        # d. return 
        if t < 5 or t%25 == 0:
            print(f'{t:3d}: p1 = {p1:12.8f} -> excess demand -> {Z1:14.8f} p2 = {p2:12.8f} -> excess demand -> {Z2:14.8f}')
        elif t == 5:
            print('   ...')
            
        t += 1    

    return p1, p2
```

Now we guess on the prices and run the function.


```python
p1 = 3
p2 = 3
kappa = 0.8
eps = 1e-8
p1, p2 = find_equilibrium(p1,p2,e1,e2,e3,betas,kappa=kappa,eps=eps)
```

      0: p1 =   3.13242998 -> excess demand -> 24830.62081187 p2 =   2.89503410 -> excess demand -> -19681.10609181
      1: p1 =   3.24953265 -> excess demand -> 21956.75110613 p2 =   2.79653311 -> excess demand -> -18468.93528892
      2: p1 =   3.35382488 -> excess demand -> 19554.79244949 p2 =   2.70436429 -> excess demand -> -17281.65490885
      3: p1 =   3.44723557 -> excess demand -> 17514.50532486 p2 =   2.61840426 -> excess demand -> -16117.50536588
      4: p1 =   3.53129068 -> excess demand -> 15760.33324137 p2 =   2.53852765 -> excess demand -> -14976.86431327
       ...
     25: p1 =   4.32935739 -> excess demand ->  3668.14907331 p2 =   1.88994781 -> excess demand ->  -737.37747432
     50: p1 =   4.69763396 -> excess demand ->  2299.14161154 p2 =   1.95010213 -> excess demand ->   778.99889458
     75: p1 =   4.96884180 -> excess demand ->  1824.74465940 p2 =   2.04851523 -> excess demand ->   679.27508163
    100: p1 =   5.18816891 -> excess demand ->  1494.51036578 p2 =   2.13043481 -> excess demand ->   558.86386466
    125: p1 =   5.36929385 -> excess demand ->  1243.02199491 p2 =   2.19816371 -> excess demand ->   464.74526696
    150: p1 =   5.52082611 -> excess demand ->  1045.38867190 p2 =   2.25480861 -> excess demand ->   390.71295581
    175: p1 =   5.64883988 -> excess demand ->   886.71796026 p2 =   2.30264612 -> excess demand ->   331.31089719
    200: p1 =   5.75780927 -> excess demand ->   757.22505816 p2 =   2.34335602 -> excess demand ->   282.85957009
    225: p1 =   5.85113097 -> excess demand ->   650.16880321 p2 =   2.37821255 -> excess demand ->   242.82156882
    250: p1 =   5.93144568 -> excess demand ->   560.73558333 p2 =   2.40820553 -> excess demand ->   209.38689340
    275: p1 =   6.00084623 -> excess demand ->   485.38651202 p2 =   2.43411891 -> excess demand ->   181.22623081
    300: p1 =   6.06101767 -> excess demand ->   421.45595138 p2 =   2.45658349 -> excess demand ->   157.33916882
    325: p1 =   6.11333460 -> excess demand ->   366.89448568 p2 =   2.47611363 -> excess demand ->   136.95713969
    350: p1 =   6.15893085 -> excess demand ->   320.09870835 p2 =   2.49313343 -> excess demand ->   119.47922356
    375: p1 =   6.19875040 -> excess demand ->   279.79517965 p2 =   2.50799584 -> excess demand ->   104.42842920
    400: p1 =   6.23358544 -> excess demand ->   244.95929392 p2 =   2.52099698 -> excess demand ->    91.42115569
    425: p1 =   6.26410534 -> excess demand ->   214.75728406 p2 =   2.53238699 -> excess demand ->    80.14538099
    450: p1 =   6.29087897 -> excess demand ->   188.50394327 p2 =   2.54237843 -> excess demand ->    70.34477323
    475: p1 =   6.31439229 -> excess demand ->   165.63125967 p2 =   2.55115281 -> excess demand ->    61.80691099
    500: p1 =   6.33506222 -> excess demand ->   145.66477716 p2 =   2.55886586 -> excess demand ->    54.35441059
    525: p1 =   6.35324786 -> excess demand ->   128.20552356 p2 =   2.56565168 -> excess demand ->    47.83814661
    550: p1 =   6.36925947 -> excess demand ->   112.91601469 p2 =   2.57162612 -> excess demand ->    42.13200407
    575: p1 =   6.38336596 -> excess demand ->    99.50928598 p2 =   2.57688956 -> excess demand ->    37.12876764
    600: p1 =   6.39580094 -> excess demand ->    87.74020338 p2 =   2.58152924 -> excess demand ->    32.73686646
    625: p1 =   6.40676783 -> excess demand ->    77.39851106 p2 =   2.58562108 -> excess demand ->    28.87777051
    650: p1 =   6.41644411 -> excess demand ->    68.30321767 p2 =   2.58923132 -> excess demand ->    25.48388912
    675: p1 =   6.42498488 -> excess demand ->    60.29802549 p2 =   2.59241785 -> excess demand ->    22.49686034
    700: p1 =   6.43252588 -> excess demand ->    53.24757962 p2 =   2.59523134 -> excess demand ->    19.86614766
    725: p1 =   6.43918608 -> excess demand ->    47.03436835 p2 =   2.59771618 -> excess demand ->    17.54788068
    750: p1 =   6.44506987 -> excess demand ->    41.55614478 p2 =   2.59991133 -> excess demand ->    15.50389093
    775: p1 =   6.45026893 -> excess demand ->    36.72376936 p2 =   2.60185101 -> excess demand ->    13.70090532
    800: p1 =   6.45486387 -> excess demand ->    32.45939446 p2 =   2.60356528 -> excess demand ->    12.10986759
    825: p1 =   6.45892559 -> excess demand ->    28.69492934 p2 =   2.60508061 -> excess demand ->    10.70536464
    850: p1 =   6.46251652 -> excess demand ->    25.37073606 p2 =   2.60642029 -> excess demand ->     9.46513932
    875: p1 =   6.46569167 -> excess demand ->    22.43451708 p2 =   2.60760485 -> excess demand ->     8.36967481
    900: p1 =   6.46849952 -> excess demand ->    19.84036280 p2 =   2.60865238 -> excess demand ->     7.40183888
    925: p1 =   6.47098281 -> excess demand ->    17.54793309 p2 =   2.60957882 -> excess demand ->     6.54657822
    950: p1 =   6.47317928 -> excess demand ->    15.52175195 p2 =   2.61039825 -> excess demand ->     5.79065504
    975: p1 =   6.47512222 -> excess demand ->    13.73059782 p2 =   2.61112309 -> excess demand ->     5.12241938
    1000: p1 =   6.47684100 -> excess demand ->    12.14697534 p2 =   2.61176431 -> excess demand ->     4.53161194
    1025: p1 =   6.47836160 -> excess demand ->    10.74665649 p2 =   2.61233160 -> excess demand ->     4.00919277
    1050: p1 =   6.47970694 -> excess demand ->     9.50828140 p2 =   2.61283349 -> excess demand ->     3.54719219
    1075: p1 =   6.48089728 -> excess demand ->     8.41301014 p2 =   2.61327757 -> excess demand ->     3.13858092
    1100: p1 =   6.48195053 -> excess demand ->     7.44421877 p2 =   2.61367049 -> excess demand ->     2.77715653
    1125: p1 =   6.48288251 -> excess demand ->     6.58723343 p2 =   2.61401818 -> excess demand ->     2.45744428
    1150: p1 =   6.48370721 -> excess demand ->     5.82909747 p2 =   2.61432584 -> excess demand ->     2.17461013
    1175: p1 =   6.48443701 -> excess demand ->     5.15836733 p2 =   2.61459810 -> excess demand ->     1.92438463
    1200: p1 =   6.48508284 -> excess demand ->     4.56493326 p2 =   2.61483904 -> excess demand ->     1.70299599
    1225: p1 =   6.48565438 -> excess demand ->     4.03986196 p2 =   2.61505226 -> excess demand ->     1.50711138
    1250: p1 =   6.48616019 -> excess demand ->     3.57525812 p2 =   2.61524095 -> excess demand ->     1.33378520
    1275: p1 =   6.48660783 -> excess demand ->     3.16414260 p2 =   2.61540795 -> excess demand ->     1.18041372
    1300: p1 =   6.48700400 -> excess demand ->     2.80034519 p2 =   2.61555574 -> excess demand ->     1.04469500
    1325: p1 =   6.48735462 -> excess demand ->     2.47841012 p2 =   2.61568655 -> excess demand ->     0.92459362
    1350: p1 =   6.48766493 -> excess demand ->     2.19351275 p2 =   2.61580231 -> excess demand ->     0.81830966
    1375: p1 =   6.48793958 -> excess demand ->     1.94138610 p2 =   2.61590477 -> excess demand ->     0.72425125
    1400: p1 =   6.48818266 -> excess demand ->     1.71825606 p2 =   2.61599545 -> excess demand ->     0.64101038
    1425: p1 =   6.48839780 -> excess demand ->     1.52078417 p2 =   2.61607571 -> excess demand ->     0.56734159
    1450: p1 =   6.48858822 -> excess demand ->     1.34601712 p2 =   2.61614675 -> excess demand ->     0.50214310
    1475: p1 =   6.48875675 -> excess demand ->     1.19134215 p2 =   2.61620963 -> excess demand ->     0.44444018
    1500: p1 =   6.48890592 -> excess demand ->     1.05444763 p2 =   2.61626527 -> excess demand ->     0.39337044
    1525: p1 =   6.48903795 -> excess demand ->     0.93328828 p2 =   2.61631453 -> excess demand ->     0.34817087
    1550: p1 =   6.48915481 -> excess demand ->     0.82605438 p2 =   2.61635812 -> excess demand ->     0.30816633
    1575: p1 =   6.48925824 -> excess demand ->     0.73114456 p2 =   2.61639671 -> excess demand ->     0.27275941
    1600: p1 =   6.48934979 -> excess demand ->     0.64714181 p2 =   2.61643086 -> excess demand ->     0.24142146
    1625: p1 =   6.48943082 -> excess demand ->     0.57279216 p2 =   2.61646109 -> excess demand ->     0.21368470
    1650: p1 =   6.48950254 -> excess demand ->     0.50698594 p2 =   2.61648785 -> excess demand ->     0.18913515
    1675: p1 =   6.48956602 -> excess demand ->     0.44874112 p2 =   2.61651153 -> excess demand ->     0.16740643
    1700: p1 =   6.48962221 -> excess demand ->     0.39718861 p2 =   2.61653249 -> excess demand ->     0.14817435
    1725: p1 =   6.48967194 -> excess demand ->     0.35155929 p2 =   2.61655104 -> excess demand ->     0.13115196
    1750: p1 =   6.48971596 -> excess demand ->     0.31117244 p2 =   2.61656747 -> excess demand ->     0.11608532
    1775: p1 =   6.48975493 -> excess demand ->     0.27542563 p2 =   2.61658200 -> excess demand ->     0.10274969
    1800: p1 =   6.48978941 -> excess demand ->     0.24378567 p2 =   2.61659487 -> excess demand ->     0.09094615
    1825: p1 =   6.48981994 -> excess demand ->     0.21578066 p2 =   2.61660625 -> excess demand ->     0.08049866
    1850: p1 =   6.48984696 -> excess demand ->     0.19099295 p2 =   2.61661633 -> excess demand ->     0.07125141
    1875: p1 =   6.48987087 -> excess demand ->     0.16905288 p2 =   2.61662526 -> excess demand ->     0.06306649
    1900: p1 =   6.48989204 -> excess demand ->     0.14963327 p2 =   2.61663315 -> excess demand ->     0.05582185
    1925: p1 =   6.48991077 -> excess demand ->     0.13244455 p2 =   2.61664014 -> excess demand ->     0.04940947
    1950: p1 =   6.48992736 -> excess demand ->     0.11723042 p2 =   2.61664633 -> excess demand ->     0.04373371
    1975: p1 =   6.48994204 -> excess demand ->     0.10376402 p2 =   2.61665180 -> excess demand ->     0.03870997
    2000: p1 =   6.48995503 -> excess demand ->     0.09184457 p2 =   2.61665665 -> excess demand ->     0.03426332
    2025: p1 =   6.48996653 -> excess demand ->     0.08129435 p2 =   2.61666094 -> excess demand ->     0.03032748
    2050: p1 =   6.48997671 -> excess demand ->     0.07195607 p2 =   2.61666474 -> excess demand ->     0.02684377
    2075: p1 =   6.48998572 -> excess demand ->     0.06369050 p2 =   2.61666810 -> excess demand ->     0.02376023
    2100: p1 =   6.48999369 -> excess demand ->     0.05637441 p2 =   2.61667108 -> excess demand ->     0.02103091
    2125: p1 =   6.49000075 -> excess demand ->     0.04989873 p2 =   2.61667371 -> excess demand ->     0.01861511
    2150: p1 =   6.49000700 -> excess demand ->     0.04416692 p2 =   2.61667604 -> excess demand ->     0.01647681
    2175: p1 =   6.49001253 -> excess demand ->     0.03909352 p2 =   2.61667810 -> excess demand ->     0.01458414
    2200: p1 =   6.49001743 -> excess demand ->     0.03460290 p2 =   2.61667993 -> excess demand ->     0.01290888
    2225: p1 =   6.49002176 -> excess demand ->     0.03062812 p2 =   2.61668155 -> excess demand ->     0.01142606
    2250: p1 =   6.49002559 -> excess demand ->     0.02710992 p2 =   2.61668298 -> excess demand ->     0.01011357
    2275: p1 =   6.49002899 -> excess demand ->     0.02399586 p2 =   2.61668424 -> excess demand ->     0.00895184
    2300: p1 =   6.49003199 -> excess demand ->     0.02123950 p2 =   2.61668536 -> excess demand ->     0.00792356
    2325: p1 =   6.49003465 -> excess demand ->     0.01879976 p2 =   2.61668636 -> excess demand ->     0.00701340
    2350: p1 =   6.49003701 -> excess demand ->     0.01664027 p2 =   2.61668723 -> excess demand ->     0.00620778
    2375: p1 =   6.49003909 -> excess demand ->     0.01472884 p2 =   2.61668801 -> excess demand ->     0.00549471
    2400: p1 =   6.49004093 -> excess demand ->     0.01303697 p2 =   2.61668870 -> excess demand ->     0.00486354
    2425: p1 =   6.49004257 -> excess demand ->     0.01153945 p2 =   2.61668931 -> excess demand ->     0.00430488
    2450: p1 =   6.49004401 -> excess demand ->     0.01021394 p2 =   2.61668985 -> excess demand ->     0.00381039
    2475: p1 =   6.49004529 -> excess demand ->     0.00904069 p2 =   2.61669032 -> excess demand ->     0.00337270
    2500: p1 =   6.49004642 -> excess demand ->     0.00800221 p2 =   2.61669075 -> excess demand ->     0.00298528
    2525: p1 =   6.49004742 -> excess demand ->     0.00708301 p2 =   2.61669112 -> excess demand ->     0.00264237
    2550: p1 =   6.49004831 -> excess demand ->     0.00626941 p2 =   2.61669145 -> excess demand ->     0.00233885
    2575: p1 =   6.49004910 -> excess demand ->     0.00554926 p2 =   2.61669174 -> excess demand ->     0.00207019
    2600: p1 =   6.49004979 -> excess demand ->     0.00491183 p2 =   2.61669200 -> excess demand ->     0.00183239
    2625: p1 =   6.49005041 -> excess demand ->     0.00434762 p2 =   2.61669223 -> excess demand ->     0.00162191
    2650: p1 =   6.49005095 -> excess demand ->     0.00384822 p2 =   2.61669244 -> excess demand ->     0.00143561
    2675: p1 =   6.49005143 -> excess demand ->     0.00340618 p2 =   2.61669262 -> excess demand ->     0.00127070
    2700: p1 =   6.49005186 -> excess demand ->     0.00301492 p2 =   2.61669277 -> excess demand ->     0.00112474
    2725: p1 =   6.49005224 -> excess demand ->     0.00266861 p2 =   2.61669292 -> excess demand ->     0.00099554
    2750: p1 =   6.49005257 -> excess demand ->     0.00236207 p2 =   2.61669304 -> excess demand ->     0.00088119
    2775: p1 =   6.49005287 -> excess demand ->     0.00209075 p2 =   2.61669315 -> excess demand ->     0.00077997
    2800: p1 =   6.49005313 -> excess demand ->     0.00185059 p2 =   2.61669325 -> excess demand ->     0.00069038
    2825: p1 =   6.49005336 -> excess demand ->     0.00163802 p2 =   2.61669333 -> excess demand ->     0.00061107
    2850: p1 =   6.49005357 -> excess demand ->     0.00144986 p2 =   2.61669341 -> excess demand ->     0.00054088
    2875: p1 =   6.49005375 -> excess demand ->     0.00128332 p2 =   2.61669348 -> excess demand ->     0.00047875
    2900: p1 =   6.49005391 -> excess demand ->     0.00113591 p2 =   2.61669354 -> excess demand ->     0.00042376
    2925: p1 =   6.49005405 -> excess demand ->     0.00100543 p2 =   2.61669359 -> excess demand ->     0.00037508
    2950: p1 =   6.49005418 -> excess demand ->     0.00088994 p2 =   2.61669364 -> excess demand ->     0.00033200
    2975: p1 =   6.49005429 -> excess demand ->     0.00078771 p2 =   2.61669368 -> excess demand ->     0.00029386
    3000: p1 =   6.49005438 -> excess demand ->     0.00069723 p2 =   2.61669372 -> excess demand ->     0.00026011


The Walras-equilibrium prices, $(p_1,p_2)$, where both excess demands are (approximately) zero have been found to be: 


```python
f'(p_1,p_2) = ({p1:.3f},{p2:.3f})'
```




    '(p_1,p_2) = (6.490,2.617)'



At last we will ensure that the excess demand is almost zero.


```python
Z1 = excess_demand_good_1_func(p1,p2,e1,e2,e3,betas)
Z2 = excess_demand_good_2_func(p1,p2,e1,e2,e3,betas)
Z3 = excess_demand_good_3_func(p1,p2,e1,e2,e3,betas)
f'Z1: {Z1:.2f}, Z2: {Z2:.2f}, Z3: {Z3:.2f}'
```




    'Z1: 0.00, Z2: 0.00, Z3: -0.01'



# Question 4
We define the utility function as presented in the text. 


```python
# The function takes the two prices, p1, p2 and the endowments, e1, e2, e3, betas and gamma.
# The income with the new prices
# The demand for each good given the new prices
# return the utility

def utility_func(p1, p2, e1, e2, e3, betas, gamma):
    
    I = p1 * e1 + p2 * e2 + e3
    x1 = beta1*I/p1
    x2 = beta2*I/p2
    x3 = beta3*I
    
    
    return (x1 ** beta1 * x2 ** beta2 * x3 ** beta3)**gamma
```

Next we calculate the mean and variance for the utility of the agents


```python
utility_vec = utility_func(p1, p2, e1, e2, e3, betas, gamma)

# Calculate mean
utility_mean = np.mean(utility_vec)

# Calculate variance
utility_var = np.var(utility_vec)

f'The distribution has the mean of {utility_mean:.3f} and the variance of {utility_var:.3f}'

```




    'The distribution has the mean of 1.010 and the variance of 0.317'



Now we will plot the distribution of utility in a histogram


```python
# Plotting the figure
fig = plt.figure(dpi=120)

ax = fig.add_subplot(1,1,1)
ax.hist(utility_vec,bins=500,density=True, alpha = 0.4, label='beta1')
ax.set_title('Distribution of utilities')
ax.set_xlabel('Utility')
ax.set_ylabel('Number of agents')
ax.set_xlim([0,6])
ax.text(3, 0.5, f'$\mu_u$: {utility_mean:.3f},  $\sigma^2_u$: {utility_var:.3f}');
```

# Question 5
We start this question by creating a new empty array for each good and then filling it up with the mean value of each endowment.


```python
# Create empty array and fill up with mean value of endowment
e1_equal = np.zeros(N)+np.mean(e1)
e2_equal = np.zeros(N)+np.mean(e2)
e3_equal = np.zeros(N)+np.mean(e3)
```

Then we can run our find_equilibrium function one more time with the endowment arrays.


```python
# Guess on prices
p1 = 3
p2 = 3

# Run function
p1, p2 = find_equilibrium(p1,p2,e1_equal,e2_equal,e3_equal,betas,kappa=kappa,eps=eps)
```

      0: p1 =   3.13240299 -> excess demand -> 24825.55975897 p2 =   2.89504860 -> excess demand -> -19678.38813939
      1: p1 =   3.24947755 -> excess demand -> 21951.48082735 p2 =   2.79656676 -> excess demand -> -18465.34392437
      2: p1 =   3.35374115 -> excess demand -> 19549.42451647 p2 =   2.70442119 -> excess demand -> -17277.29411287
      3: p1 =   3.44712313 -> excess demand -> 17509.12132286 p2 =   2.61848798 -> excess demand -> -16112.47732727
      4: p1 =   3.53114976 -> excess demand -> 15754.99444829 p2 =   2.53864120 -> excess demand -> -14971.27125771
       ...
     25: p1 =   4.32881196 -> excess demand ->  3665.82708308 p2 =   1.89065163 -> excess demand ->  -734.78080912
     50: p1 =   4.69682912 -> excess demand ->  2297.33952155 p2 =   1.95095356 -> excess demand ->   779.53750330
     75: p1 =   4.96780468 -> excess demand ->  1823.05317155 p2 =   2.04941887 -> excess demand ->   679.57023695
    100: p1 =   5.18691322 -> excess demand ->  1492.91876997 p2 =   2.13136828 -> excess demand ->   559.02987191
    125: p1 =   5.36783268 -> excess demand ->  1241.52683984 p2 =   2.19911271 -> excess demand ->   464.82072065
    150: p1 =   5.51917205 -> excess demand ->  1043.98617748 p2 =   2.25576291 -> excess demand ->   390.72314748
    175: p1 =   5.64700501 -> excess demand ->   885.40419047 p2 =   2.30359834 -> excess demand ->   331.27379011
    200: p1 =   5.75580516 -> excess demand ->   755.99602954 p2 =   2.34430079 -> excess demand ->   282.78830357
    225: p1 =   5.84896863 -> excess demand ->   649.02053716 p2 =   2.37914604 -> excess demand ->   242.72594953
    250: p1 =   5.92913560 -> excess demand ->   559.66413687 p2 =   2.40912503 -> excess demand ->   209.27434732
    275: p1 =   5.99839840 -> excess demand ->   484.38799846 p2 =   2.43502256 -> excess demand ->   181.10243703
    300: p1 =   6.05844153 -> excess demand ->   420.52655557 p2 =   2.45747012 -> excess demand ->   157.20849654
    325: p1 =   6.11063911 -> excess demand ->   366.03047617 p2 =   2.47698256 -> excess demand ->   136.82295853
    350: p1 =   6.15612449 -> excess demand ->   319.29644750 p2 =   2.49398437 -> excess demand ->   119.34412945
    375: p1 =   6.19584114 -> excess demand ->   279.05113221 p2 =   2.50882883 -> excess demand ->   104.29441254
    400: p1 =   6.23058082 -> excess demand ->   244.27003433 p2 =   2.52181226 -> excess demand ->    91.28972878
    425: p1 =   6.26101242 -> excess demand ->   214.11950259 p2 =   2.53318500 -> excess demand ->    80.01767626
    450: p1 =   6.28770439 -> excess demand ->   187.91445113 p2 =   2.54315972 -> excess demand ->    70.22161994
    475: p1 =   6.31114228 -> excess demand ->   165.08699297 p2 =   2.55191803 -> excess demand ->    61.68889594
    500: p1 =   6.33174261 -> excess demand ->   145.16279991 p2 =   2.55961573 -> excess demand ->    54.24192648
    525: p1 =   6.34986409 -> excess demand ->   127.74302968 p2 =   2.56638695 -> excess demand ->    47.73143108
    550: p1 =   6.36581661 -> excess demand ->   112.49032901 p2 =   2.57234756 -> excess demand ->    42.03117126
    575: p1 =   6.37986874 -> excess demand ->    99.11786447 p2 =   2.57759797 -> excess demand ->    37.03383396
    600: p1 =   6.39225377 -> excess demand ->    87.38063248 p2 =   2.58222540 -> excess demand ->    32.64777158
    625: p1 =   6.40317479 -> excess demand ->    77.06850625 p2 =   2.58630575 -> excess demand ->    28.79439463
    650: p1 =   6.41280900 -> excess demand ->    68.00062145 p2 =   2.58990526 -> excess demand ->    25.40606713
    675: p1 =   6.42131120 -> excess demand ->    60.02080470 p2 =   2.59308179 -> excess demand ->    22.42439343
    700: p1 =   6.42881688 -> excess demand ->    52.99382226 p2 =   2.59588598 -> excess demand ->    19.79881281
    725: p1 =   6.43544477 -> excess demand ->    46.80227995 p2 =   2.59836219 -> excess demand ->    17.48543844
    750: p1 =   6.44129902 -> excess demand ->    41.34404447 p2 =   2.60054934 -> excess demand ->    15.44609178
    775: p1 =   6.44647110 -> excess demand ->    36.53008556 p2 =   2.60248161 -> excess demand ->    13.64749474
    800: p1 =   6.45104140 -> excess demand ->    32.28266044 p2 =   2.60418905 -> excess demand ->    12.06059022
    825: p1 =   6.45508065 -> excess demand ->    28.53377851 p2 =   2.60569808 -> excess demand ->    10.65996752
    850: p1 =   6.45865110 -> excess demand ->    25.22389729 p2 =   2.60703197 -> excess demand ->     9.42337443
    875: p1 =   6.46180759 -> excess demand ->    22.30080988 p2 =   2.60821120 -> excess demand ->     8.33130112
    900: p1 =   6.46459846 -> excess demand ->    19.71869251 p2 =   2.60925383 -> excess demand ->     7.36662376
    925: p1 =   6.46706631 -> excess demand ->    17.43728612 p2 =   2.61017578 -> excess demand ->     6.51429858
    950: p1 =   6.46924874 -> excess demand ->    15.42119109 p2 =   2.61099110 -> excess demand ->     5.76109802
    975: p1 =   6.47117891 -> excess demand ->    13.63925765 p2 =   2.61171218 -> excess demand ->     5.09538288
    1000: p1 =   6.47288611 -> excess demand ->    12.06405771 p2 =   2.61234996 -> excess demand ->     4.50690482
    1025: p1 =   6.47439620 -> excess demand ->    10.67142625 p2 =   2.61291410 -> excess demand ->     3.98663493
    1050: p1 =   6.47573200 -> excess demand ->     9.44006218 p2 =   2.61341312 -> excess demand ->     3.52661449
    1075: p1 =   6.47691369 -> excess demand ->     8.35118038 p2 =   2.61385458 -> excess demand ->     3.11982498
    1100: p1 =   6.47795911 -> excess demand ->     7.38820788 p2 =   2.61424512 -> excess demand ->     2.76007444
    1125: p1 =   6.47888399 -> excess demand ->     6.53651815 p2 =   2.61459064 -> excess demand ->     2.44189824
    1150: p1 =   6.47970227 -> excess demand ->     5.78319844 p2 =   2.61489633 -> excess demand ->     2.16047206
    1175: p1 =   6.48042626 -> excess demand ->     5.11684578 p2 =   2.61516680 -> excess demand ->     1.91153572
    1200: p1 =   6.48106683 -> excess demand ->     4.52738798 p2 =   2.61540610 -> excess demand ->     1.69132625
    1225: p1 =   6.48163362 -> excess demand ->     4.00592640 p2 =   2.61561784 -> excess demand ->     1.49651914
    1250: p1 =   6.48213513 -> excess demand ->     3.54459764 p2 =   2.61580519 -> excess demand ->     1.32417666
    1275: p1 =   6.48257889 -> excess demand ->     3.13645198 p2 =   2.61597097 -> excess demand ->     1.17170246
    1300: p1 =   6.48297156 -> excess demand ->     2.77534624 p2 =   2.61611766 -> excess demand ->     1.03680150
    1325: p1 =   6.48331902 -> excess demand ->     2.45584946 p2 =   2.61624746 -> excess demand ->     0.91744484
    1350: p1 =   6.48362648 -> excess demand ->     2.17315978 p2 =   2.61636233 -> excess demand ->     0.81183857
    1375: p1 =   6.48389856 -> excess demand ->     1.92303109 p2 =   2.61646396 -> excess demand ->     0.71839640
    1400: p1 =   6.48413931 -> excess demand ->     1.70170837 p2 =   2.61655391 -> excess demand ->     0.63571554
    1425: p1 =   6.48435236 -> excess demand ->     1.50587065 p2 =   2.61663350 -> excess demand ->     0.56255530
    1450: p1 =   6.48454090 -> excess demand ->     1.33258058 p2 =   2.61670393 -> excess demand ->     0.49781836
    1475: p1 =   6.48470774 -> excess demand ->     1.17923997 p2 =   2.61676625 -> excess demand ->     0.44053408
    1500: p1 =   6.48485538 -> excess demand ->     1.04355049 p2 =   2.61682141 -> excess demand ->     0.38984385
    1525: p1 =   6.48498603 -> excess demand ->     0.92347897 p2 =   2.61687022 -> excess demand ->     0.34498813
    1550: p1 =   6.48510165 -> excess demand ->     0.81722673 p2 =   2.61691341 -> excess demand ->     0.30529496
    1575: p1 =   6.48520396 -> excess demand ->     0.72320246 p2 =   2.61695163 -> excess demand ->     0.27016986
    1600: p1 =   6.48529451 -> excess demand ->     0.63999827 p2 =   2.61698546 -> excess demand ->     0.23908688
    1625: p1 =   6.48537464 -> excess demand ->     0.56636850 p2 =   2.61701539 -> excess demand ->     0.21158067
    1650: p1 =   6.48544555 -> excess demand ->     0.50121103 p2 =   2.61704188 -> excess demand ->     0.18723950
    1675: p1 =   6.48550830 -> excess demand ->     0.44355067 p2 =   2.61706532 -> excess demand ->     0.16569906
    1700: p1 =   6.48556383 -> excess demand ->     0.39252455 p2 =   2.61708607 -> excess demand ->     0.14663701
    1725: p1 =   6.48561298 -> excess demand ->     0.34736916 p2 =   2.61710443 -> excess demand ->     0.12976811
    1750: p1 =   6.48565647 -> excess demand ->     0.30740891 p2 =   2.61712068 -> excess demand ->     0.11483999
    1775: p1 =   6.48569496 -> excess demand ->     0.27204598 p2 =   2.61713506 -> excess demand ->     0.10162931
    1800: p1 =   6.48572902 -> excess demand ->     0.24075136 p2 =   2.61714778 -> excess demand ->     0.08993845
    1825: p1 =   6.48575916 -> excess demand ->     0.21305696 p2 =   2.61715904 -> excess demand ->     0.07959254
    1850: p1 =   6.48578584 -> excess demand ->     0.18854854 p2 =   2.61716900 -> excess demand ->     0.07043683
    1875: p1 =   6.48580944 -> excess demand ->     0.16685953 p2 =   2.61717782 -> excess demand ->     0.06233438
    1900: p1 =   6.48583033 -> excess demand ->     0.14766556 p2 =   2.61718563 -> excess demand ->     0.05516401
    1925: p1 =   6.48584882 -> excess demand ->     0.13067959 p2 =   2.61719253 -> excess demand ->     0.04881849
    1950: p1 =   6.48586518 -> excess demand ->     0.11564759 p2 =   2.61719865 -> excess demand ->     0.04320292
    1975: p1 =   6.48587966 -> excess demand ->     0.10234477 p2 =   2.61720406 -> excess demand ->     0.03823334
    2000: p1 =   6.48589248 -> excess demand ->     0.09057220 p2 =   2.61720884 -> excess demand ->     0.03383541
    2025: p1 =   6.48590382 -> excess demand ->     0.08015385 p2 =   2.61721308 -> excess demand ->     0.02994339
    2050: p1 =   6.48591385 -> excess demand ->     0.07093394 p2 =   2.61721683 -> excess demand ->     0.02649907
    2075: p1 =   6.48592273 -> excess demand ->     0.06277459 p2 =   2.61722015 -> excess demand ->     0.02345095
    2100: p1 =   6.48593059 -> excess demand ->     0.05555381 p2 =   2.61722308 -> excess demand ->     0.02075345
    2125: p1 =   6.48593755 -> excess demand ->     0.04916363 p2 =   2.61722568 -> excess demand ->     0.01836625
    2150: p1 =   6.48594370 -> excess demand ->     0.04350850 p2 =   2.61722798 -> excess demand ->     0.01625364
    2175: p1 =   6.48594915 -> excess demand ->     0.03850387 p2 =   2.61723001 -> excess demand ->     0.01438404
    2200: p1 =   6.48595397 -> excess demand ->     0.03407491 p2 =   2.61723182 -> excess demand ->     0.01272950
    2225: p1 =   6.48595824 -> excess demand ->     0.03015540 p2 =   2.61723341 -> excess demand ->     0.01126527
    2250: p1 =   6.48596201 -> excess demand ->     0.02668675 p2 =   2.61723482 -> excess demand ->     0.00996947
    2275: p1 =   6.48596535 -> excess demand ->     0.02361708 p2 =   2.61723607 -> excess demand ->     0.00882272
    2300: p1 =   6.48596831 -> excess demand ->     0.02090051 p2 =   2.61723717 -> excess demand ->     0.00780788
    2325: p1 =   6.48597093 -> excess demand ->     0.01849641 p2 =   2.61723815 -> excess demand ->     0.00690978
    2350: p1 =   6.48597324 -> excess demand ->     0.01636885 p2 =   2.61723902 -> excess demand ->     0.00611497
    2375: p1 =   6.48597529 -> excess demand ->     0.01448601 p2 =   2.61723978 -> excess demand ->     0.00541160
    2400: p1 =   6.48597711 -> excess demand ->     0.01281975 p2 =   2.61724046 -> excess demand ->     0.00478912
    2425: p1 =   6.48597871 -> excess demand ->     0.01134515 p2 =   2.61724106 -> excess demand ->     0.00423825
    2450: p1 =   6.48598013 -> excess demand ->     0.01004017 p2 =   2.61724159 -> excess demand ->     0.00375075
    2475: p1 =   6.48598139 -> excess demand ->     0.00888530 p2 =   2.61724206 -> excess demand ->     0.00331932
    2500: p1 =   6.48598250 -> excess demand ->     0.00786327 p2 =   2.61724247 -> excess demand ->     0.00293751
    2525: p1 =   6.48598349 -> excess demand ->     0.00695879 p2 =   2.61724284 -> excess demand ->     0.00259962
    2550: p1 =   6.48598436 -> excess demand ->     0.00615835 p2 =   2.61724317 -> excess demand ->     0.00230060
    2575: p1 =   6.48598513 -> excess demand ->     0.00544999 p2 =   2.61724346 -> excess demand ->     0.00203597
    2600: p1 =   6.48598581 -> excess demand ->     0.00482310 p2 =   2.61724371 -> excess demand ->     0.00180179
    2625: p1 =   6.48598641 -> excess demand ->     0.00426832 p2 =   2.61724394 -> excess demand ->     0.00159453
    2650: p1 =   6.48598695 -> excess demand ->     0.00377736 p2 =   2.61724414 -> excess demand ->     0.00141112
    2675: p1 =   6.48598742 -> excess demand ->     0.00334287 p2 =   2.61724431 -> excess demand ->     0.00124881
    2700: p1 =   6.48598784 -> excess demand ->     0.00295835 p2 =   2.61724447 -> excess demand ->     0.00110516
    2725: p1 =   6.48598821 -> excess demand ->     0.00261807 p2 =   2.61724461 -> excess demand ->     0.00097804
    2750: p1 =   6.48598854 -> excess demand ->     0.00231693 p2 =   2.61724473 -> excess demand ->     0.00086554
    2775: p1 =   6.48598883 -> excess demand ->     0.00205042 p2 =   2.61724484 -> excess demand ->     0.00076598
    2800: p1 =   6.48598909 -> excess demand ->     0.00181457 p2 =   2.61724493 -> excess demand ->     0.00067788
    2825: p1 =   6.48598931 -> excess demand ->     0.00160585 p2 =   2.61724502 -> excess demand ->     0.00059990
    2850: p1 =   6.48598951 -> excess demand ->     0.00142114 p2 =   2.61724509 -> excess demand ->     0.00053090
    2875: p1 =   6.48598969 -> excess demand ->     0.00125767 p2 =   2.61724516 -> excess demand ->     0.00046983
    2900: p1 =   6.48598985 -> excess demand ->     0.00111301 p2 =   2.61724522 -> excess demand ->     0.00041579
    2925: p1 =   6.48598999 -> excess demand ->     0.00098498 p2 =   2.61724527 -> excess demand ->     0.00036796
    2950: p1 =   6.48599011 -> excess demand ->     0.00087169 p2 =   2.61724532 -> excess demand ->     0.00032564
    2975: p1 =   6.48599022 -> excess demand ->     0.00077142 p2 =   2.61724536 -> excess demand ->     0.00028818
    3000: p1 =   6.48599031 -> excess demand ->     0.00068269 p2 =   2.61724539 -> excess demand ->     0.00025503


The Walras-equilibrium prices, $(p_1,p_2)$, where both excess demands are (approximately) zero have been found to be: 


```python
f'(p_1,p_2) = ({p1:.3f},{p2:.3f})'
```




    '(p_1,p_2) = (6.486,2.617)'



We see that there is a small change in p1, but beside that th equilibrium prices are almost identical.

At last we will create an interactive figure for changing values of gamma, which will show how the two distributions behave and how their mean and variance are affected.


```python
def utility_fig(gamma):
    
    # Define old utility vector and calculate mean and variance
    utility_vec = utility_func(p1, p2, e1, e2, e3, betas, gamma)
    utility_mean = np.mean(utility_vec)
    utility_var = np.var(utility_vec)
    
     # Define new utility vector and calculate mean and variance
    utility_vec_equal= utility_func(p1, p2, e1_equal, e1_equal, e1_equal, betas, gamma)
    utility_mean_equal = np.mean(utility_vec_equal)
    utility_var_equal = np.var(utility_vec_equal)
    
    # Plot the figures
    fig = plt.figure(figsize=(15,10))
   
    ax1 = fig.add_subplot(2,2,1)
    ax1.hist(utility_vec, bins=500, color ="#6ab2f2")
    ax1.set_title('Distribution of utilities')
    ax1.set_xlabel('Utility')
    ax1.set_ylabel('Number of agents')
    plt.axis([-0.1, 5, 0, 1600])
        
    ax2 = fig.add_subplot(2,2,2)
    ax2.hist(utility_vec_equal, bins=500, color ="#f7d56f")
    ax2.set_title('Distribution of utilities with equal endowments')
    ax2.set_xlabel('Utility')
    ax2.set_ylabel('Number of agents')

    # Plot text with mean and variance
    plt.text(-4, 900, f'$\mu$: {utility_mean:.3f} $\sigma^2$: {utility_var:.3f}')
    plt.text(2.0, 900, f'$\mu$: {utility_mean_equal:.3f} $\sigma^2$: {utility_var_equal:.3f}')
    plt.axis([-0.1, 5, 0, 1600])
    
    plt.show()
# Create slider
interactive_plot = interactive(utility_fig, 
                               gamma=widgets.FloatSlider(description="$\gamma$", 
                                                       min=0.0,
                                                       max=2.0,
                                                       step=0.1,
                                                       value=0.8,
                                                       continuous_update=True) 
                              )

interactive_plot
```


    interactive(children=(FloatSlider(value=0.8, description='$\\gamma$', max=2.0), Output()), _dom_classes=('widg…


First of all we notice how concentrated around the mean the utitilty is for the second plot, due to the variance being much smaller. We also notice how the utility is never below 1. When we increase the value of gamma, we notice how the variance of the first plot increases by a much higher ratio than the second plot. Having the same equal endowments means the change in utility is only due to small changes in preferences.

e59e1d79e889ef119a5dd96c5be36e9019e6230b
ipynb
Jupyter Notebook
examproject/exam_2019.ipynb
NumEconCopenhagen/projects-2019-pickles
7f4d66612bf8f575b745a2c1c32477938e3dcbb6
MIT
examproject/exam_2019.ipynb
NumEconCopenhagen/projects-2019-pickles
7f4d66612bf8f575b745a2c1c32477938e3dcbb6
MIT
2019-04-08T17:31:57.000Z
2019-05-14T18:47:13.000Z
examproject/exam_2019.ipynb
NumEconCopenhagen/projects-2019-pickles
7f4d66612bf8f575b745a2c1c32477938e3dcbb6
MIT
2019-05-14T08:10:26.000Z
2019-12-09T09:29:09.000Z
Qwen/Qwen-72B
1. YES
2. YES