Option pricing with the Black-Scholes model. Code not running

Question

(*initialisation of variables*)
n = 10;(*expected value greater than 3*)
r = 0.03;
σ = 0.2;
sMax = 2;
sMin = -2;

(*define Δx.*)
Δx = (sMax - sMin)/(n - 1);

(*adding vector matrix... any assignment of v value should be 
included here*)
vecMatrix = Array[v, n];

(*reusable internal variables*)
p = 0.5 (σ^2)/(Δx^2);
q = (2 r - (σ^2))/Δx;

x = (2 p + r) - 0.75 q;
y = 5 p + q;
a = p - 0.25 q;
b = -2 p - r;
c = p + 0.25 q;

(*introducing matC for one time initialisation*)
matC = SparseArray[{Band[{1, 1}, {n, n}, n - 1] -> x, 
    Band[{1, 2}, {n, n}, {n - 1, n - 3}] -> y, 
    Band[{2, 1}, {n - 1, n}] -> a, Band[{2, 2}, {n - 1, n}] -> b, 
    Band[{2, 3}, {n - 1, n}] -> c}, {n, n}];
 MatrixForm [matC]*
 MatrixForm[vecMatrix]

Answer 1

I've done a lot of option pricing with MMA, but I can't relate your code to plain old Black-Scholes:

d1[spot_,strike_,ir_,div_,vol_, T_] = (Log[spot/strike] + (ir-div+vol^2/2) T) / (vol T^0.5]);  
d2[spot_,strike_,ir_,div_,vol_, T_] = (Log[spot/strike] + (ir-div-vol^2/2) T) / (vol T^0.5]);
N[z_] = (1 + Erf[z/Sqrt[2]])/2 ; (*cumulative normal density*)

bsCall[spot_,strike_,ir_,div_,vol_,T_] :=
   spot*E^(-div*T)*N[d1[spot,strike,ir,div,vol, T]] - 
   strike*E^(-r*T)*N[d2spot,strike,ir,div,vol,T]]

where

ir = interest rate  
div=dividends paid  
vol=volatility of the stock as measured by the standard deviation of its price
T = time to expiry  
Log[] = the natural logarithm

(I added this simple stuff to help others understand your question.) I don't see where the Black-Scholes part comes into your code? It would help educate others if you put your question in a context, for example, black-scholes option pricing is easy if you have one option, but if you have 100,000 at varying strikes, etc. But your question seems to really be about matrix algebra and sparce arrays.

Answer 2

You have a syntax error in the last line. I suspect you mean:

matC.vecMatrix // MatrixForm

which yields:

0.19875 v[1]+0.55125 v[2]
0.09 v[1]-0.2325 v[2]+0.1125 v[3]
0.09 v[2]-0.2325 v[3]+0.1125 v[4]
0.09 v[3]-0.2325 v[4]+0.1125 v[5]
0.09 v[4]-0.2325 v[5]+0.1125 v[6]
0.09 v[5]-0.2325 v[6]+0.1125 v[7]
0.09 v[6]-0.2325 v[7]+0.1125 v[8]
0.09 v[7]-0.2325 v[8]+0.1125 v[9]
0.09 v[8]-0.2325 v[9]+0.1125 v[10]
0.55125 v[9]+0.19875 v[10]

As an aside, comments that just say (* define \[Delta] x *) aren't that helpful to the next person looking at the code. I would recommend using commenting to explain what the variables are for.

As a further comment on your code, all your variables are global to your session, so the comment "reusable internal variables" is misleading. If you want to scope variables as internal to a particular calculation, you need to use scoping constructs like Module, Block and With. You might find this question helpful.

Answer 3

Write down the Black-Scholes equation in the following way

TraditionalForm[BS = {-(r*c[t, s]) + r*s*D[c[t, s], {s}] + (1/2)*s^2*\[Sigma]^2*D[c[t, s], 
{s, 2}] + D[c[t, s], {t}] == 0, c[T, s] == Max[s - k, 0]}]

which will spit out the Latex-like form that, in turn, can be solved analytically as follows

DSolve[BS, c[t, s], {t, s}][[1]]
(dsol = c[t, s] /. DSolve[BS, c[t, s], {t, s}][[1]]) // TraditionalForm

