Move kolmogorov tutorial to advanced and add dependencies.

tutorials/advanced/03-kolmogorov_equations.jmd

@@ -0,0 +1,134 @@
+---
+title:  Kolmogorov Backward Equations
+author: Ashutosh Bharambe
+---
+
+```julia
+using Flux, StochasticDiffEq
+using NeuralNetDiffEq
+using Plots
+using CuArrays
+using CUDAnative
+```
+## Introduction on Backward Kolmogorov Equations
+
+The backward Kolmogorov Equation deals with a terminal condtion.
+The one dimensional backward kolmogorov equation that we are going to deal with is of the form :
+
+$$
+  \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} ,\hspace{0.5cm} p(T , x) = \varphi(x)
+$$
+for all $ t \in{ [0 , T] } $ and for all $ x \in R^d $
+
+#### The Black Scholes Model
+
+The Black-Scholes Model governs the price evolution of the European put or call option. In the below equation V is the price of some derivative , S is the Stock Price , r is the risk free interest
+rate and σ the volatility of the stock returns. The payoff at a time T is known to us. And this makes it a terminal PDE. In case of an European put option the PDE is:
+$$
+  \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  ,\hspace{0.5cm} V(T , S) =  max\{\mathcal{K} - S , 0 \}
+$$
+for all $ t \in{ [0 , T] } $ and for all $ S \in R^d $
+
+In order to make the above equation in the form of the Backward - Kolmogorov PDE we should substitute
+
+$$
+  V(S , t) = e^{r(t-T)}p(S , t)
+$$
+and thus we get
+$$
+  e^{r(t-T)}\frac{\partial p}{\partial t} + re^{r(t-T)}p(S , t)  = -\mu(x)\frac{\partial p}{\partial x}e^{r(t-T)} - \frac{1}{2}{\sigma^2}(x)\frac{\partial^2 p}{\partial x^2}e^{r(t-T)}
+  + re^{r(t-T)}p(S , t)
+$$
+And the terminal condition
+$$
+  p(S , T) = max\{ \mathcal{K} - x , 0 \}
+$$
+We will train our model and the model itself will be the solution of the equation
+## Defining the problem and the solver
+We should start defining the terminal condition for our equation:
+```julia
+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
+```
+Now we shall define the problem :
+We will define the σ and μ by comparing it to the orignal equation. The xspan is the span of initial stock prices.
+```julia
+d = 1
+r = 0.02
+sigma = 0.4
+xspan = (80.00 , 115.0)
+tspan = (0.0 , 1.0)
+σ(du , u , p , t) = du .= sigma.*u
+μ(du , u , p , t) = du .= r.*u
+prob = KolmogorovPDEProblem(μ , σ , phi , xspan , tspan, d)
+```
+Now once we have defined our problem it is necessary to define the parameters for the solver.
+```julia
+sdealg = EM()
+ensemblealg = EnsembleThreads()
+dt = 0.01
+dx = 0.001
+trajectories = 100000
+```
+
+Now lets define our model m and the optimiser
+```julia
+m = Chain(Dense(d, 8, leakyrelu),Dense(8, 16, leakyrelu),Dense(16 , 8 , leakyrelu) , Dense(8 , 1))
+use_gpu = false
+if CUDAnative.functional() == true
+  m = fmap(CuArrays.cu , m)
+  use_gpu = true
+end
+opt = Flux.ADAM(0.01)
+```
+And then finally call the solver
+```julia
+@time sol = solve(prob, NeuralNetDiffEq.NNKolmogorov(m, opt, sdealg, ensemblealg), verbose = true, dt = dt,
+            dx = dx , trajectories = trajectories , abstol=1e-6, maxiters = 4200 , use_gpu = use_gpu)
+```
+## Analyzing the solution
+Now let us find a Monte-Carlo Solution and plot the both:
+```julia
+monte_carlo_sol = []
+x_out = collect(85:5.00:110.00)
+for x in x_out
+  u₀= [x]
+  g_val(du , u , p , t) = du .= 0.4.*u
+  f_val(du , u , p , t) = du .= 0.02.*u
+  dt = 0.01
+  tspan = (0.0,1.0)
+  prob = SDEProblem(f_val,g_val,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)
+  mean_phi = sum(phi(s))/length(phi(s))
+  global monte_carlo_sol = push!(monte_carlo_sol , mean_phi)
+end
+
+```
+
+##Plotting the Solutions
+We should reshape the inputs and outputs to make it compatible with our model. This is the most important part. The algorithm gives a distributed function over all initial prices in the xspan.
+```julia
+x_model = reshape(x_out, 1 , size(x_out)[1])
+if use_gpu == true
+  m = fmap(cpu , m)
+end
+y_out = m(x_model)
+y_out = reshape(y_out , 6 , 1)
+```
+And now finally we can plot the solutions
+```julia
+plot(x_out , y_out , lw = 3 ,  xaxis="Initial Stock Price", yaxis="Payoff" , label = "NNKolmogorov")
+plot!(x_out , monte_carlo_sol , lw = 3 ,  xaxis="Initial Stock Price", yaxis="Payoff" ,label = "Monte Carlo Solutions")
+
+```

tutorials/advanced/Project.toml

@@ -5,27 +5,30 @@ CUDAnative = "be33ccc6-a3ff-5ff2-a52e-74243cff1e17"
 CuArrays = "3a865a2d-5b23-5a0f-bc46-62713ec82fae"
 DiffEqOperators = "9fdde737-9c7f-55bf-ade8-46b3f136cc48"
 DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
+Flux = "587475ba-b771-5e3f-ad9e-33799f191a9c"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
+NeuralNetDiffEq = "8faf48c0-8b73-11e9-0e63-2155955bfa4d"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SparseDiffTools = "47a9eef4-7e08-11e9-0b38-333d64bd3804"
 SparsityDetection = "684fba80-ace3-11e9-3d08-3bc7ed6f96df"
+StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0"
 Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4"
 
 [compat]
 AlgebraicMultigrid = "0.3"
 BenchmarkTools = "0.5"
+CUDAnative = "3.1"
+CuArrays = "2.2"
 DiffEqOperators = "4.10"
 DifferentialEquations = "6.14"
-Sundials = "4.2"
-SparsityDetection = "0.3"
-OrdinaryDiffEq = "5.41"
 ModelingToolkit = "3.10"
-CuArrays = "2.2"
-CUDAnative = "3.1"
+NLsolve = "4.4"
+OrdinaryDiffEq = "5.41"
 Plots = "1.4"
 SparseDiffTools = "1.9"
-NLsolve = "4.4"
+SparsityDetection = "0.3"
+Sundials = "4.2"