Notebook

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 [ ]: