@@ -330,7 +330,31 @@ function
330
330
331
331
In our setting, we have the following key result
332
332
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
334
358
335
359
The intuition is similar to the intuition for the Bellman equation, which was
336
360
provided after :eq: `fpb30 `
@@ -524,35 +548,39 @@ We use an interpolation function from the
524
548
`interpolation.py package <https://github.com/EconForge/interpolation.py >`_
525
549
because it comes in handy later when we want to just-in-time compile our code
526
550
551
+ This library can be installed via `pip ` with the following command
552
+
527
553
.. code-block :: python3
528
554
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
535
558
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
539
564
540
- c_grid = np.linspace(0, 1, 6)
541
565
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
544
569
545
- f_grid = np.linspace(0, 1, 150)
570
+ def Af(x):
571
+ return interp(c_grid, f(c_grid), x)
546
572
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)
554
575
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()
556
584
557
585
Another advantage of piecewise linear interpolation is that it preserves useful shape properties such as monotonicity and concavity / convexity
558
586
@@ -594,33 +622,48 @@ Here's a function that generates a Bellman operator using linear interpolation
594
622
.. code-block :: python3
595
623
596
624
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
+ """
599
631
600
632
f, u = og.f, og.u
601
633
y_grid, shocks = og.y_grid, og.shocks
602
634
603
635
@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
606
641
w_func = lambda x: interp(y_grid, w, x)
607
642
return u(c) + β * np.mean(w_func(f(y - c) * shocks))
608
643
609
644
@njit(parallel=parallel_flag)
610
645
def T(w):
646
+ """
647
+ The Bellman operator
648
+ """
611
649
w_new = np.empty_like(w)
612
650
for i in prange(len(y_grid)):
613
651
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]
615
654
w_new[i] = w_max
616
655
return w_new
617
-
656
+
618
657
@njit
619
658
def get_greedy(v):
659
+ """
660
+ Computes the v-greedy policy of a given function v
661
+ """
620
662
σ = np.empty_like(v)
621
663
for i in range(len(y_grid)):
622
664
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]
624
667
σ[i] = c_max
625
668
return σ
626
669
@@ -684,11 +727,15 @@ We will define functions to compute the closed form solutions to check our answe
684
727
.. code-block :: python3
685
728
686
729
def σ_star(y, α, β):
687
- # True optimal policy
730
+ """
731
+ True optimal policy
732
+ """
688
733
return (1 - α * β) * y
689
734
690
735
def v_star(y, α, β, μ):
691
- # True value function
736
+ """
737
+ True value function
738
+ """
692
739
c1 = np.log(1 - α * β) / (1 - β)
693
740
c2 = (μ + α * np.log(α * β)) / (1 - α)
694
741
c3 = 1 / (1 - β)
@@ -712,10 +759,27 @@ We first need to define a jitted version of the production function
712
759
713
760
@njit
714
761
def f(k):
715
- # Production function
762
+ """
763
+ Cobb-Douglas production function
764
+ """
716
765
return k**α
717
766
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
+
718
776
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
+
719
783
T, get_greedy = bellman_function_factory(og)
720
784
721
785
@@ -851,7 +915,7 @@ The Policy Function
851
915
single: Optimal Growth; Policy Function
852
916
853
917
To 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
855
919
from the optimal wage rate
856
920
857
921
The 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:`
887
951
The next figure shows a simulation of 100 elements of this sequence for three different discount factors (and hence three different policies)
888
952
889
953
.. figure :: /_static/figures/solution_og_ex2.png
890
- :scale: 100 %
954
+ :scale: 60 %
891
955
892
956
In each sequence, the initial condition is :math: `y_0 = 0.1 `
893
957
0 commit comments