which gives you the analytical solution of the Black-Scholes model

$\frac{1}{2} e^{-r T} \left(\frac{s e^{r T} \left(\text{erfc} \left(2 k \log +\left(2 r+\sigma ^2\right) (t-T)-2 (s \log )\right)\right)}{2 \sqrt{2} \sigma \sqrt{T-t}}-\frac{k e^{r t} \left(\text{erfc} \left(2 k \log +\left(2 r-\sigma ^2\right) (t-T)-2 (s \log )\right)\right)}{2 \sqrt{2} \sigma \sqrt{T-t}}\right)$

Mind that dsol is our analytical solution here.

For our conditions we will use the obtained above solution to get the price of the call option, for example

NumberForm[dsol /. {t -> 0, s -> 100, k -> 100., \[Sigma] -> 0.2, T -> 1,r -> 0.05}, 16]

resulting

10.45058357218558

To make the plot for the above conditions

t = 0; k = 100; \[Sigma] = 0.2; T = 1; r = 0.05;
Plot3D[1/
2 E^(-r T) (-E^(r t)
   k Erfc[((t - T) (2 r - [Sigma]^2) + 2 Log[k] - 2 Log[s])/(
  2 Sqrt[2] Sqrt[-t + T] [Sigma])] + 
E^(r T) s Erfc[((t - T) (2 r + [Sigma]^2) + 2 Log[k] - 
   2 Log[s])/(2 Sqrt[2] Sqrt[-t + T] [Sigma])]),
    {t, 0.01, 1}, {s, 0, 100}, PlotRange -> All, 
AxesLabel -> {"T-t", "S", "C"}, LabelStyle -> Larger, 
ImageSize -> 600, ColorFunction -> Function[{x, y, z}, Hue[y*28/29]]]

If you need numerical solutions and explanations, look at https://lucynowacki.github.io/blog/method-of-lines-for-black-scholes-explicit/index.html or if you need code in Mathematica, touch me.

You can also solve the BS numerically. The simplest is the Explicit Finite Difference method as follows

K=100.; S0=100.; r=0.05; sig=0.2; T=12/12.; M=300; Nt=4000; 
Smax=220.;
n=Range[1,M-1];
dt=N[T/Nt];
ds=Smax/M//N;
U=ConstantArray[0.,{M+1,Nt+1}];
Sgrid = Range[0, Smax, ds];
Tgrid=Range[T,0,-dt];
(*boundary conditions*)
U[[1,All]]=0.;
U[[-1,All]]=Smax- K E^(-r*Tgrid);
U[[All,-1]]= Max[#-K,0]& /@ Sgrid; (*terminal condition*)
(*coefficients for finite difference matrix*)
\[Alpha] = 0.5*dt*(sig^2*n^2-r*n);
\[Beta] = 1-dt*(sig^2*n^2+r);
\[Gamma] = 0.5*dt*(sig^2*n^2+r*n);
(*set up the finite difference matrix*)
L = Table[KroneckerDelta[n,m](1-dt* 
(sig^2*n^2+r))+KroneckerDelta[n+1,m](0.5*dt* 
(sig^2*n^2+r*n))+KroneckerDelta[n-1,m](0.5*dt*(sig^2*n^2-r*n)), 
{n,1,M-1}, {m,1,M-1}];
Max[Abs[Eigenvalues[L]]] (*must be less than 1 for stability*)
offsetConstants = {\[Alpha][[1]], \[Gamma][[-1]]}; (* offsets to 
accomodate boundary conditions*)
(*to populate finite difference matrix with numerical solution*)
Do[ U[[2;;-2,it]] = L . U[[2;;-2, it+1]];
U[[{2,-2},it]] = U[[{2,-2},it]]+ 
offsetConstants*U[[{1,-1},it+1]],
{it,Nt,1,-1}]
Print["The call price is: 
",Interpolation[Thread[{Sgrid,U[[All,1]]}],S0]]

To plot the solution matrix

ListPlot3D[Transpose[U],ColorFunction->Hue, PlotRange-> 
{{0,100},All,{0,15}},DataRange->{{0,Smax},{0,1}},AspectRatio->1, 
LabelStyle -> Larger]

Note that here the solution matrix is transposed; however, it is irrelevant.