## 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()