@@ -330,7 +330,31 @@ function
330330
331331In our setting, we have the following key result
332332
333- * A feasible consumption policy is optimal if and only if it is :math: `v^*`-greedy
333+ * A feasible consumption policy is optimal if and only import numpy as np
334+ import matplotlib.pyplot as plt
335+ from interpolation import interp
336+ from numba import njit, prange
337+ from quantecon.optimize.scalar_maximization import brent_max
338+
339+
340+ def f(x):
341+ y1 = 2 * np.cos(6 * x) + np.sin(14 * x)
342+ return y1 + 2.5
343+
344+ def Af(x):
345+ return interp(c_grid, f(c_grid), x)
346+
347+ c_grid = np.linspace(0, 1, 6)
348+ f_grid = np.linspace(0, 1, 150)
349+
350+ fig, ax = plt.subplots(figsize=(10, 6))
351+
352+ ax.plot(f_grid, f(f_grid), 'b-', label='true function')
353+ ax.plot(f_grid, Af(f_grid), 'g-', label='linear approximation')
354+ ax.vlines(c_grid, c_grid * 0, f(c_grid), linestyle='dashed', alpha=0.5)
355+
356+ ax.set(xlim=(0, 1), ylim=(0, 6))
357+ plt.show() it is :math: `v^*`-greedy
334358
335359The intuition is similar to the intuition for the Bellman equation, which was
336360provided after :eq: `fpb30 `
@@ -524,35 +548,39 @@ We use an interpolation function from the
524548`interpolation.py package <https://github.com/EconForge/interpolation.py >`_
525549because it comes in handy later when we want to just-in-time compile our code
526550
551+ This library can be installed via `pip ` with the following command
552+
527553.. code-block :: python3
528554
529- import numpy as np
530- import matplotlib.pyplot as plt
531- from interpolation import interp
532- from numba import njit, prange
533- from quantecon.optimize.scalar_maximization import brent_max
534-
555+ !pip install interpolation
556+
557+ .. code-block :: python3
535558
536- def f(x):
537- y1 = 2 * np.cos(6 * x) + np.sin(14 * x)
538- return y1 + 2.5
559+ import numpy as np
560+ import matplotlib.pyplot as plt
561+ from interpolation import interp
562+ from numba import njit, prange
563+ from quantecon.optimize.scalar_maximization import brent_max
539564
540- c_grid = np.linspace(0, 1, 6)
541565
542- def Af(x):
543- return interp(c_grid, f(c_grid), x)
566+ def f(x):
567+ y1 = 2 * np.cos(6 * x) + np.sin(14 * x)
568+ return y1 + 2.5
544569
545- f_grid = np.linspace(0, 1, 150)
570+ def Af(x):
571+ return interp(c_grid, f(c_grid), x)
546572
547- fig, ax = plt.subplots(figsize=(10, 6))
548- ax.plot(f_grid, f(f_grid), 'b-', lw=2, alpha=0.8, label='true function')
549- ax.plot(f_grid, Af(f_grid), 'g-', lw=2, alpha=0.8,
550- label='linear approximation')
551- ax.vlines(c_grid, c_grid * 0, f(c_grid), linestyle='dashed', alpha=0.5)
552- ax.legend(loc='upper center')
553- ax.set(xlim=(0, 1), ylim=(0, 6))
573+ c_grid = np.linspace(0, 1, 6)
574+ f_grid = np.linspace(0, 1, 150)
554575
555- plt.show()
576+ fig, ax = plt.subplots(figsize=(10, 6))
577+
578+ ax.plot(f_grid, f(f_grid), 'b-', label='true function')
579+ ax.plot(f_grid, Af(f_grid), 'g-', label='linear approximation')
580+ ax.vlines(c_grid, c_grid * 0, f(c_grid), linestyle='dashed', alpha=0.5)
581+
582+ ax.set(xlim=(0, 1), ylim=(0, 6))
583+ plt.show()
556584
557585 Another advantage of piecewise linear interpolation is that it preserves useful shape properties such as monotonicity and concavity / convexity
558586
@@ -594,33 +622,48 @@ Here's a function that generates a Bellman operator using linear interpolation
594622.. code-block :: python3
595623
596624 def bellman_function_factory(og, parallel_flag=True):
597-
598- '''og is an instance of the OptimalGrowthModel'''
625+ """
626+ A function factory for building the Bellman operator, as well as
627+ a function that computes greedy policies.
628+
629+ Here og is an instance of OptimalGrowthModel.
630+ """
599631
600632 f, u = og.f, og.u
601633 y_grid, shocks = og.y_grid, og.shocks
602634
603635 @njit
604- def objective(c, w, y):
605- # Right hand side of Bellman equation
636+ def objective(c, w, y):
637+ """
638+ The right hand side of the Bellman equation
639+ """
640+ # First turn w into a function via interpolation
606641 w_func = lambda x: interp(y_grid, w, x)
607642 return u(c) + β * np.mean(w_func(f(y - c) * shocks))
608643
609644 @njit(parallel=parallel_flag)
610645 def T(w):
646+ """
647+ The Bellman operator
648+ """
611649 w_new = np.empty_like(w)
612650 for i in prange(len(y_grid)):
613651 y = y_grid[i]
614- w_max = brent_max(objective, 1e-10, y, args=(w, y))[1] # Solve for optimal w at y
652+ # Solve for optimal w at y
653+ w_max = brent_max(objective, 1e-10, y, args=(w, y))[1]
615654 w_new[i] = w_max
616655 return w_new
617-
656+
618657 @njit
619658 def get_greedy(v):
659+ """
660+ Computes the v-greedy policy of a given function v
661+ """
620662 σ = np.empty_like(v)
621663 for i in range(len(y_grid)):
622664 y = y_grid[i]
623- c_max = brent_max(objective, 1e-10, y, args=(v, y))[0] # Solve for optimal c at y
665+ # Solve for optimal c at y
666+ c_max = brent_max(objective, 1e-10, y, args=(v, y))[0]
624667 σ[i] = c_max
625668 return σ
626669
@@ -684,11 +727,15 @@ We will define functions to compute the closed form solutions to check our answe
684727.. code-block :: python3
685728
686729 def σ_star(y, α, β):
687- # True optimal policy
730+ """
731+ True optimal policy
732+ """
688733 return (1 - α * β) * y
689734
690735 def v_star(y, α, β, μ):
691- # True value function
736+ """
737+ True value function
738+ """
692739 c1 = np.log(1 - α * β) / (1 - β)
693740 c2 = (μ + α * np.log(α * β)) / (1 - α)
694741 c3 = 1 / (1 - β)
@@ -712,10 +759,27 @@ We first need to define a jitted version of the production function
712759
713760 @njit
714761 def f(k):
715- # Production function
762+ """
763+ Cobb-Douglas production function
764+ """
716765 return k**α
717766
767+ og = OptimalGrowthModel(f=f, u=np.log)
768+ T, get_greedy = bellman_function_factory(og)
769+
770+ Now we will create an instance of the model and assign it to the variable `og `
771+
772+ This instance will use the Cobb-Douglas production function and log utility
773+
774+ .. code-block :: python3
775+
718776 og = OptimalGrowthModel(f=f, u=np.log)
777+
778+ We will use `og ` to generate the Bellman operator and a function that computes
779+ greedy policies
780+
781+ .. code-block :: python3
782+
719783 T, get_greedy = bellman_function_factory(og)
720784
721785
@@ -851,7 +915,7 @@ The Policy Function
851915 single: Optimal Growth; Policy Function
852916
853917To compute an approximate optimal policy, we will use the second function
854- return from `bellman_function_factory ` that backs out the optimal policy
918+ returned from `bellman_function_factory ` that backs out the optimal policy
855919from the optimal wage rate
856920
857921The next figure compares the result to the exact solution, which, as mentioned
@@ -887,7 +951,7 @@ Once an optimal consumption policy :math:`\sigma` is given, income follows :eq:`
887951The next figure shows a simulation of 100 elements of this sequence for three different discount factors (and hence three different policies)
888952
889953.. figure :: /_static/figures/solution_og_ex2.png
890- :scale: 100 %
954+ :scale: 60 %
891955
892956In each sequence, the initial condition is :math: `y_0 = 0.1 `
893957
0 commit comments