upload materials

spring2024/math/codes/career_choice.py

@@ -0,0 +1,122 @@
+# kimura and yasui 2007
+import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
+from scipy.optimize import fsolve
+
+# latex font
+plt.rc('text', usetex=True)
+plt.rc('font', family='serif')
+
+# parameters
+alpha = 0.33
+tau = 0.6
+gamma = 0.6
+z = 0.2
+b = 0.1
+A = 4.5
+
+# functions
+
+
+def theta_():
+    term1 = (1-tau)**( (1/(alpha*(1-gamma))) - 1  ) / (1-gamma) 
+    term2 = ( (1-alpha)/b) ** (1/alpha)
+    return term1 * term2
+
+theta = theta_()
+k_bar = 1/theta
+
+def k1_(k):
+    if k < k_bar:
+        C1 = A * z * (1-gamma)/gamma 
+        C2 = 1 / (1-tau*theta*k)
+        B1 = (1-alpha) * ((1-tau)**(1-alpha)) * (theta**(1-alpha)) / ((1-gamma)**alpha)
+        return C1 * C2 * (B1 * k + (1-theta*k)*b) 
+    else:
+        C1 = A * z * (1-gamma)/gamma
+        C2 = (1-alpha) / (((1-tau)**alpha) * ((1-gamma)**alpha))
+        return C1 * C2 * (k**alpha)
+
+# finding  m and phi given k
+
+
+def phi_(k):
+    phi_crude = theta * k
+    if phi_crude < 1:
+        return phi_crude
+    else:
+        return 1
+
+
+def m_(k):
+    return (1 - tau * phi_(k)) * gamma / z
+
+k_vals = np.linspace(0, 1, 100)
+k1_vals = np.zeros(100)
+
+for i in range(100):
+    k1_vals[i] = k1_(k_vals[i])
+
+# find the steady states
+def k1L_(k):
+    C1 = A * z * (1-gamma)/gamma 
+    C2 = 1 / (1-tau*theta*k)
+    B1 = (1-alpha) * ((1-tau)**(1-alpha)) * (theta**(1-alpha)) / ((1-gamma)**alpha)
+    return C1 * C2 * (B1 * k + (1-theta*k)*b) - k 
+
+def k1H_(k):
+    C1 = A * z * (1-gamma)/gamma
+    C2 = (1-alpha) / (((1-tau)**alpha) * ((1-gamma)**alpha))
+    return C1 * C2 * (k**alpha) - k
+
+
+
+
+# steady states
+kL = fsolve(k1L_, 0.1)
+kL_unstable = fsolve(k1L_, 0.5)
+kH = fsolve(k1H_, k_bar)
+
+mL = m_(kL)
+phiL = phi_(kL)
+mH = m_(kH)
+phiH = phi_(kH)
+
+plt.figure(figsize=(4, 3))
+# plot the function
+plt.plot(k_vals, k1_vals)
+# plot the 45 degree line
+plt.plot(k_vals, k_vals, linestyle=':', color='grey')
+# plot the steady states
+plt.scatter(kL, kL, color='blue')
+plt.scatter(kH, kH, color='red')
+# plot the unstable steady state
+plt.scatter(kL_unstable, kL_unstable, color='orange')
+# plot the threshold
+plt.vlines(k_bar, 0, k1_(k_bar), linestyle='--', color='grey')
+
+# vertical and horizontal lines
+plt.vlines(kL, 0, kL, linestyle='--', color='grey')
+plt.vlines(kH, 0, kH, linestyle='--', color='grey')
+plt.vlines(kL_unstable, 0, kL_unstable, linestyle='--', color='grey')
+
+# some cosmetics
+plt.xlabel('$k_t$', fontsize=15, loc='right')
+plt.ylabel('$k_{t+1}$', fontsize=15, rotation=0, loc='top')
+# disable the top and right frame
+ax = plt.gca()
+ax.spines['top'].set_visible(False)
+ax.spines['right'].set_visible(False)
+plt.ylim(0, 1)
+ax.margins(0)
+# annotate xticks and yticks
+plt.xticks([kL[0], kH[0], kL_unstable[0], k_bar], [
+           r'$k_L$', r'$k_H$', r'$k_U$', r'$\bar{k}$'], fontsize=12)
+# print the steady states m and phi
+print('kL:', kL)
+print('mL:', mL)
+print('phiL:', phiL)
+print('kH:', kH)
+print('mH:', mH)
+print('phiH:', phiH)

