```python
import random 
import sympy as sm
import numpy as np
import matplotlib.pyplot as plt
import ipywidgets as widgets
from scipy import optimize

#setting nice printing and seed
sm.init_printing(use_unicode=True)
random.seed(12848174)
```


```python
rho = sm.symbols('rho')
beta = sm.symbols('beta')
gamma = sm.symbols('gamma')
w = sm.symbols('w')
b = sm.symbols('b')
delta = sm.symbols('Delta')
h1 = sm.symbols('h_1')
h2 = sm.symbols('h_2')
c1 = sm.symbols('c1')
c2 = sm.symbols('c2')
l1 = sm.symbols('l1')
l2 = sm.symbols('l2')
exp = sm.symbols('E[v_2(h_2)]')
u = sm.symbols("u")
```


```python
u2 = (c2**(1-rho))/(1-rho)
dis_u2 = gamma*l2
total_u2 = u2 - dis_u2
h2_f = (h1 + l1)
c2_f = w*h2*l2
total_u2
```


```python
total_u2_sub2 = total_u2.subs(c2,c2_f)
#foc2 = sm.diff(total_u2_sub2,l2)
#foc2
total_u2_sub2
```


```python
#we should not solve for l2 rigth? this is just a dummy of work or not
#sol2 = sm.solve(sm.Eq(foc2,0),l2)[0]
#sol2
h2_sol = sm.solve(sm.Eq(total_u2_sub2,u),h2)[0]
h2_sol
```


```python
def utility_work(h2,l2,w,rho,gamma):
    solution = (((h2*l2*w)**(1-rho))/(1-rho))-gamma*l2
    return solution

def utility_notwork(benefits,rho,gamma,l2):
 
    b = ((benefits**(1-rho))/(1-rho))-gamma*l2
    return b
def hum_cap_u(l2,gamma,b,w,rho):
    # i use h2_sol but replace u with the "cutoff"-value, that is the value of unemployment benefits
    nw_u = utility_notwork(b,rho,gamma,0)
    hum_cap = ((-gamma*l2*rho+gamma*l2-rho*nw_u+nw_u)**-(1/(rho-1)))/(l2*w)
    
    return hum_cap
```


```python
rho = 2
beta = 0.96
gamma = 0.1
w = 2
b = 1
delta = 0.1
```


```python
#I calculatet the humam capital given b
hum_cap_u(1,gamma,b,w,rho)
```


```python
def plot_func(h2,w,rho,gamma,b):
    
    work = utility_work(h2,1,w,rho,gamma)
    not_work = utility_notwork(b,rho,gamma,0)
    cutoff = hum_cap_u(1,gamma,b,w,rho)
    
    plt.plot(h2,work,label = "Utility of working")
    plt.plot(h2,not_work,color="red",label="Utility of not working")
    plt.legend()
    
    plt.xlabel("Human capital")
    plt.ylabel("Utility")

    plt.show
    
    print("The worker will supply labour if human capital is above",round(cutoff[0],2), "and no labour if it is below" )
    return 
```


```python
h_vec = np.linspace(0.1,1.5,100) 
b_vec = np.linspace(1,1,100)
plot_func(h_vec,w,rho,gamma,b_vec)
```


```python
#Another way of showing this is using a binary plot. This requires us to define 
#each level om human capital as either a capital where the worker provies labour
# or a level where she does not - the cutoff is as shown before
cutoff = hum_cap_u(1,gamma,b,w,rho)
print(round(cutoff,2),"human capital")
```

    0.56 human capital



```python
#We define a list
willwork = list(range(100))

#And fill it with either success (1) if the human capital i above the cutoff or
#or falier (0) if it is on the cutoff or below 
for i in range(len(h_vec)): 
    if h_vec[i] > cutoff : 
        willwork[i] = 1
    else : 
        willwork[i] = 0
```


```python
cutoff_vec = np.linspace(cutoff,cutoff,100)
plt.plot(h_vec,willwork)
plt.plot(cutoff_vec,willwork,linestyle="dashed",color="red",label=("Cutoff human capital level"))
plt.legend()
plt.xlabel("Human capital")
plt.ylabel("Labour supply")

plt.show
```

some text here


```python
#the labour supply as a function of human captal is then:
def labour_supply(hum_cap,w,rho,gamma,b):
    cutoff = hum_cap_u(1,gamma,b,w,rho)
    
    if hum_cap > cutoff :
        print("with a human capital level of",hum_cap,"the worker wil supply 1 labour (l = 1)")
        l = 1
    else :
        print("with a human capital level of",hum_cap,"the worker will supply 0 labour (l=0)")
        l = 0
    return l
```


