<jupyter_start><jupyter_text>Value Function Iteration  Log-Utility<jupyter_code>from numpy import log, linspace, zeros, interp, sqrt
from scipy.interpolate import PchipInterpolator
import matplotlib.pyplot as plt

# Problem Params
initial_resource = 100;
beta = 0.9;
utility = lambda x: log(x)

# Hyperparams
grid = 50;
maxIter = 50;
min_state =1e-5;
exp = 2 # oversampling

# Arrays
state = linspace(min_state**(1/exp),initial_resource**(1/exp),grid)**exp ;
value = zeros((2, grid));
policy = zeros(grid);
current = iter = 1; 
value[current, :] = utility(state)

while (iter <= maxIter):
  iter = iter + 1;
  previous = current; 
  current = 1 - (previous == 1);
  for i in range(grid):
    min_control = min_state/10
    old = PchipInterpolator(state, value[previous, :], axis=0,extrapolate=True)
    def w(x):
      #oldval = interp(state[i] - x, state, value[previous, :])   
      return - utility(x) - beta * old(state[i] - x)
    res = minimize(w, min_control, bounds= ([(min_control, state[i])]))
    newpol, newval = res.x, -res.fun
    value[current, i] = newval
    policy[i] = newpol

value_actual = lambda x: (1/(1-beta))*log((1-beta)*x) + beta*log(beta)/((1-beta)**2)

# Verification
plt.figure(1)
plt.plot(value[current, :], label = 'Numerical' ) # numerical
plt.plot(value_actual(state), label = 'Theoretical') # theoretical
plt.title('Value Function')
plt.legend()
plt.figure(2)
plt.plot(state, policy, label = 'Numerical' ) # numerical
plt.plot(state, (1-beta)*state, label = 'Theoretical') # theoretical
plt.title('Optimal Policy Rule')
plt.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>CARA Utility<jupyter_code>from numpy import log, linspace, zeros, interp
from scipy.interpolate import PchipInterpolator
import matplotlib.pyplot as plt

# Problem Params
initial_resource = 100;
beta = 0.9;rho = 0.75
utility = lambda x: x**(1-rho)/(1-rho)

# Hyperparams
grid = 100;
maxIter = 100;
min_state =1e-5;
exp = 2 # oversampling

# Arrays
state = linspace(min_state**(1/exp),initial_resource**(1/exp),grid)**exp ;
value = zeros((2, grid));
policy = zeros(grid);
current = iter = 1; 
value[current, :] = utility(state)

while (iter <= maxIter):
  iter = iter + 1;
  previous = current; 
  current = 1 - (previous == 1);
  for i in range(grid):
    min_control = min_state/10
    old = PchipInterpolator(state, value[previous, :], axis=0,extrapolate=True)
    def w(x):
      #oldval = interp(state[i] - x, state, value[previous, :])   
      return - utility(x) - beta * old(state[i] - x)
    res = minimize(w, min_control, bounds= ([(min_control, state[i])]))
    newpol, newval = res.x, -res.fun
    value[current, i] = newval
    policy[i] = newpol

# Verification
plt.figure(1)
plt.plot(state, value[current, :], label = 'Numerical' ) # numerical
plt.plot(state, (1-beta**(1/rho))**(-rho)*utility(state), label = 'Analytical') # numerical
plt.title('Value Function')
plt.legend()
plt.figure(2)
plt.plot(state, policy, label = 'Numerical' ) # numerical
plt.plot(state, (1-beta**(1/rho))*state, label = 'Analytical' ) # numerical
plt.title('Optimal Policy Rule')
plt.legend()
plt.show()

from numpy import log, linspace, zeros, interp, sqrt
from scipy.interpolate import PchipInterpolator
import matplotlib.pyplot as plt

# Problem Params
initial_resource = 100;
beta = 0.9;
utility = lambda x: log(x)

# Hyperparams
grid = 50;
maxIter = 50;
min_state =1e-5;
exp = 2 # oversampling

# Arrays
state = linspace(min_state**(1/exp),initial_resource**(1/exp),grid)**exp ;
value = zeros((2, grid));
policy = zeros(grid);
current = iter = 1; 
value[current, :] = utility(state)

while (iter <= maxIter):
  iter = iter + 1;
  previous = current; 
  current = 1 - (previous == 1);
  for i in range(grid):
    min_control = min_state/10
    old = PchipInterpolator(state, value[previous, :], axis=0,extrapolate=True)
    def w(x):
      #oldval = interp(state[i] - x, state, value[previous, :])   
      return - utility(x) - beta * old(state[i] - x)
    res = minimize(w, min_control, bounds= ([(min_control, state[i])]))
    newpol, newval = res.x, -res.fun
    value[current, i] = newval
    policy[i] = newpol

value_actual = lambda x: (1/(1-beta))*log((1-beta)*x) + beta*log(beta)/((1-beta)**2)

# Verification
plt.figure(1)
plt.plot(value[current, :], label = 'Numerical' ) # numerical
plt.plot(value_actual(state), label = 'Theoretical') # theoretical
plt.title('Value Function')
plt.legend()
plt.figure(2)
plt.plot(state, policy, label = 'Numerical' ) # numerical
plt.plot(state, (1-beta)*state, label = 'Theoretical') # theoretical
plt.title('Optimal Policy Rule')
plt.legend()
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1664516129032258, 'fraction_numerical': 0.02043010752688172, 'mean_word_length': 4.583433373349339, 'pattern_counts': {'":': 0, '<': 12, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '2', 'dir': 'jupyter-structured-clean-dedup', 'id': '106061', 'license': 'MIT', 'n_tokens_mistral': 1802, 'n_tokens_neox': 1675, 'n_words': 396, 'path': 'CakeEatingDP.ipynb', 'repo_name': 'pranjalrawat007/QuantitativeMacroeconomics', 'size': '4650'}