+using Flux, StochasticDiffEq +using NeuralNetDiffEq +using Plots +using LinearAlgebra ++ + + +
The one dimensional backward kolmogorov equation that we are going to deal with is of the form :
+\[ +\begin{equation} + \frac{\partial p}{\partial t} = \mu(x)\frac{\partial p}{\partial x} + \frac{1}{2}{\sigma^2}(x)\frac{\partial^2 p}{\partial x^2} +\end{equation} +\]
+Consider a Black Scholes Partial Differential Equation of the form :
+\[ +\begin{equation} + \frac{\partial V}{\partial t} + rS\frac{\partial V}{\partial S} + \frac{1}{2}{\sigma^2}{S^2}\frac{\partial^2 V}{\partial S^2} -rV = 0 +\end{equation} +\]
+In order to make the above equation in the form of the Backward - Kolmogorov PDE we should substitute
+\[ +\begin{equation} + V(S , t) = e^{rt}p(S , t) +\end{equation} +\begin{equation} + p(S , T) = max\{ x - \mathcal{K} , 0 \} +\end{equation} +\]
+We should start defining the terminal condition for our equation:
+ + ++function phi(xi) + y = Float64[] + K = 100 + for x in eachcol(xi) + val = max(K - maximum(x) , 0.00) + y = push!(y , val) + end + y = reshape(y , 1 , size(y)[1] ) + return y +end ++ + +
+phi (generic function with 1 method) ++ + +
Now we shall define the problem : We will start with a simple one dimensional Black-Scholes PDE and then move to the higher dimensions We will define the σ and μ by comparing it to the orignal equation.
+ + ++d = 1 +r = 0.06 +sigma = 0.4 +xspan = (95.00 , 110.0) +tspan = (0.0 , 1.0) +σ(du , u , p , t) = du .= sigma.*u +μ(du , u , p , t) = du .= r.*u ++ + +
+μ (generic function with 1 method) ++ + +
Now once we have defined our problem it is necessary to define the parameters for the solver.
+ + ++sdealg = EM() +ensemblealg = EnsembleThreads() +dt = 0.01 +dx = 0.01 +trajectories = 100000 ++ + +
+100000 ++ + +
Now lets define our model m and the optimiser
+ + ++m = Chain(Dense(d, 32, leakyrelu),Dense(32, 128, leakyrelu),Dense(128 , 32 , leakyrelu) , Dense(32 , 1)) +m = fmap(cu , m) ++ + +
+ERROR: UndefVarError: cu not defined ++ + + +
+prob = KolmogorovPDEProblem(μ , σ , phi , xspan , tspan, d) +opt = Flux.ADAM(0.01) ++ + +
+ADAM(0.01, (0.9, 0.999), IdDict{Any,Any}())
+
+
+
+And then finally call the solver
+ + ++sol = solve(prob, NeuralNetDiffEq.NNKolmogorov(m, opt, sdealg, ensemblealg), verbose = true, dt = dt, + dx = dx , trajectories = trajectories , abstol=1e-6, maxiters = 3000) ++ + +
+Current loss is: 561.711498088225 +Current loss is: 386.22229456357377 +Current loss is: 296.80721579035924 +Current loss is: 324.38432342387557 +Current loss is: 349.2953855668709 +Current loss is: 335.98255012439193 +Current loss is: 311.27680726367424 +Current loss is: 295.37787544376033 +Current loss is: 302.71222054144835 +Current loss is: 312.29963874523617 +Current loss is: 310.0594473214611 +Current loss is: 302.37353839562803 +Current loss is: 296.45023439229044 +Current loss is: 295.15561561986084 +Current loss is: 297.2961876000978 +Current loss is: 299.91116127358544 +Current loss is: 300.614684811973 +Current loss is: 299.1410869807606 +Current loss is: 296.78057049421614 +Current loss is: 295.25791255616923 +Current loss is: 295.19178606538486 +Current loss is: 296.2725870709754 +Current loss is: 297.12846342631974 +Current loss is: 297.0123166653831 +Current loss is: 296.17366586047854 +Current loss is: 295.32546177836065 +Current loss is: 295.0267334278429 +Current loss is: 295.346872950253 +Current loss is: 295.73072278145463 +Current loss is: 295.7533136410972 +Current loss is: 295.5129969768086 +Current loss is: 295.23461197689846 +Current loss is: 295.05682093896877 ++ +
+ERROR: InterruptException: ++ + +
Now let us find a Monte-Carlo Solution and plot the both:
+ + ++monte_carlo_sol = [] +x_out = collect(98:1.00:110.00) +for x in x_out + u₀= [x] + g_val(du , u , p , t) = 0.4.*u + f_val(du , u , p , t) = du .= 0.04.*u + dt = 0.01 + tspan = (0.0,1.0) + prob = SDEProblem(f,g,u₀,tspan) + output_func(sol,i) = (sol[end],false) + ensembleprob_val = EnsembleProblem(prob , output_func = output_func ) + sim_val = solve(ensembleprob_val, EM(), EnsembleThreads() , dt=0.01, trajectories=100000,adaptive=false) + s = reduce(hcat , sim_val.u) + global monte_carlo_sol = push!(monte_carlo_sol , mean(phi(s))) +end ++ + +
+ERROR: UndefVarError: f not defined ++ + + +
+m = fmap(cpu , m) +x_model = reshape(x_out, 1 , size(x_out)[1]) +y_out = m(x_model) +y_out = reshape(y_out , 13 , 1) ++ + +
+13×1 Array{Float32,2}:
+ 11.862513
+ 11.9830265
+ 12.103537
+ 12.224047
+ 12.3445635
+ 12.46507
+ 12.585586
+ 12.706097
+ 12.826607
+ 12.947119
+ 13.067632
+ 13.188141
+ 13.308652
+
+
+
+Plotting the solutions
+ + ++plot(x_out , y_out , lw = 3 , xaxis="Initial Stock Price", yaxis="Payoff") +plot!(x_out , monte_carlo_sol , lw = 3 , xaxis="Initial Stock Price", yaxis="Payoff") ++ + +
+ERROR: BoundsError: attempt to access 0-element Array{Float64,1} at index [1:13]
+
+
+
+
+