Merge pull request #81 from QuantEcon/issue49

removed Dierckx: replaced Spline1D with CubicSpline from DataInterpolations

Author: jbrightuniverse

lectures/Manifest.toml

@@ -352,6 +352,12 @@ git-tree-sha1 = "55be907a531471de062f321147006c04b8c1e75f"
 uuid = "1313f7d8-7da2-5740-9ea0-a2ca25f37964"
 version = "0.8.0"
 
+[[DataInterpolations]]
+deps = ["ChainRulesCore", "LinearAlgebra", "Optim", "RecipesBase", "RecursiveArrayTools", "Reexport"]
+git-tree-sha1 = "7d235e13a5de66ffc0cb85446bcd3edb0556f992"
+uuid = "82cc6244-b520-54b8-b5a6-8a565e85f1d0"
+version = "3.6.0"
+
 [[DataStructures]]
 deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
 git-tree-sha1 = "4437b64df1e0adccc3e5d1adbc3ac741095e4677"
@@ -377,18 +383,6 @@ uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"
 deps = ["Mmap"]
 uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"
 
-[[Dierckx]]
-deps = ["Dierckx_jll"]
-git-tree-sha1 = "5fefbe52e9a6e55b8f87cb89352d469bd3a3a090"
-uuid = "39dd38d3-220a-591b-8e3c-4c3a8c710a94"
-version = "0.5.1"
-
-[[Dierckx_jll]]
-deps = ["CompilerSupportLibraries_jll", "Libdl", "Pkg"]
-git-tree-sha1 = "a580560f526f6fc6973e8bad2b036514a4e3b013"
-uuid = "cd4c43a9-7502-52ba-aa6d-59fb2a88580b"
-version = "0.0.1+0"
-
 [[DiffEqBase]]
 deps = ["ArrayInterface", "ChainRulesCore", "DataStructures", "DocStringExtensions", "FastBroadcast", "FunctionWrappers", "IterativeSolvers", "LabelledArrays", "LinearAlgebra", "Logging", "MuladdMacro", "NonlinearSolve", "Parameters", "Printf", "RecursiveArrayTools", "RecursiveFactorization", "Reexport", "Requires", "SciMLBase", "Setfield", "SparseArrays", "StaticArrays", "Statistics", "SuiteSparse", "ZygoteRules"]
 git-tree-sha1 = "9d312bb0b7c8ace440a71c64330cf1bea0ade0c8"

lectures/Project.toml

@@ -14,7 +14,7 @@ CategoricalArrays = "324d7699-5711-5eae-9e2f-1d82baa6b597"
 Compat = "34da2185-b29b-5c13-b0c7-acf172513d20"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 DataFramesMeta = "1313f7d8-7da2-5740-9ea0-a2ca25f37964"
-Dierckx = "39dd38d3-220a-591b-8e3c-4c3a8c710a94"
+DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0"
 DiffEqJump = "c894b116-72e5-5b58-be3c-e6d8d4ac2b12"
 DiffEqOperators = "9fdde737-9c7f-55bf-ade8-46b3f136cc48"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"

lectures/dynamic_programming_squared/amss.md