spring2024/math/codes/childcare.py

@@ -0,0 +1,56 @@
+import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
+
+# latex font
+plt.rc('text', usetex=True)
+plt.rc('font', family='serif')
+
+
+# parameters
+delta = 0.6
+psi = 0.1
+s_bar = 1
+phi = 0.4
+w_m = 1
+p_s = 0.5
+
+# functions
+
+
+def sthres_():
+    return (w_m * phi - psi) / ((p_s + w_m)*phi)
+
+
+def s_(w_f):
+    if w_f < p_s:
+        return 0
+    else:
+        return s_bar
+
+
+def n_(w_f):
+    s = s_(w_f)
+    return (delta/(1+delta)) * (w_m + w_f) / (psi + (s*p_s + (1-s)*w_f)*phi)
+
+
+w_f_vals = np.linspace(0, w_m, 100)
+n_vals = np.zeros(100)
+for i in range(100):
+    n_vals[i] = n_(w_f_vals[i])
+
+plt.figure(figsize=(4, 3))
+plt.plot(w_f_vals, n_vals)
+plt.xlabel('$w_f$', fontsize=15, loc='right')
+plt.ylabel('$n$', fontsize=15, rotation=0, loc='top')
+plt.scatter(sthres_(), n_(p_s), color='red')
+plt.vlines(p_s, 0, n_(p_s), linestyle='--', color='grey')
+plt.hlines(n_(p_s), 0, p_s, linestyle='--', color='grey')
+plt.ylim(1.4, 4)
+ax.margins(0)
+# disable the top and right frame
+ax = plt.gca()
+ax.spines['top'].set_visible(False)
+ax.spines['right'].set_visible(False)
+
+plt.show()

spring2024/math/codes/fertility_edu.py

@@ -0,0 +1,68 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+
+def w_threshold(theta, phi, eta):
+    return theta / (phi * eta)
+
+
+def e_(theta, phi, eta, w, w_thres):
+    if w > w_thres:
+        return (eta * phi * w - theta) / (1-eta)
+    else:
+        return 0
+
+
+def n_(phi, theta, eta, gamma, w, w_thres):
+    if w > w_thres:
+        return (1-eta) * gamma * w / ((1+gamma) * (phi*w - theta))
+    else:
+        return gamma / (phi*(1+gamma))
+
+
+theta = 51.61
+phi = 0.039
+eta = 0.572
+gamma = 0.1
+w_thres = w_threshold(theta, phi, eta)
+w_thres
+
+# latex
+plt.rc('text', usetex=True)
+plt.rc('font', family='serif')
+
+wvals = np.linspace(w_thres-2000, w_thres + 10000, 1000)
+evals = [e_(theta, phi, eta, w, w_thres) for w in wvals]
+nvals = [n_(phi, theta, eta, gamma, w, w_thres) for w in wvals]
+
+plt.figure(figsize=(4, 3))
+plt.plot(wvals, evals)
+plt.xticks([w_thres], [r'$\theta/(\eta \phi)$'], fontsize=14)
+plt.yticks([0])
+plt.xlabel('w', fontsize=14, loc='right')
+plt.ylabel('e', fontsize=14, loc='top', rotation=0)
+# remove top and right frames
+ax = plt.gca()
+ax.spines['top'].set_visible(False)
+ax.spines['right'].set_visible(False)
+
+N_max = gamma/(phi*(1+gamma))
+N_min = gamma*(1-eta)/(phi*(1+gamma))
+print(N_max, N_min)
+
+plt.figure(figsize=(4, 3))
+plt.plot(wvals, nvals)
+plt.xticks([w_thres], [r'$\theta/(\eta \phi)$'], fontsize=14)
+# plt.yticks([])
+plt.xlabel('w', fontsize=14, loc='right')
+plt.ylabel('n', fontsize=14, loc='top', rotation=0)
+plt.ylim(N_min-0.2, N_max+0.2)
+# remove top and right frames
+ax = plt.gca()
+ax.spines['top'].set_visible(False)
+ax.spines['right'].set_visible(False)
+plt.vlines(w_thres, N_max, 0, linestyle='--', color='grey')
+plt.hlines(N_min, w_thres, wvals[-1], linestyle='--', color='black')
+plt.yticks([N_min, N_max], [
+           r'$\frac{\gamma(1-\eta)}{\phi(1+\gamma)}$', r'$\frac{\gamma}{\phi(1+\gamma)}$'], fontsize=14)
+plt.show()

