# ContinuousDP: Stochastic Optimal Growth Model¶

In [1]:
using QuantEcon
using BasisMatrices
using ContinuousDPs
using PyPlot
using Random

In [2]:
alpha = 0.4
beta = 0.96
mu = 0
sigma = 0.1;

In [3]:
f(s, x) = log(x)
g(s, x, e) = (s - x)^alpha * e;

In [4]:
shock_size = 250
shocks = exp.(mu .+ sigma * randn(shock_size))
weights = fill(1/shock_size, shock_size);

In [5]:
grid_max = 4.
n = 30
s_min, s_max = 1e-5, grid_max
basis = Basis(ChebParams(n, s_min, s_max))

Out[5]:
1 dimensional Basis on the hypercube formed by (1.0e-5,) × (4.0,).
Basis families are Cheb

In [6]:
x_lb(s) = s_min
x_ub(s) = s;

In [7]:
ab = alpha * beta
c1 = log(1 - ab) / (1 - beta)
c2 = (mu + alpha * log(ab)) / (1 - alpha)
c3 = 1 / (1 - beta)
c4 = 1 / (1 - ab)

# True optimal policy
c_star(y) = (1 - alpha * beta) * y

# True value function
v_star(y) = c1 + c2 * (c3 - c4) + c4 * log(y);

In [8]:
cdp = ContinuousDP(f, g, beta, shocks, weights, x_lb, x_ub, basis);

In [9]:
@code_warntype ContinuousDP(f, g, beta, shocks, weights, x_lb, x_ub, basis)

Body::ContinuousDP{1,Array{Float64,1},Array{Float64,1},typeof(f),typeof(g),typeof(x_lb),typeof(x_ub)}
│  60 1 ─ %1 = invoke ContinuousDPs.Interp(_9::Basis{1,Tuple{ChebParams{Float64}}})::ContinuousDPs.Interp{1,Array{Float64,1},Array{Float64,2},LinearAlgebra.LU{Float64,Array{Float64,2}}}
│╻╷ Type61 │   %2 = %new(ContinuousDP{1,Array{Float64,1},Array{Float64,1},typeof(f),typeof(g),typeof(x_lb),typeof(x_ub)}, f, g, discount, shocks, weights, x_lb, x_ub, %1)::ContinuousDP{1,Array{Float64,1},Array{Float64,1},typeof(f),typeof(g),typeof(x_lb),typeof(x_ub)}
│  62 └──      return %2


## First test¶

In [10]:
C_star = cdp.interp.Phi \ v_star.(cdp.interp.S)
Tv = Array{Float64}(undef, cdp.interp.length)
C = copy(C_star)
bellman_operator!(cdp, C, Tv);

In [11]:
@code_warntype bellman_operator!(cdp, C, Tv)

Body::Array{Float64,1}
│╻ getproperty198 1 ─ %1 = (Base.getfield)(cdp, :interp)::ContinuousDPs.Interp{1,Array{Float64,1},TM,TL} where TL<:LinearAlgebra.Factorization where TM<:(AbstractArray{T,2} where T)
││    │   %2 = (Base.getfield)(%1, :S)::Array{Float64,1}
│     │   %3 = invoke ContinuousDPs.s_wise_max!(_2::ContinuousDP{1,Array{Float64,1},Array{Float64,1},typeof(f),typeof(g),typeof(x_lb),typeof(x_ub)}, %2::Array{Float64,1}, _3::Array{Float64,1}, _4::Array{Float64,1})::Array{Float64,1}
│╻ getproperty199 │   %4 = (Base.getfield)(cdp, :interp)::ContinuousDPs.Interp{1,Array{Float64,1},TM,TL} where TL<:LinearAlgebra.Factorization where TM<:(AbstractArray{T,2} where T)
││    │   %5 = (Base.getfield)(%4, :Phi_lu)::LinearAlgebra.Factorization
│     │        (ContinuousDPs.ldiv!)(C, %5, %3)
│ 200 └──      return C

In [12]:
grid_size = 200
grid_y = collect(range(s_min, stop=s_max, length=grid_size))
V_approx = funeval(C, cdp.interp.basis, grid_y);

