From eea3b1f9d5f78d31d8c887e8164d8646175e1fb6 Mon Sep 17 00:00:00 2001 From: ashutosh-b-b Date: Tue, 23 Jan 2024 16:00:53 +0530 Subject: [PATCH 01/37] feat: Add NNStopping for Optimal Stopping Problems --- src/NNStopping.jl | 112 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 112 insertions(+) create mode 100644 src/NNStopping.jl diff --git a/src/NNStopping.jl b/src/NNStopping.jl new file mode 100644 index 0000000..6609bca --- /dev/null +++ b/src/NNStopping.jl @@ -0,0 +1,112 @@ +struct NNStopping{M, O} <: HighDimPDEAlgorithm + m::M + opt::O + + function NNStopping(m, opt = Flux.Optimise.Adam(0.1)) + new{typeof(m), typeof(opt)}(m, opt) + end +end + +function DiffEqBase.solve(prob::PIDEProblem, + pdealg::NNStopping, + sdealg; + verbose = false, + maxiters = 300, + trajectories = 100, + dt = eltype(prob.tspan)(0), + ensemblealg = EnsembleThreads(), + kwargs...) + + tspan = prob.tspan + sigma = prob.g + μ = prob.f + g = prob.kwargs.data.g + u0 = prob.u0 + ts = tspan[1]:dt:tspan[2] + T = tspan[2] + + m = alg.m + opt = alg.opt + + + prob = SDEProblem(μ,sigma,u0,tspan) + ensembleprob = EnsembleProblem(prob) + sim = solve(ensembleprob, sdealg, ensemblealg, dt=dt,trajectories=trajectories,adaptive=false) + payoff = [] + times = [] + iter = 0 + # for u in sim.u + # un = [] + # function Un(n , X ) + # if size(un)[1] >= n + # return un[n] + # else + # if(n == 1) + # ans = first(m(X[1])[1]) + # un = [ans] + # return ans + # else + # ans = max(first(m(X[n])[n]) , n + 1 - size(ts)[1])*(1 - sum(Un(i , X ) for i in 1:n-1)) + # un = vcat( un , ans) + # return ans + # end + # end + # end + + function Un!(un_cache, n, X) + !isnan(un_cache[n]) && return un_cache[n] + if n == 1 + un_cache[n] = first(m(X[1])[1]) + else + un_cache[n] = max(first(m(X[n])[n]) , n + 1 - length(ts))*(1 - sum(un_cache[1:n-1])) + end + return un_cache[n] + end + + + function loss() + reward = 0.00 + map(sim.u) do X + un_cache = fill(NaN, length(ts)) + reward += mapreduce((x, t) -> Un!(un_cache, i, X)*g(t,X), +, X, ts) + end + + return 10000 - reward + end + + dataset = Iterators.repeated(() , maxiters) + + callback = function () + l = loss() + println("Current loss is: $l") + end + + Flux.train!(loss, Flux.params(m), dataset, opt; cb = callback) + + Usum = 0 + ti = 0 + Xt = sim.u[1].u + un_cache = fill(NaN, N) + for i in 1:N + Un = Un!(un_cache, i , Xt) + Usum = Usum + Un + if Usum >= 1 - Un + ti = i + break + end + end + for u in sim.u + X = u.u + price = g(ts[ti] , X[ti]) + payoff = vcat(payoff , price) + times = vcat(times, ti) + iter = iter + 1 + # println("SUM : $sump") + # println("TIME : $ti") + end + sum(payoff)/size(payoff)[1] + + + + +end \ No newline at end of file From cfd5b454806f04fa0cd2947f30222fda733c2b0b Mon Sep 17 00:00:00 2001 From: ashutosh-b-b Date: Tue, 23 Jan 2024 16:01:23 +0530 Subject: [PATCH 02/37] test: add tests for NNStopping --- test/NNStopping.jl | 49 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 49 insertions(+) create mode 100644 test/NNStopping.jl diff --git a/test/NNStopping.jl b/test/NNStopping.jl new file mode 100644 index 0000000..e46a97a --- /dev/null +++ b/test/NNStopping.jl @@ -0,0 +1,49 @@ +using Test, Flux , StochasticDiffEq , LinearAlgebra +println("Optimal Stopping Time Test") +using NeuralPDE + +using Random +Random.seed!(100) + +d = 1 +r = 0.04f0 +beta = 0.2f0 +T = 1 +u0 = fill(80.00 , d , 1) +f(du,u,p,t) = (du .= r*u) +sigma(du,u,p,t) = (du .= Diagonal(beta*u)) +tspan = (0.0 , 1.0) +N = 50 +dt = tspan[2]/49 +K = 100.00 +function g(t , x) + return exp(-r*t)*(max(K - maximum(x) , 0)) +end + +prob = SDEProblem(f , sigma , u0 , tspan ; g = g) +opt = Flux.ADAM(0.1) +m = Chain(Dense(d , 5, tanh), Dense(5, 16 , tanh) , Dense(16 , N ), softmax) +sol = solve(prob, NeuralPDE.NNStopping(m, opt), verbose = true, dt = dt, sde_alg = EM() + abstol=1e-6, maxiters = 20 , trajectories = 200) + +##Analytical Binomial Tree approach for American Options +function BinomialTreeAM1D(S0 , N , r , beta) + V = zeros(N+1) + dT = T/N + u = exp(beta*sqrt(dT)) + d = 1/u + S_T = [S0*(u^j)* (d^(N-j)) for j in 0:N] + a = exp(r*dT) + p = (a - d)/(u - d) + q = 1.0 - p + V = [max(K - x , 0) for x in S_T] + for i in N-1:-1:0 + V[1:end-1] = exp(-r*dT).*(p*V[2:end] + q*V[1:end-1]) + S_T = S_T*u + V = [max(K - S_T[i] , V[i]) for i in 1:size(S_T)[1]] + end + return V[1] +end +real_sol = BinomialTreeAM1D(u0[1] , N , r , beta) +error = abs(sol - real_sol) +@test error < 0.5 \ No newline at end of file From 2c516a3dd7c2a87e8b2fbfc38038bd0f52462723 Mon Sep 17 00:00:00 2001 From: ashutosh-b-b Date: Sat, 27 Jan 2024 00:33:25 +0530 Subject: [PATCH 03/37] feat: add `NNStopping` for OptimalStopping Problem Refactors old code from NeuralPDE --- src/NNStopping.jl | 169 ++++++++++++++++++++-------------------------- 1 file changed, 75 insertions(+), 94 deletions(-) diff --git a/src/NNStopping.jl b/src/NNStopping.jl index 6609bca..0aaa54b 100644 --- a/src/NNStopping.jl +++ b/src/NNStopping.jl @@ -1,112 +1,93 @@ -struct NNStopping{M, O} <: HighDimPDEAlgorithm +struct NNStopping{M, O} m::M opt::O +end - function NNStopping(m, opt = Flux.Optimise.Adam(0.1)) - new{typeof(m), typeof(opt)}(m, opt) - end +struct NNStoppingModelArray{M} + ms::M end -function DiffEqBase.solve(prob::PIDEProblem, - pdealg::NNStopping, - sdealg; - verbose = false, - maxiters = 300, - trajectories = 100, - dt = eltype(prob.tspan)(0), - ensemblealg = EnsembleThreads(), - kwargs...) +Flux.@functor NNStoppingModelArray +function (model::NNStoppingModelArray)(X, G) + XG = cat(X, reshape(G, 1, size(G)...), dims = 1) + broadcast((x, m) -> m(x), eachslice(XG, dims = 2)[2:end], model.ms) +end + +function DiffEqBase.solve(prob::SDEProblem, + pdealg::NNStopping, + sdealg; + verbose = false, + maxiters = 300, + trajectories = 100, + dt = eltype(prob.tspan)(0), + ensemblealg = EnsembleThreads(), + kwargs...) + g = prob.kwargs[:payoff] + + sde_prob = SDEProblem(prob.f, prob.u0, prob.tspan) + ensemble_prob = EnsembleProblem(sde_prob) + sim = solve(ensemble_prob, + sdealg, + trajectories = trajectories, + dt = dt, + adaptive = false) + + m = NNStoppingModelArray(pdealg.m) + opt = pdealg.opt tspan = prob.tspan - sigma = prob.g - μ = prob.f - g = prob.kwargs.data.g - u0 = prob.u0 ts = tspan[1]:dt:tspan[2] - T = tspan[2] - - m = alg.m - opt = alg.opt - - - prob = SDEProblem(μ,sigma,u0,tspan) - ensembleprob = EnsembleProblem(prob) - sim = solve(ensembleprob, sdealg, ensemblealg, dt=dt,trajectories=trajectories,adaptive=false) - payoff = [] - times = [] - iter = 0 - # for u in sim.u - # un = [] - # function Un(n , X ) - # if size(un)[1] >= n - # return un[n] - # else - # if(n == 1) - # ans = first(m(X[1])[1]) - # un = [ans] - # return ans - # else - # ans = max(first(m(X[n])[n]) , n + 1 - size(ts)[1])*(1 - sum(Un(i , X ) for i in 1:n-1)) - # un = vcat( un , ans) - # return ans - # end - # end - # end - - function Un!(un_cache, n, X) - !isnan(un_cache[n]) && return un_cache[n] - if n == 1 - un_cache[n] = first(m(X[1])[1]) - else - un_cache[n] = max(first(m(X[n])[n]) , n + 1 - length(ts))*(1 - sum(un_cache[1:n-1])) + N = length(ts) - 1 + + G = reduce(hcat, map(u -> map(i -> g(u.t[i], u.u[i]), 1:(N + 1)), sim)) + + function nn_loss(m) + preds = m(Array(sim), G) + preds = reduce(vcat, preds) + un_gs = map(eachcol(preds), eachcol(G)) do pred, g_ + local u_sum = pred[1] + uns = map(1:N) do i + i == 1 && return pred[i] + # i == N && return (1 - u_sum) + res = pred[i] * (1 - u_sum) + u_sum += res + return res + end + + return sum(uns .* g_[2:end]) end - return un_cache[n] + return -1 * mean(un_gs) end - - function loss() - reward = 0.00 - map(sim.u) do X - un_cache = fill(NaN, length(ts)) - reward += mapreduce((x, t) -> Un!(un_cache, i, X)*g(t,X), +, X, ts) + opt_state = Flux.setup(opt, m) + for epoch in 1:maxiters + # sim = solve(ensemble_prob, EM(), trajectories = M, dt = dt) + gs = Flux.gradient(m) do model + nn_loss(model) end - - return 10000 - reward + Flux.update!(opt_state, m, gs[1]) + l = nn_loss(m) + verbose && @info "Current Epoch : $epoch Current Loss :$l" end - dataset = Iterators.repeated(() , maxiters) - - callback = function () - l = loss() - println("Current loss is: $l") + final_preds = m(Array(sim), G) + final_preds = reduce(vcat, final_preds) + un_gs = map(eachcol(final_preds), eachcol(G)) do pred, g_ + local u_sum = pred[1] + uns = map(1:N) do i + i == 1 && return pred[i] + res = pred[i] * (1 - u_sum) + u_sum += res + return res + end + uns end - Flux.train!(loss, Flux.params(m), dataset, opt; cb = callback) - - Usum = 0 - ti = 0 - Xt = sim.u[1].u - un_cache = fill(NaN, N) - for i in 1:N - Un = Un!(un_cache, i , Xt) - Usum = Usum + Un - if Usum >= 1 - Un - ti = i - break - end + ns = map(un_gs) do un + argmax(cumsum(un) + un .>= 1) end - for u in sim.u - X = u.u - price = g(ts[ti] , X[ti]) - payoff = vcat(payoff , price) - times = vcat(times, ti) - iter = iter + 1 - # println("SUM : $sump") - # println("TIME : $ti") - end - sum(payoff)/size(payoff)[1] - - - -end \ No newline at end of file + tss = getindex.(Ref(ts), ns .+ 1) + payoff = mean(map((t, u) -> g(t, u(t)), tss, sim)) + return (payoff = payoff, stopping_time = tss) +end From 1fe1d3f166c251ca21b5d04a723ec074cf98ef52 Mon Sep 17 00:00:00 2001 From: ashutosh-b-b Date: Sat, 27 Jan 2024 00:34:11 +0530 Subject: [PATCH 04/37] chore: export `NNStopping` --- src/HighDimPDE.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/HighDimPDE.jl b/src/HighDimPDE.jl index 3980c14..912a543 100644 --- a/src/HighDimPDE.jl +++ b/src/HighDimPDE.jl @@ -149,8 +149,9 @@ include("DeepSplitting.jl") include("DeepBSDE.jl") include("DeepBSDE_Han.jl") include("MLP.jl") +include("NNStopping.jl") -export PIDEProblem, PIDEProblem, PIDESolution, DeepSplitting, DeepBSDE, MLP +export PIDEProblem, PIDEProblem, PIDESolution, DeepSplitting, DeepBSDE, MLP, NNStopping export NormalSampling, UniformSampling, NoSampling, solve end From 3faa879c8ee1f4bcfb8204ab7463788ca45a474c Mon Sep 17 00:00:00 2001 From: ashutosh-b-b Date: Sat, 27 Jan 2024 00:34:31 +0530 Subject: [PATCH 05/37] test: add tests for NNStopping - adds test for a Bermudan max call example from paper - adds test for a Put Basket in Duprie's local volatility model --- test/NNStopping.jl | 101 +++++++++++++++++++++++++++------------------ 1 file changed, 60 insertions(+), 41 deletions(-) diff --git a/test/NNStopping.jl b/test/NNStopping.jl index e46a97a..875b776 100644 --- a/test/NNStopping.jl +++ b/test/NNStopping.jl @@ -1,49 +1,68 @@ -using Test, Flux , StochasticDiffEq , LinearAlgebra +using Test, Flux, StochasticDiffEq, LinearAlgebra println("Optimal Stopping Time Test") -using NeuralPDE +using HighDimPDE using Random -Random.seed!(100) - -d = 1 -r = 0.04f0 -beta = 0.2f0 -T = 1 -u0 = fill(80.00 , d , 1) -f(du,u,p,t) = (du .= r*u) -sigma(du,u,p,t) = (du .= Diagonal(beta*u)) -tspan = (0.0 , 1.0) -N = 50 -dt = tspan[2]/49 +Random.seed!(10) +# Bermudan Max-Call Standard Example +d = 3 +r = 0.05 +beta = 0.2 +T = 3 +u0 = fill(90.0, d) +delta = 0.1 +f(du, u, p, t) = du .= (r - delta) * u +sigma(du, u, p, t) = du .= beta * u +tspan = (0.0, T) +N = 9 +dt = T / (N) K = 100.00 -function g(t , x) - return exp(-r*t)*(max(K - maximum(x) , 0)) +function g(t, x) + return exp(-r * t) * (max(maximum(x) - K, 0)) end -prob = SDEProblem(f , sigma , u0 , tspan ; g = g) -opt = Flux.ADAM(0.1) -m = Chain(Dense(d , 5, tanh), Dense(5, 16 , tanh) , Dense(16 , N ), softmax) -sol = solve(prob, NeuralPDE.NNStopping(m, opt), verbose = true, dt = dt, sde_alg = EM() - abstol=1e-6, maxiters = 20 , trajectories = 200) +prob = SDEProblem(f, sigma, u0, tspan; payoff = g) +models = [Chain(Dense(d + 1, 32, tanh), BatchNorm(32, tanh), Dense(32, 1, sigmoid)) + for i in 1:N] +opt = Flux.Optimisers.Adam(0.01) +alg = NNStopping(models, opt) +sol = solve(prob, alg, SRIW1(); dt = dt, trajectories = 1000, maxiters = 1000, verbose = true) +analytical_sol = 11.278 # Ref [1] +@test abs(analytical_sol - sol.payoff) < 0.5 + +# Put basket option in Dupire’s local volatility model +d = 5 +r = 0.05 +beta = 0.2 +T = 1.0 +u0 = fill(100.0, d) +delta = 0.1 +f(u, p, t) = (r - delta) * u + +function B(t, x) + x * 0.6 * exp(-0.05 * sqrt(t)) * + (1.2 - exp(-0.1 * t - 0.001 * (exp(-r * t * x - first(u0))^2))) +end +sigma(u, p, t) = B.(Ref(t), u) -##Analytical Binomial Tree approach for American Options -function BinomialTreeAM1D(S0 , N , r , beta) - V = zeros(N+1) - dT = T/N - u = exp(beta*sqrt(dT)) - d = 1/u - S_T = [S0*(u^j)* (d^(N-j)) for j in 0:N] - a = exp(r*dT) - p = (a - d)/(u - d) - q = 1.0 - p - V = [max(K - x , 0) for x in S_T] - for i in N-1:-1:0 - V[1:end-1] = exp(-r*dT).*(p*V[2:end] + q*V[1:end-1]) - S_T = S_T*u - V = [max(K - S_T[i] , V[i]) for i in 1:size(S_T)[1]] - end - return V[1] +tspan = (0.0f0, T) +N = 10 +dt = T / (N) +K = 100.00 +function g(t, x) + return exp(-r * t) * (max(K - mean(x), 0)) end -real_sol = BinomialTreeAM1D(u0[1] , N , r , beta) -error = abs(sol - real_sol) -@test error < 0.5 \ No newline at end of file +prob = SDEProblem(f, sigma, u0, tspan; payoff = g) +models = [Chain(Dense(d + 1, 32, tanh), BatchNorm(32, tanh), Dense(32, 1, sigmoid)) + for i in 1:N] +opt = Flux.Optimisers.Adam(0.01) +alg = NNStopping(models, opt) +sol = solve(prob, alg, SRIW1(); dt = dt, trajectories = 1000, maxiters = 500, verbose = true) +analytical_sol = 6.301 # Ref [2] +@test abs(analytical_sol - sol.payoff) < 0.5 + +### References for Analytical Payoffs: +#= +[1] Section 4.4.1: Solving high-dimensional optimal stopping problems using deep learning (https://arxiv.org/pdf/1908.01602.pdf) +[2] Section 4.4.3: Solving high-dimensional optimal stopping problems using deep learning (https://arxiv.org/pdf/1908.01602.pdf) +=# From d346a80c113fd1f7b6b96a397567037d818cfdd0 Mon Sep 17 00:00:00 2001 From: ashutosh-b-b Date: Sat, 27 Jan 2024 01:35:54 +0530 Subject: [PATCH 06/37] docs: add docstring for `NNStopping` --- src/NNStopping.jl | 40 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 40 insertions(+) diff --git a/src/NNStopping.jl b/src/NNStopping.jl index 0aaa54b..df53dad 100644 --- a/src/NNStopping.jl +++ b/src/NNStopping.jl @@ -1,3 +1,43 @@ +""" +```julia +NNStopping(models, opt) +``` + +[Deep Optimal Stopping](https://arxiv.org/pdf/1908.01602.pdf), S Becker, P Cheridito, A Jentzen3, and T Welti. + +## Arguments +- `models::Vector{Flux.Chain}`: A vector of Flux.Chain where each model corresponds to a specific timestep from the timespan (tspan). The overall length of the vector should be `length(timesteps) - 1`. +- `opt`: the optimization algorithm to be used to optimize the neural networks. Defaults to `ADAM(0.1)`. + +## Example +d = 3 # Number of assets in the stock +r = 0.05 # interest rate +beta = 0.2 # volatility +T = 3 # maturity +u0 = fill(90.0, d) # initial stock value +delta = 0.1 # delta +f(du, u, p, t) = du .= (r - delta) * u # drift +sigma(du, u, p, t) = du .= beta * u # diffusion +tspan = (0.0, T) +N = 9 # discretization parameter +dt = T / (N) +K = 100.00 # strike price + +# payoff function +function g(t, x) + return exp(-r * t) * (max(maximum(x) - K, 0)) +end + +prob = SDEProblem(f, sigma, u0, tspan; payoff = g) +models = [Chain(Dense(d + 1, 32, tanh), BatchNorm(32, tanh), Dense(32, 1, sigmoid)) + for i in 1:N] + +opt = Flux.Optimisers.Adam(0.01) +alg = NNStopping(models, opt) + +sol = solve(prob, alg, SRIW1(); dt = dt, trajectories = 1000, maxiters = 1000, verbose = true) +``` +""" struct NNStopping{M, O} m::M opt::O From f5c4fd5fcedd18713f07de6027443d1ba6796e7b Mon Sep 17 00:00:00 2001 From: ashutosh-b-b Date: Sat, 27 Jan 2024 01:36:18 +0530 Subject: [PATCH 07/37] test: enable tests --- test/runtests.jl | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index 4f28da3..6aafef9 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -2,10 +2,11 @@ using Test: include using HighDimPDE, Test @testset "HighDimPDE" begin - include("reflect.jl") - include("MLP.jl") - include("DeepSplitting.jl") - include("DeepBSDE.jl") - include("DeepBSDE_Han.jl") - include("MCSample.jl") + # include("reflect.jl") + # include("MLP.jl") + # include("DeepSplitting.jl") + # include("DeepBSDE.jl") + # include("DeepBSDE_Han.jl") + # include("MCSample.jl") + include("NNStopping.jl") end From ab340cf9003e5d69812c2272f049721ff564b417 Mon Sep 17 00:00:00 2001 From: ashutosh-b-b Date: Sat, 27 Jan 2024 01:36:40 +0530 Subject: [PATCH 08/37] chore: fix formatting --- test/NNStopping.jl | 19 ++++++++++++++++--- 1 file changed, 16 insertions(+), 3 deletions(-) diff --git a/test/NNStopping.jl b/test/NNStopping.jl index 875b776..b923ebe 100644 --- a/test/NNStopping.jl +++ b/test/NNStopping.jl @@ -1,9 +1,10 @@ using Test, Flux, StochasticDiffEq, LinearAlgebra println("Optimal Stopping Time Test") using HighDimPDE +using Statistics using Random -Random.seed!(10) +Random.seed!(101) # Bermudan Max-Call Standard Example d = 3 r = 0.05 @@ -26,7 +27,13 @@ models = [Chain(Dense(d + 1, 32, tanh), BatchNorm(32, tanh), Dense(32, 1, sigmoi for i in 1:N] opt = Flux.Optimisers.Adam(0.01) alg = NNStopping(models, opt) -sol = solve(prob, alg, SRIW1(); dt = dt, trajectories = 1000, maxiters = 1000, verbose = true) +sol = solve(prob, + alg, + SRIW1(); + dt = dt, + trajectories = 1000, + maxiters = 1000, + verbose = true) analytical_sol = 11.278 # Ref [1] @test abs(analytical_sol - sol.payoff) < 0.5 @@ -57,7 +64,13 @@ models = [Chain(Dense(d + 1, 32, tanh), BatchNorm(32, tanh), Dense(32, 1, sigmoi for i in 1:N] opt = Flux.Optimisers.Adam(0.01) alg = NNStopping(models, opt) -sol = solve(prob, alg, SRIW1(); dt = dt, trajectories = 1000, maxiters = 500, verbose = true) +sol = solve(prob, + alg, + SRIW1(); + dt = dt, + trajectories = 1000, + maxiters = 500, + verbose = true) analytical_sol = 6.301 # Ref [2] @test abs(analytical_sol - sol.payoff) < 0.5 From e670a5a66d684a025b01c048ecd7402d8d8c4db9 Mon Sep 17 00:00:00 2001 From: ashutosh-b-b Date: Tue, 30 Jan 2024 13:17:33 +0530 Subject: [PATCH 09/37] docs: add docs for NNStopping --- docs/src/NNStopping.md | 29 +++++++++++++++++++++++++++++ 1 file changed, 29 insertions(+) create mode 100644 docs/src/NNStopping.md diff --git a/docs/src/NNStopping.md b/docs/src/NNStopping.md new file mode 100644 index 0000000..a5fb052 --- /dev/null +++ b/docs/src/NNStopping.md @@ -0,0 +1,29 @@ +# [The `NNStopping` algorithm](@id nn_stopping) + +```@autodocs +Modules = [HighDimPDE] +Pages = ["NNStopping.jl"] +``` +## The general idea 💡 + +Similar to DeepSplitting and DeepBSDE, NNStopping evaluates the PDE as a Stochastic Differential Equation. Consider an Obstacle PDE of the form: +```math + max\lbrace\partial_t u(t,x) + \mu(t, x) \nabla_x u(t,x) + \frac{1}{2} \sigma^2(t, x) \Delta_x u(t,x) , g(t,x) - u(t,x)\rbrace +``` + +Such PDEs are commonly used as representations for the dynamics of stock prices that can be exercised before maturity, such as American Options. + +Using the Feynman-Kac formula, the underlying SDE will be: + +```math +dX_{t}=\mu(X,t)dt + \sigma(X,t)\ dW_{t}^{Q} +``` + +The payoff of the option would then be: + +```math +sup\lbrace\mathbb{E}[g(\tau, X_\tau)]\rbrace +``` +Where τ is the stopping (exercising) time. The goal is to retrive both the optimal exercising strategy (τ) and the payoff. + +We approximate each stopping decision with a neural network architecture, inorder to maximise the expected payoff. From 71ff37eaeef8511038eecf1cb6f9a4f9125ffb45 Mon Sep 17 00:00:00 2001 From: ashutosh-b-b Date: Tue, 30 Jan 2024 13:22:13 +0530 Subject: [PATCH 10/37] docs: add tutorial for NNStopping --- docs/src/tutorials/nnstopping.md | 50 ++++++++++++++++++++++++++++++++ 1 file changed, 50 insertions(+) create mode 100644 docs/src/tutorials/nnstopping.md diff --git a/docs/src/tutorials/nnstopping.md b/docs/src/tutorials/nnstopping.md new file mode 100644 index 0000000..80c874b --- /dev/null +++ b/docs/src/tutorials/nnstopping.md @@ -0,0 +1,50 @@ +## Solving for optimal strategy and expected payoff of a Bermudan Max-Call option + +We will calculate optimal strategy for Bermudan Max-Call option with following drift, diffusion and payoff: +```math +μ(x) =(r − δ) x, σ(x) = β diag(x1, ... , xd),\\ +g(t, x) = e^{-rt}max\lbrace max\lbrace x1, ... , xd \rbrace − K, 0\rbrace +``` +We define the parameters, drift function and the diffusion function for the dynamics of the option. +```julia +d = 3 # Number of assets in the stock +r = 0.05 # interest rate +beta = 0.2 # volatility +T = 3 # maturity +u0 = fill(90.0, d) # initial stock value +delta = 0.1 # delta +f(du, u, p, t) = du .= (r - delta) * u # drift +sigma(du, u, p, t) = du .= beta * u # diffusion +tspan = (0.0, T) +N = 9 # discretization parameter +dt = T / (N) +K = 100.00 # strike price + +# payoff function +function g(t, x) + return exp(-r * t) * (max(maximum(x) - K, 0)) +end + +``` +We then define a `SDEProblem`: +```julia +prob = SDEProblem(f, sigma, u0, tspan; payoff = g) +``` +!!! note + We provide the payoff function with a keyword argument `payoff` + +And now we define our models: +```julia +models = [Chain(Dense(d + 1, 32, tanh), BatchNorm(32, tanh), Dense(32, 1, sigmoid)) + for i in 1:N] +``` +!!! note + The number of models should be equal to the time discritization. + +And finally we define our optimizer and algorithm, and call `solve`: +```julia +opt = Flux.Optimisers.Adam(0.01) +alg = NNStopping(models, opt) + +sol = solve(prob, alg, SRIW1(); dt = dt, trajectories = 1000, maxiters = 1000, verbose = true) +``` From de266d26818e634b3a82ca0e3b4719dd2b3e7ded Mon Sep 17 00:00:00 2001 From: ashutosh-b-b Date: Tue, 30 Jan 2024 13:23:40 +0530 Subject: [PATCH 11/37] docs: add tutorial and page to pages.jl --- docs/pages.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/docs/pages.jl b/docs/pages.jl index c3b4f2a..c953803 100644 --- a/docs/pages.jl +++ b/docs/pages.jl @@ -3,11 +3,13 @@ pages = [ "Getting started" => "getting_started.md", "Solver Algorithms" => ["MLP.md", "DeepSplitting.md", - "DeepBSDE.md"], + "DeepBSDE.md", + "NNStopping.md"], "Tutorials" => [ "tutorials/deepsplitting.md", "tutorials/deepbsde.md", "tutorials/mlp.md", + "tutorials/nnstopping.md" ], "Feynman Kac formula" => "Feynman_Kac.md", ] From a9351234ef2e8fe31211df8cd2a7196112e64756 Mon Sep 17 00:00:00 2001 From: ashutosh-b-b Date: Wed, 31 Jan 2024 01:26:37 +0530 Subject: [PATCH 12/37] feat: add dispatch on `PIDEProblem` for non linear term 0 This dipatch will cater to defining optimal stopping problems and kolmogorov PDE problems --- src/HighDimPDE.jl | 55 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 55 insertions(+) diff --git a/src/HighDimPDE.jl b/src/HighDimPDE.jl index 912a543..f4e64b5 100644 --- a/src/HighDimPDE.jl +++ b/src/HighDimPDE.jl @@ -133,6 +133,61 @@ function PIDESolution(x0, ts, losses, usols, ufuns, limits = nothing) limits) end +""" +Defines a Partial Differential Equation without the non linear term `f` : +## Arguments : +* `μ` : drift function, of the form `μ(x, p, t)`. +* `σ` : diffusion function `σ(x, p, t)`. +* `x0`: initial condition for the `x` +* `tspan`: timespan of the problem. +* `g` : terminal condition, of the form `g(x)`. Defaults to `nothing` + +## Keyword Arguments: +* `payoff`: While defining an optimal stopping problem provide the `payoff` function `pay(x,t)`` +* `xspan`: The domain of system state. This can be a tuple of floats for single dimension, and a vector of tuples for multiple dimensions. Where each tuple corresponds to a dimension of state vector. +* `noise_rate_prototype` : Incase of a non diagonal noise, the prototype of `dx` in `σ` +""" +function PIDEProblem(μ, + σ, + x0, + tspan, + g = nothing; + xspan = nothing, + payoff = nothing, + p = nothing, + x0_sample = NoSampling(), + kwargs...) + kwargs = merge(NamedTuple(kwargs), (xspan = xspan, payoff = payoff)) + + @assert !isnothing(x0)||!isnothing(xspan) "Both `x0` and `xspan` cannot be nothing" + xspan = !isnothing(xspan) && !isa(xspan, AbstractVector) ? [xspan] : xspan + + # For Optimal stopping problem + u_est = isnothing(x0) && !isnothing(xspan) ? g(first.(xspan)) : payoff(x0, tspan[1]) + + PIDEProblem{typeof(u_est), + typeof(g), + Nothing, + typeof(μ), + typeof(σ), + typeof(x0), + eltype(tspan), + typeof(p), + typeof(x0_sample), + Nothing, + typeof(kwargs)}(u_est, + g, + nothing, + μ, + σ, + x0, + tspan, + p, + x0_sample, + nothing, + kwargs) +end + Base.summary(prob::PIDESolution) = string(nameof(typeof(prob))) function Base.show(io::IO, A::PIDESolution) From 71fe64c35aae5e0a14762ff9b2ae2f5ba46d80d2 Mon Sep 17 00:00:00 2001 From: ashutosh-b-b Date: Wed, 31 Jan 2024 01:27:34 +0530 Subject: [PATCH 13/37] fix: update NNStopping to work with PIDEProblem --- src/NNStopping.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/NNStopping.jl b/src/NNStopping.jl index df53dad..f8ff8c1 100644 --- a/src/NNStopping.jl +++ b/src/NNStopping.jl @@ -24,11 +24,11 @@ dt = T / (N) K = 100.00 # strike price # payoff function -function g(t, x) +function g(x, t) return exp(-r * t) * (max(maximum(x) - K, 0)) end -prob = SDEProblem(f, sigma, u0, tspan; payoff = g) +prob = PIDEProblem(f, sigma, u0, tspan; payoff = g) models = [Chain(Dense(d + 1, 32, tanh), BatchNorm(32, tanh), Dense(32, 1, sigmoid)) for i in 1:N] @@ -54,7 +54,7 @@ function (model::NNStoppingModelArray)(X, G) broadcast((x, m) -> m(x), eachslice(XG, dims = 2)[2:end], model.ms) end -function DiffEqBase.solve(prob::SDEProblem, +function DiffEqBase.solve(prob::PIDEProblem, pdealg::NNStopping, sdealg; verbose = false, @@ -65,7 +65,7 @@ function DiffEqBase.solve(prob::SDEProblem, kwargs...) g = prob.kwargs[:payoff] - sde_prob = SDEProblem(prob.f, prob.u0, prob.tspan) + sde_prob = SDEProblem(prob.μ, prob.σ, prob.x, prob.tspan) ensemble_prob = EnsembleProblem(sde_prob) sim = solve(ensemble_prob, sdealg, @@ -79,7 +79,7 @@ function DiffEqBase.solve(prob::SDEProblem, ts = tspan[1]:dt:tspan[2] N = length(ts) - 1 - G = reduce(hcat, map(u -> map(i -> g(u.t[i], u.u[i]), 1:(N + 1)), sim)) + G = reduce(hcat, map(u -> map(i -> g(u.u[i], u.t[i]), 1:(N + 1)), sim)) function nn_loss(m) preds = m(Array(sim), G) @@ -128,6 +128,6 @@ function DiffEqBase.solve(prob::SDEProblem, end tss = getindex.(Ref(ts), ns .+ 1) - payoff = mean(map((t, u) -> g(t, u(t)), tss, sim)) + payoff = mean(map((t, u) -> g(u(t), t), tss, sim)) return (payoff = payoff, stopping_time = tss) end From d9d9f4c0dae1400f6ff38ff5c4ea3f8353675c6b Mon Sep 17 00:00:00 2001 From: ashutosh-b-b Date: Wed, 31 Jan 2024 01:27:55 +0530 Subject: [PATCH 14/37] test: update tests --- test/NNStopping.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/NNStopping.jl b/test/NNStopping.jl index b923ebe..5657eea 100644 --- a/test/NNStopping.jl +++ b/test/NNStopping.jl @@ -22,7 +22,7 @@ function g(t, x) return exp(-r * t) * (max(maximum(x) - K, 0)) end -prob = SDEProblem(f, sigma, u0, tspan; payoff = g) +prob = PIDEProblem(f, sigma, u0, tspan; payoff = g) models = [Chain(Dense(d + 1, 32, tanh), BatchNorm(32, tanh), Dense(32, 1, sigmoid)) for i in 1:N] opt = Flux.Optimisers.Adam(0.01) From 425ab1516a54ec396b0b3181ca2ec28bb3a76510 Mon Sep 17 00:00:00 2001 From: ashutosh-b-b Date: Wed, 31 Jan 2024 01:28:33 +0530 Subject: [PATCH 15/37] docs: update docs --- docs/src/NNStopping.md | 2 +- docs/src/tutorials/nnstopping.md | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/docs/src/NNStopping.md b/docs/src/NNStopping.md index a5fb052..d9db1a0 100644 --- a/docs/src/NNStopping.md +++ b/docs/src/NNStopping.md @@ -22,7 +22,7 @@ dX_{t}=\mu(X,t)dt + \sigma(X,t)\ dW_{t}^{Q} The payoff of the option would then be: ```math -sup\lbrace\mathbb{E}[g(\tau, X_\tau)]\rbrace +sup\lbrace\mathbb{E}[g(X_\tau, \tau)]\rbrace ``` Where τ is the stopping (exercising) time. The goal is to retrive both the optimal exercising strategy (τ) and the payoff. diff --git a/docs/src/tutorials/nnstopping.md b/docs/src/tutorials/nnstopping.md index 80c874b..54fb78a 100644 --- a/docs/src/tutorials/nnstopping.md +++ b/docs/src/tutorials/nnstopping.md @@ -21,14 +21,14 @@ dt = T / (N) K = 100.00 # strike price # payoff function -function g(t, x) +function g(x, t) return exp(-r * t) * (max(maximum(x) - K, 0)) end ``` -We then define a `SDEProblem`: +We then define a `PIDEProblem` with no non linear term: ```julia -prob = SDEProblem(f, sigma, u0, tspan; payoff = g) +prob = PIDEProblem(f, sigma, u0, tspan; payoff = g) ``` !!! note We provide the payoff function with a keyword argument `payoff` From 0ada89c375b45d42de9b1e3a40954353473cf967 Mon Sep 17 00:00:00 2001 From: ashutosh-b-b Date: Wed, 31 Jan 2024 01:28:53 +0530 Subject: [PATCH 16/37] chore: fix formatting in pages.jl --- docs/pages.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/pages.jl b/docs/pages.jl index c953803..ce1a4bf 100644 --- a/docs/pages.jl +++ b/docs/pages.jl @@ -3,13 +3,13 @@ pages = [ "Getting started" => "getting_started.md", "Solver Algorithms" => ["MLP.md", "DeepSplitting.md", - "DeepBSDE.md", + "DeepBSDE.md", "NNStopping.md"], "Tutorials" => [ "tutorials/deepsplitting.md", "tutorials/deepbsde.md", "tutorials/mlp.md", - "tutorials/nnstopping.md" + "tutorials/nnstopping.md", ], "Feynman Kac formula" => "Feynman_Kac.md", ] From ab4be5a1031f77accddbb86c7bbfffa10434b6be Mon Sep 17 00:00:00 2001 From: ashutosh-b-b Date: Wed, 31 Jan 2024 01:31:04 +0530 Subject: [PATCH 17/37] docs: fix spelling errors --- docs/src/NNStopping.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/NNStopping.md b/docs/src/NNStopping.md index d9db1a0..5fb4e37 100644 --- a/docs/src/NNStopping.md +++ b/docs/src/NNStopping.md @@ -24,6 +24,6 @@ The payoff of the option would then be: ```math sup\lbrace\mathbb{E}[g(X_\tau, \tau)]\rbrace ``` -Where τ is the stopping (exercising) time. The goal is to retrive both the optimal exercising strategy (τ) and the payoff. +Where τ is the stopping (exercising) time. The goal is to retrieve both the optimal exercising strategy (τ) and the payoff. We approximate each stopping decision with a neural network architecture, inorder to maximise the expected payoff. From 31a64f178681579334b05aee00d3e16d6d6e35c6 Mon Sep 17 00:00:00 2001 From: ashutosh-b-b Date: Wed, 31 Jan 2024 01:58:44 +0530 Subject: [PATCH 18/37] test: update tests --- test/NNStopping.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/NNStopping.jl b/test/NNStopping.jl index 5657eea..ca2a667 100644 --- a/test/NNStopping.jl +++ b/test/NNStopping.jl @@ -18,7 +18,7 @@ tspan = (0.0, T) N = 9 dt = T / (N) K = 100.00 -function g(t, x) +function g(x, t) return exp(-r * t) * (max(maximum(x) - K, 0)) end @@ -56,10 +56,10 @@ tspan = (0.0f0, T) N = 10 dt = T / (N) K = 100.00 -function g(t, x) +function g(x, t) return exp(-r * t) * (max(K - mean(x), 0)) end -prob = SDEProblem(f, sigma, u0, tspan; payoff = g) +prob = PIDEProblem(f, sigma, u0, tspan; payoff = g) models = [Chain(Dense(d + 1, 32, tanh), BatchNorm(32, tanh), Dense(32, 1, sigmoid)) for i in 1:N] opt = Flux.Optimisers.Adam(0.01) From a698b96ae8f98ede6d7a85dc9aa3dddc7e38d2dd Mon Sep 17 00:00:00 2001 From: ashutosh-b-b Date: Thu, 1 Feb 2024 01:17:09 +0530 Subject: [PATCH 19/37] test: Increase maxiters in DeepSplitting heat eqn tests --- test/DeepSplitting.jl | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/test/DeepSplitting.jl b/test/DeepSplitting.jl index e401d22..aad8597 100644 --- a/test/DeepSplitting.jl +++ b/test/DeepSplitting.jl @@ -99,7 +99,7 @@ end sol = solve(prob, alg, dt, verbose = false, use_cuda = use_cuda, - maxiters = 1500, + maxiters = 1800, batch_size = batch_size, cuda_device = cuda_device) xs = x0_sample(repeat(x0, 1, batch_size)) @@ -235,7 +235,7 @@ end Dense(hls, 1)) # Neural network used by the scheme opt = Flux.Optimise.Adam(1e-3) #optimiser - alg = DeepSplitting(nn, opt = opt ) + alg = DeepSplitting(nn, opt = opt) X0 = fill(0.0f0, d) # initial point g(X) = 1.0f0 ./ (2.0f0 .+ 4.0f-1 * sum(X .^ 2, dims = 1)) # initial condition @@ -282,7 +282,7 @@ end Dense(hls, 1)) # Neural network used by the scheme opt = Flux.Optimise.Adam(1e-2) #optimiser - alg = DeepSplitting(nn, opt = opt ) + alg = DeepSplitting(nn, opt = opt) X0 = fill(0.0f0, d) # initial point g(X) = exp.(-0.25f0 * sum(X .^ 2, dims = 1)) # initial condition @@ -330,7 +330,7 @@ if false Dense(hls, 1)) # Neural network used by the scheme opt = Flux.Optimise.Adam(1e-3) #optimiser - alg = DeepSplitting(nn, opt = opt ) + alg = DeepSplitting(nn, opt = opt) X0 = repeat([1.0f0, 0.5f0], div(d, 2)) # initial point g(X) = sum(X .^ 2, dims = 1) # initial condition @@ -382,7 +382,7 @@ if false Dense(hls, 1)) # Neural network used by the scheme opt = Flux.Optimise.Adam(1e-3) #optimiser - alg = DeepSplitting(nn, opt = opt ) + alg = DeepSplitting(nn, opt = opt) X0 = fill(0.0f0, d) # initial point g(X) = log.(5.0f-1 .+ 5.0f-1 * sum(X .^ 2, dims = 1)) # initial condition @@ -431,7 +431,7 @@ end Dense(hls, 1)) # Neural network used by the scheme opt = Flux.Optimise.Adam() - alg = DeepSplitting(nn, opt = opt, λs = [1e-2,1e-3] ) + alg = DeepSplitting(nn, opt = opt, λs = [1e-2, 1e-3]) X0 = fill(100.0f0, d) # initial point g(X) = minimum(X, dims = 1) # initial condition @@ -520,7 +520,7 @@ end Dense(hls, 1, x -> x^2)) # positive function opt = Flux.Optimise.Adam(1e-2)#optimiser - alg = DeepSplitting(nn_batch, K=K, opt = opt, mc_sample = x0_sample) + alg = DeepSplitting(nn_batch, K = K, opt = opt, mc_sample = x0_sample) function g(x) Float32((2 * π)^(-d / 2)) * ss0^(-Float32(d) * 5.0f-1) * @@ -576,7 +576,7 @@ end Dense(hls, 1)) # Neural network used by the scheme opt = Flux.Optimise.Adam(1e-2) #optimiser - alg = DeepSplitting(nn, K=K, opt = opt, mc_sample = UniformSampling(-∂, ∂) ) + alg = DeepSplitting(nn, K = K, opt = opt, mc_sample = UniformSampling(-∂, ∂)) x0 = fill(0.0f0, d) # initial point g(X) = exp.(-0.25f0 * sum(X .^ 2, dims = 1)) # initial condition From 901470e18b57ed68978cd5fa1806aa84e4b29c72 Mon Sep 17 00:00:00 2001 From: ashutosh-b-b Date: Thu, 1 Feb 2024 12:29:12 +0530 Subject: [PATCH 20/37] feat: update PIDEProblem definition to accomodate all kinds of equations --- src/HighDimPDE.jl | 188 +++++++++++++++++++++------------------------- 1 file changed, 87 insertions(+), 101 deletions(-) diff --git a/src/HighDimPDE.jl b/src/HighDimPDE.jl index f4e64b5..9fa8e4f 100644 --- a/src/HighDimPDE.jl +++ b/src/HighDimPDE.jl @@ -28,30 +28,8 @@ function Base.show(io::IO, A::AbstractPDEProblem) show(io, A.tspan) end -""" -$(SIGNATURES) - -Defines a Partial Integro Differential Problem, of the form -```math -\\begin{aligned} - \\frac{du}{dt} &= \\tfrac{1}{2} \\text{Tr}(\\sigma \\sigma^T) \\Delta u(x, t) + \\mu \\nabla u(x, t) \\\\ - &\\quad + \\int f(x, y, u(x, t), u(y, t), ( \\nabla_x u )(x, t), ( \\nabla_x u )(y, t), p, t) dy, -\\end{aligned} -``` -with `` u(x,0) = g(x)``. - -## Arguments +include("MCSample.jl") -* `g` : initial condition, of the form `g(x, p, t)`. -* `f` : nonlinear function, of the form `f(x, y, u(x, t), u(y, t), ∇u(x, t), ∇u(y, t), p, t)`. -* `μ` : drift function, of the form `μ(x, p, t)`. -* `σ` : diffusion function `σ(x, p, t)`. -* `x`: point where `u(x,t)` is approximated. Is required even in the case where `x0_sample` is provided. Determines the dimensionality of the PDE. -* `tspan`: timespan of the problem. -* `p`: the parameter vector. -* `x0_sample` : sampling method for `x0`. Can be `UniformSampling(a,b)`, `NormalSampling(σ_sampling, shifted)`, or `NoSampling` (by default). If `NoSampling`, only solution at the single point `x` is evaluated. -* `neumann_bc`: if provided, Neumann boundary conditions on the hypercube `neumann_bc[1] × neumann_bc[2]`. -""" struct PIDEProblem{uType, G, F, Mu, Sigma, xType, tType, P, UD, NBC, K} <: DiffEqBase.AbstractODEProblem{uType, tType, false} u0::uType @@ -67,43 +45,107 @@ struct PIDEProblem{uType, G, F, Mu, Sigma, xType, tType, P, UD, NBC, K} <: kwargs::K end -function PIDEProblem(g, f, μ, σ, x::Vector{X}, tspan; - p = nothing, - x0_sample = NoSampling(), - neumann_bc::NBC = nothing, - kwargs...) where {X <: AbstractFloat, NBC <: Union{Nothing, AbstractVector}} - @assert eltype(tspan)<:AbstractFloat "`tspan` should be a tuple of Float" - - isnothing(neumann_bc) ? nothing : @assert eltype(eltype(neumann_bc)) <: eltype(x) - @assert eltype(g(x))==eltype(x) "Type of `g(x)` must match type of x" - try - @assert(eltype(f(x, x, g(x), g(x), x, x, p, tspan[1]))==eltype(x), - "Type of non linear function `f(x)` must type of x") - catch e - if isa(e, MethodError) - @assert(eltype(f(x, eltype(x)(0.0), x, p, tspan[1]))==eltype(x), +""" +$(SIGNATURES) + +Defines one of the following equations: +- Partial Integro Differential Equation + * f -> f(x, y, u(x, t), u(y, t), ∇u(x, t), ∇u(y, t), p, t) +- Semilinear Parabolic Partial Differential Equation + * f -> f(X, u, σᵀ∇u, p, t) +- Kolmogorov Differential Equation + * f -> `nothing` + * x0 -> nothing, xspan must be provided. +- Obstacle Partial Differential Equation + * f -> `nothing` + * g -> `nothing` + * discounted payoff function provided. + +## Arguments + +* `μ` : drift function, of the form `μ(x, p, t)`. +* `σ` : diffusion function `σ(x, p, t)`. +* `x`: point where `u(x,t)` is approximated. Is required even in the case where `x0_sample` is provided. Determines the dimensionality of the PDE. +* `tspan`: timespan of the problem. +* `g` : initial condition, of the form `g(x, p, t)`. +* `f` : when defining PIDE : nonlinear function, of the form `f(x, y, u(x, t), u(y, t), ∇u(x, t), ∇u(y, t), p, t)`. When defining Semilinear PDE: `f(X, u, σᵀ∇u, p, t)` + +## Optional Arguments +* `p`: the parameter vector. +* `x0_sample` : sampling method for `x0`. Can be `UniformSampling(a,b)`, `NormalSampling(σ_sampling, shifted)`, or `NoSampling` (by default). If `NoSampling`, only solution at the single point `x` is evaluated. +* `neumann_bc`: if provided, Neumann boundary conditions on the hypercube `neumann_bc[1] × neumann_bc[2]`. +* `xspan`: The domain of the independent variable `x` +* `payoff`: The discounted payoff function. Required when solving for optimal stopping problem (Obstacle PDEs). +""" +function PIDEProblem(μ, + σ, + x0::Union{Nothing, AbstractArray}, + tspan::TF, + g = nothing, + f = nothing; + p::Union{Nothing, AbstractVector} = nothing, + xspan::Union{Nothing, TF, AbstractVector{<:TF}} = nothing, + x0_sample::Union{Nothing, AbstractSampling} = NoSampling(), + neumann_bc::Union{Nothing, AbstractVector} = nothing, + payoff = nothing, + kw...) where {TF <: Tuple{AbstractFloat, AbstractFloat}} + + # Check the Initial Condition Function returns correct types. + isnothing(g) && @assert !isnothing(payoff) "Either of `g` or `payoff` must be provided." + + isnothing(neumann_bc) ? nothing : @assert eltype(eltype(neumann_bc)) <: eltype(x0) + + @assert !isnothing(x0)||!isnothing(xspan) "Either of `x0` or `xspan` must be provided." + + # Check if the Non Linear Function `f` returns correct values. + if !isnothing(f) + try + @assert(eltype(f(x0, x0, g(x0), g(x0), x0, x0, p, tspan[1]))==eltype(x0), "Type of non linear function `f(x)` must type of x") - else - throw(e) + catch e + if isa(e, MethodError) + @assert(eltype(f(x0, eltype(x0)(0.0), x0, p, tspan[1]))==eltype(x0), + "Type of non linear function `f(x)` must type of x") + else + throw(e) + end end end - PIDEProblem{typeof(g(x)), + # Wrap kwargs : + kw = NamedTuple(kw) + prob_kw = (p = p, xspan = xspan, payoff = payoff) + kwargs = merge(prob_kw, kw) + + # If xspan isa Tuple, then convert it as a Vector{Tuple} with single element + xspan = isa(xspan, Tuple) ? [xspan] : xspan + + # if `x0` is not provided, pick up the lower-bound of `xspan`. + x0 = isnothing(x0) ? first.(xspan) : x0 + + # Initial Condition + u0 = if haskey(kw, :p_prototype) + u0 = g(x0, kw.p_prototype.p_phi) + else + !isnothing(g) ? g(x0) : payoff(x0, 0.0) + end + @assert eltype(u0)==eltype(x0) "Type of `g(x)` must match the Type of x" + PIDEProblem{typeof(u0), typeof(g), typeof(f), typeof(μ), typeof(σ), - typeof(x), + typeof(x0), eltype(tspan), typeof(p), typeof(x0_sample), typeof(neumann_bc), - typeof(kwargs)}(g(x), + typeof(kwargs)}(u0, g, f, μ, σ, - x, + x0, tspan, p, x0_sample, @@ -133,61 +175,6 @@ function PIDESolution(x0, ts, losses, usols, ufuns, limits = nothing) limits) end -""" -Defines a Partial Differential Equation without the non linear term `f` : -## Arguments : -* `μ` : drift function, of the form `μ(x, p, t)`. -* `σ` : diffusion function `σ(x, p, t)`. -* `x0`: initial condition for the `x` -* `tspan`: timespan of the problem. -* `g` : terminal condition, of the form `g(x)`. Defaults to `nothing` - -## Keyword Arguments: -* `payoff`: While defining an optimal stopping problem provide the `payoff` function `pay(x,t)`` -* `xspan`: The domain of system state. This can be a tuple of floats for single dimension, and a vector of tuples for multiple dimensions. Where each tuple corresponds to a dimension of state vector. -* `noise_rate_prototype` : Incase of a non diagonal noise, the prototype of `dx` in `σ` -""" -function PIDEProblem(μ, - σ, - x0, - tspan, - g = nothing; - xspan = nothing, - payoff = nothing, - p = nothing, - x0_sample = NoSampling(), - kwargs...) - kwargs = merge(NamedTuple(kwargs), (xspan = xspan, payoff = payoff)) - - @assert !isnothing(x0)||!isnothing(xspan) "Both `x0` and `xspan` cannot be nothing" - xspan = !isnothing(xspan) && !isa(xspan, AbstractVector) ? [xspan] : xspan - - # For Optimal stopping problem - u_est = isnothing(x0) && !isnothing(xspan) ? g(first.(xspan)) : payoff(x0, tspan[1]) - - PIDEProblem{typeof(u_est), - typeof(g), - Nothing, - typeof(μ), - typeof(σ), - typeof(x0), - eltype(tspan), - typeof(p), - typeof(x0_sample), - Nothing, - typeof(kwargs)}(u_est, - g, - nothing, - μ, - σ, - x0, - tspan, - p, - x0_sample, - nothing, - kwargs) -end - Base.summary(prob::PIDESolution) = string(nameof(typeof(prob))) function Base.show(io::IO, A::PIDESolution) @@ -198,7 +185,6 @@ function Base.show(io::IO, A::PIDESolution) show(io, A.us) end -include("MCSample.jl") include("reflect.jl") include("DeepSplitting.jl") include("DeepBSDE.jl") @@ -206,7 +192,7 @@ include("DeepBSDE_Han.jl") include("MLP.jl") include("NNStopping.jl") -export PIDEProblem, PIDEProblem, PIDESolution, DeepSplitting, DeepBSDE, MLP, NNStopping +export PIDEProblem, PIDESolution, DeepSplitting, DeepBSDE, MLP, NNStopping export NormalSampling, UniformSampling, NoSampling, solve end From 130cd20175c063405f7ffb62362181e19fe9aac4 Mon Sep 17 00:00:00 2001 From: ashutosh-b-b Date: Thu, 1 Feb 2024 12:30:43 +0530 Subject: [PATCH 21/37] test: add tests for PIDEProblem --- test/PIDEProblem.jl | 205 ++++++++++++++++++++++++++++++++++++++++++++ test/runtests.jl | 1 + 2 files changed, 206 insertions(+) create mode 100644 test/PIDEProblem.jl diff --git a/test/PIDEProblem.jl b/test/PIDEProblem.jl new file mode 100644 index 0000000..0baeb2e --- /dev/null +++ b/test/PIDEProblem.jl @@ -0,0 +1,205 @@ +using HighDimPDE + +@testset "PIDEs" begin + for d in [1, 10] + x0 = fill(8.0f0, d) # initial points + tspan = (0.0f0, 5.0f0) + + μ_f(X, p, t) = zero(X) #Vector d x 1 + σ_f(X, p, t) = Diagonal(ones(Float32, d)) |> Matrix + g(X) = sum(X .^ 2) # terminal condition + f_nonlinear(y, z, v_y, v_z, ∇v_y, ∇v_z, p, t) = 0.0f0 .* v_y + + prob = PIDEProblem(μ_f, σ_f, x0, tspan, g, f_nonlinear) + + @test prob.f == f_nonlinear + @test prob.g == g + @test prob.μ == μ_f + @test prob.σ == σ_f + @test prob.x == x0 + @test prob.u0 == g(x0) + @test prob.p == nothing + @test prob.tspan == tspan + + @testset "With neumann_bc" begin + neumann_bc = [-1 * fill(5.0f-1, d), fill(5.0f-1, d)] + prob = PIDEProblem(μ_f, σ_f, x0, tspan, g, f_nonlinear; neumann_bc) + + @test prob.f == f_nonlinear + @test prob.g == g + @test prob.μ == μ_f + @test prob.σ == σ_f + @test prob.x == x0 + @test prob.u0 == g(x0) + @test prob.p == nothing + @test prob.tspan == tspan + @test prob.neumann_bc == neumann_bc + end + + @testset "With x0_sample" begin + x0_sample = UniformSampling(fill(-5.0f-1, d), fill(5.0f-1, d)) + + prob = PIDEProblem(μ_f, σ_f, x0, tspan, g, f_nonlinear; x0_sample) + + @test prob.f == f_nonlinear + @test prob.g == g + @test prob.μ == μ_f + @test prob.σ == σ_f + @test prob.x == x0 + @test prob.u0 == g(x0) + @test prob.p == nothing + @test prob.tspan == tspan + @test prob.x0_sample == x0_sample + end + + @testset "With x0_sample and neumann_bc" begin + x0_sample = UniformSampling(fill(-5.0f-1, d), fill(5.0f-1, d)) + neumann_bc = [-1 * fill(5.0f-1, d), fill(5.0f-1, d)] + + prob = PIDEProblem(μ_f, σ_f, x0, tspan, g, f_nonlinear; x0_sample, neumann_bc) + + @test prob.f == f_nonlinear + @test prob.g == g + @test prob.μ == μ_f + @test prob.σ == σ_f + @test prob.x == x0 + @test prob.u0 == g(x0) + @test prob.p == nothing + @test prob.tspan == tspan + @test prob.x0_sample == x0_sample + @test prob.neumann_bc == neumann_bc + end + end +end + +@testset "Semilinear Parabolic PDEs" begin + for d in [1, 10] + x0 = fill(8.0f0, d) # initial points + tspan = (0.0f0, 5.0f0) + + μ_f(X, p, t) = zero(X) #Vector d x 1 + σ_f(X, p, t) = Diagonal(ones(Float32, d)) |> Matrix + g(X) = sum(X .^ 2) # terminal condition + f_semilinear(X, u, σᵀ∇u, p, t) = Float32(0.0) + + prob = PIDEProblem(μ_f, σ_f, x0, tspan, g, f_semilinear) + + @test prob.f == f_semilinear + @test prob.g == g + @test prob.μ == μ_f + @test prob.σ == σ_f + @test prob.x == x0 + @test prob.u0 == g(x0) + @test prob.p == nothing + @test prob.tspan == tspan + end +end + +@testset "Obstacle PDEs : Optimal Stopping Problems" begin + for d in [1, 10] + r = 0.05 + beta = 0.2 + T = 3.0 + u0 = fill(90.0, d) + delta = 0.1 + mu(du, u, p, t) = du .= (r - delta) * u + sigma(du, u, p, t) = du .= beta * u + tspan = (0.0, T) + N = 9 + dt = T / (N) + K = 100.00 + function payoff(x, t) + return exp(-r * t) * (max(maximum(x) - K, 0)) + end + + prob = PIDEProblem(mu, sigma, u0, tspan; payoff = payoff) + + @test prob.f == nothing + @test prob.g == nothing + @test prob.μ == mu + @test prob.σ == sigma + @test prob.x == u0 + @test prob.u0 == payoff(x0, 0.0) + @test prob.p == nothing + @test prob.tspan == tspan + @test prob.kwargs.payoff == payoff + end +end + +@testset "Kolmogorov PDEs" begin + for d in [1, 10] + xspan = d == 1 ? (-6.0, 6.0) : [(-6.0, 6.0) for _ in 1:d] + tspan = (0.0, 1.0) + σ(u, p, t) = 0.5 * u + μ(u, p, t) = 0.5 * 0.25 * u + + function g(x) + 1.77 .* x .- 0.015 .* x .^ 3 + end + + prob = PIDEProblem(μ, σ, nothing, tspan, g; xspan) + + @test prob.f == nothing + @test prob.g == g + @test prob.μ == μ + @test prob.σ == σ + + if d == 1 + @test prob.u0 == g([xspan[1]]) + @test prob.x == [xspan[1]] + else + @test prob.u0 == g(first.(xspan)) + @test prob.x == first.(xspan) + end + @test prob.p == nothing + @test prob.tspan == tspan + @test prob.kwargs.xspan == xspan + end +end + +@testset "Kolmogorov PDEs - Parametric Family" begin + for d in [1, 10] + γ_mu_prototype = nothing + γ_sigma_prototype = zeros(d, d, 1) + γ_phi_prototype = nothing + + tspan = (0.00, 1.00) + xspan = d == 1 ? (0.00, 3.00) : [(0.00, 3.00) for _ in 1:d] + + function g(x, p_phi) + x .^ 2 + end + + sigma(dx, x, p_sigma, t) = dx .= p_sigma[:, :, 1] + mu(dx, x, p_mu, t) = dx .= 0.00 + + p_domain = (p_sigma = (0.00, 2.00), p_mu = nothing, p_phi = nothing) + p_prototype = (p_sigma = γ_sigma_prototype, + p_mu = γ_mu_prototype, + p_phi = γ_phi_prototype) + + prob = PIDEProblem(mu, sigma, nothing, tspan, g; + xspan, + p_domain = p_domain, + p_prototype = p_prototype) + + @test prob.f == nothing + @test prob.g == g + @test prob.μ == mu + @test prob.σ == sigma + + if d == 1 + @test prob.u0 == g([xspan[1]], p_prototype.p_phi) + @test prob.x == [xspan[1]] + else + @test prob.u0 == g(first.(xspan), p_prototype.p_phi) + @test prob.x == first.(xspan) + end + + @test prob.p == nothing + @test prob.tspan == tspan + @test prob.kwargs.xspan == xspan + @test prob.kwargs.p_domain == p_domain + @test prob.kwargs.p_prototype == p_prototype + end +end diff --git a/test/runtests.jl b/test/runtests.jl index 0c54ab4..d4829d4 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,6 +1,7 @@ using SafeTestsets, Test @testset "HighDimPDE" begin @time @safetestset "Quality Assurance" include("qa.jl") + @time @safetestset "PIDEProblem" include("PIDEProblem.jl") @time @safetestset "reflect" include("reflect.jl") @time @safetestset "reflect" include("reflect.jl") @time @safetestset "MLP" include("MLP.jl") From cd1d23b766a4c19765565edcaf95cf09feda14ba Mon Sep 17 00:00:00 2001 From: ashutosh-b-b Date: Thu, 1 Feb 2024 13:32:28 +0530 Subject: [PATCH 22/37] chore: update all tests with the new definition of PIDEProblem --- test/DeepBSDE.jl | 12 ++++++------ test/DeepBSDE_Han.jl | 12 ++++++------ test/DeepSplitting.jl | 22 +++++++++++----------- test/MLP.jl | 18 +++++++++--------- test/NNStopping.jl | 2 +- test/PIDEProblem.jl | 4 ++-- 6 files changed, 35 insertions(+), 35 deletions(-) diff --git a/test/DeepBSDE.jl b/test/DeepBSDE.jl index dfffefb..2f1a8d9 100644 --- a/test/DeepBSDE.jl +++ b/test/DeepBSDE.jl @@ -27,7 +27,7 @@ end f(X, u, σᵀ∇u, p, t) = Float32(0.0) μ_f(X, p, t) = zero(X) #Vector d x 1 σ_f(X, p, t) = Diagonal(ones(Float32, d)) |> Matrix #Matrix d x d - prob = PIDEProblem(g, f, μ_f, σ_f, x0, tspan) + prob = PIDEProblem(μ_f, σ_f, x0, tspan, g, f) hls = 10 + d #hidden layer size opt = Flux.Optimise.Adam(0.005) #optimizer @@ -69,7 +69,7 @@ end f(X, u, σᵀ∇u, p, t) = Float32(0.0) μ_f(X, p, t) = zero(X) #Vector d x 1 σ_f(X, p, t) = Diagonal(ones(Float32, d)) |> Matrix #Matrix d x d - prob = PIDEProblem(g, f, μ_f, σ_f, x0, tspan) + prob = PIDEProblem(μ_f, σ_f, x0, tspan, g, f) hls = 10 + d #hidden layer size opt = Flux.Optimise.Adam(0.005) #optimizer @@ -113,7 +113,7 @@ end g(X) = sum(X .^ 2) μ_f(X, p, t) = zero(X) #Vector d x 1 σ_f(X, p, t) = Diagonal(sigma * X) |> Matrix #Matrix d x d - prob = PIDEProblem(g, f, μ_f, σ_f, x0, tspan) + prob = PIDEProblem(μ_f, σ_f, x0, tspan, g, f) hls = 10 + d #hide layer size opt = Flux.Optimise.Adam(0.001) @@ -153,7 +153,7 @@ end f(X, u, σᵀ∇u, p, t) = u .- u .^ 3 μ_f(X, p, t) = zero(X) #Vector d x 1 σ_f(X, p, t) = Diagonal(ones(Float32, d)) |> Matrix #Matrix d x d - prob = PIDEProblem(g, f, μ_f, σ_f, x0, tspan) + prob = PIDEProblem(μ_f, σ_f, x0, tspan, g, f) hls = 20 + d #hidden layer size opt = Flux.Optimise.Adam(5^-3) #optimizer @@ -195,7 +195,7 @@ end f(X, u, σᵀ∇u, p, t) = -λ * sum(σᵀ∇u .^ 2) μ_f(X, p, t) = zero(X) #Vector d x 1 λ σ_f(X, p, t) = Diagonal(sqrt(2.0f0) * ones(Float32, d)) #Matrix d x d - prob = PIDEProblem(g, f, μ_f, σ_f, x0, tspan) + prob = PIDEProblem(μ_f, σ_f, x0, tspan, g, f) # TODO: This is a very large neural networks which size must be reduced. hls = 256 #hidden layer size @@ -266,7 +266,7 @@ end μ_f(X, p, t) = µc * X #Vector d x 1 σ_f(X, p, t) = σc * Diagonal(X) #Matrix d x d - prob = PIDEProblem(g, f, μ_f, σ_f, x0, tspan) + prob = PIDEProblem(μ_f, σ_f, x0, tspan, g, f) hls = 256 #hidden layer size opt = Flux.Optimise.Adam(0.008) #optimizer diff --git a/test/DeepBSDE_Han.jl b/test/DeepBSDE_Han.jl index a3d859b..ddafd5e 100644 --- a/test/DeepBSDE_Han.jl +++ b/test/DeepBSDE_Han.jl @@ -28,7 +28,7 @@ end f(X, u, σᵀ∇u, p, t) = 0.0f0 # function from solved equation μ_f(X, p, t) = 0.0 σ_f(X, p, t) = 1.0 - prob = PIDEProblem(g, f, μ_f, σ_f, x0, tspan) + prob = PIDEProblem(μ_f, σ_f, x0, tspan, g, f) hls = 10 + d #hidden layer size opt = Flux.Optimise.Adam(0.005) #optimizer @@ -72,7 +72,7 @@ end f(X, u, σᵀ∇u, p, t) = 0.0f0 μ_f(X, p, t) = 0.0 σ_f(X, p, t) = 1.0 - prob = PIDEProblem(g, f, μ_f, σ_f, x0, tspan) + prob = PIDEProblem(μ_f, σ_f, x0, tspan, g, f) hls = 10 + d #hidden layer size opt = Flux.Optimise.Adam(0.005) #optimizer @@ -117,7 +117,7 @@ end g(X) = sum(X .^ 2) μ_f(X, p, t) = 0.0 σ_f(X, p, t) = Diagonal(sigma * X) - prob = PIDEProblem(g, f, μ_f, σ_f, x0, tspan) + prob = PIDEProblem(μ_f, σ_f, x0, tspan, g, f) hls = 10 + d #hide layer size opt = Flux.Optimise.Adam(0.001) @@ -159,7 +159,7 @@ end f(X, u, σᵀ∇u, p, t) = u .- u .^ 3 μ_f(X, p, t) = 0.0 σ_f(X, p, t) = 1.0 - prob = PIDEProblem(g, f, μ_f, σ_f, x0, tspan) + prob = PIDEProblem(μ_f, σ_f, x0, tspan, g, f) hls = 10 + d #hidden layer size opt = Flux.Optimise.Adam(5^-4) #optimizer @@ -205,7 +205,7 @@ end f(X, u, σᵀ∇u, p, t) = -λ * sum(σᵀ∇u .^ 2) μ_f(X, p, t) = 0.0 σ_f(X, p, t) = sqrt(2) - prob = PIDEProblem(g, f, μ_f, σ_f, x0, tspan) + prob = PIDEProblem(μ_f, σ_f, x0, tspan, g, f) hls = 12 + d #hidden layer size opt = Flux.Optimise.Adam(0.03) #optimizer @@ -278,7 +278,7 @@ end σc = 0.2f0 μ_f(X, p, t) = µc * X #Vector d x 1 σ_f(X, p, t) = σc * Diagonal(X) #Matrix d x d - prob = PIDEProblem(g, f, μ_f, σ_f, x0, tspan) + prob = PIDEProblem(μ_f, σ_f, x0, tspan, g, f) hls = 10 + d #hidden layer size opt = Flux.Optimise.Adam(0.008) #optimizer diff --git a/test/DeepSplitting.jl b/test/DeepSplitting.jl index aad8597..41f67d4 100644 --- a/test/DeepSplitting.jl +++ b/test/DeepSplitting.jl @@ -50,7 +50,7 @@ end f(y, z, v_y, v_z, ∇v_y, ∇v_z, p, t) = 0.0f0 .* v_y # defining the problem - prob = PIDEProblem(g, f, μ, σ, x0, tspan) + prob = PIDEProblem(μ, σ, x0, tspan, g, f) # solving sol = solve(prob, alg, dt, verbose = true, @@ -94,7 +94,7 @@ end f(y, z, v_y, v_z, ∇v_y, ∇v_z, p, t) = 0.0f0 .* v_y #TODO: this fix is not nice # defining the problem - prob = PIDEProblem(g, f, μ, σ, x0, tspan, x0_sample = x0_sample) + prob = PIDEProblem(μ, σ, x0, tspan, g, f; x0_sample = x0_sample) # solving sol = solve(prob, alg, dt, verbose = false, @@ -140,7 +140,7 @@ end f(y, z, v_y, v_z, ∇v_y, ∇v_z, p, t) = 0.0f0 .* v_y #TODO: this fix is not nice # defining the problem - prob = PIDEProblem(g, f, μ, σ, x0, tspan, + prob = PIDEProblem(μ, σ, x0, tspan, g, f; x0_sample = x0_sample, neumann_bc = [-∂, ∂]) # solving @@ -198,7 +198,7 @@ end f(y, z, v_y, v_z, ∇v_y, ∇v_z, p, t) = r * v_y #TODO: this fix is not nice # defining the problem - prob = PIDEProblem(g, f, μ, σ, x0, tspan, + prob = PIDEProblem(μ, σ, x0, tspan, g, f; x0_sample = x0_sample) # solving sol = solve(prob, alg, dt, @@ -243,7 +243,7 @@ end f(y, z, v_y, v_z, ∇v_y, ∇v_z, p, t) = -a.(v_y) # nonlocal nonlinear function # defining the problem - prob = PIDEProblem(g, f, μ, σ, X0, tspan) + prob = PIDEProblem(μ, σ, X0, tspan, g, f) # solving @time sol = solve(prob, alg, @@ -290,7 +290,7 @@ end f(y, z, v_y, v_z, ∇v_y, ∇v_z, p, t) = a.(v_y) # nonlocal nonlinear function # defining the problem - prob = PIDEProblem(g, f, μ, σ, X0, tspan, neumann_bc = [-∂, ∂]) + prob = PIDEProblem(μ, σ, X0, tspan, g, f; neumann_bc = [-∂, ∂]) # solving @time sol = solve(prob, alg, @@ -337,7 +337,7 @@ if false f(y, z, v_y, v_z, ∇v_y, ∇v_z, p, t) = r * (v_y .- sum(y .* ∇v_y, dims = 1)) # defining the problem - prob = PIDEProblem(g, f, μ, σ, X0, tspan) + prob = PIDEProblem(μ, σ, X0, tspan, g, f) # solving @time xs, ts, sol = solve(prob, alg, @@ -389,7 +389,7 @@ if false f(y, z, v_y, v_z, ∇v_y, ∇v_z, p, t) = λ * sum(∇v_y .^ 2, dims = 1) # defining the problem - prob = PIDEProblem(g, f, μ, σ, X0, tspan) + prob = PIDEProblem(μ, σ, X0, tspang, f) # solving @time sol = solve(prob, alg, @@ -453,7 +453,7 @@ end f(y, z, v_y, v_z, ∇v_y, ∇v_z, p, t) = -(1.0f0 - δ) * Q.(v_y) .* v_y .- R * v_y # defining the problem - prob = PIDEProblem(g, f, μ, σ, X0, tspan) + prob = PIDEProblem(μ, σ, X0, tspan, g, f) # solving @time sol = solve(prob, alg, @@ -533,7 +533,7 @@ end end # nonlocal nonlinear part of the # defining the problem - prob = PIDEProblem(g, f, μ, σ, x0, tspan, + prob = PIDEProblem(μ, σ, x0, tspan, g, f; x0_sample = x0_sample) # solving sol = solve(prob, @@ -584,7 +584,7 @@ end f(y, z, v_y, v_z, ∇v_y, ∇v_z, p, t) = a.(v_y) .- a.(v_z) #.* Float32(π^(d/2)) * σ_sampling^d .* exp.(sum(z.^2, dims = 1) / σ_sampling^2) # nonlocal nonlinear part of the # defining the problem - prob = PIDEProblem(g, f, μ, σ, x0, tspan, neumann_bc = [-∂, ∂]) + prob = PIDEProblem(μ, σ, x0, tspan, g, f; neumann_bc = [-∂, ∂]) # solving @time sol = solve(prob, alg, diff --git a/test/MLP.jl b/test/MLP.jl index 9dbf4fa..6f8f8d0 100644 --- a/test/MLP.jl +++ b/test/MLP.jl @@ -33,7 +33,7 @@ end f(y, z, v_y, v_z, ∇v_y, ∇v_z, p, t) = 0e0 .* v_y #TODO: this fix is not nice # defining the problem - prob = PIDEProblem(g, f, μ, σ, x0, tspan) + prob = PIDEProblem(μ, σ, x0, tspan, g, f) # solving sol = solve(prob, alg, multithreading = false) u1 = sol.us[end] @@ -64,7 +64,7 @@ end f(y, z, v_y, v_z, ∇v_y, ∇v_z, p, t) = 0e0 .* v_y #TODO: this fix is not nice # defining the problem - prob = PIDEProblem(g, f, μ, σ, x0, tspan) + prob = PIDEProblem(μ, σ, x0, tspan, g, f) # solving sol = solve(prob, alg, multithreading = true) u1 = sol.us[end] @@ -90,7 +90,7 @@ end f(y, z, v_y, v_z, ∇v_y, ∇v_z, p, t) = 0e0 .* v_y #TODO: this fix is not nice # defining the problem - prob = PIDEProblem(g, f, μ, σ, x0, tspan, neumann_bc = neumann_bc) + prob = PIDEProblem(μ, σ, x0, tspan, g, f, neumann_bc = neumann_bc) # solving sol = solve(prob, alg) push!(u1s, sol.us[end]) @@ -120,7 +120,7 @@ end f(y, z, v_y, v_z, ∇v_y, ∇v_z, p, t) = r * v_y #TODO: this fix is not nice # defining the problem - prob = PIDEProblem(g, f, μ, σ, x, tspan) + prob = PIDEProblem(μ, σ, x, tspan, g, f) # solving sol = solve(prob, alg) u1 = sol.us[end] @@ -146,7 +146,7 @@ end f(y, z, v_y, v_z, ∇v_y, ∇v_z, p, t) = -a.(v_y) # defining the problem - prob = PIDEProblem(g, f, μ, σ, X0, tspan) + prob = PIDEProblem(μ, σ, X0, tspan, g, f) # solving sol = solve(prob, alg) u1 = sol.us[end] @@ -174,7 +174,7 @@ end f(y, z, v_y, v_z, ∇v_y, ∇v_z, p, t) = a.(v_y) # nonlocal nonlinear part of the # defining the problem - prob = PIDEProblem(g, f, μ, σ, X0, tspan, neumann_bc = neumann_bc) + prob = PIDEProblem(μ, σ, X0, tspan, g, f; neumann_bc = neumann_bc) # solving sol = solve(prob, alg) u1 = sol.us[end] @@ -213,7 +213,7 @@ end f(y, z, v_y, v_z, ∇v_y, ∇v_z, p, t) = -(1.0f0 - δ) * Q.(v_y) .* v_y .- R * v_y # defining the problem - prob = PIDEProblem(g, f, μ, σ, X0, tspan) + prob = PIDEProblem(μ, σ, X0, tspan, g, f) # solving sol = solve(prob, alg) @@ -268,7 +268,7 @@ end end # nonlocal nonlinear part of the # defining the problem - prob = PIDEProblem(g, f, μ, σ, x, tspan) + prob = PIDEProblem(μ, σ, x, tspan, g, f) # solving sol = solve(prob, alg) u1 = sol.us[end] @@ -298,7 +298,7 @@ end f(y, z, v_y, v_z, ∇v_y, ∇v_z, p, t) = a.(v_y) .- a.(v_z) #.* Float32(π^(d/2)) * σ_sampling^d .* exp.(sum(z.^2, dims = 1) / σ_sampling^2) # nonlocal nonlinear part of the # defining the problem - prob = PIDEProblem(g, f, μ, σ, x, tspan, neumann_bc = neumann_bc) + prob = PIDEProblem(μ, σ, x, tspan, g, f, neumann_bc = neumann_bc) # solving sol = solve(prob, alg) push!(u1s, sol.us[end]) diff --git a/test/NNStopping.jl b/test/NNStopping.jl index ca2a667..a48ffa1 100644 --- a/test/NNStopping.jl +++ b/test/NNStopping.jl @@ -9,7 +9,7 @@ Random.seed!(101) d = 3 r = 0.05 beta = 0.2 -T = 3 +T = 3.0 u0 = fill(90.0, d) delta = 0.1 f(du, u, p, t) = du .= (r - delta) * u diff --git a/test/PIDEProblem.jl b/test/PIDEProblem.jl index 0baeb2e..1f35722 100644 --- a/test/PIDEProblem.jl +++ b/test/PIDEProblem.jl @@ -1,4 +1,4 @@ -using HighDimPDE +using HighDimPDE @testset "PIDEs" begin for d in [1, 10] @@ -119,7 +119,7 @@ end @test prob.μ == mu @test prob.σ == sigma @test prob.x == u0 - @test prob.u0 == payoff(x0, 0.0) + @test prob.u0 == payoff(u0, 0.0) @test prob.p == nothing @test prob.tspan == tspan @test prob.kwargs.payoff == payoff From f5ed8e41deb7f459d8fd6fa31988e4fbddc8d916 Mon Sep 17 00:00:00 2001 From: ashutosh-b-b Date: Thu, 1 Feb 2024 13:33:03 +0530 Subject: [PATCH 23/37] docs: update all docs with the new definition of PIDEProblem --- docs/src/DeepSplitting.md | 4 ++-- docs/src/getting_started.md | 6 +++--- docs/src/tutorials/deepbsde.md | 6 +++--- docs/src/tutorials/deepsplitting.md | 4 +--- docs/src/tutorials/mlp.md | 5 ++--- paper/paper.md | 6 +++--- src/DeepBSDE.jl | 2 +- 7 files changed, 15 insertions(+), 18 deletions(-) diff --git a/docs/src/DeepSplitting.md b/docs/src/DeepSplitting.md index c041454..bd14ffd 100644 --- a/docs/src/DeepSplitting.md +++ b/docs/src/DeepSplitting.md @@ -62,14 +62,14 @@ In `HighDimPDE.jl` the right parameter combination $\theta$ is found by iterativ `DeepSplitting` allows obtaining $u(t,x)$ on a single point $x \in \Omega$ with the keyword $x$. ```julia -prob = PIDEProblem(g, f, μ, σ, x, tspan) +prob = PIDEProblem(μ, σ, x, tspan, g, f,) ``` ### Hypercube Yet more generally, one wants to solve Eq. (1) on a $d$-dimensional cube $[a,b]^d$. This is offered by `HighDimPDE.jl` with the keyword `x0_sample`. ```julia -prob = PIDEProblem(g, f, μ, σ, x, tspan, x0_sample = x0_sample) +prob = PIDEProblem(μ, σ, x, tspan, g, f; x0_sample = x0_sample) ``` Internally, this is handled by assigning a random variable as the initial point of the particles, i.e. ```math diff --git a/docs/src/getting_started.md b/docs/src/getting_started.md index d68a10d..8cbe510 100644 --- a/docs/src/getting_started.md +++ b/docs/src/getting_started.md @@ -33,7 +33,7 @@ g(x) = exp(-sum(x.^2)) # initial condition μ(x, p, t) = 0.0 # advection coefficients σ(x, p, t) = 0.1 # diffusion coefficients f(x, y, v_x, v_y, ∇v_x, ∇v_y, p, t) = max(0.0, v_x) * (1 - max(0.0, v_x)) # nonlocal nonlinear part of the -prob = PIDEProblem(g, f, μ, σ, x0, tspan) # defining the problem +prob = PIDEProblem(μ, σ, x0, tspan, g, f) # defining the problem ## Definition of the algorithm alg = MLP() # defining the algorithm. We use the Multi Level Picard algorithm @@ -62,7 +62,7 @@ g(x) = exp( -sum(x.^2) ) # initial condition σ(x, p, t) = 0.1 # diffusion coefficients mc_sample = UniformSampling(fill(-5f-1, d), fill(5f-1, d)) f(x, y, v_x, v_y, ∇v_x, ∇v_y, p, t) = max(0.0, v_x) * (1 - max(0.0, v_y)) -prob = PIDEProblem(g, f, μ, σ, x0, tspan) # defining x0_sample is sufficient to implement Neumann boundary conditions +prob = PIDEProblem(μ, σ, x0, tspan, g, f) # defining x0_sample is sufficient to implement Neumann boundary conditions ## Definition of the algorithm alg = MLP(mc_sample = mc_sample) @@ -87,7 +87,7 @@ g(x) = exp.(-sum(x.^2, dims=1)) # initial condition σ(x, p, t) = 0.1f0 # diffusion coefficients x0_sample = UniformSampling(fill(-5f-1, d), fill(5f-1, d)) f(x, y, v_x, v_y, ∇v_x, ∇v_y, p, t) = v_x .* (1f0 .- v_y) -prob = PIDEProblem(g, f, μ, σ, x0, tspan, +prob = PIDEProblem(μ, σ, x0, tspan, g, f; x0_sample = x0_sample) ## Definition of the neural network to use diff --git a/docs/src/tutorials/deepbsde.md b/docs/src/tutorials/deepbsde.md index 6f00113..76bac9f 100644 --- a/docs/src/tutorials/deepbsde.md +++ b/docs/src/tutorials/deepbsde.md @@ -19,7 +19,7 @@ g(X) = log(0.5f0 + 0.5f0 * sum(X.^2)) f(X,u,σᵀ∇u,p,t) = -λ * sum(σᵀ∇u.^2) μ_f(X,p,t) = zero(X) # Vector d x 1 λ σ_f(X,p,t) = Diagonal(sqrt(2.0f0) * ones(Float32, d)) # Matrix d x d -prob = PIDEProblem(g, f, μ_f, σ_f, X0, tspan) +prob = PIDEProblem(μ_f, σ_f, X0, tspan, g, f) hls = 10 + d # hidden layer size opt = Optimisers.Adam(0.01) # optimizer # sub-neural network approximating solutions at the desired point @@ -75,7 +75,7 @@ g(X) = log(0.5f0 + 0.5f0*sum(X.^2)) f(X,u,σᵀ∇u,p,t) = -λ*sum(σᵀ∇u.^2) μ_f(X,p,t) = zero(X) #Vector d x 1 λ σ_f(X,p,t) = Diagonal(sqrt(2.0f0)*ones(Float32,d)) #Matrix d x d -prob = PIDEProblem(g, f, μ_f, σ_f, X0, tspan) +prob = PIDEProblem(μ_f, σ_f, X0, tspan, g, f) ``` #### Define the Solver Algorithm @@ -135,7 +135,7 @@ f(X,u,σᵀ∇u,p,t) = r * (u - sum(X.*σᵀ∇u)) g(X) = sum(X.^2) μ_f(X,p,t) = zero(X) #Vector d x 1 σ_f(X,p,t) = Diagonal(sigma*X) #Matrix d x d -prob = PIDEProblem(g, f, μ_f, σ_f, X0, tspan) +prob = PIDEProblem(μ_f, σ_f, X0, tspan, g, f) ``` As described in the API docs, we now need to define our `NNPDENS` algorithm diff --git a/docs/src/tutorials/deepsplitting.md b/docs/src/tutorials/deepsplitting.md index b04b208..0c94410 100644 --- a/docs/src/tutorials/deepsplitting.md +++ b/docs/src/tutorials/deepsplitting.md @@ -22,9 +22,7 @@ g(x) = exp.(- sum(x.^2, dims=1) ) # initial condition σ(x, p, t) = 0.1f0 # diffusion coefficients x0_sample = UniformSampling(fill(-5f-1, d), fill(5f-1, d)) f(x, y, v_x, v_y, ∇v_x, ∇v_y, p, t) = v_x .* (1f0 .- v_y) -prob = PIDEProblem(g, f, μ, - σ, x0, tspan, - x0_sample = x0_sample) +prob = prob = PIDEProblem(μ, σ, x0, tspan, g, f; x0_sample = x0_sample) ## Definition of the neural network to use using Flux # needed to define the neural network diff --git a/docs/src/tutorials/mlp.md b/docs/src/tutorials/mlp.md index a5042bf..6156571 100644 --- a/docs/src/tutorials/mlp.md +++ b/docs/src/tutorials/mlp.md @@ -17,7 +17,7 @@ g(x) = exp(- sum(x.^2) ) # initial condition μ(x, p, t) = 0.0 # advection coefficients σ(x, p, t) = 0.1 # diffusion coefficients f(x, y, v_x, v_y, ∇v_x, ∇v_y, p, t) = max(0.0, v_x) * (1 - max(0.0, v_x)) # nonlocal nonlinear part of the -prob = PIDEProblem(g, f, μ, σ, x0, tspan) # defining the problem +prob = PIDEProblem(μ, σ, x0, tspan, g, f) # defining the problem ## Definition of the algorithm alg = MLP() # defining the algorithm. We use the Multi Level Picard algorithm @@ -44,8 +44,7 @@ g(x) = exp( -sum(x.^2) ) # initial condition σ(x, p, t) = 0.1 # diffusion coefficients mc_sample = UniformSampling(fill(-5f-1, d), fill(5f-1, d)) f(x, y, v_x, v_y, ∇v_x, ∇v_y, t) = max(0.0, v_x) * (1 - max(0.0, v_y)) -prob = PIDEProblem(g, f, μ, - σ, x0, tspan) # defining x0_sample is sufficient to implement Neumann boundary conditions +prob = PIDEProblem(μ, σ, x0, tspan, g, f) # defining x0_sample is sufficient to implement Neumann boundary conditions ## Definition of the algorithm alg = MLP(mc_sample = mc_sample ) diff --git a/paper/paper.md b/paper/paper.md index ce631cc..bc66b02 100644 --- a/paper/paper.md +++ b/paper/paper.md @@ -127,7 +127,7 @@ vol = prod(x0_sample[2] - x0_sample[1]) f(y, z, v_y, v_z, p, t) = max.(v_y, 0f0) .* (m(y) .- vol * max.(v_z, 0f0) .* m(z)) # nonlocal nonlinear part of the # defining the problem -prob = PIDEProblem(g, f, μ, σ, tspan, +prob = PIDEProblem(μ, σ, tspan, g, f, x0_sample = x0_sample ) # solving @@ -162,8 +162,8 @@ g(x) = exp( -sum(x.^2) ) # initial condition σ(x, p, t) = 0.1 # diffusion coefficients x0_sample = [-1/2, 1/2] f(x, y, v_x, v_y, ∇v_x, ∇v_y, t) = max(0.0, v_x) * (1 - max(0.0, v_y)) -prob = PIDEProblem(g, f, μ, - σ, x0, tspan, +prob = PIDEProblem(μ, + σ, x0, tspan, g, f, x0_sample = x0_sample) # defining x0_sample is sufficient to implement Neumann boundary conditions ## Definition of the algorithm diff --git a/src/DeepBSDE.jl b/src/DeepBSDE.jl index 8b19909..8844e80 100644 --- a/src/DeepBSDE.jl +++ b/src/DeepBSDE.jl @@ -26,7 +26,7 @@ f(X,u,σᵀ∇u,p,t) = r * (u - sum(X.*σᵀ∇u)) g(X) = sum(X.^2) μ_f(X,p,t) = zero(X) #Vector d x 1 σ_f(X,p,t) = Diagonal(sigma*X) #Matrix d x d -prob = PIDEProblem(g, f, μ_f, σ_f, x0, tspan) +prob = PIDEProblem(μ_f, σ_f, x0, tspan, g, f) hls = 10 + d #hidden layer size opt = Flux.Optimise.Adam(0.001) From 9656d8ed1c932779e3878f9496cdcd7426d032d3 Mon Sep 17 00:00:00 2001 From: ashutosh-b-b Date: Thu, 1 Feb 2024 17:56:26 +0530 Subject: [PATCH 24/37] test: add Random.seed! to MLP tests --- test/MLP.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/test/MLP.jl b/test/MLP.jl index 6f8f8d0..3a09e0a 100644 --- a/test/MLP.jl +++ b/test/MLP.jl @@ -3,7 +3,8 @@ using Random using Test using Statistics # using Revise - +using Random +Random.seed!(100) #relative error l2 function rel_error_l2(u, uanal) if abs(uanal) >= 10 * eps(eltype(uanal)) From 32a66cc510723038291ad7f1203699b769c8fc4f Mon Sep 17 00:00:00 2001 From: ashutosh-b-b Date: Thu, 1 Feb 2024 20:06:57 +0530 Subject: [PATCH 25/37] test: add Random.seed! to DeepSplitting --- test/DeepSplitting.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/test/DeepSplitting.jl b/test/DeepSplitting.jl index 41f67d4..69fb32b 100644 --- a/test/DeepSplitting.jl +++ b/test/DeepSplitting.jl @@ -4,6 +4,9 @@ using Test using Flux using Statistics using CUDA +using Random +Random.seed!(100) + if CUDA.functional() use_cuda = true cuda_device = 0 From 518a2f1602ac6aa6541431604abb48f4ef1150cd Mon Sep 17 00:00:00 2001 From: Ashutosh Bharambe Date: Sun, 4 Feb 2024 00:53:50 +0530 Subject: [PATCH 26/37] feat: add `ParabolicPDEProblem` --- src/HighDimPDE.jl | 115 ++++++++++++++++++++++++++++++++++++---------- 1 file changed, 92 insertions(+), 23 deletions(-) diff --git a/src/HighDimPDE.jl b/src/HighDimPDE.jl index 9fa8e4f..82f7335 100644 --- a/src/HighDimPDE.jl +++ b/src/HighDimPDE.jl @@ -48,9 +48,89 @@ end """ $(SIGNATURES) -Defines one of the following equations: -- Partial Integro Differential Equation - * f -> f(x, y, u(x, t), u(y, t), ∇u(x, t), ∇u(y, t), p, t) +Defines a Partial Integro Differential Problem, of the form +```math +\\begin{aligned} + \\frac{du}{dt} &= \\tfrac{1}{2} \\text{Tr}(\\sigma \\sigma^T) \\Delta u(x, t) + \\mu \\nabla u(x, t) \\\\ + &\\quad + \\int f(x, y, u(x, t), u(y, t), ( \\nabla_x u )(x, t), ( \\nabla_x u )(y, t), p, t) dy, +\\end{aligned} +``` +with `` u(x,0) = g(x)``. + +## Arguments + +* `g` : initial condition, of the form `g(x, p, t)`. +* `f` : nonlinear function, of the form `f(x, y, u(x, t), u(y, t), ∇u(x, t), ∇u(y, t), p, t)`. +* `μ` : drift function, of the form `μ(x, p, t)`. +* `σ` : diffusion function `σ(x, p, t)`. +* `x`: point where `u(x,t)` is approximated. Is required even in the case where `x0_sample` is provided. Determines the dimensionality of the PDE. +* `tspan`: timespan of the problem. +* `p`: the parameter vector. +* `x0_sample` : sampling method for `x0`. Can be `UniformSampling(a,b)`, `NormalSampling(σ_sampling, shifted)`, or `NoSampling` (by default). If `NoSampling`, only solution at the single point `x` is evaluated. +* `neumann_bc`: if provided, Neumann boundary conditions on the hypercube `neumann_bc[1] × neumann_bc[2]`. +""" +function PIDEProblem(μ, + σ, + x0::Union{Nothing, AbstractArray}, + tspan::TF, + g, + f; + p::Union{Nothing, AbstractVector} = nothing, + x0_sample::Union{Nothing, AbstractSampling} = NoSampling(), + neumann_bc::Union{Nothing, AbstractVector} = nothing, + kw...) where {TF <: Tuple{AbstractFloat, AbstractFloat}} + + + isnothing(neumann_bc) ? nothing : @assert eltype(eltype(neumann_bc)) <: eltype(x0) + + @assert(eltype(f(x0, x0, g(x0), g(x0), x0, x0, p, tspan[1]))==eltype(x0), + "Type returned by non linear function `f` must match the type of `x0`") + + @assert eltype(g(x0))==eltype(x0) "Type of `g(x)` must match the Type of x" + + PIDEProblem{typeof(g(x0)), + typeof(g), + typeof(f), + typeof(μ), + typeof(σ), + typeof(x0), + eltype(tspan), + typeof(p), + typeof(x0_sample), + typeof(neumann_bc), + typeof(kw)}(g(x0), + g, + f, + μ, + σ, + x0, + tspan, + p, + x0_sample, + neumann_bc, + kw) + +end + +struct ParabolicPDEProblem{uType, G, F, Mu, Sigma, xType, tType, P, UD, NBC, K} <: + DiffEqBase.AbstractODEProblem{uType, tType, false} + u0::uType + g::G # initial condition + f::F # nonlinear part + μ::Mu + σ::Sigma + x::xType + tspan::Tuple{tType, tType} + p::P + x0_sample::UD # the domain of u to be solved + neumann_bc::NBC # neumann boundary conditions + kwargs::K +end + +""" +$(SIGNATURES) + +Defines a Parabolic Partial Differential Equation of the form: - Semilinear Parabolic Partial Differential Equation * f -> f(X, u, σᵀ∇u, p, t) - Kolmogorov Differential Equation @@ -68,7 +148,7 @@ Defines one of the following equations: * `x`: point where `u(x,t)` is approximated. Is required even in the case where `x0_sample` is provided. Determines the dimensionality of the PDE. * `tspan`: timespan of the problem. * `g` : initial condition, of the form `g(x, p, t)`. -* `f` : when defining PIDE : nonlinear function, of the form `f(x, y, u(x, t), u(y, t), ∇u(x, t), ∇u(y, t), p, t)`. When defining Semilinear PDE: `f(X, u, σᵀ∇u, p, t)` +* `f` : nonlinear function, of the form `f(X, u, σᵀ∇u, p, t)` ## Optional Arguments * `p`: the parameter vector. @@ -77,12 +157,12 @@ Defines one of the following equations: * `xspan`: The domain of the independent variable `x` * `payoff`: The discounted payoff function. Required when solving for optimal stopping problem (Obstacle PDEs). """ -function PIDEProblem(μ, +function ParabolicPDEProblem(μ, σ, x0::Union{Nothing, AbstractArray}, - tspan::TF, + tspan::TF; g = nothing, - f = nothing; + f = nothing, p::Union{Nothing, AbstractVector} = nothing, xspan::Union{Nothing, TF, AbstractVector{<:TF}} = nothing, x0_sample::Union{Nothing, AbstractSampling} = NoSampling(), @@ -97,24 +177,12 @@ function PIDEProblem(μ, @assert !isnothing(x0)||!isnothing(xspan) "Either of `x0` or `xspan` must be provided." - # Check if the Non Linear Function `f` returns correct values. - if !isnothing(f) - try - @assert(eltype(f(x0, x0, g(x0), g(x0), x0, x0, p, tspan[1]))==eltype(x0), - "Type of non linear function `f(x)` must type of x") - catch e - if isa(e, MethodError) - @assert(eltype(f(x0, eltype(x0)(0.0), x0, p, tspan[1]))==eltype(x0), + !isnothing(f) && @assert(eltype(f(x0, eltype(x0)(0.0), x0, p, tspan[1]))==eltype(x0), "Type of non linear function `f(x)` must type of x") - else - throw(e) - end - end - end # Wrap kwargs : kw = NamedTuple(kw) - prob_kw = (p = p, xspan = xspan, payoff = payoff) + prob_kw = (xspan = xspan, payoff = payoff) kwargs = merge(prob_kw, kw) # If xspan isa Tuple, then convert it as a Vector{Tuple} with single element @@ -130,7 +198,8 @@ function PIDEProblem(μ, !isnothing(g) ? g(x0) : payoff(x0, 0.0) end @assert eltype(u0)==eltype(x0) "Type of `g(x)` must match the Type of x" - PIDEProblem{typeof(u0), + + ParabolicPDEProblem{typeof(u0), typeof(g), typeof(f), typeof(μ), @@ -192,7 +261,7 @@ include("DeepBSDE_Han.jl") include("MLP.jl") include("NNStopping.jl") -export PIDEProblem, PIDESolution, DeepSplitting, DeepBSDE, MLP, NNStopping +export PIDEProblem, ParabolicPDEProblem, PIDESolution, DeepSplitting, DeepBSDE, MLP, NNStopping export NormalSampling, UniformSampling, NoSampling, solve end From dc66a82c37a5f3333c8f02f7640bfb82538e0d8f Mon Sep 17 00:00:00 2001 From: Ashutosh Bharambe Date: Sun, 4 Feb 2024 00:54:34 +0530 Subject: [PATCH 27/37] fix: update `solve` dispatches for ParabolicPDEProblem --- src/DeepBSDE.jl | 2 +- src/DeepBSDE_Han.jl | 2 +- src/DeepSplitting.jl | 10 ++++++++-- src/MLP.jl | 9 +++++++-- src/NNStopping.jl | 2 +- 5 files changed, 18 insertions(+), 7 deletions(-) diff --git a/src/DeepBSDE.jl b/src/DeepBSDE.jl index 8844e80..56d17f9 100644 --- a/src/DeepBSDE.jl +++ b/src/DeepBSDE.jl @@ -75,7 +75,7 @@ Returns a `PIDESolution` object. [Deep Primal-Dual algorithm for BSDEs](https://papers.ssrn.com/sol3/papers.cfm?abstract_id=3071506). - Extra keyword arguments passed to `solve` will be further passed to the SDE solver. """ -function DiffEqBase.solve(prob::PIDEProblem, +function DiffEqBase.solve(prob::ParabolicPDEProblem, pdealg::DeepBSDE, sdealg; verbose = false, diff --git a/src/DeepBSDE_Han.jl b/src/DeepBSDE_Han.jl index cc925c3..eb785b7 100644 --- a/src/DeepBSDE_Han.jl +++ b/src/DeepBSDE_Han.jl @@ -1,5 +1,5 @@ # called whenever sdealg is not specified. -function DiffEqBase.solve(prob::PIDEProblem, +function DiffEqBase.solve(prob::ParabolicPDEProblem, alg::DeepBSDE; dt, abstol = 1.0f-6, diff --git a/src/DeepSplitting.jl b/src/DeepSplitting.jl index 939b20e..9541c52 100644 --- a/src/DeepSplitting.jl +++ b/src/DeepSplitting.jl @@ -64,7 +64,7 @@ Returns a `PIDESolution` object. - `use_cuda` : set to `true` to use CUDA. - `cuda_device` : integer, to set the CUDA device used in the training, if `use_cuda == true`. """ -function DiffEqBase.solve(prob::PIDEProblem, +function DiffEqBase.solve(prob::Union{PIDEProblem, ParabolicPDEProblem}, alg::DeepSplitting, dt; batch_size = 1, @@ -98,7 +98,13 @@ function DiffEqBase.solve(prob::PIDEProblem, K = alg.K opt = alg.opt λs = alg.λs - g, f, μ, σ, p = prob.g, prob.f, prob.μ, prob.σ, prob.p + g, μ, σ, p = prob.g, prob.μ, prob.σ, prob.p + + f = if isa(prob, ParabolicPDEProblem) + (y, z, v_y, v_z, ∇v_y, ∇v_z, p, t) -> prob.f(y, v_y, ∇v_y, p, t ) + else + prob.f + end T = eltype(x0) # neural network model diff --git a/src/MLP.jl b/src/MLP.jl index 1e68a7b..a6e7e62 100644 --- a/src/MLP.jl +++ b/src/MLP.jl @@ -29,7 +29,7 @@ Returns a `PIDESolution` object. * `multithreading` : if `true`, distributes the job over all the threads available. * `verbose`: print information over the iterations. """ -function DiffEqBase.solve(prob::PIDEProblem, +function DiffEqBase.solve(prob::Union{PIDEProblem, ParabolicPDEProblem}, alg::MLP; multithreading = true, verbose = false,) @@ -41,7 +41,12 @@ function DiffEqBase.solve(prob::PIDEProblem, M = alg.M L = alg.L mc_sample! = alg.mc_sample! - g, f = prob.g, prob.f + g= prob.g + f = if isa(prob, ParabolicPDEProblem) + (y, z, v_y, v_z, ∇v_y, ∇v_z, p, t) -> prob.f(y, v_y, ∇v_y, p, t ) + else + prob.f + end # errors prob.x0_sample isa NoSampling ? nothing : diff --git a/src/NNStopping.jl b/src/NNStopping.jl index f8ff8c1..39d1793 100644 --- a/src/NNStopping.jl +++ b/src/NNStopping.jl @@ -54,7 +54,7 @@ function (model::NNStoppingModelArray)(X, G) broadcast((x, m) -> m(x), eachslice(XG, dims = 2)[2:end], model.ms) end -function DiffEqBase.solve(prob::PIDEProblem, +function DiffEqBase.solve(prob::ParabolicPDEProblem, pdealg::NNStopping, sdealg; verbose = false, From ee934ade7f79af1a7c4d5086b0077cfba2b66cd2 Mon Sep 17 00:00:00 2001 From: Ashutosh Bharambe Date: Sun, 4 Feb 2024 00:55:11 +0530 Subject: [PATCH 28/37] test: update tests --- test/DeepBSDE.jl | 12 +++---- test/DeepBSDE_Han.jl | 14 ++++---- test/DeepSplitting.jl | 36 +++++++++---------- test/MLP.jl | 28 +++++++-------- test/NNStopping.jl | 4 +-- .../{PIDEProblem.jl => ProblemConstructor.jl} | 8 ++--- 6 files changed, 51 insertions(+), 51 deletions(-) rename test/{PIDEProblem.jl => ProblemConstructor.jl} (95%) diff --git a/test/DeepBSDE.jl b/test/DeepBSDE.jl index 2f1a8d9..5684153 100644 --- a/test/DeepBSDE.jl +++ b/test/DeepBSDE.jl @@ -27,7 +27,7 @@ end f(X, u, σᵀ∇u, p, t) = Float32(0.0) μ_f(X, p, t) = zero(X) #Vector d x 1 σ_f(X, p, t) = Diagonal(ones(Float32, d)) |> Matrix #Matrix d x d - prob = PIDEProblem(μ_f, σ_f, x0, tspan, g, f) + prob = ParabolicPDEProblem(μ_f, σ_f, x0, tspan, g, f) hls = 10 + d #hidden layer size opt = Flux.Optimise.Adam(0.005) #optimizer @@ -69,7 +69,7 @@ end f(X, u, σᵀ∇u, p, t) = Float32(0.0) μ_f(X, p, t) = zero(X) #Vector d x 1 σ_f(X, p, t) = Diagonal(ones(Float32, d)) |> Matrix #Matrix d x d - prob = PIDEProblem(μ_f, σ_f, x0, tspan, g, f) + prob = ParabolicPDEProblem(μ_f, σ_f, x0, tspan, g, f) hls = 10 + d #hidden layer size opt = Flux.Optimise.Adam(0.005) #optimizer @@ -113,7 +113,7 @@ end g(X) = sum(X .^ 2) μ_f(X, p, t) = zero(X) #Vector d x 1 σ_f(X, p, t) = Diagonal(sigma * X) |> Matrix #Matrix d x d - prob = PIDEProblem(μ_f, σ_f, x0, tspan, g, f) + prob = ParabolicPDEProblem(μ_f, σ_f, x0, tspan, g, f) hls = 10 + d #hide layer size opt = Flux.Optimise.Adam(0.001) @@ -153,7 +153,7 @@ end f(X, u, σᵀ∇u, p, t) = u .- u .^ 3 μ_f(X, p, t) = zero(X) #Vector d x 1 σ_f(X, p, t) = Diagonal(ones(Float32, d)) |> Matrix #Matrix d x d - prob = PIDEProblem(μ_f, σ_f, x0, tspan, g, f) + prob = ParabolicPDEProblem(μ_f, σ_f, x0, tspan, g, f) hls = 20 + d #hidden layer size opt = Flux.Optimise.Adam(5^-3) #optimizer @@ -195,7 +195,7 @@ end f(X, u, σᵀ∇u, p, t) = -λ * sum(σᵀ∇u .^ 2) μ_f(X, p, t) = zero(X) #Vector d x 1 λ σ_f(X, p, t) = Diagonal(sqrt(2.0f0) * ones(Float32, d)) #Matrix d x d - prob = PIDEProblem(μ_f, σ_f, x0, tspan, g, f) + prob = ParabolicPDEProblem(μ_f, σ_f, x0, tspan, g, f) # TODO: This is a very large neural networks which size must be reduced. hls = 256 #hidden layer size @@ -266,7 +266,7 @@ end μ_f(X, p, t) = µc * X #Vector d x 1 σ_f(X, p, t) = σc * Diagonal(X) #Matrix d x d - prob = PIDEProblem(μ_f, σ_f, x0, tspan, g, f) + prob = ParabolicPDEProblem(μ_f, σ_f, x0, tspan, g, f) hls = 256 #hidden layer size opt = Flux.Optimise.Adam(0.008) #optimizer diff --git a/test/DeepBSDE_Han.jl b/test/DeepBSDE_Han.jl index ddafd5e..f915f64 100644 --- a/test/DeepBSDE_Han.jl +++ b/test/DeepBSDE_Han.jl @@ -28,7 +28,7 @@ end f(X, u, σᵀ∇u, p, t) = 0.0f0 # function from solved equation μ_f(X, p, t) = 0.0 σ_f(X, p, t) = 1.0 - prob = PIDEProblem(μ_f, σ_f, x0, tspan, g, f) + prob = ParabolicPDEProblem(μ_f, σ_f, x0, tspan, g, f) hls = 10 + d #hidden layer size opt = Flux.Optimise.Adam(0.005) #optimizer @@ -72,7 +72,7 @@ end f(X, u, σᵀ∇u, p, t) = 0.0f0 μ_f(X, p, t) = 0.0 σ_f(X, p, t) = 1.0 - prob = PIDEProblem(μ_f, σ_f, x0, tspan, g, f) + prob = ParabolicPDEProblem(μ_f, σ_f, x0, tspan, g, f) hls = 10 + d #hidden layer size opt = Flux.Optimise.Adam(0.005) #optimizer @@ -117,7 +117,7 @@ end g(X) = sum(X .^ 2) μ_f(X, p, t) = 0.0 σ_f(X, p, t) = Diagonal(sigma * X) - prob = PIDEProblem(μ_f, σ_f, x0, tspan, g, f) + prob = ParabolicPDEProblem(μ_f, σ_f, x0, tspan, g, f) hls = 10 + d #hide layer size opt = Flux.Optimise.Adam(0.001) @@ -159,7 +159,7 @@ end f(X, u, σᵀ∇u, p, t) = u .- u .^ 3 μ_f(X, p, t) = 0.0 σ_f(X, p, t) = 1.0 - prob = PIDEProblem(μ_f, σ_f, x0, tspan, g, f) + prob = ParabolicPDEProblem(μ_f, σ_f, x0, tspan, g, f) hls = 10 + d #hidden layer size opt = Flux.Optimise.Adam(5^-4) #optimizer @@ -205,7 +205,7 @@ end f(X, u, σᵀ∇u, p, t) = -λ * sum(σᵀ∇u .^ 2) μ_f(X, p, t) = 0.0 σ_f(X, p, t) = sqrt(2) - prob = PIDEProblem(μ_f, σ_f, x0, tspan, g, f) + prob = ParabolicPDEProblem(μ_f, σ_f, x0, tspan, g, f) hls = 12 + d #hidden layer size opt = Flux.Optimise.Adam(0.03) #optimizer @@ -278,7 +278,7 @@ end σc = 0.2f0 μ_f(X, p, t) = µc * X #Vector d x 1 σ_f(X, p, t) = σc * Diagonal(X) #Matrix d x d - prob = PIDEProblem(μ_f, σ_f, x0, tspan, g, f) + prob = ParabolicPDEProblem(μ_f, σ_f, x0, tspan, g, f) hls = 10 + d #hidden layer size opt = Flux.Optimise.Adam(0.008) #optimizer @@ -321,7 +321,7 @@ end # σ_f(X,p,t) = 1.0 # u_domain = -500:0.1:500 # A = -2:0.01:2 -# prob = PIDEProblem(g, f, μ_f, σ_f, x0, tspan ;A = A , x0_sample = u_domain) +# prob = ParabolicPDEProblem(g, f, μ_f, σ_f, x0, tspan ;A = A , x0_sample = u_domain) # hls = 10 + d # hidden layer size # opt = Flux.Optimise.Adam(0.005) # optimizer diff --git a/test/DeepSplitting.jl b/test/DeepSplitting.jl index 69fb32b..fd8b393 100644 --- a/test/DeepSplitting.jl +++ b/test/DeepSplitting.jl @@ -50,10 +50,10 @@ end opt = Flux.Optimise.Adam(0.01) #optimiser alg = DeepSplitting(nn, opt = opt) - f(y, z, v_y, v_z, ∇v_y, ∇v_z, p, t) = 0.0f0 .* v_y + f(y, v_y, ∇v_y, p, t) = 0.0f0 .* v_y # defining the problem - prob = PIDEProblem(μ, σ, x0, tspan, g, f) + prob = ParabolicPDEProblem(μ, σ, x0, tspan; g, f) # solving sol = solve(prob, alg, dt, verbose = true, @@ -94,10 +94,10 @@ end opt = Flux.Optimise.Adam(0.01) #optimiser alg = DeepSplitting(nn, opt = opt) - f(y, z, v_y, v_z, ∇v_y, ∇v_z, p, t) = 0.0f0 .* v_y #TODO: this fix is not nice + f(y, v_y, ∇v_y, p, t) = 0.0f0 .* v_y #TODO: this fix is not nice # defining the problem - prob = PIDEProblem(μ, σ, x0, tspan, g, f; x0_sample = x0_sample) + prob = ParabolicPDEProblem(μ, σ, x0, tspan; g, f, x0_sample = x0_sample) # solving sol = solve(prob, alg, dt, verbose = false, @@ -140,10 +140,10 @@ end opt = Flux.Optimise.Adam(0.01) #optimiser alg = DeepSplitting(nn, opt = opt) - f(y, z, v_y, v_z, ∇v_y, ∇v_z, p, t) = 0.0f0 .* v_y #TODO: this fix is not nice + f(y, v_y, ∇v_y, p, t) = 0.0f0 .* v_y #TODO: this fix is not nice # defining the problem - prob = PIDEProblem(μ, σ, x0, tspan, g, f; + prob = ParabolicPDEProblem(μ, σ, x0, tspan; g, f, x0_sample = x0_sample, neumann_bc = [-∂, ∂]) # solving @@ -198,10 +198,10 @@ end opt = Flux.Optimise.Adam(0.01) #optimiser alg = DeepSplitting(nn, opt = opt) - f(y, z, v_y, v_z, ∇v_y, ∇v_z, p, t) = r * v_y #TODO: this fix is not nice + f(y, v_y, ∇v_y, p, t) = r * v_y #TODO: this fix is not nice # defining the problem - prob = PIDEProblem(μ, σ, x0, tspan, g, f; + prob = ParabolicPDEProblem(μ, σ, x0, tspan; g, f, x0_sample = x0_sample) # solving sol = solve(prob, alg, dt, @@ -243,10 +243,10 @@ end X0 = fill(0.0f0, d) # initial point g(X) = 1.0f0 ./ (2.0f0 .+ 4.0f-1 * sum(X .^ 2, dims = 1)) # initial condition a(u) = u - u^3 - f(y, z, v_y, v_z, ∇v_y, ∇v_z, p, t) = -a.(v_y) # nonlocal nonlinear function + f(y, v_y, ∇v_y, p, t) = -a.(v_y) # nonlocal nonlinear function # defining the problem - prob = PIDEProblem(μ, σ, X0, tspan, g, f) + prob = ParabolicPDEProblem(μ, σ, X0, tspan; g, f) # solving @time sol = solve(prob, alg, @@ -290,10 +290,10 @@ end X0 = fill(0.0f0, d) # initial point g(X) = exp.(-0.25f0 * sum(X .^ 2, dims = 1)) # initial condition a(u) = u - u^3 - f(y, z, v_y, v_z, ∇v_y, ∇v_z, p, t) = a.(v_y) # nonlocal nonlinear function + f(y, v_y, ∇v_y, p, t) = a.(v_y) # nonlocal nonlinear function # defining the problem - prob = PIDEProblem(μ, σ, X0, tspan, g, f; neumann_bc = [-∂, ∂]) + prob = ParabolicPDEProblem(μ, σ, X0, tspan; g, f, neumann_bc = [-∂, ∂]) # solving @time sol = solve(prob, alg, @@ -337,10 +337,10 @@ if false X0 = repeat([1.0f0, 0.5f0], div(d, 2)) # initial point g(X) = sum(X .^ 2, dims = 1) # initial condition - f(y, z, v_y, v_z, ∇v_y, ∇v_z, p, t) = r * (v_y .- sum(y .* ∇v_y, dims = 1)) + f(y, v_y, ∇v_y, p, t) = r * (v_y .- sum(y .* ∇v_y, dims = 1)) # defining the problem - prob = PIDEProblem(μ, σ, X0, tspan, g, f) + prob = ParabolicPDEProblem(μ, σ, X0, tspan; g, f) # solving @time xs, ts, sol = solve(prob, alg, @@ -389,10 +389,10 @@ if false X0 = fill(0.0f0, d) # initial point g(X) = log.(5.0f-1 .+ 5.0f-1 * sum(X .^ 2, dims = 1)) # initial condition - f(y, z, v_y, v_z, ∇v_y, ∇v_z, p, t) = λ * sum(∇v_y .^ 2, dims = 1) + f(y, v_y, ∇v_y, p, t) = λ * sum(∇v_y .^ 2, dims = 1) # defining the problem - prob = PIDEProblem(μ, σ, X0, tspang, f) + prob = ParabolicPDEProblem(μ, σ, X0, tspan; g, f) # solving @time sol = solve(prob, alg, @@ -453,10 +453,10 @@ end µc = 0.02f0 σc = 0.2f0 - f(y, z, v_y, v_z, ∇v_y, ∇v_z, p, t) = -(1.0f0 - δ) * Q.(v_y) .* v_y .- R * v_y + f(y, v_y, ∇v_y, p, t) = -(1.0f0 - δ) * Q.(v_y) .* v_y .- R * v_y # defining the problem - prob = PIDEProblem(μ, σ, X0, tspan, g, f) + prob = ParabolicPDEProblem(μ, σ, X0, tspan; g, f) # solving @time sol = solve(prob, alg, diff --git a/test/MLP.jl b/test/MLP.jl index 3a09e0a..c242e5c 100644 --- a/test/MLP.jl +++ b/test/MLP.jl @@ -31,10 +31,10 @@ end alg = MLP() - f(y, z, v_y, v_z, ∇v_y, ∇v_z, p, t) = 0e0 .* v_y #TODO: this fix is not nice + f(y, v_y, ∇v_y, p, t) = 0e0 .* v_y #TODO: this fix is not nice # defining the problem - prob = PIDEProblem(μ, σ, x0, tspan, g, f) + prob = ParabolicPDEProblem(μ, σ, x0, tspan; g, f) # solving sol = solve(prob, alg, multithreading = false) u1 = sol.us[end] @@ -62,10 +62,10 @@ end alg = MLP() - f(y, z, v_y, v_z, ∇v_y, ∇v_z, p, t) = 0e0 .* v_y #TODO: this fix is not nice + f(y, v_y, ∇v_y, p, t) = 0e0 .* v_y #TODO: this fix is not nice # defining the problem - prob = PIDEProblem(μ, σ, x0, tspan, g, f) + prob = ParabolicPDEProblem(μ, σ, x0, tspan; g, f) # solving sol = solve(prob, alg, multithreading = true) u1 = sol.us[end] @@ -88,10 +88,10 @@ end # d = 10 x0 = fill(3e-1, d) alg = MLP() - f(y, z, v_y, v_z, ∇v_y, ∇v_z, p, t) = 0e0 .* v_y #TODO: this fix is not nice + f(y, v_y, ∇v_y, p, t) = 0e0 .* v_y #TODO: this fix is not nice # defining the problem - prob = PIDEProblem(μ, σ, x0, tspan, g, f, neumann_bc = neumann_bc) + prob = ParabolicPDEProblem(μ, σ, x0, tspan; g, f, neumann_bc = neumann_bc) # solving sol = solve(prob, alg) push!(u1s, sol.us[end]) @@ -118,10 +118,10 @@ end alg = MLP() - f(y, z, v_y, v_z, ∇v_y, ∇v_z, p, t) = r * v_y #TODO: this fix is not nice + f(y, v_y, ∇v_y, p, t) = r * v_y #TODO: this fix is not nice # defining the problem - prob = PIDEProblem(μ, σ, x, tspan, g, f) + prob = ParabolicPDEProblem(μ, σ, x, tspan; g, f) # solving sol = solve(prob, alg) u1 = sol.us[end] @@ -144,10 +144,10 @@ end X0 = fill(0.0f0, d) # initial point g(X) = 1.0f0 ./ (2.0f0 .+ 4.0f-1 * sum(X .^ 2)) # initial condition a(u) = u - u^3 - f(y, z, v_y, v_z, ∇v_y, ∇v_z, p, t) = -a.(v_y) + f(y, v_y, ∇v_y, p, t) = -a.(v_y) # defining the problem - prob = PIDEProblem(μ, σ, X0, tspan, g, f) + prob = ParabolicPDEProblem(μ, σ, X0, tspan; g, f) # solving sol = solve(prob, alg) u1 = sol.us[end] @@ -172,10 +172,10 @@ end X0 = fill(0.0f0, d) # initial point g(X) = exp.(-0.25f0 * sum(X .^ 2)) # initial condition a(u) = u - u^3 - f(y, z, v_y, v_z, ∇v_y, ∇v_z, p, t) = a.(v_y) # nonlocal nonlinear part of the + f(y, v_y, ∇v_y, p, t) = a.(v_y) # nonlocal nonlinear part of the # defining the problem - prob = PIDEProblem(μ, σ, X0, tspan, g, f; neumann_bc = neumann_bc) + prob = ParabolicPDEProblem(μ, σ, X0, tspan; g, f, neumann_bc = neumann_bc) # solving sol = solve(prob, alg) u1 = sol.us[end] @@ -211,10 +211,10 @@ end µc = 0.02f0 σc = 0.2f0 - f(y, z, v_y, v_z, ∇v_y, ∇v_z, p, t) = -(1.0f0 - δ) * Q.(v_y) .* v_y .- R * v_y + f(y, v_y, ∇v_y, p, t) = -(1.0f0 - δ) * Q.(v_y) .* v_y .- R * v_y # defining the problem - prob = PIDEProblem(μ, σ, X0, tspan, g, f) + prob = ParabolicPDEProblem(μ, σ, X0, tspan; g, f) # solving sol = solve(prob, alg) diff --git a/test/NNStopping.jl b/test/NNStopping.jl index a48ffa1..8858e30 100644 --- a/test/NNStopping.jl +++ b/test/NNStopping.jl @@ -22,7 +22,7 @@ function g(x, t) return exp(-r * t) * (max(maximum(x) - K, 0)) end -prob = PIDEProblem(f, sigma, u0, tspan; payoff = g) +prob = ParabolicPDEProblem(f, sigma, u0, tspan; payoff = g) models = [Chain(Dense(d + 1, 32, tanh), BatchNorm(32, tanh), Dense(32, 1, sigmoid)) for i in 1:N] opt = Flux.Optimisers.Adam(0.01) @@ -59,7 +59,7 @@ K = 100.00 function g(x, t) return exp(-r * t) * (max(K - mean(x), 0)) end -prob = PIDEProblem(f, sigma, u0, tspan; payoff = g) +prob = ParabolicPDEProblem(f, sigma, u0, tspan; payoff = g) models = [Chain(Dense(d + 1, 32, tanh), BatchNorm(32, tanh), Dense(32, 1, sigmoid)) for i in 1:N] opt = Flux.Optimisers.Adam(0.01) diff --git a/test/PIDEProblem.jl b/test/ProblemConstructor.jl similarity index 95% rename from test/PIDEProblem.jl rename to test/ProblemConstructor.jl index 1f35722..f555f90 100644 --- a/test/PIDEProblem.jl +++ b/test/ProblemConstructor.jl @@ -82,7 +82,7 @@ end g(X) = sum(X .^ 2) # terminal condition f_semilinear(X, u, σᵀ∇u, p, t) = Float32(0.0) - prob = PIDEProblem(μ_f, σ_f, x0, tspan, g, f_semilinear) + prob = ParabolicPDEProblem(μ_f, σ_f, x0, tspan; g, f = f_semilinear) @test prob.f == f_semilinear @test prob.g == g @@ -112,7 +112,7 @@ end return exp(-r * t) * (max(maximum(x) - K, 0)) end - prob = PIDEProblem(mu, sigma, u0, tspan; payoff = payoff) + prob = ParabolicPDEProblem(mu, sigma, u0, tspan; payoff = payoff) @test prob.f == nothing @test prob.g == nothing @@ -137,7 +137,7 @@ end 1.77 .* x .- 0.015 .* x .^ 3 end - prob = PIDEProblem(μ, σ, nothing, tspan, g; xspan) + prob = ParabolicPDEProblem(μ, σ, nothing, tspan; g, xspan) @test prob.f == nothing @test prob.g == g @@ -178,7 +178,7 @@ end p_mu = γ_mu_prototype, p_phi = γ_phi_prototype) - prob = PIDEProblem(mu, sigma, nothing, tspan, g; + prob = ParabolicPDEProblem(mu, sigma, nothing, tspan; g, xspan, p_domain = p_domain, p_prototype = p_prototype) From 3fc048a900e8bbd71632d1149e7a2a0429e76328 Mon Sep 17 00:00:00 2001 From: Ashutosh Bharambe Date: Sun, 4 Feb 2024 01:25:30 +0530 Subject: [PATCH 29/37] test: make chain f64 in NNStopping tests --- test/NNStopping.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/NNStopping.jl b/test/NNStopping.jl index 8858e30..91ece07 100644 --- a/test/NNStopping.jl +++ b/test/NNStopping.jl @@ -23,7 +23,7 @@ function g(x, t) end prob = ParabolicPDEProblem(f, sigma, u0, tspan; payoff = g) -models = [Chain(Dense(d + 1, 32, tanh), BatchNorm(32, tanh), Dense(32, 1, sigmoid)) +models = [Chain(Dense(d + 1, 32, tanh), BatchNorm(32, tanh), Dense(32, 1, sigmoid)) |> f64 for i in 1:N] opt = Flux.Optimisers.Adam(0.01) alg = NNStopping(models, opt) From ce172b4a24521ab582c4b3472e52b498215d2de1 Mon Sep 17 00:00:00 2001 From: ashutosh-b-b Date: Sun, 4 Feb 2024 09:11:24 +0530 Subject: [PATCH 30/37] test: fix bug in runtests --- test/runtests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index d4829d4..016c8fc 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,7 +1,7 @@ using SafeTestsets, Test @testset "HighDimPDE" begin @time @safetestset "Quality Assurance" include("qa.jl") - @time @safetestset "PIDEProblem" include("PIDEProblem.jl") + @time @safetestset "ProblemConstructors" include("ProblemConstructor.jl") @time @safetestset "reflect" include("reflect.jl") @time @safetestset "reflect" include("reflect.jl") @time @safetestset "MLP" include("MLP.jl") From 9ba145589fa6766522f0d142686d2ac491fdb030 Mon Sep 17 00:00:00 2001 From: ashutosh-b-b Date: Sun, 4 Feb 2024 11:27:59 +0530 Subject: [PATCH 31/37] test: increase error tolerance --- test/NNStopping.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/NNStopping.jl b/test/NNStopping.jl index 91ece07..6bceb95 100644 --- a/test/NNStopping.jl +++ b/test/NNStopping.jl @@ -35,7 +35,7 @@ sol = solve(prob, maxiters = 1000, verbose = true) analytical_sol = 11.278 # Ref [1] -@test abs(analytical_sol - sol.payoff) < 0.5 +@test abs(analytical_sol - sol.payoff) < 0.7 # Put basket option in Dupire’s local volatility model d = 5 From 125182dd89c3ed31d378d88c010fd82b8de5811d Mon Sep 17 00:00:00 2001 From: Ashutosh Bharambe Date: Fri, 9 Feb 2024 15:25:30 +0530 Subject: [PATCH 32/37] docs: update docs for MLP and DeepSplitting --- docs/src/tutorials/deepsplitting.md | 2 +- docs/src/tutorials/mlp.md | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/src/tutorials/deepsplitting.md b/docs/src/tutorials/deepsplitting.md index 0c94410..60ac77c 100644 --- a/docs/src/tutorials/deepsplitting.md +++ b/docs/src/tutorials/deepsplitting.md @@ -22,7 +22,7 @@ g(x) = exp.(- sum(x.^2, dims=1) ) # initial condition σ(x, p, t) = 0.1f0 # diffusion coefficients x0_sample = UniformSampling(fill(-5f-1, d), fill(5f-1, d)) f(x, y, v_x, v_y, ∇v_x, ∇v_y, p, t) = v_x .* (1f0 .- v_y) -prob = prob = PIDEProblem(μ, σ, x0, tspan, g, f; x0_sample = x0_sample) +prob = PIDEProblem(μ, σ, x0, tspan, g, f; x0_sample = x0_sample) ## Definition of the neural network to use using Flux # needed to define the neural network diff --git a/docs/src/tutorials/mlp.md b/docs/src/tutorials/mlp.md index 6156571..3774808 100644 --- a/docs/src/tutorials/mlp.md +++ b/docs/src/tutorials/mlp.md @@ -16,8 +16,8 @@ x0 = fill(0.,d) # initial point g(x) = exp(- sum(x.^2) ) # initial condition μ(x, p, t) = 0.0 # advection coefficients σ(x, p, t) = 0.1 # diffusion coefficients -f(x, y, v_x, v_y, ∇v_x, ∇v_y, p, t) = max(0.0, v_x) * (1 - max(0.0, v_x)) # nonlocal nonlinear part of the -prob = PIDEProblem(μ, σ, x0, tspan, g, f) # defining the problem +f(x, v_x, ∇v_x, p, t) = max(0.0, v_x) * (1 - max(0.0, v_x)) # nonlocal nonlinear part of the +prob = ParabolicPDEProblem(μ, σ, x0, tspan, g, f) # defining the problem ## Definition of the algorithm alg = MLP() # defining the algorithm. We use the Multi Level Picard algorithm From 953476130eb99c1cdf8871ca05311d0a9d4dc3ea Mon Sep 17 00:00:00 2001 From: ashutosh-b-b Date: Fri, 9 Feb 2024 18:57:08 +0530 Subject: [PATCH 33/37] docs: add eqn for ParabolicPDEProblem --- src/HighDimPDE.jl | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/src/HighDimPDE.jl b/src/HighDimPDE.jl index 82f7335..769c8c3 100644 --- a/src/HighDimPDE.jl +++ b/src/HighDimPDE.jl @@ -131,6 +131,13 @@ end $(SIGNATURES) Defines a Parabolic Partial Differential Equation of the form: +```math +\\begin{aligned} + \\frac{du}{dt} &= \\tfrac{1}{2} \\text{Tr}(\\sigma \\sigma^T) \\Delta u(x, t) + \\mu \\nabla u(x, t) \\\\ + &\\quad + f(x, u(x, t), ( \\nabla_x u )(x, t), p, t) +\\end{aligned} +``` + - Semilinear Parabolic Partial Differential Equation * f -> f(X, u, σᵀ∇u, p, t) - Kolmogorov Differential Equation From a7358f7dbd32110b5482478679bd0c8df9e67466 Mon Sep 17 00:00:00 2001 From: ashutosh-b-b Date: Fri, 9 Feb 2024 18:57:37 +0530 Subject: [PATCH 34/37] docs: add section on problems and add notes on defining problems --- docs/pages.jl | 1 + docs/src/index.md | 15 ++++++++++++++- docs/src/problems.md | 8 ++++++++ docs/src/tutorials/nnstopping.md | 2 ++ 4 files changed, 25 insertions(+), 1 deletion(-) create mode 100644 docs/src/problems.md diff --git a/docs/pages.jl b/docs/pages.jl index ce1a4bf..5a0650e 100644 --- a/docs/pages.jl +++ b/docs/pages.jl @@ -1,6 +1,7 @@ pages = [ "Home" => "index.md", "Getting started" => "getting_started.md", + "Problems" => "problems.md", "Solver Algorithms" => ["MLP.md", "DeepSplitting.md", "DeepBSDE.md", diff --git a/docs/src/index.md b/docs/src/index.md index 49f6a9d..32c70e0 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -1,8 +1,9 @@ # HighDimPDE.jl -**HighDimPDE.jl** is a Julia package to **solve Highly Dimensional non-linear, non-local PDEs** of the form +**HighDimPDE.jl** is a Julia package to **solve Highly Dimensional non-linear, non-local PDEs** of the forms: +1. Partial Integro Differential Equations: ```math \begin{aligned} (\partial_t u)(t,x) &= \int_{\Omega} f\big(t,x,{\bf x}, u(t,x),u(t,{\bf x}), ( \nabla_x u )(t,x ),( \nabla_x u )(t,{\bf x} ) \big) \, d{\bf x} \\ @@ -12,6 +13,18 @@ where $u \colon [0,T] \times \Omega \to \R$, $\Omega \subseteq \R^d$ is subject to initial and boundary conditions, and where $d$ is large. +2. Parabolic Partial Differential Equations: +```math +\begin{aligned} + (\partial_t u)(t,x) &= f\big(t,x, u(t,x), ( \nabla_x u )(t,x )\big) + + \big\langle \mu(t,x), ( \nabla_x u )( t,x ) \big\rangle + \tfrac{1}{2} \text{Trace} \big(\sigma(t,x) [ \sigma(t,x) ]^* ( \text{Hess}_x u)(t, x ) \big). +\end{aligned} +``` + +where $u \colon [0,T] \times \Omega \to \R$, $\Omega \subseteq \R^d$ is subject to initial and boundary conditions, and where $d$ is large. + +!!! note + The difference between the two problems is that in Partial Integro Differential Equations, the integrand is integrated over **x**, while in Parabolic Integro Differential Equations, the function `f` is just evaluated for `x`. **HighDimPDE.jl** implements solver algorithms that break down the curse of dimensionality, including diff --git a/docs/src/problems.md b/docs/src/problems.md new file mode 100644 index 0000000..a10555c --- /dev/null +++ b/docs/src/problems.md @@ -0,0 +1,8 @@ +```@docs +PIDEProblem +ParabolicPDEProblem +``` + +!!! note + While choosing to define a PDE using `PIDEProblem`, not that the function being integrated `f` is a function of `f(x, y, v_x, v_y, ∇v_x, ∇v_y)` out of which `y` is the integrating variable and `x` is constant throughout the integration. + If a PDE has no integral and the non linear term `f` is just evaluated as `f(x, v_x, ∇v_x)` then we suggest using `ParabolicPDEProblem` \ No newline at end of file diff --git a/docs/src/tutorials/nnstopping.md b/docs/src/tutorials/nnstopping.md index 54fb78a..25edad2 100644 --- a/docs/src/tutorials/nnstopping.md +++ b/docs/src/tutorials/nnstopping.md @@ -1,3 +1,5 @@ +# `NNStopping` + ## Solving for optimal strategy and expected payoff of a Bermudan Max-Call option We will calculate optimal strategy for Bermudan Max-Call option with following drift, diffusion and payoff: From 53cc92baf63822df7a279e721a9414d6b8f99338 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Fri, 9 Feb 2024 08:39:10 -0500 Subject: [PATCH 35/37] Update docs/src/problems.md --- docs/src/problems.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/problems.md b/docs/src/problems.md index a10555c..6ae5eeb 100644 --- a/docs/src/problems.md +++ b/docs/src/problems.md @@ -4,5 +4,5 @@ ParabolicPDEProblem ``` !!! note - While choosing to define a PDE using `PIDEProblem`, not that the function being integrated `f` is a function of `f(x, y, v_x, v_y, ∇v_x, ∇v_y)` out of which `y` is the integrating variable and `x` is constant throughout the integration. + While choosing to define a PDE using `PIDEProblem`, note that the function being integrated `f` is a function of `f(x, y, v_x, v_y, ∇v_x, ∇v_y)` out of which `y` is the integrating variable and `x` is constant throughout the integration. If a PDE has no integral and the non linear term `f` is just evaluated as `f(x, v_x, ∇v_x)` then we suggest using `ParabolicPDEProblem` \ No newline at end of file From 0f1c148a521272394dd12fb8fa8e32d1ea9b5b49 Mon Sep 17 00:00:00 2001 From: ashutosh-b-b Date: Fri, 9 Feb 2024 19:53:38 +0530 Subject: [PATCH 36/37] docs: make docstrings have TYPEDSIGNATURES --- src/DeepBSDE.jl | 4 +++- src/DeepBSDE_Han.jl | 11 +++++++++++ src/DeepSplitting.jl | 2 +- src/MLP.jl | 2 +- src/NNStopping.jl | 16 ++++++++++++++++ 5 files changed, 32 insertions(+), 3 deletions(-) diff --git a/src/DeepBSDE.jl b/src/DeepBSDE.jl index 56d17f9..c421b21 100644 --- a/src/DeepBSDE.jl +++ b/src/DeepBSDE.jl @@ -59,7 +59,7 @@ end DeepBSDE(u0, σᵀ∇u; opt = Flux.Optimise.Adam(0.1)) = DeepBSDE(u0, σᵀ∇u, opt) """ -$(SIGNATURES) +$(TYPEDSIGNATURES) Returns a `PIDESolution` object. @@ -73,6 +73,8 @@ Returns a `PIDESolution` object. [DifferentialEquations.jl doc](https://diffeq.sciml.ai/stable/solvers/sde_solve/). - `limits`: if `true`, upper and lower limits will be calculated, based on [Deep Primal-Dual algorithm for BSDEs](https://papers.ssrn.com/sol3/papers.cfm?abstract_id=3071506). +- `maxiters`: The number of training epochs. Defaults to `300` +- `trajectories`: The number of trajectories simulated for training. Defaults to `100` - Extra keyword arguments passed to `solve` will be further passed to the SDE solver. """ function DiffEqBase.solve(prob::ParabolicPDEProblem, diff --git a/src/DeepBSDE_Han.jl b/src/DeepBSDE_Han.jl index eb785b7..75744c8 100644 --- a/src/DeepBSDE_Han.jl +++ b/src/DeepBSDE_Han.jl @@ -1,4 +1,15 @@ # called whenever sdealg is not specified. +""" +$(TYPEDSIGNATURES) + +Returns a `PIDESolution` object. + +# Arguments: +- `maxiters`: The number of training epochs. Defaults to `300` +- `trajectories`: The number of trajectories simulated for training. Defaults to `100` + +To use [SDE Algorithms](https://diffeq.sciml.ai/stable/solvers/sde_solve/) use [`DeepBSDE`](@ref) +""" function DiffEqBase.solve(prob::ParabolicPDEProblem, alg::DeepBSDE; dt, diff --git a/src/DeepSplitting.jl b/src/DeepSplitting.jl index 9541c52..91ca483 100644 --- a/src/DeepSplitting.jl +++ b/src/DeepSplitting.jl @@ -51,7 +51,7 @@ function DeepSplitting(nn; end """ -$(SIGNATURES) +$(TYPEDSIGNATURES) Returns a `PIDESolution` object. diff --git a/src/MLP.jl b/src/MLP.jl index a6e7e62..c640649 100644 --- a/src/MLP.jl +++ b/src/MLP.jl @@ -21,7 +21,7 @@ end MLP(; M = 4, L = 4, K = 10, mc_sample = NoSampling()) = MLP(M, L, K, mc_sample) """ -$(SIGNATURES) +$(TYPEDSIGNATURES) Returns a `PIDESolution` object. diff --git a/src/NNStopping.jl b/src/NNStopping.jl index 39d1793..0f49f25 100644 --- a/src/NNStopping.jl +++ b/src/NNStopping.jl @@ -54,6 +54,22 @@ function (model::NNStoppingModelArray)(X, G) broadcast((x, m) -> m(x), eachslice(XG, dims = 2)[2:end], model.ms) end +""" +$(TYPEDSIGNATURES) + +Returns a NamedTuple with `payoff` and `stopping_time` + +Arguments: +- `sdealg`: a SDE solver from [DifferentialEquations.jl](https://diffeq.sciml.ai/stable/solvers/sde_solve/). + If not provided, the plain vanilla [DeepBSDE](https://arxiv.org/abs/1707.02568) method will be applied. + If provided, the SDE associated with the PDE problem will be solved relying on + methods from DifferentialEquations.jl, using [Ensemble solves](https://diffeq.sciml.ai/stable/features/ensemble/) + via `sdealg`. Check the available `sdealg` on the + [DifferentialEquations.jl doc](https://diffeq.sciml.ai/stable/solvers/sde_solve/). +- `maxiters`: The number of training epochs. Defaults to `300` +- `trajectories`: The number of trajectories simulated for training. Defaults to `100` +- Extra keyword arguments passed to `solve` will be further passed to the SDE solver. +""" function DiffEqBase.solve(prob::ParabolicPDEProblem, pdealg::NNStopping, sdealg; From 66372d38755269e88789351da6d1c2d04276d8e0 Mon Sep 17 00:00:00 2001 From: ashutosh-b-b Date: Fri, 9 Feb 2024 19:54:19 +0530 Subject: [PATCH 37/37] docs: explicilty add section for supported problems for each solver --- docs/src/DeepBSDE.md | 3 +++ docs/src/DeepSplitting.md | 4 ++++ docs/src/MLP.md | 4 ++++ docs/src/NNStopping.md | 3 +++ 4 files changed, 14 insertions(+) diff --git a/docs/src/DeepBSDE.md b/docs/src/DeepBSDE.md index 95e4f8d..428f4c3 100644 --- a/docs/src/DeepBSDE.md +++ b/docs/src/DeepBSDE.md @@ -1,5 +1,8 @@ # [The `DeepBSDE` algorithm](@id deepbsde) +### Problems Supported: +1. [`ParabolicPDEProblem`](@ref) + ```@autodocs Modules = [HighDimPDE] Pages = ["DeepBSDE.jl", "DeepBSDE_Han.jl"] diff --git a/docs/src/DeepSplitting.md b/docs/src/DeepSplitting.md index bd14ffd..2639036 100644 --- a/docs/src/DeepSplitting.md +++ b/docs/src/DeepSplitting.md @@ -1,5 +1,9 @@ # [The `DeepSplitting` algorithm](@id deepsplitting) +### Problems Supported: +1. [`PIDEProblem`](@ref) +2. [`ParabolicPDEProblem`](@ref) + ```@autodocs Modules = [HighDimPDE] Pages = ["DeepSplitting.jl"] diff --git a/docs/src/MLP.md b/docs/src/MLP.md index 61e6bee..3e603da 100644 --- a/docs/src/MLP.md +++ b/docs/src/MLP.md @@ -1,5 +1,9 @@ # [The `MLP` algorithm](@id mlp) +### Problems Supported: +1. [`PIDEProblem`](@ref) +2. [`ParabolicPDEProblem`](@ref) + ```@autodocs Modules = [HighDimPDE] Pages = ["MLP.jl"] diff --git a/docs/src/NNStopping.md b/docs/src/NNStopping.md index 5fb4e37..df773dc 100644 --- a/docs/src/NNStopping.md +++ b/docs/src/NNStopping.md @@ -1,5 +1,8 @@ # [The `NNStopping` algorithm](@id nn_stopping) +### Problems Supported: +1. [`ParabolicPDEProblem`](@ref) + ```@autodocs Modules = [HighDimPDE] Pages = ["NNStopping.jl"]