spring2024/math/codes/life_expectancy.py

@@ -0,0 +1,81 @@
+# library
+import numpy as np
+import matplotlib.pyplot as plt
+from scipy.optimize import fsolve
+
+# parameters
+A = 5
+b = 0.2
+z = 0.12
+beta = 0.99**35
+gamma = 0.4
+sigma = 0.45
+alpha = 0.35
+
+# key param
+x1 = 0.6
+x2 = 0.7
+
+# functions
+def eta_(x):
+    return (1 + gamma + beta*x) / (1 + beta*x)
+
+def theta_(x):
+    return ((1-alpha)/b)**(1/alpha) * (1-sigma)**((eta_(x)/alpha) - 1) * eta_(x)
+
+def kbar_(x):
+    return ((b/(1-alpha))**(1/alpha)) * ((1-sigma)**(1 - (eta_(x))/alpha)) / eta_(x)
+
+def k1_(k, x):
+    if k < kbar_(x):
+        C = A * beta * z * x / (gamma * (1 - sigma * theta_(x) * k))
+        B1 = (((1-alpha)/(b**(1-alpha)))**(1/alpha)) * ((1-sigma)**((1-alpha)*eta_(x)/alpha)) * eta_(x) 
+        return C * (B1 * k + (1 - theta_(x)*k)*b)
+    else:
+        C = A * beta * z * x * (1-alpha) / (gamma*((1-sigma)**alpha))
+        return C * ((eta_(x))**alpha) * (k**alpha)
+
+k_vals = np.linspace(0.01, 1, 100)
+k1_vals = np.zeros(100)
+k2_vals = np.zeros(100)
+
+for i in range(100):
+    k1_vals[i] = k1_(k_vals[i], x1)
+    k2_vals[i] = k1_(k_vals[i], x2)
+
+# find the steady states
+def k1L_(k, x):
+    C = A * beta * z * x / (gamma * (1 - sigma * theta_(x) * k))
+    B1 = (((1-alpha)/(b**(1-alpha)))**(1/alpha)) * ((1-sigma)**((1-alpha)*eta_(x)/alpha)) * eta_(x) 
+    return C * (B1 * k + (1 - theta_(x)*k)*b) - k
+
+def k1H_(k, x):
+    C = A * beta * z * x * (1-alpha) / (gamma*((1-sigma)**alpha))
+    return C * ((eta_(x))**alpha) * (k**alpha) - k
+
+# solve for the steady states
+kL1 = fsolve(k1L_, kbar_(x1), args=(x1))
+kH1 = fsolve(k1H_, kbar_(x1), args=(x1))
+kH2 = fsolve(k1H_, kbar_(x2), args=(x2))
+kL2 = fsolve(k1L_, kbar_(x2), args=(x2))
+
+
+plt.figure(figsize=(4, 3))
+# plot the function
+plt.plot(k_vals, k1_vals, label=r'$x=0.6$')
+plt.plot(k_vals, k2_vals, label=r'$x=0.7$')
+# plot the 45 degree line
+plt.plot(k_vals, k_vals, linestyle=':', color='grey')
+# plot the steady states
+#plt.scatter(kL, kL, color='blue')
+plt.scatter(kH1, kH1, color='red')
+plt.scatter(kH2, kH2, color='orange')
+# delete top and right frame
+ax = plt.gca()
+ax.spines['top'].set_visible(False)
+ax.spines['right'].set_visible(False)
+ax.margins(0)
+plt.legend()
+plt.xlabel('$k_t$', fontsize=15, loc='right')
+plt.ylabel('$k_{t+1}$', fontsize=15, rotation=0, loc='top')
+plt.show()