Solving cass-koopmans planning model with python value function iteration

Quoted from Forex.com.bd-Editorial External question — Economics Stack Exchange Author: dueun Source score (net votes, not local likes): 6 Original post: https://economics.stackexchange.com/questions/51160 License: CC BY-SA 4.0 — https://creativecommons.org/licenses/by-sa/4.0/ Adaptation: HTML converted to plain text; contact email addresses removed. edit: sorry, I used stack exchange first time, so I didn't know I can type python code. you can now copy and run the code. I'm practicing solving economic model using python programming. I refer to https://python.quantecon.org/optgrowth.html and https://python.quantecon.org/optgrowth_fast.html . I wanna solve Cass-Koopmans planning problem with bellman equation. There is only one agent who governs all resource allocations in the economy. She produces a (single) good through a production function: $$ y_{t} = f(k_{t}) \tag{1} $$ The good can be used for consumption and investment. In this exercise, we introduce a partial capital depreciation, i.e. capital depreciates at the rate $0<\delta<1$ each period. Hence the resource constraint for the agent is: $$ k_{t+1} + c_t = y_t + (1-\delta)k_t \tag{2} $$ The agent wants to maximize $$ \mathbb E \left[ \sum_{t = 0}^{\infty} \beta^t u(c_t) \right] \tag{3} $$ subject to (2), where $ \beta \in (0, 1) $ is a discount factor. The model is almost identical to Cass-Koopmans planning model. In order to solve this model, we will apply value function iteration algorithm. Hence, we first reformulate the maximization problem into a Bellman equation. We can write the value function for the agent's utility maximization problem in the form of Bellman equation . In this exercise we will use capital ( $k$ ) as state variable instead of output ( $y$ ). The value function is defined as follows: $$ v(k) = \max_{0 \leq c \leq y + (1-\delta)k} \left\{ u(c) + \beta v(k') \right\} \tag{4} $$ subject to $$ k' = f(k) + (1-\delta)k - c $$ This formulation takes consumption ( $c$ ) as control variable. For the computation below we will use the next period capital ( $k'$ ) as control variable. Then (4) can be rewritten as follows: $$ v(k) = \max_{0 \leq k' \leq y + (1-\delta)k} \left\{ u\left(f(k) + (1-\delta)k - k'\right) + \beta v(k') \right\} \tag{5} $$ Essentially we converted the constrained maximization problem in (4) into the unconstrained maximization problem in (5). So I made my python code: import numpy as np import matplotlib.pyplot as plt from numba import njit, float64 from numba.experimental import jitclass from quantecon.optimize.scalar_maximization import brent_max from interpolation import interp opt_growth_data = [('α', float64), ('β', float64), ('γ', float64), ('δ', float64), ('grid',float64[:])] @jitclass(opt_growth_data) class OptimalGrowth_VI: def __init__(self, α=0.4, β=0.96, γ=2.0, δ=0.1, grid_max=10, grid_size=500): self.α, self.β, self.γ, self.δ = α, β, γ, δ self.grid = np.linspace(0.1, grid_max, grid_size) def f(self, k): return k**self.α def u(self, c): return c**(1 - self.γ) / (1 - self.γ) def objective(self, k, kp, v_array): f, u, β, δ = self.f, self.u, self.β, self.δ v = lambda x: interp(self.grid, v_array, x) return u(f(k)+(1-δ)*k-kp) + β*v(kp) @njit def T(v, og_VI): v_greedy = np.empty_like(v) v_new = np.empty_like(v) for i in range(len(og_VI.grid)): k = og_VI.grid[i] lower = 1e-10 upper = og_VI.f(k) + (1-og_VI.δ)*k result = brent_max(og_VI.objective, lower, upper, args=(k,v)) v_greedy[i], v_new[i] = result[0], result[1] return v_greedy, v_new def solve_model_VI(og_VI, tol=1e-4, max_iter=1000, print_skip=20): v = og_VI.grid error = tol+1 i=0 while i < max_iter and error > tol: v_greedy, v_new = T(v, og_VI) error = np.max(np.abs(v - v_new)) i += 1 if i % print_skip == 0: print(f"Error at iteration {i} is {error}.") v = v_new if i == max_iter: print("Failed to converge!") if i < max_iter: print(f"\nConverged in {i} iterations.") return v_greedy, v_new og_VI = OptimalGrowth_VI() v_greedy, v_solution = solve_model_VI(og_VI) plt.plot(og_VI.grid, v_greedy) I thought that I made my code correctly, but it produced weird optimal k's graph. Could you review my code and tell me what's wrong with my implementation? It's a nuisance, but any help would be really appreciated.

Cancel quote