Add infinite case (#14)

Co-authored-by: Michał Dobranowski <mdbrnowski@gmail.com>

example.jl

@@ -31,8 +31,6 @@ best_allocation = solve(my_model, 20, 0.2)
 # ╔═╡ e527fc57-60e0-44dc-9943-bc9095a2166f
 plot_allocation(best_allocation)
 
-
-
 # ╔═╡ 69167b33-034d-480b-9a06-193c11c61c49
 kss = RamseyGrowthModel.steady_state_K(my_model)
 
@@ -42,6 +40,12 @@ kss_best_allocation = solve(my_model, 100, kss)
 # ╔═╡ 47943756-3c8d-4241-ad2a-e02e55965685
 plot_allocation(kss_best_allocation, kss)
 
+# ╔═╡ 97b50f22-b883-47fe-a4f1-7c3d8d435d82
+infinite_best_allocation = solve(my_model, Inf, 3.1)
+
+# ╔═╡ e67e0207-8d63-4137-9d38-9f60a10091a5
+plot_allocation(infinite_best_allocation, kss)
+
 # ╔═╡ 00000000-0000-0000-0000-000000000001
 PLUTO_PROJECT_TOML_CONTENTS = """
 [deps]
@@ -1185,9 +1189,10 @@ version = "1.4.1+1"
 # ╠═500630d6-97ff-4266-ab70-6f3889097314
 # ╠═7354f6ed-fb29-4247-90c7-7718a3a395ff
 # ╠═e527fc57-60e0-44dc-9943-bc9095a2166f
-# ╠═41792cc7-80d5-4689-ab54-d94e76eb05c7
 # ╠═69167b33-034d-480b-9a06-193c11c61c49
 # ╠═77398e1c-55a3-4ae7-868a-cdb7cdf2a8f2
 # ╠═47943756-3c8d-4241-ad2a-e02e55965685
+# ╠═97b50f22-b883-47fe-a4f1-7c3d8d435d82
+# ╠═e67e0207-8d63-4137-9d38-9f60a10091a5
 # ╟─00000000-0000-0000-0000-000000000001
 # ╟─00000000-0000-0000-0000-000000000002

src/RamseyGrowthModel.jl

@@ -85,7 +85,7 @@ end
 
 function steady_state_K(model::GrowthModel)::Real
     f′(k::Real) = ForwardDiff.derivative(model.f, k)
-    find_zero(x -> f′(x) - 1/model.β + 1 - model.δ, (0, Inf64))
+    find_zero(x -> f′(x) - 1 / model.β + 1 - model.δ, (0, Inf64))
 end
 
 """
@@ -105,11 +105,20 @@ solve(model::GrowthModel, T::Integer, K₀::Real; kwargs...)
 
 The algorithm uses a binary search; if you want, you can override the default maximum number of iterations (`max_iter=1000`) or error tolerance (`tol=K₀/1e6`).
 """
-function solve(model::GrowthModel, T::Integer, K₀::Real; tol::Real=K₀ / 1e6, max_iter::Integer=1000)::DataFrame
+function solve(model::GrowthModel, T::Union{Integer,typeof(Inf)}, K₀::Real; tol::Real=K₀ / 1e6, max_iter::Integer=1000)::DataFrame
     K₀ > 0 || throw(DomainError("Initial capital `K₀` should be positive."))
     T > 0 || throw(DomainError("Time horizon `T` should be positive."))
+    T isa Integer || isinf(T) || throw(DomainError("Time horizon `T` should be an Integer or `Inf`."))
     C_low, C_high = 0, model.f(K₀)
 
+    K_ter = 0
+
+    if isinf(T)
+        @info "T specified as infinity; setting T to 100 for calculation and plotting purposes."
+        T = 100
+        K_ter = steady_state_K(model)
+    end
+
     if last(shooting(model, T + 1, K₀, C_high)).K > 0
         @info "Your initial capital is very high; the best allocation has been found using the highest possible initial consumption."
         return shooting(model, T, K₀, C_high)
@@ -119,13 +128,16 @@ function solve(model::GrowthModel, T::Integer, K₀::Real; tol::Real=K₀ / 1e6,
         C_mid = (C_low + C_high) / 2
         allocation = shooting(model, T, K₀, C_mid)
         next_K = next_K_C(model, last(allocation).K, last(allocation).C)[1]
-        if isnan(next_K)
-            C_high = C_mid
-        elseif next_K > tol
-            C_low = C_mid
-        else
+
+        error = next_K - K_ter
+
+        if abs(error) < tol
             @info "The best allocation has been found after $iter iterations."
             return allocation
+        elseif error > 0
+            C_low = C_mid
+        else
+            C_high = C_mid
         end
     end
 

test/runtests.jl

@@ -50,6 +50,8 @@ using DataFrames
         @test_throws DomainError solve(model, 20, 0)
         @test_throws DomainError solve(model, -1, K₀)
         @test_throws DomainError solve(model, 0, K₀)
+        @test_throws DomainError solve(model, 3.14, K₀)
+        @test_throws DomainError solve(model, -Inf, K₀)
 
         # check if both keyword arguments work
         @test (
@@ -81,6 +83,13 @@ using DataFrames
                 @test_logs (:info, info_pattern) solve(model_tech, 20, K₀)  # better techonology
             ) isa DataFrame
         end
+
+        @testset "infinite time" begin
+            model = GrowthModel(0.95, 0.02, 2.0, 0.3, 1)
+
+            @test_logs (:info, r"T specified as infinity;") (:info, info_pattern) solve(model, Inf, 3)
+            @test_logs (:info, r"T specified as infinity;") (:info, r"capital is very high") solve(model, Inf, 100)
+        end
     end
 
 end