The Black-Scholes partial differential equation is the foundation of much of modern finance. Suppose you have a European-style call option (right to buy) for a certain stock, to be exercised at a time $T$ in the future. Then, as a function of the stock price $s$, the value $v(t,s)$ of the option is governed by $$ v_t = -\frac{\sigma^2}{2}s^2v_{ss} -rsv_s + rv, \quad s>0, \quad t > 0, $$ subject to $$ v(0, T) = 0, \quad v(t, s)\rightarrow s \text{ as } s\rightarrow \infty. $$
It's much more convenient for us to work on a finite domain, so we will truncate and replace the asymptotic condition with a Neumann condition.
d = [0 500]; s = chebfun('s', d); sigma = 0.45; r = 0.03; A = chebop(@(s,v) -sigma^2/2*s.^2.*diff(v,2) - r*s.*diff(v) + r*v, d); A.lbc = 0; A.rbc = @(v) diff(v) - 1; % replaces v->s as s->infinity
Because the PDE is linear, we can solve it using an operator exponential (i.e., $C_0$ semigroup) via expm
. However, that function requires homogeneous boundary conditions, in the form $Bv=q$.
There is a very general trick to remove the nonhomoegeneous conditions, by finding a particular solution. Suppose $u(s)$ satisfies $Au=0$, $Bu=q$. Then $(v-u)_t = v_t=Av = A(v-u)$, and $B(u-v)=0$. Hence one can propagate $w=v-u$ with homogeneous conditions in combination with the solution of a boundary value problem for $u$.
u = A\0; % u satisfies u_t = A*u = 0, B*u = q. A.rbc = 0; % Homogenous BCs.
At the final time, the value of the option is the positive part of the difference between the strike price (for us, 50) and the stock price.
vT = max( 0, s-50 ); clf plot(vT, [40 60]), hold on, ylim([-.5 14]) wT = vT - u; % Adjusted variable - satisfies w_t = A*w, B*w = 0.
We can now evaluate at any time directly (not by marching) using the operator exponential implemented by expm
. Advance in time to $t = 0.1, 0.2, ..., 0.5$ using expm
:
for t = 0.1:0.1:0.5 w = expm(A, -t, wT); v = w + u; plot(v, [40 60], 'k') end
The value of this option if $s=55$, six months before the maturity date, is
v(55)
ans = 9.849887661936435
Note that the initial condition for $w$ is piecewise defined. Thus, w
will be piecewise at all times in Chebfun.
w
w = chebfun column (1 smooth piece) interval length endpoint values [ 0, 5e+02] 97 5e-11 7.9e-12 Epslevel = 4.510522e-13. Vscale = 4.925018e+01.
However, the solution is actually smooth.
wss = diff(w,2); jump2 = wss(50+100*eps) - wss(50-100*eps)
jump2 = -2.428612866367530e-17