@@ -649,7 +649,7 @@ function get_policies_time1(T::BellmanEquation,
             init[i] = lb[i]
         end
     end
-    (minf, minx, ret) = optimize(opt, init)
+    (minf, minx, ret) = NLopt.optimize(opt, init)
     T.z0[i_x, s] = vcat(minx[1], minx[1] + G[s], minx[2:end])
     return vcat(-minf, T.z0[i_x, s])
 end
@@ -687,7 +687,7 @@ function get_policies_time0(T::BellmanEquation,
             init[i] = lb[i]
         end
     end
-    (minf, minx, ret) = optimize(opt, init)
+    (minf, minx, ret) = NLopt.optimize(opt, init)
     return vcat(-minf, vcat(minx[1], minx[1]+G[s0], minx[2:end]))
 end
 ```
@@ -961,8 +961,8 @@ assets, returning any excess revenues to the household as nonnegative lump sum t
 The recursive formulation is implemented as follows
 
 ```{code-cell} julia
-using Dierckx
 
+using DataInterpolations
 
 mutable struct BellmanEquation_Recursive{TP <: Model, TI <: Integer, TR <: Real}
     model::TP
@@ -1036,19 +1036,15 @@ function solve_time1_bellman(model::Model{TR}, μgrid::AbstractArray) where {TR
         xprimes = repeat(x, 1, S)
         xgrid[s_, :] = x
         for sprime = 1:S
-            splc = Spline1D(x[end:-1:1], c[:, sprime][end:-1:1], k=3)
-            spln = Spline1D(x[end:-1:1], n[:, sprime][end:-1:1], k=3)
-            splx = Spline1D(x[end:-1:1], xprimes[:, sprime][end:-1:1], k=3)
+            splc = CubicSpline(c[:, sprime][end:-1:1], x[end:-1:1])
+            spln = CubicSpline(n[:, sprime][end:-1:1], x[end:-1:1])
+            splx = CubicSpline(xprimes[:, sprime][end:-1:1], x[end:-1:1])
             cf[s_, sprime] = y -> splc(y)
             nf[s_, sprime] = y -> spln(y)
             xprimef[s_, sprime] = y -> splx(y)
-            # cf[s_, sprime] = LinInterp(x[end:-1:1], c[:, sprime][end:-1:1])
-            # nf[s_, sprime] = LinInterp(x[end:-1:1], n[:, sprime][end:-1:1])
-            # xprimef[s_, sprime] = LinInterp(x[end:-1:1], xprimes[:, sprime][end:-1:1])
         end
-        splV = Spline1D(x[end:-1:1], V[end:-1:1], k=3)
+        splV = CubicSpline(V[end:-1:1], x[end:-1:1])
         Vf[s_] = y -> splV(y)
-        # Vf[s_] = LinInterp(x[end:-1:1], V[end:-1:1])
     end
 
     policies = [cf, nf, xprimef]
@@ -1094,14 +1090,13 @@ function fit_policy_function(T::BellmanEquation_Recursive,
         for (i_x, x) in enumerate(xgrid)
             PFvec[:, i_x] = PF(i_x, x, s_)
         end
-        splV = Spline1D(xgrid, PFvec[1,:], k=3)
+        splV = CubicSpline(PFvec[1,:], xgrid)
         Vf[s_] = y -> splV(y)
-        # Vf[s_] = LinInterp(xgrid, PFvec[1, :])
         for sprime=1:S
-            splc = Spline1D(xgrid, PFvec[1 + sprime, :], k=3)
-            spln = Spline1D(xgrid, PFvec[1 + S + sprime, :], k=3)
-            splxprime = Spline1D(xgrid, PFvec[1 + 2S + sprime, :], k=3)
-            splTT = Spline1D(xgrid, PFvec[1 + 3S + sprime, :], k=3)
+            splc = CubicSpline(PFvec[1 + sprime, :], xgrid)
+            spln = CubicSpline(PFvec[1 + S + sprime, :], xgrid)
+            splxprime = CubicSpline(PFvec[1 + 2S + sprime, :], xgrid)
+            splTT = CubicSpline(PFvec[1 + 3S + sprime, :], xgrid)
             cf[s_, sprime] = y -> splc(y)
             nf[s_, sprime] = y -> spln(y)
             xprimef[s_, sprime] = y -> splxprime(y)
@@ -1281,7 +1276,7 @@ function get_policies_time1(T::BellmanEquation_Recursive,
     ftol_rel!(opt, 1e-8)
     ftol_abs!(opt, 1e-8)
 
-    (minf, minx, ret) = optimize(opt, init)
+    (minf, minx, ret) = NLopt.optimize(opt, init)
 
     if ret != :SUCCESS && ret != :ROUNDOFF_LIMITED && ret != :MAXEVAL_REACHED &&
         ret != :FTOL_REACHED && ret != :MAXTIME_REACHED
@@ -1356,7 +1351,7 @@ function get_policies_time0(T::BellmanEquation_Recursive,
     maxeval!(opt, 100000000)
     maxtime!(opt, 30)
 
-    (minf, minx, ret) = optimize(opt, init)
+    (minf, minx, ret) = NLopt.optimize(opt, init)
 
     if ret != :SUCCESS && ret != :ROUNDOFF_LIMITED && ret != :MAXEVAL_REACHED &&
         ret != :FTOL_REACHED
@@ -1681,12 +1676,12 @@ p
 tags: [remove-cell]
 ---
 @testset begin
-  #test sim_seq_long_plot[50, 3] ≈ 0.3951985593686047
-  #test sim_bel_long_plot[50, 3] ≈ 0.05684753244006188 atol = 1e-2
-  #test sim_seq_long_plot[100, 4] ≈ 0.340233842670859
-  #test sim_bel_long_plot[100, 4] ≈ 0.2093423366870517 atol = 1e-3
-  #test sim_seq_long_plot[200, 2] ≈ 0.5839693539786998
-  #test sim_bel_long_plot[200, 2] ≈ 0.6324036099550768 atol = 1e-3
+  @test sim_seq_long_plot[50, 3] ≈ 0.3951985593686047
+  @test sim_bel_long_plot[50, 3] ≈ 0.05684753244006188 atol = 1e-2
+  @test sim_seq_long_plot[100, 4] ≈ 0.340233842670859
+  @test sim_bel_long_plot[100, 4] ≈ 0.2093423366870517 atol = 1e-3
+  @test sim_seq_long_plot[200, 2] ≈ 0.5839693539786998
+  @test sim_bel_long_plot[200, 2] ≈ 0.6324036099550768 atol = 1e-3
 end
 ```