```python
human_capital = widgets.FloatSlider(min=0.1,max=1.5,value=0.1,step=0.01)
widgets.interact(labour_supply,hum_cap=human_capital,w=widgets.fixed(w),rho=widgets.fixed(rho),gamma=widgets.fixed(gamma),b=widgets.fixed(b))
```


    interactive(children=(FloatSlider(value=0.1, description='hum_cap', max=1.5, min=0.1, step=0.01), Output()), _…





    <function __main__.labour_supply(hum_cap, w, rho, gamma, b)>


Interactive solution to question 1

```python

```

## Question 2


```python
rho = sm.symbols('rho')
beta = sm.symbols('beta')
gamma = sm.symbols('gamma')
w = sm.symbols('w')
b = sm.symbols('b')
delta = sm.symbols('Delta')
h1 = sm.symbols('h_1')
h2 = sm.symbols('h_2')
c1 = sm.symbols('c1')
c2 = sm.symbols('c2')
l1 = sm.symbols('l1')
l2 = sm.symbols('l2')
exp = sm.symbols('E[v_2(h_2)]')
u = sm.symbols("u")
```


```python
u1 = (c1**(1-rho))/(1-rho)
dis_u1 = gamma*l1
total_u1 = u1 - dis_u1 + beta*exp
total_u1
```


```python
# I define funcitons for period 2 
def utility_2(h1,w=1,l1=1,rho=2,gamma=0.1,beta=0.96,delta=0.1,b=1):
    if l1 == 1:
        h2 = h1+1 +(0.5*delta)
        c = w*h1*1
        
        u_p2 = utility_work(h2,1,w,rho,gamma)
    
        solution = ((c**(1-rho))/(1-rho))-gamma*1 +beta*u_p2
    else :
        h2 = h1 + 0 + (0.5*delta)
        c = b
        
        u_nw_p2 = utility_notwork(b,rho,gamma,0)
        solution = ((c**(1-rho))/(1-rho))-gamma*0 +beta*u_nw_p2
    
    return solution



```


```python
h1 = 0.1
utility_2(h1,w,1,rho,gamma,beta,delta,b)
```


```python
rho = 2
beta = 0.96
gamma = 0.1
w = 2
b = 1
delta = 0.1
h_vec = np.linspace(0.1,1.5,100) 

```


```python
b_vec = np.linspace(1,1,100)
x = utility_2(h_vec,w,1,rho,gamma,beta,delta,b)
y = utility_2(h_vec,w,0,rho,gamma,beta,delta,b_vec)

```


```python
#The utility of not working is 
y[0]
```


```python
#The cutof point can be found using a numerical optimizer 

rho = 2
beta = 0.96
gamma = 0.1
w = 2
b = 1
delta = 0.1

#because it is negative values we minimize (and not -minimize)
def objective(h):
    return utility_2(h,w,1,rho,gamma,beta,delta,b)
def constraint(h):
    return (1.96 + utility_2(h,w,1,rho,gamma,beta,delta,b))
```


```python
cons = ({"type": "ineq","fun":constraint})
h_init = 0.1

```


```python
sol = optimize.minimize(objective,h_init,method="SLSQP",constraints=cons)
sol
```




         fun: -1.9599999941585289
         jac: array([4.28598835])
     message: 'Optimization terminated successfully.'
        nfev: 21
         nit: 7
        njev: 7
      status: 0
     success: True
           x: array([0.35172508])




```python
#The equilibrium value is 
sol.x[0]
```


```python
cutoff = np.linspace(sol.x[0],sol.x[0],100)
plt.plot(h_vec,x, label = "Utility of working",alpha= 0.8)
plt.plot(h_vec,y, color="green", label = "Utility of not working",alpha= 0.8)
plt.plot(cutoff,x,color = "red", linestyle = "dashed",label="Labour cutoff",alpha= 0.4)
plt.scatter(cutoff[0],y[0],color="purple",alpha = 1)
plt.xlabel("Human capital level")
plt.ylabel("Utility")
plt.legend()
plt.show()
```


```python
willwork_p1 = list(range(100))

#And fill it with either success (1) if the human capital i above the cutoff or
#or falier (0) if it is on the cutoff or below 
for i in range(len(h_vec)): 
    if utility_2(h_vec[i],w,1,rho,gamma,beta,delta,1) > utility_2(h_vec[i],w,0,rho,gamma,beta,delta,1) : 
        willwork_p1[i] = 1
    else : 
        willwork_p1[i] = 0
```


```python
plt.plot(h_vec,willwork_p1)
plt.plot(cutoff,np.linspace(0,1,100),linestyle = "dashed",color="red")
```

This means that the worker will supply all her labour (1) if her human capital level i above 0.35, but 0 if her human capital level is below this. 

## Question 3


```python

```


```python
#itinerary 
## make st of w's with 0.1,0.2 and so on til 0.9 
### calculate utility for those 
#### see if any of them are higher that b

```

We will examine if workers will ever work if their wage is lower than the unemployment benefits. In this case the unemployment benefits are given as 1. It is theoretically posible that a worker would chose to work, if the acumulation of human capital can more than make up for the loss in utility. 


```python
#We define a list of different wage levels to examnie 
wages = np.linspace(0,0.99,100)
#We have chosen these for the sake of interpertation 
```


```python
rho = 2
beta = 0.96
gamma = 0.1
b = 1
delta = 0.1
h_vec = np.linspace(0.1,1.5,100)

```


```python
#This for loop test if the utility gained from working, even if wages are lover than b, is higher than -1.96 - the utility of
#not working in both periods 
for wage in wages: 
    name = "wage_" + str(wage)
    cutoff_w = utility_2(h1 = h_vec, w = wage, l1 = 1, rho = 2, gamma = 0.1, beta= 0.96, delta = 0.1, b = 1)
    if max(cutoff_w) > -1.96:
        print(name) 
        #This will be the lowest wage where it is profitable to work 
        break
```

    wage_0.6


    C:\Users\augus\Anaconda3\lib\site-packages\ipykernel_launcher.py:2: RuntimeWarning: divide by zero encountered in reciprocal
      
    C:\Users\augus\Anaconda3\lib\site-packages\ipykernel_launcher.py:9: RuntimeWarning: divide by zero encountered in reciprocal
      if __name__ == '__main__':



```python
#It tursout that it is utility profitable for the worker to work even if the wage is 0.6 while the unemployment benefits are 1
## Now we will show this 
```


```python
cutoff = np.linspace(sol.x[0],sol.x[0],100)
plt.plot(h_vec,x, label = "Wage = 2",alpha= 0.6)
plt.plot(h_vec,y, color = "red", linestyle = "dashed",alpha= 0.5, label = "Utility of not working")
plt.plot(h_vec,cutoff_w, color = "green", label = "Wage = 0.6")
plt.xlabel("Human capital level")
plt.ylabel("Utility")
plt.ylim(-5,0)
plt.legend()
plt.show()
```


```python
#It is evident that at very high levels of human capital, it can be rational to work, even if the wage is lower than the 
#unemployment benefits. From our loop before we know that it is the lowest wage where this is true is 0.6. Therefore we 
# use optimization to find the human capital level of at which i makes sense to work if your wage is 0.6. 

```


```python

```


```python
rho = 2
beta = 0.96
gamma = 0.1
w = 0.6
b = 1
delta = 0.1

#because it is negative values we minimize (and not -minimize)
def objective(h):
    return utility_2(h,w,1,rho,gamma,beta,delta,b)
def constraint(h):
    return (1.96 + utility_2(h,w,1,rho,gamma,beta,delta,b))
```


```python
cons = ({"type": "ineq","fun":constraint})
h_init = 0.1

```


```python
sol = optimize.minimize(objective,h_init,method="SLSQP",constraints=cons)
sol
```




         fun: -1.960000000389105
         jac: array([1.01749761])
     message: 'Optimization terminated successfully.'
        nfev: 27
         nit: 9
        njev: 9
      status: 0
     success: True
           x: array([1.47461354])




```python
sol.x[0]
```


```python
#I test that there is no wage value below 0.6 where the worker would work 
w = 0.59

#because it is negative values we minimize (and not -minimize)
def objective(h):
    return utility_2(h,w,1,rho,gamma,beta,delta,b)
def constraint(h):
    return (1.96 + utility_2(h,w,1,rho,gamma,beta,delta,b))
cons = ({"type": "ineq","fun":constraint})
h_init = 0.1
sol = optimize.minimize(objective,h_init,method="SLSQP",constraints=cons)
if sol.x[0] > 1.5:
    print("No level of human capital where the worker would work, at human capital level of",str(w))
else: 
    print("The worker til work if she has a human capital level of", str(round(sol.x[0],2)))
```

    No level of human capital where the worker would work, at human capital level of 0.59



```python
#So the worker will work at a wage of 0.6, if she has a human capital level og 1.47 or higher 
```


```python
#Define wages to plot 
wages = [0.6,0.7,0.8,0.9,0.99]
#And the parameters again
rho = 2
beta = 0.96
gamma = 0.1
b = 1
delta = 0.1 
h_vec = np.linspace(0.1,1.5,100)
```


```python
for wage in wages: 
    wage_vec = utility_2(h1 = h_vec, w = wage, l1 = 1, rho = 2, gamma = 0.1, beta= 0.96, delta = 0.1, b = 1)
    plt.plot(h_vec,wage_vec, label = "Wage ="+str(wage))

plt.plot(h_vec,x, label = "Wage = 2",alpha= 0.6)
plt.plot(h_vec,y, color = "red", linestyle = "dashed",alpha= 0.5, label = "Utility of not working")
plt.xlabel("Human capital level")
plt.ylabel("Utility")
plt.ylim(-5,-0.5)
plt.legend()
plt.show()
```


```python
#The respective "cutoff"-levels of human capital can be calculated using numerical optimazation 
```


```python
for wage in wages:
    w = wage 
    def objective(h):
        return utility_2(h,w,1,rho,gamma,beta,delta,b)
    def constraint(h):
        return (1.96 + utility_2(h,w,1,rho,gamma,beta,delta,b))
    cons = ({"type": "ineq","fun":constraint})
    h_init = 0.1
    
    sol = optimize.minimize(objective,h_init,method="SLSQP",constraints=cons)
    print("Human capital level =",round(sol.x[0],2),"with wages =", wage)
```

    Human capital level = 1.47 with wages = 0.6
    Human capital level = 1.23 with wages = 0.7
    Human capital level = 1.05 with wages = 0.8
    Human capital level = 0.91 with wages = 0.9
    Human capital level = 0.81 with wages = 0.99



```python
#We could make this interactiv - it will take some time though 
```


```python
#Explain why people chose to work even if the wage is lowwer than b
```

## Assigmnet 2 


```python
def corr_func_opt(phi,v_m=0, s_m=0):

    np.random.seed(1)
    pi_m = 0
    yn = 0
    y_m = 0
    xt = np.random.normal(loc=0,scale=pms['sigma_x'],size=1000)
    ct = np.random.normal(loc=0,scale=pms['sigma_c'],size=1000)

    
    list_output = list(range(1000))
    list_inflation = list(range(1000))
    for i in range(1000):
        if i == 0 :
            
            y = equili_output_ar(v_m,pi_m,s_m,y_m,ct[i],xt[i])
            list_output[i] = y
            yn = y
            p = equili_pi_phi(y_m,yn,pi_m,s_m,ct[i],phi)
            list_inflation[i] = p  
            pi_m = p
            y_m = y
            v_m = pms["delta"] * v_m + xt[i]
            s_m = pms["omega"] * s_m + ct [i] 
        else :
            y = equili_output_ar(v_m,pi_m,s_m,y_m,ct[i],xt[i])
            
            list_output[i] = y
            yn = y
            p = equili_pi_phi(y_m,yn,pi_m,s_m,ct[i],phi)
            list_inflation[i] = p  
            pi_m = p
            y_m = y
            v_m = pms["delta"] * v_m + xt[i]
            s_m = pms["omega"] * s_m + ct[i]
    
    correlation = np.corrcoef(list_inflation,list_output)[0,1]

    return correlation

```


```python
v_m = 0 
s_m = 0
def objective(phi):
    return corr_func_opt(phi, v_m, s_m)
def constraint(phi):
    return (0.31 - corr_func_opt(phi, v_m, s_m))
```


```python

```


```python
constraint(0.3)
```


```python
cons = ({"type": "ineq","fun":constraint})
phi_init = 0.1

```


```python
sol = optimize.minimize(objective,phi_init,method="SLSQP",constraints=cons)
sol
```


```python

```

3262705d6593d7323b11f59c238cdec9aa5f573e
ipynb
Jupyter Notebook
examproject/August-checkpoint.ipynb
NumEconCopenhagen/projects-2019-al
a5966ce338b1aab9dc52f1b05a1f7a54bd1b25f9
MIT
examproject/August-checkpoint.ipynb
NumEconCopenhagen/projects-2019-al
a5966ce338b1aab9dc52f1b05a1f7a54bd1b25f9
MIT
2019-04-10T07:46:39.000Z
2019-05-14T16:40:44.000Z
examproject/August-checkpoint.ipynb
NumEconCopenhagen/projects-2019-al
a5966ce338b1aab9dc52f1b05a1f7a54bd1b25f9
MIT
2019-04-08T14:04:53.000Z
2019-05-13T08:06:52.000Z
Qwen/Qwen-72B
1. YES
2. YES