diff --git a/docs/pages.jl b/docs/pages.jl index c3b4f2a..5a0650e 100644 --- a/docs/pages.jl +++ b/docs/pages.jl @@ -1,13 +1,16 @@ pages = [ "Home" => "index.md", "Getting started" => "getting_started.md", + "Problems" => "problems.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", ] 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 c041454..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"] @@ -62,14 +66,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/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 new file mode 100644 index 0000000..df773dc --- /dev/null +++ b/docs/src/NNStopping.md @@ -0,0 +1,32 @@ +# [The `NNStopping` algorithm](@id nn_stopping) + +### Problems Supported: +1. [`ParabolicPDEProblem`](@ref) + +```@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(X_\tau, \tau)]\rbrace +``` +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. 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/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..6ae5eeb --- /dev/null +++ b/docs/src/problems.md @@ -0,0 +1,8 @@ +```@docs +PIDEProblem +ParabolicPDEProblem +``` + +!!! note + 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 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..60ac77c 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 = 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..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(g, f, μ, σ, x0, tspan) # 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 @@ -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/docs/src/tutorials/nnstopping.md b/docs/src/tutorials/nnstopping.md new file mode 100644 index 0000000..25edad2 --- /dev/null +++ b/docs/src/tutorials/nnstopping.md @@ -0,0 +1,52 @@ +# `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: +```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(x, t) + return exp(-r * t) * (max(maximum(x) - K, 0)) +end + +``` +We then define a `PIDEProblem` with no non linear term: +```julia +prob = PIDEProblem(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) +``` 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..c421b21 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) @@ -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,9 +73,11 @@ 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::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..75744c8 100644 --- a/src/DeepBSDE_Han.jl +++ b/src/DeepBSDE_Han.jl @@ -1,5 +1,16 @@ # called whenever sdealg is not specified. -function DiffEqBase.solve(prob::PIDEProblem, +""" +$(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, abstol = 1.0f-6, diff --git a/src/DeepSplitting.jl b/src/DeepSplitting.jl index 939b20e..91ca483 100644 --- a/src/DeepSplitting.jl +++ b/src/DeepSplitting.jl @@ -51,7 +51,7 @@ function DeepSplitting(nn; end """ -$(SIGNATURES) +$(TYPEDSIGNATURES) Returns a `PIDESolution` object. @@ -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/HighDimPDE.jl b/src/HighDimPDE.jl index 23e6d74..5eb3568 100644 --- a/src/HighDimPDE.jl +++ b/src/HighDimPDE.jl @@ -28,6 +28,23 @@ function Base.show(io::IO, A::AbstractPDEProblem) show(io, A.tspan) end +include("MCSample.jl") + +struct PIDEProblem{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) @@ -52,58 +69,159 @@ with `` u(x,0) = g(x)``. * `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 - 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 +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 -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), - "Type of non linear function `f(x)` must type of x") - else - throw(e) - 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 + * 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` : nonlinear function, of the form `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 ParabolicPDEProblem(μ, + σ, + 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." + + !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") + + # Wrap kwargs : + kw = NamedTuple(kw) + 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 + 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(g(x)), + ParabolicPDEProblem{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, @@ -143,14 +261,14 @@ function Base.show(io::IO, A::PIDESolution) show(io, A.us) end -include("MCSample.jl") include("reflect.jl") 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, ParabolicPDEProblem, PIDESolution, DeepSplitting, DeepBSDE, MLP, NNStopping export NormalSampling, UniformSampling, NoSampling, solve end diff --git a/src/MLP.jl b/src/MLP.jl index 1e68a7b..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. @@ -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 new file mode 100644 index 0000000..0f49f25 --- /dev/null +++ b/src/NNStopping.jl @@ -0,0 +1,149 @@ +""" +```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(x, t) + return exp(-r * t) * (max(maximum(x) - K, 0)) +end + +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) +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 +end + +struct NNStoppingModelArray{M} + ms::M +end + +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 + +""" +$(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; + verbose = false, + maxiters = 300, + trajectories = 100, + dt = eltype(prob.tspan)(0), + ensemblealg = EnsembleThreads(), + kwargs...) + g = prob.kwargs[:payoff] + + sde_prob = SDEProblem(prob.μ, prob.σ, prob.x, 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 + ts = tspan[1]:dt:tspan[2] + N = length(ts) - 1 + + 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) + 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 -1 * mean(un_gs) + end + + 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 + Flux.update!(opt_state, m, gs[1]) + l = nn_loss(m) + verbose && @info "Current Epoch : $epoch Current Loss :$l" + end + + 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 + + ns = map(un_gs) do un + argmax(cumsum(un) + un .>= 1) + end + + tss = getindex.(Ref(ts), ns .+ 1) + payoff = mean(map((t, u) -> g(u(t), t), tss, sim)) + return (payoff = payoff, stopping_time = tss) +end diff --git a/test/DeepBSDE.jl b/test/DeepBSDE.jl index dfffefb..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(g, f, μ_f, σ_f, x0, tspan) + 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(g, f, μ_f, σ_f, x0, tspan) + 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(g, f, μ_f, σ_f, x0, tspan) + 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(g, f, μ_f, σ_f, x0, tspan) + 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(g, f, μ_f, σ_f, x0, tspan) + 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(g, f, μ_f, σ_f, x0, tspan) + 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 a3d859b..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(g, f, μ_f, σ_f, x0, tspan) + 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(g, f, μ_f, σ_f, x0, tspan) + 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(g, f, μ_f, σ_f, x0, tspan) + 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(g, f, μ_f, σ_f, x0, tspan) + 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(g, f, μ_f, σ_f, x0, tspan) + 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(g, f, μ_f, σ_f, x0, tspan) + 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 e401d22..fd8b393 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 @@ -47,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(g, f, μ, σ, x0, tspan) + prob = ParabolicPDEProblem(μ, σ, x0, tspan; g, f) # solving sol = solve(prob, alg, dt, verbose = true, @@ -91,15 +94,15 @@ 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(g, f, μ, σ, x0, tspan, x0_sample = x0_sample) + prob = ParabolicPDEProblem(μ, σ, x0, tspan; g, f, x0_sample = x0_sample) # solving 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)) @@ -137,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(g, f, μ, σ, x0, tspan, + prob = ParabolicPDEProblem(μ, σ, x0, tspan; g, f, x0_sample = x0_sample, neumann_bc = [-∂, ∂]) # solving @@ -195,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(g, f, μ, σ, x0, tspan, + prob = ParabolicPDEProblem(μ, σ, x0, tspan; g, f, x0_sample = x0_sample) # solving sol = solve(prob, alg, dt, @@ -235,15 +238,15 @@ 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 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(g, f, μ, σ, X0, tspan) + prob = ParabolicPDEProblem(μ, σ, X0, tspan; g, f) # solving @time sol = solve(prob, alg, @@ -282,15 +285,15 @@ 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 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(g, f, μ, σ, X0, tspan, neumann_bc = [-∂, ∂]) + prob = ParabolicPDEProblem(μ, σ, X0, tspan; g, f, neumann_bc = [-∂, ∂]) # solving @time sol = solve(prob, alg, @@ -330,14 +333,14 @@ 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 - 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(g, f, μ, σ, X0, tspan) + prob = ParabolicPDEProblem(μ, σ, X0, tspan; g, f) # solving @time xs, ts, sol = solve(prob, alg, @@ -382,14 +385,14 @@ 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 - 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(g, f, μ, σ, X0, tspan) + prob = ParabolicPDEProblem(μ, σ, X0, tspan; g, f) # solving @time sol = solve(prob, alg, @@ -431,7 +434,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 @@ -450,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(g, f, μ, σ, X0, tspan) + prob = ParabolicPDEProblem(μ, σ, X0, tspan; g, f) # solving @time sol = solve(prob, alg, @@ -520,7 +523,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) * @@ -533,7 +536,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, @@ -576,7 +579,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 @@ -584,7 +587,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..c242e5c 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)) @@ -30,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(g, f, μ, σ, x0, tspan) + prob = ParabolicPDEProblem(μ, σ, x0, tspan; g, f) # solving sol = solve(prob, alg, multithreading = false) u1 = sol.us[end] @@ -61,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(g, f, μ, σ, x0, tspan) + prob = ParabolicPDEProblem(μ, σ, x0, tspan; g, f) # solving sol = solve(prob, alg, multithreading = true) u1 = sol.us[end] @@ -87,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(g, f, μ, σ, x0, tspan, neumann_bc = neumann_bc) + prob = ParabolicPDEProblem(μ, σ, x0, tspan; g, f, neumann_bc = neumann_bc) # solving sol = solve(prob, alg) push!(u1s, sol.us[end]) @@ -117,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(g, f, μ, σ, x, tspan) + prob = ParabolicPDEProblem(μ, σ, x, tspan; g, f) # solving sol = solve(prob, alg) u1 = sol.us[end] @@ -143,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(g, f, μ, σ, X0, tspan) + prob = ParabolicPDEProblem(μ, σ, X0, tspan; g, f) # solving sol = solve(prob, alg) u1 = sol.us[end] @@ -171,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(g, f, μ, σ, X0, tspan, neumann_bc = neumann_bc) + prob = ParabolicPDEProblem(μ, σ, X0, tspan; g, f, neumann_bc = neumann_bc) # solving sol = solve(prob, alg) u1 = sol.us[end] @@ -210,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(g, f, μ, σ, X0, tspan) + prob = ParabolicPDEProblem(μ, σ, X0, tspan; g, f) # solving sol = solve(prob, alg) @@ -268,7 +269,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 +299,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 new file mode 100644 index 0000000..6bceb95 --- /dev/null +++ b/test/NNStopping.jl @@ -0,0 +1,81 @@ +using Test, Flux, StochasticDiffEq, LinearAlgebra +println("Optimal Stopping Time Test") +using HighDimPDE +using Statistics + +using Random +Random.seed!(101) +# Bermudan Max-Call Standard Example +d = 3 +r = 0.05 +beta = 0.2 +T = 3.0 +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(x, t) + return exp(-r * t) * (max(maximum(x) - K, 0)) +end + +prob = ParabolicPDEProblem(f, sigma, u0, tspan; payoff = g) +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) +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.7 + +# 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) + +tspan = (0.0f0, T) +N = 10 +dt = T / (N) +K = 100.00 +function g(x, t) + return exp(-r * t) * (max(K - mean(x), 0)) +end +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) +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) +=# diff --git a/test/ProblemConstructor.jl b/test/ProblemConstructor.jl new file mode 100644 index 0000000..f555f90 --- /dev/null +++ b/test/ProblemConstructor.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 = ParabolicPDEProblem(μ_f, σ_f, x0, tspan; g, f = 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 = ParabolicPDEProblem(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(u0, 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 = ParabolicPDEProblem(μ, σ, 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 = ParabolicPDEProblem(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 145ffe2..016c8fc 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,9 +1,11 @@ using SafeTestsets, Test @testset "HighDimPDE" begin @time @safetestset "Quality Assurance" include("qa.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") @time @safetestset "Deep Splitting" include("DeepSplitting.jl") @time @safetestset "MC Sample" include("MCSample.jl") + @time @safetestset "NNStopping" include("NNStopping.jl") end