In [13]:
fig, ax = subplots(figsize=(9, 5))
ax[:set_ylim](-35, -24)
ax[:plot](grid_y, V_approx, lw=2, alpha=0.6, label=L"$Tv^*$")
ax[:plot](grid_y, v_star.(grid_y), lw=2, alpha=0.6, label=L"$v^*$")
ax[:legend](loc="lower right")
show()

In [14]:
@time bellman_operator!(cdp, C, Tv)
@time bellman_operator!(cdp, C, Tv)
@time bellman_operator!(cdp, C, Tv);

  0.106662 seconds (1.21 M allocations: 54.788 MiB, 19.19% gc time)
0.083827 seconds (1.13 M allocations: 51.089 MiB, 11.62% gc time)
0.091252 seconds (1.20 M allocations: 54.406 MiB, 13.58% gc time)

In [15]:
s = 2.
@code_warntype ContinuousDPs._s_wise_max(cdp, s, C)

Body::Tuple{Float64,Float64}
│╻       getproperty151 1 ── %1  = (Base.getfield)(cdp, :shocks)::Array{Float64,1}
│╻       size    │    %2  = (Base.arraysize)(%1, 1)::Int64
││╻       Type    │    %3  = $(Expr(:foreigncall, :(:jl_alloc_array_2d), Array{Float64,2}, svec(Any, Int64, Int64), :(:ccall), 3, Array{Float64,2}, :(%2), 1, 1, :(%2)))::Array{Float64,2} │ 152 │ %4 = %new(getfield(ContinuousDPs, Symbol("#objective#3")){ContinuousDP{1,Array{Float64,1},Array{Float64,1},typeof(f),typeof(g),typeof(x_lb),typeof(x_ub)},Float64,Array{Float64,1},Array{Float64,2}}, cdp, s, C, %3)::getfield(ContinuousDPs, Symbol("#objective#3")){ContinuousDP{1,Array{Float64,1},Array{Float64,1},typeof(f),typeof(g),typeof(x_lb),typeof(x_ub)},Float64,Array{Float64,1},Array{Float64,2}} │ 161 │ %5 = Optim.optimize::Core.Compiler.Const(Optim.optimize, false) │╻ x_lb │ %6 = Main.s_min::Any │ │ %7 = (isa)(%6, Float64)::Bool │ └─── goto #12 if not %7 │ 2 ── %9 = π (%6, Float64) ││╻ sqrt │ %10 = (Base.Math.sqrt_llvm)(2.220446049250313e-16)::Float64 ││╻ #optimize#78 │ %11 = π (1, Int64) │││ │ %12 = π (%11, Int64) │││ │ %13 = π (false, Bool) │││ │ %14 = π (%12, Int64) │││ │ invoke Core.kwfunc(Optim.optimize::Any) │││ │ %16 = Optim.optimize::typeof(Optim.optimize) │││╻╷╷╷╷ #optimize │ %17 = (Base.slt_int)(0, 1)::Bool ││││┃│││ isempty └─── goto #4 if not %17 │││││┃││ iterate 3 ── goto #5 ││││││┃│ iterate 4 ── invoke Base.getindex(()::Tuple{}, 1::Int64) │││││││┃ iterate └───$(Expr(:unreachable))
│││││││     5 ┄─       goto #6
│││││╻       iterate    6 ──       goto #7
│││││       7 ──       goto #8
││││        8 ── %25 = invoke Optim.:(#optimize#71)(%10::Float64, 2.220446049250313e-16::Float64, 1000::Int64, false::Bool, %13::Bool, nothing::Nothing, %14::Int64, false::Bool, %16::typeof(Optim.optimize), %4::getfield(ContinuousDPs, Symbol("#objective#3")){ContinuousDP{1,Array{Float64,1},Array{Float64,1},typeof(f),typeof(g),typeof(x_lb),typeof(x_ub)},Float64,Array{Float64,1},Array{Float64,2}}, %9::Float64, _3::Float64, $(QuoteNode(Optim.Brent()))::Optim.Brent)::Optim.UnivariateOptimizationResults{Float64,Float64,_1,_2,Float64,Optim.Brent} where _2 where _1 ││││ └─── goto #9 │││ 9 ── goto #10 ││ 10 ─ goto #11 │ 11 ─ goto #13 │ 12 ─ %30 = (%5)(%4, %6, s)::Optim.UnivariateOptimizationResults{Float64,Float64,_1,_2,Float64,Optim.Brent} where _2 where _1 │ └─── goto #13 │ 13 ┄ %32 = φ (#11 => %25, #12 => %30)::Optim.UnivariateOptimizationResults{Float64,Float64,_1,_2,Float64,Optim.Brent} where _2 where _1 │╻ getproperty162 │ %33 = (Base.getfield)(%32, :minimum)::Any │ │ (Core.typeassert)(%33, ContinuousDPs.Float64) │ │ %35 = π (%33, Float64) │╻ - │ %36 = (Base.neg_float)(%35)::Float64 │╻ getproperty163 │ %37 = (Base.getfield)(%32, :minimizer)::Any │ │ (Core.typeassert)(%37, ContinuousDPs.Float64) │ │ %39 = π (%37, Float64) │ 164 │ %40 = (Core.tuple)(%36, %39)::Tuple{Float64,Float64} │ └─── return %40  In [16]: @time ContinuousDPs._s_wise_max(cdp, s, C) @time ContinuousDPs._s_wise_max(cdp, s, C) @time ContinuousDPs._s_wise_max(cdp, s, C);   0.002056 seconds (28.31 k allocations: 1.281 MiB) 0.002404 seconds (28.29 k allocations: 1.278 MiB) 0.002289 seconds (28.29 k allocations: 1.278 MiB)  In [17]: v_init_func(s) = 5 * log(s) w = v_init_func.(grid_y) n = 35 fig, ax = subplots(figsize=(9, 6)) ax[:set_ylim](-50, 10) ax[:set_xlim](minimum(grid_y), maximum(grid_y)) lb = "initial condition" jet = ColorMap("jet") ax[:plot](grid_y, w, color=jet(0), lw=2, alpha=0.6, label=lb) S = cdp.interp.S V = v_init_func.(S) for i in 1:n C = cdp.interp.Phi \ V bellman_operator!(cdp, C, V) w = funeval(C, cdp.interp.basis, grid_y) ax[:plot](grid_y, w, color=jet(i / n), lw=2, alpha=0.6) end lb = "true value function" ax[:plot](grid_y, v_star.(grid_y), "k-", lw=2, alpha=0.8, label=lb) ax[:legend](loc="lower right") show()  ## Solve by policy iteration¶ In [18]: res = solve(cdp);  Compute iterate 6 with error 1.0658141036401503e-13 Converged in 6 steps  In [19]: @code_warntype solve(cdp)  Body::ContinuousDPs.CDPSolveResult{PFI,1,Array{Float64,1},Array{Float64,1}} │ 296 1 ─ %1 = ContinuousDPs.PFI::Core.Compiler.Const(PFI, false) │╻ solve │ %2 = (Base.Math.sqrt_llvm)(2.220446049250313e-16)::Float64 ││ │ %3 = invoke ContinuousDPs.:(#solve#7)(%2::Float64, 500::Int64, 2::Int64, 50::Int64, _1::Function, _2::ContinuousDP{1,Array{Float64,1},Array{Float64,1},typeof(f),typeof(g),typeof(x_lb),typeof(x_ub)}, %1::Type{PFI})::ContinuousDPs.CDPSolveResult{PFI,1,Array{Float64,1},Array{Float64,1}} │ └── return %3  In [20]: @time res = solve(cdp) @time res = solve(cdp) @time res = solve(cdp);  Compute iterate 6 with error 1.0658141036401503e-13 Converged in 6 steps 0.909546 seconds (11.46 M allocations: 571.261 MiB, 11.57% gc time) Compute iterate 6 with error 1.0658141036401503e-13 Converged in 6 steps 0.914953 seconds (11.46 M allocations: 571.270 MiB, 11.24% gc time) Compute iterate 6 with error 1.0658141036401503e-13 Converged in 6 steps 0.910241 seconds (11.46 M allocations: 571.270 MiB, 11.05% gc time)  In [21]: set_eval_nodes!(res, grid_y);  In [22]: @code_warntype set_eval_nodes!(res, grid_y)  Body::ContinuousDPs.CDPSolveResult{PFI,1,Array{Float64,1},Array{Float64,1}} │╻╷ set_eval_nodes!137 1 ─ %1 = (Base.getfield)(s_nodes_coord, 1, true)::Array{Float64,1} ││╻ setproperty! │ (Base.setfield!)(res, :eval_nodes, %1) ││╻ setproperty! │ (Base.setfield!)(res, :eval_nodes_coord, s_nodes_coord) │╻ set_eval_nodes! │ %4 = invoke ContinuousDPs.evaluate!(_2::ContinuousDPs.CDPSolveResult{PFI,1,Array{Float64,1},Array{Float64,1}})::ContinuousDPs.CDPSolveResult{PFI,1,Array{Float64,1},Array{Float64,1}} │ └── return %4  In [23]: @time set_eval_nodes!(res, grid_y) @time set_eval_nodes!(res, grid_y) @time set_eval_nodes!(res, grid_y);   0.575157 seconds (7.74 M allocations: 349.491 MiB, 11.62% gc time) 0.573510 seconds (7.74 M allocations: 349.491 MiB, 10.07% gc time) 0.592195 seconds (7.74 M allocations: 349.491 MiB, 9.38% gc time)  In [24]: fig, ax = subplots(figsize=(9, 5)) ax[:set_ylim](-35, -24) ax[:plot](grid_y, res.V, lw=2, alpha=0.6, label="approximate value function") ax[:plot](grid_y, v_star.(grid_y), lw=2, alpha=0.6, label="true value function") ax[:legend](loc="lower right") show()  In [25]: fig, ax = subplots(figsize=(9, 5)) ax[:plot](grid_y, res.X, lw=2, alpha=0.6, label="approximate policy function") ax[:plot](grid_y, c_star.(grid_y), lw=2, alpha=0.6, label="true policy function") ax[:legend](loc="lower right") show()  In [26]: fig, ax = subplots(figsize=(9, 5)) ax[:plot](grid_y, res.resid, lw=2, alpha=0.6, label="residual") show()  ## Simulate the controlled Markov process¶ In [27]: s_init = 0.1 ts_length = 100 y = simulate(res, s_init, ts_length) fig, ax = subplots(figsize=(9, 6)) ax[:plot](1:ts_length, y, lw=2, alpha=0.6, label="beta =$(cdp.discount)" )
ax[:legend](loc="lower right")
show()

In [28]:
fig, ax = subplots(figsize=(9, 6))

for beta in (0.9, 0.94, 0.98)
cdp.discount = beta
res = solve(cdp, verbose=0)
set_eval_nodes!(res, grid_y)
y = simulate(res, s_init, ts_length)
ax[:plot](1:ts_length, y, lw=2, alpha=0.6, label="beta = $(cdp.discount)" ) end ax[:legend](loc="lower right") show()  In [29]: @code_warntype simulate!(MersenneTwister(0), Array{Float64}(undef, ts_length), res, s_init)  Body::Array{Float64,1} │╻ size335 1 ── %1 = (Base.arraysize)(s_path, 1)::Int64 │╻ getproperty336 │ %2 = (Base.getfield)(res, :cdp)::ContinuousDP{1,Array{Float64,1},Array{Float64,1},Tf,Tg,Tlb,Tub} where Tub<:Function where Tlb<:Function where Tg<:Function where Tf<:Function ││ │ %3 = (Base.getfield)(%2, :weights)::Array{Float64,1} │╻ cumsum │ invoke Core.kwfunc(Base.cumsum::Any) ││ │ %5 = Base.cumsum::typeof(cumsum) ││╻╷╷╷╷ #cumsum │ %6 = (Base.slt_int)(0, 1)::Bool │││┃│││ isempty └─── goto #3 if not %6 ││││┃││ iterate 2 ── goto #4 │││││┃│ iterate 3 ── invoke Base.getindex(()::Tuple{}, 1::Int64) ││││││┃ iterate └───$(Expr(:unreachable))
││││││        4 ┄─        goto #5
││││╻         iterate    5 ──        goto #6
││││          6 ──        goto #7
│││           7 ── %14  = invoke Base.:(#cumsum#572)(1::Int64, %5::Function, %3::Array{Float64,1})::Array{Float64,1}
│││           └───        goto #8
││            8 ──        goto #9
│╻         -337 9 ── %17  = (Base.sub_int)(%1, 1)::Int64
││╻╷╷╷      rand    │    %18  = $(Expr(:foreigncall, :(:jl_alloc_array_1d), Array{Float64,1}, svec(Any, Int64), :(:ccall), 2, Array{Float64,1}, :(%17), :(%17)))::Array{Float64,1} │││╻╷ rand! │ %19 = (Base.arraylen)(%18)::Int64 ││││╻ rand! │ %20 = (Base.mul_int)(8, %19)::Int64 │││││╻ _rand! │ %21 = (Base.arraylen)(%18)::Int64 ││││││╻ * │ %22 = (Base.mul_int)(8, %21)::Int64 ││││││╻ <= │ %23 = (Base.sle_int)(%20, %22)::Bool ││││││ └─── goto #11 if not %23 ││││││ 10 ─ goto #12 │ 11 ─ nothing ││││││ 12 ┄ %27 = φ (#10 => true, #11 => false)::Bool ││││││ └─── goto #14 if not %27 ││││││╻ macro expansion 13 ─ %29 =$(Expr(:gc_preserve_begin, :(%18)))
│││││││╻╷        pointer    │    %30  = $(Expr(:foreigncall, :(:jl_array_ptr), Ptr{Float64}, svec(Any), :(:ccall), 1, :(%18)))::Ptr{Float64} │││││││╻ Type │ %31 = %new(Random.UnsafeView{Float64}, %30, %19)::Random.UnsafeView{Float64} │││││││ │ invoke Random.rand!(_2::MersenneTwister, %31::Random.UnsafeView{Float64},$(QuoteNode(Random.SamplerTrivial{Random.CloseOpen01{Float64},Float64}(Random.CloseOpen01{Float64}())))::Random.SamplerTrivial{Random.CloseOpen01{Float64},Float64})
│││││││       │           $(Expr(:gc_preserve_end, :(%29))) │││││││ └─── goto #15 ││││││╻ Type 14 ─ %35 = %new(Core.AssertionError, "sizeof(Float64) * n64 <= sizeof(T) * length(A) && isbitstype(T)")::AssertionError ││││││ │ (Base.throw)(%35) ││││││ └───$(Expr(:unreachable))
│││││         15 ┄        goto #16
││││          16 ─        goto #17
│││           17 ─        goto #18
││            18 ─        goto #19
│╻         -338 19 ─ %42  = (Base.sub_int)(%1, 1)::Int64
││╻         Type    │    %43  = $(Expr(:foreigncall, :(:jl_alloc_array_1d), Array{Int64,1}, svec(Any, Int64), :(:ccall), 2, Array{Int64,1}, :(%42), :(%42)))::Array{Int64,1} │╻ -339 │ %44 = (Base.sub_int)(%1, 1)::Int64 ││╻╷╷╷ Type │ %45 = (Base.sle_int)(1, %44)::Bool │││╻ unitrange_last │ (Base.sub_int)(%44, 1) ││││ │ %47 = (Base.ifelse)(%45, %44, 0)::Int64 ││╻╷╷ isempty │ %48 = (Base.slt_int)(%47, 1)::Bool ││ └─── goto #21 if not %48 ││ 20 ─ goto #22 ││ 21 ─ goto #22 │ 22 ┄ %52 = φ (#20 => true, #21 => false)::Bool │ │ %53 = φ (#21 => 1)::Int64 │ │ %54 = φ (#21 => 1)::Int64 │ │ %55 = (Base.not_int)(%52)::Bool │ └─── goto #38 if not %55 │ 23 ┄ %57 = φ (#22 => %53, #37 => %92)::Int64 │ │ %58 = φ (#22 => %54, #37 => %93)::Int64 │╻ getindex340 │ %59 = (Base.arrayref)(true, %18, %57)::Float64 ││╻╷╷╷╷ #searchsortedlast#5 │ %60 = (Base.arraysize)(%14, 1)::Int64 │││╻╷╷╷ searchsortedlast │ %61 = (Base.slt_int)(%60, 0)::Bool ││││┃│││││ axes │ %62 = (Base.ifelse)(%61, 0, %60)::Int64 ││││╻╷ searchsortedlast └─── %63 = (Base.add_int)(%62, 1)::Int64 ││││╻ searchsortedlast 24 ┄ %64 = φ (#23 => 0, #28 => %78)::Int64 │││││ │ %65 = φ (#23 => %63, #28 => %79)::Int64 │││││╻ - │ %66 = (Base.sub_int)(%65, 1)::Int64 │││││╻ < │ %67 = (Base.slt_int)(%64, %66)::Bool │││││ └─── goto #29 if not %67 │││││╻ + 25 ─ %69 = (Base.add_int)(%64, %65)::Int64 ││││││╻ >>> │ %70 = (Base.lshr_int)(%69, 0x0000000000000001)::Int64 ││││││╻ << │ %71 = (Base.shl_int)(%69, 0xffffffffffffffff)::Int64 ││││││ │ %72 = (Base.ifelse)(true, %70, %71)::Int64 │││││╻ getindex │ %73 = (Base.arrayref)(false, %14, %72)::Float64 ││││││╻ isless │ %74 = (Base.fpislt)(%59, %73)::Bool │││││ └─── goto #27 if not %74 │││││ 26 ─ goto #28 │ 27 ─ nothing │││││ 28 ┄ %78 = φ (#26 => %64, #27 => %72)::Int64 │││││ │ %79 = φ (#26 => %72, #27 => %65)::Int64 │││││ └─── goto #24 │││││ 29 ─ goto #30 ││││ 30 ─ goto #31 │││ 31 ─ goto #32 ││ 32 ─ goto #33 │╻ + 33 ─ %85 = (Base.add_int)(%64, 1)::Int64 │╻ setindex! │ (Base.arrayset)(true, %43, %85, %57) ││╻ == │ %87 = (%58 === %47)::Bool ││ └─── goto #35 if not %87 ││ 34 ─ goto #36 ││╻ + 35 ─ %90 = (Base.add_int)(%58, 1)::Int64 │╻ iterate └─── goto #36 │ 36 ┄ %92 = φ (#35 => %90)::Int64 │ │ %93 = φ (#35 => %90)::Int64 │ │ %94 = φ (#34 => true, #35 => false)::Bool │ │ %95 = (Base.not_int)(%94)::Bool │ └─── goto #38 if not %95 │ 37 ─ goto #23 │╻ getproperty343 38 ─ %98 = (Base.getfield)(res, :eval_nodes_coord)::Tuple{Array{Float64,1}} ││╻ getindex │ %99 = (Base.getfield)(%98, 1, true)::Array{Float64,1} ││╻ Type │ %100 = invoke LinParams{Array{Float64,1}}(%99::Array{Float64,1}, 0::Int64)::LinParams{Array{Float64,1}} ││ │ %101 = (Core.tuple)(%100)::Tuple{LinParams{Array{Float64,1}}} ││╻╷╷ _Basis │ %102 = %new(Basis{1,Tuple{LinParams{Array{Float64,1}}}}, %101)::Basis{1,Tuple{LinParams{Array{Float64,1}}}} │╻ getproperty344 │ %103 = (Base.getfield)(res, :X)::Array{Float64,1} ││╻ Type │ %104 = invoke BasisMatrices.BasisMatrix(BasisMatrices.Nothing::Type{Nothing}, %102::Basis{1,Tuple{LinParams{Array{Float64,1}}}},$(QuoteNode(Tensor()))::Tensor)::BasisMatrix{Tensor,SparseArrays.SparseMatrixCSC{Float64,Int64}}
││╻         Type    │    %105 = invoke BasisMatrices.get_coefs(%102::Basis{1,Tuple{LinParams{Array{Float64,1}}}}, %104::BasisMatrix{Tensor,SparseArrays.SparseMatrixCSC{Float64,Int64}}, %103::Array{Float64,1})::Array{Float64,1}
│╻╷        axes346 │    %106 = (Base.arraysize)(s_path, 1)::Int64
││╻╷╷╷      map    │    %107 = (Base.slt_int)(%106, 0)::Bool
│││┃││       Type    │           (Base.ifelse)(%107, 0, %106)
│╻         getproperty347 │    %109 = (Base.getfield)(res, :cdp)::ContinuousDP{1,Array{Float64,1},Array{Float64,1},Tf,Tg,Tlb,Tub} where Tub<:Function where Tlb<:Function where Tg<:Function where Tf<:Function
││            │    %110 = (Base.getfield)(%109, :shocks)::Array{Float64,1}
││╻         size    │    %111 = (Base.arraysize)(%110, 1)::Int64
│││╻╷╷╷      Type    │    %112 = (Base.slt_int)(%111, 0)::Bool
││││┃│        Type    │           (Base.ifelse)(%112, 0, %111)
│╻         setindex!348 │           (Base.arrayset)(true, s_path, s_init, 1)
│╻         -349 │    %115 = (Base.sub_int)(%1, 1)::Int64
││╻╷╷╷      Type    │    %116 = (Base.sle_int)(1, %115)::Bool
│││╻         unitrange_last    │           (Base.sub_int)(%115, 1)
││││          │    %118 = (Base.ifelse)(%116, %115, 0)::Int64
││╻╷╷       isempty    │    %119 = (Base.slt_int)(%118, 1)::Bool
││            └───        goto #40 if not %119
││            39 ─        goto #41
││            40 ─        goto #41
│             41 ┄ %123 = φ (#39 => true, #40 => false)::Bool
│             │    %124 = φ (#40 => 1)::Int64
│             │    %125 = φ (#40 => 1)::Int64
│             │    %126 = (Base.not_int)(%123)::Bool
│             └───        goto #47 if not %126
│             42 ┄ %128 = φ (#41 => %124, #46 => %149)::Int64
│             │    %129 = φ (#41 => %125, #46 => %150)::Int64
│╻         getindex350 │    %130 = (Base.arrayref)(true, s_path, %128)::Float64
│╻╷╷╷╷╷╷   Interpoland351 │    %131 = \$(Expr(:foreigncall, :(:jl_alloc_array_2d), Array{Float64,2}, svec(Any, Int64, Int64), :(:ccall), 3, Array{Float64,2}, 1, 1, 1, 1))::Array{Float64,2}
││┃│││      Interpoland    │    %132 = invoke Base.fill!(%131::Array{Float64,2}, %130::Float64)::Array{Float64,2}
│││┃         funeval    │    %133 = invoke BasisMatrices.funeval(%105::Array{Float64,1}, %102::Basis{1,Tuple{LinParams{Array{Float64,1}}}}, %132::Array{Float64,2}, 0::Int64)::Array{Float64,2}
││││╻         getindex    │    %134 = (Base.arrayref)(true, %133, 1)::Float64
│╻         getindex352 │    %135 = (Base.arrayref)(true, %43, %128)::Int64
│╻         getproperty    │    %136 = (Base.getfield)(res, :cdp)::ContinuousDP{1,Array{Float64,1},Array{Float64,1},Tf,Tg,Tlb,Tub} where Tub<:Function where Tlb<:Function where Tg<:Function where Tf<:Function
││            │    %137 = (Base.getfield)(%136, :shocks)::Array{Float64,1}
│╻         getindex    │    %138 = (Base.arrayref)(true, %137, %135)::Float64
│╻         +353 │    %139 = (Base.add_int)(%128, 1)::Int64
│╻         getproperty    │    %140 = (Base.getfield)(res, :cdp)::ContinuousDP{1,Array{Float64,1},Array{Float64,1},Tf,Tg,Tlb,Tub} where Tub<:Function where Tlb<:Function where Tg<:Function where Tf<:Function
││            │    %141 = (Base.getfield)(%140, :g)::Function
│             │    %142 = (%141)(%130, %134, %138)::Any
│             │           (Base.setindex!)(s_path, %142, %139)
││╻         ==    │    %144 = (%129 === %118)::Bool
││            └───        goto #44 if not %144
││            43 ─        goto #45
││╻         +    44 ─ %147 = (Base.add_int)(%129, 1)::Int64
│╻         iterate    └───        goto #45
│             45 ┄ %149 = φ (#44 => %147)::Int64
│             │    %150 = φ (#44 => %147)::Int64
│             │    %151 = φ (#43 => true, #44 => false)::Bool
│             │    %152 = (Base.not_int)(%151)::Bool
│             └───        goto #47 if not %152
│             46 ─        goto #42
│         356 47 ─        return s_path

In [30]:
@time simulate!(MersenneTwister(0), Array{Float64}(undef, ts_length), res, s_init)
@time simulate!(MersenneTwister(0), Array{Float64}(undef, ts_length), res, s_init)
@time simulate!(MersenneTwister(0), Array{Float64}(undef, ts_length), res, s_init);

  0.000365 seconds (3.56 k allocations: 291.109 KiB)
0.000315 seconds (3.55 k allocations: 289.953 KiB)
0.000329 seconds (3.55 k allocations: 289.953 KiB)

In [ ]: