# Black Scholes PDE solved by method of lines

I've been trying to solve the Black Scholes partial differential equation (BSPDE) by transforming it to the heat equation (HE) form and using method of lines (discretizing spatial argument and leaving time continuous) to change the problem to a system of ODEs. I then solved those ODEs using NDSolve and now I'm left with a bunch of functions which, put together, make the solution of the transformed BSPDE.

The problem is that I need to transform the solution back to the original arguments and I have no idea wheter that's possible or how to do it.

Any help would be greatly appreciated.

PS: I am aware that that is not how the solution to the BSPDE should look. I'm working on that too ;)

The code so far:

S0 =  10;
T = 1;
sigma = 0.3;
r = 0.05;
K = 10;
n = 100;

Subscript[h, n] = 10/n ;

U[t_] = Table[Subscript[u, i][t], {i, 0, n}] ;

equations =
Join[{0},
ListCorrelate[{1, -2, 1}/Subscript[h, n]^2,
U[t], {1, 2}, {Subscript[u, n - 1][t]}]]];

initialcondition =
Table[Max[0,
Exp[i Subscript[h, n] (r/sigma^2 + 0.5)] *
Exp[i Subscript[h, n] (r/sigma^2 - 0.5)]], {i, 0, n}]];

lines = NDSolve[{equations, initialcondition}, U[t], {t, 0, 4}]

ParametricPlot3D[
Evaluate[
Table[{i Subscript[h, n], t,
First[Subscript[u, i][t] /. lines]}, {i, 0, n}]], {t, 0, 4},
PlotRange -> All, AxesLabel -> {"x", "t", "u"}]

• Welcome to Mathematica.SE! I suggest the following: 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Read the faq! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign!
– user9660
May 2, 2015 at 18:40
• NDSolve can use the method of lines by default. May 2, 2015 at 20:09
• @2012rcampion thank you for your answer. yes, I was aware that the article you posted exists, but I couldn't find the answer to my problem there. I can get the solution, but I don't know how to return all the transformations I made to get the solution with the original arguments. I'm sorry if I didn't make the problem clear enough in the first post. The article is helpful, but does not solve my issue. May 2, 2015 at 20:39
• If you can post your code so far we might be able to help you... May 2, 2015 at 20:40
• what I have so far actually follows the solution presented in the article you linked in your first comment, but I will edit it in if it helps May 2, 2015 at 21:05

First, the analytical solution

truVal =
FinancialDerivative[{"European", "Call"}, {"StrikePrice" -> 100.00,
"Expiration" -> 1},
{"InterestRate" -> 0.05, "Volatility" -> 0.2,
"CurrentPrice" -> 100}]


gives

Out[6]= 10.4502


To apply MOL to the Black-Scholes model, first, we have to build a space-discredited grid. Time is continuous in the MOL (but its discretisation on the grid is still necessary); however, time is not considered discrete per se in MOL. As an integrator for the explicit method, we can use Runge-Kutta-3 with variance reduction or the Taylor method. Naturally, we have applied initial and boundary conditions. For a discrete grid, we can use the NDSolveFiniteDifferenceDerivative Mathematica method. The function below for explicit Black Scholes is based on the RK-TVD-3 integrator below. If you need an implicit way, use the Trapezoidal or Hermite scheme. Function eats consecutively: strike KK_, min grid value SSmin_, max grid value SSmax_, interest rate rr_, time of contract TT_, current price of the stock SS0_, space(price) discretisation NNs_, time discretisation NNt_, difference order for grid orda_ as follows:

BSEx[100.,10.,230.,0.2,0.05, 12/12.,100.,500,5*10^4,2][[4]]


Function spits out consecutively: matrix of prices U, error, Courant number(max eigenvector), price of an option, 16-digit price of the option, price discretisation interval ds, time discretisation interval dt. Thus to get the price of the option use [[4]] next to the function definition as above.

There is a full implementation below:

BSEx[KK_, SSmin_, SSmax_, ssigma_, rr_, TT_, SS0_, NNs_, NNt_, orda_]:= Module[{K=KK,
Smin=SSmin, Smax=SSmax, sig=ssigma, r=rr, T=TT, S0=SS0, Ns=NNs, Nt=NNt, it, m},
exact=FinancialDerivative[{"European","Call"},{"StrikePrice"->K,"Expiration"->T},
{"InterestRate"->r,"Volatility"->sig,"CurrentPrice"->S0}]; (*needed to compute the error*)
Subscript[Dn, nn_][Xgrid_]:=NDSolveFiniteDifferenceDerivative[Derivative[nn], Xgrid,
"DifferenceOrder"-> orda]@ "DifferentiationMatrix";
dt=N[T/Nt]; (*time step_size*)
ds=(Smax-Smin)/(Ns-1); (*price step_size*)
U=ConstantArray[0.,{Ns,Nt}]; (*Grid of prices which will be populated during solving*)
Sgrid = Range[Smin, Smax, ds]; (*span of prices on which equation will be solved*)
Len[Sgrid];
tgrid=T-Range[0,T,dt]; (*span of times*)
U[[All,1]]= Max[#-K,0]& /@ Sgrid; (*Initial condition*)
upper = Table[0., {m, 0, Nt+1}]; (*Left Boundary condition*)
lower = Table[Smax- K E^(-r*m*dt), {m,0,Nt+1}]; (*Right Boundary condition*)
II = IdentityMatrix[Ns];
D1=Subscript[Dn, 1][Sgrid]; (*finite difference matrix for the first derivative*)
D2=Subscript[Dn, 2][Sgrid]; (*finite difference matrix for the second derivative*)
L = 0.5*sig^2*Sgrid^2*D2+r*Sgrid*D1-r*II; (*linear operator*)

M=DeveloperToPackedArray[M]; (*faster*)
Courant=Max[Abs[Eigenvalues[M]]]; (*must be less than 1 for stability*)

Do[
U1 =  U[[All,it]]+dt*L . U[[All,it]]; (*RK-TVD-3*)
U2= 3/4 U[[All,it]]+1/4 U1+dt*1/4 L . U1;
U[[All,it+1]]= 1/3 U[[All,it]]+2/3 U2+dt*2/3*L . U2,
{it,1,Nt-1}];
price= Interpolation[Thread[{Sgrid,U[[All,-1]]}],S0, InterpolationOrder->3];(* price must be interpolated*)
{U,Abs[exact-price], Courant, price, NumberForm[price,16], ds, dt}] (*results*)


To execute

BSEx[100.,10.,230.,0.2,0.05, 12/12.,100.,500,5*10^4,2][[4]]


gives

Out[51]= 10.4502


For better estimation increase grid, but it will need more time and memory.

To make a sexy plot:

U=BSEx[100.,10.,200.,0.2,0.05, 12/12.,100.,20,2*10^1,2][[1]];
ListPointPlot3D[Transpose[U],ColorFunction->Function[{x,y,z},Hue[x*58/29]],
PlotRange->{{10,100},All,{0,10.5}},DataRange->{{10,200},{0,1}}, ImageSize->700,
AxesLabel->{"x", "t", "U"}, PlotStyle-> PointSize[0.005],BoxRatios->{3,3,1.5},Background->Black,
PlotLabel->Style["MOL for BS", Pink,15],BoxStyle->Directive[Dashed],LabelStyle->{Pink,18} ]
`

If you would like know how to implement implicit method look at https://lucynowacki.github.io/blog/method-of-lines-for-black-scholes-implicit/index.html or contact me.