## initialize
alpha = 0.3
delta = 0.015
rho = 0.05
epsilon = 0.5
def f(k):
return k**alpha
def f_prime(k):
return alpha*k**(alpha - 1)
k_star = ((delta + rho)/alpha) ** (1/(alpha - 1))
c_star = f(k_star) - delta * k_star
print(k_star, c_star)
from scipy.integrate import solve_ivp
from scipy.optimize import brentq
import matplotlib.pyplot as plt
import numpy as np
sigma=1/epsilon
k_ss=(1/alpha*(delta+rho))**(1/(alpha-1))
c_ss=f(k_ss)-delta*k_ss
print(k_ss, c_ss)
def system(t,y):
c,k=y[0],y[1]
dc=epsilon*c*(f_prime(k)-delta-rho)
dk=f(k)-delta*k-c
return [dc,dk]
def hit_0(t,y):
return np.min(y)-1e-4
hit_0.terminal=True
def reach_ss(t,y):
return (y[0]-c_ss)**2+(y[1]-k_ss)**2-1e-3
reach_ss.terminal=True
k0= k_ss / 2
T=100
threshold = 1e-3
K=np.linspace(0,60,200)
fig,ax=plt.subplots()
ax.axvline(x=k_ss,c='black',lw=1)
ax.axvline(x=k0,ls='--',c='orange')
ax.plot(K,f(K)-delta*K,c='black',lw=1)
ax.scatter(k_ss,c_ss,c='black')
ax.set_xlabel('k')
ax.set_ylabel('c')
ax.set_ylim(0,3)
ax.set_xlim(0,40)
fig.show()
def solve(T):
threshhold = 1e-3
guess_low = 0
guess_high = 12
solved = False
i = 0
while solved is False:
guess = (guess_low + guess_high) / 2
sol = solve_ivp(system,[0,T],y0=(guess,k0),max_step=1e-1,events=hit_0)
c_sol = sol.y[0]
k_sol = sol.y[1]
# print(len(c_sol), len(k_sol))
if len(k_sol) < T*10+2 - 1:
pad = T*10 + 2 - len(k_sol)
k_sol = list(k_sol) + [0]*pad
# print(k_sol)
k_end = k_sol[-1]
k_before_end = k_sol[-2]\
if abs(k_end - 0) < threshhold and abs(k_before_end - 0) > threshhold:
solved = True
elif abs(k_end - 0) > threshhold:
guess_low = guess
else:
guess_high = guess
i += 1
if i % 100 == 0:
print(guess)
return k_sol, c_sol
k_sol_5, c_sol_5 = solve(5)
k_sol_20, c_sol_20 = solve(20)
k_sol_50, c_sol_50 = solve(50)
K=np.linspace(0,60,200)
fig,ax=plt.subplots()
ax.scatter(k_ss,c_ss,c='black')
ax.axvline(x=k_ss,c='black',lw=1)
ax.axvline(x=k0,ls='--',c='orange')
ax.plot(K,f(K)-delta*K,c='black',lw=1)
ax.plot(k_sol_50,c_sol_50, label = "T = 50")
ax.plot(k_sol_20,c_sol_20, label = "T = 20")
ax.plot(k_sol_5, c_sol_5, label = "T = 5")
ax.legend()
ax.set_xlabel('k')
ax.set_ylabel('c')
ax.set_ylim(0,4)
ax.set_xlim(0,15)
fig.show()