Presented by Arnav Sood and Jesse Perla

This notebook is a computational appendix to the earlier notebook on optimal stopping problems, and their formulations as linear complementarity problems (LCPs).

As before, we run some setup code

# if you are running locally for the first time, run ] add InstantiateFromURL using InstantiateFromURL activate_github_path("QuantEcon/SimpleDifferentialOperators.jl", path = "docs/examples", activate = true)

Updating registry at `~/.julia/registries/General` Updating git-repo `https://github.com/JuliaRegistries/General.git` piling project...

using SimpleDifferentialOperators using Plots, LinearAlgebra, Suppressor, BenchmarkTools

And define a set of model objects

μ = -0.1 σ = 0.1 ρ = 0.05 γ = 0.5 ρ = 0.05 S₀ = 20.0 u(x) = x^γ # (13) S(x) = S₀ # (15) @assert μ <= 0 # otherwise we need to swap the L₁₋bc operator

We can also define the boundary conditions and LCP matrices at this stage

M = 300 x̄ = range(0.0, 5.0, length = (M+2)) x = interiornodes(x̄) bc = (Reflecting(), Reflecting()) # assume x_t reflects at 0.0 and 5.0 L = I * ρ - μ*L₁₋bc(x̄, bc) - σ^2 / 2 * L₂bc(x̄, bc) # (2) and (14) q = -u.(x) + L*S.(x) # (9) w(z) = L*z + q; # (10)

This is the first solver we tried. To install it, simply run

# install if required # ] add PATHSolver

And to use it

using PATHSolver lb = zeros(M) ub = 300*ones(M) # Need to have upper bounds for the z. code, z, w_val = @btime @suppress solveLCP($w, $lb, $ub) # Solves (12)

5.461 ms (356 allocations: 2.44 MiB)

@assert code == :Solved # otherwise, an error in convergence

We can also try the MCP solver. We expect this to be slower, since we forgo specializations for the linear case.

code_, z_, w_val_ = @btime @suppress solveMCP($w, $lb, $ub) # Solves (12)

23.856 ms (1944 allocations: 33.55 MiB)

@assert code_ == :Solved # otherwise, an error in convergence

As a sanity check, we can plot

v = z + S.(x) plot(x, v, title = "Value Function for Optimal Stopping", legend = false, ylabel = "v(x)", xlabel = "x")

**Note:** As mentioned, for anything other than a rough estimate, the Julia code you're benchmarking should be inside of a function.

There is also the NLsolve.jl package, which offers a `mcpsolve()`

method for mixed complementarity problems.

# install if required # ] add NLsolve

And run

using NLsolve r = @btime @suppress mcpsolve($w, $lb, $ub, zeros(300), # initial condition reformulation = :smooth, # uses a so-called "Fischer function" to smooth out the problem autodiff = :forward, inplace = false, ftol = 1e-12);

203.362 ms (9519 allocations: 51.24 MiB)

@assert converged(r) == true # otherwise, an error in convergence

And the plot

v = r.zero + S.(x) plot(x, v, title = "Value Function for Optimal Stopping", legend = false, ylabel = "v(x)", xlabel = "x")

We can also use JuMP, Julia's DSL for numerical programming.

The solver is the Interior Point OPTimizer (IPOPT).

# install if required # ] add JuMP Ipopt

The setup here is a bit different, since we're using the quadratic formulation.

using JuMP, Ipopt m = Model(with_optimizer(Ipopt.Optimizer, tol=1e-12)); @variable(m, z[1:300]); @constraint(m, z .>= 0); @constraint(m, L*z + q .>= 0); @objective(m, Min, z'*L*z + z'*q); @suppress optimize!(m); termination_status(m)

LOCALLY_SOLVED::TerminationStatusCode = 4

And the plot

v = value.(z) + S.(x) plot(x, v, title = "Value Function for Optimal Stopping", legend = false, ylabel = "v(x)", xlabel = "x")

Note that we didn't compute a benchmark. As `optimize!()`

is in-place (that is, it modifies its arguments), we need to ensure that runs are actually independent.

Here's one approach.

function ipoptFunc() m = Model(with_optimizer(Ipopt.Optimizer, tol=1e-12)); @variable(m, z[1:300]); @constraint(m, z .>= 0); @constraint(m, L*z + q .>= 0); @objective(m, Min, z'*L*z + z'*q); @suppress optimize!(m); end @benchmark ipoptFunc()

BenchmarkTools.Trial: memory estimate: 3.61 MiB allocs estimate: 57669 -------------- minimum time: 103.205 ms (0.00% GC) median time: 121.977 ms (0.00% GC) mean time: 123.279 ms (0.91% GC) maximum time: 138.994 ms (5.32% GC) -------------- samples: 41 evals/sample: 1

In other words, we're benchmarking the process of building and solving a model from scratch.

Note that if you want more granular performance data, you can use the `@benchmark`

macro as above.

The flexibility of JuMP is that, after formulating the problem, swapping out solvers is fairly trivial.

Here's the same code, using the Oxford OSQP solver (Operator Splitting for Quadratic Programs).

As you can see, only one line needed to change.

# install if required # ] add OSQP

And run

using OSQP m = Model(with_optimizer(OSQP.MathOptInterfaceOSQP.Optimizer)) @variable(m, z[1:300]); @constraint(m, z .>= 0); @constraint(m, L*z + q .>= 0); @objective(m, Min, z'*L*z + z'*q); @suppress optimize!(m); @show termination_status(m)

termination_status(m) = ALMOST_OPTIMAL::TerminationStatusCode = 7 ALMOST_OPTIMAL::TerminationStatusCode = 7

And the plot

v = value.(z) + S.(x) plot(x, v, title = "Value Function for Optimal Stopping", legend = false, ylabel = "v(x)", xlabel = "x")

We can also write benchmarking code as above

function osqpFunc() m = Model(with_optimizer(OSQP.MathOptInterfaceOSQP.Optimizer)) @variable(m, z[1:300]); @constraint(m, z .>= 0); @constraint(m, L*z + q .>= 0); @objective(m, Min, z'*L*z + z'*q); @suppress optimize!(m); end @benchmark osqpFunc()

BenchmarkTools.Trial: memory estimate: 4.14 MiB allocs estimate: 61518 -------------- minimum time: 48.092 ms (0.00% GC) median time: 55.280 ms (0.00% GC) mean time: 56.038 ms (2.47% GC) maximum time: 69.130 ms (0.00% GC) -------------- samples: 90 evals/sample: 1

The outcomes of this notebook are:

Expository code for solving an LCP-formulated stopping problem using a variety of Julia tools.

The finding that using code optimized for the linear case leads to noticeable speed gains, and that

`PATHSolver.jl`

is a good bet for small problems like this one.Note: The size is crucial for this kind of analysis.

`JuMP`

(for example) is designed with massive problems in mind, and as such might be doing work that's costly here, but vital for the larger case.