Problem formulation and model

We model the price processes of two stocks in the two-dimensional Black-Scholes model [1] as solutions of the following stochastic differential equations:

$$ dS_1/S_1 = (r - 0.5 \sigma_1^2) dt + \sigma_1 dW^1_t $$

$$ dS_2/S_2 = (r - 0.5 \sigma_2^2) dt + \sigma_2 dW^2_t, $$

where $W^1_t$ and $W^2_t$ denote Brownian motions/Wiener processes with correlation $\rho$.

We are interested in pricing spread options. These are important instruments on several markets such as on equity, fixed income, foreign exchange, commodities, or energy markets, see [2]. In particular, when prices of spread options are used to fit model parameters, fast and accurate pricing methods are called for.

Payoff of a spread option with strike price $K$ and maturity $T$ is

$$ \mbox{payoff} = \max \{ S_1(T) - S_2(T) - K, 0 \}. $$

Using risk-neutral valuation, the price of the spread option is given by the following expectation:

$$ spread_{price} = \exp(-r T) E[\max \{ S_1(T) - S_2(T) - K, 0\}]. $$

In order to compute the price of the spread option in the Black-Scholes model, the expected value, a two-dimensional integral, needs to be evaluated. This is considerably more time-demanding than pricing one-dimensional plain vanilla options with the famous and explicit Black-Scholes formula. For model calibration (parameter fitting) one determines the parameters $\sigma_1$ and $\sigma_2$ from plain vanilla options that only involve either stock 1 or stock 2. The spread option prices then can be used to determine the parameter $\rho$.

During the model calibration, for fixed $\sigma_1$ and $\sigma_2$, interest rate $r$ and initial stock prices $S_1(0)$ and $S_2(0)$, the spread option prices need to be evaluated for a set of maturities $T$ and strikes $K$, and different parameters $\rho$.

Our approach is to interpret $spread_{price}$ as a function of the three model parameters $T$, $K$ and $\rho$ and use 3D Chebyshev interpolation of these parameters [3].

We first specify the fixed parameters

S01 = 1; % initial value of stock 1
S02 = 1; % initial value of stock 2
r = 0;   % interest rate
x01 = log(S01);
x02 = log(S02);
sig1 = 0.3; % model parameter sigma_1
sig2 = 0.3; % model parameter sigma_2
%
% Parameter domain for T x K x rho:
%
T_min = 0.3;
T_max = 2;
K_min = 0.3;
K_max = 2;
rho_min = -0.9;
rho_max = 0.9;
dom = [T_min, T_max, K_min, K_max, rho_min, rho_max];
%
% Definition of problem dependent functions:
%
payoff = @(K, x1, x2) max(exp(x01+x1) - exp(x02+x2)-K, 0);
normdensity_2d = @(sig1, sig2, rho, x1, x2, mu1, mu2) (1./(sqrt((sig1.^2) .* ...
    (sig2.^2).*(1-rho.^2)).*2.*pi)).*exp(-0.5*(1./((sig1.^2).*(sig2.^2).*...
    (1-rho.^2))).*(((sig2.^2).*(x1-mu1).^2)-2.*rho.*sig1.*sig2 .* ...
    (x1-mu1).*(x2-mu2)+(sig1.^2).*(x2-mu2).^2));
integrand = @(T, K, rho, x1, x2) payoff(K,x1, x2) .* ...
    normdensity_2d(sig1.*sqrt(T), sig2.*sqrt(T), rho, x1, x2, (r-0.5.*sig1.^2).*T, (r-0.5.*sig2.^2).*T);
%
% Derivation of two-dimensional integral/expectation via 2D Clenshaw-Curtis quadrature:
% (We choose a grid of 150x150 quadrature points. For higher accuracy,
% more points need to be chosen.)
%
ff1 = @(T, K, rho, x1, x2) integrand(T, K, rho, x1, x2);
[xx, wx] = chebpts(150, [-10, 10]);
[yy, wy] = chebpts(150, [-10, 10]);
[xx2, yy2] = meshgrid(xx, yy);
pricetmp2d = @(T, K, rho) wy * feval(ff1, T, K, rho, xx2, yy2) * wx';

price = @(T, K, rho) exp(-r.*T) .* pricetmp2d(T, K, rho);

price is a 3D function that returns the spread option price based on the three parameters $T$, $K$ and $\rho$.

Calling the Chebfun3 constructor

Here, we set the tolerance to $10^{-5}$ (This is done for speed. With the default tolerance, this example would take a factor of 10 times longer. See Subsection 18.10 of the Chebfun Guide.) The run-time of the constructor scales linearly with the run-time of the bivariate integration. The longer this step takes, the higher the expected efficiency gain of the interpolation.

tic
chebPrice = chebfun3(price, 'eps', 1e-5, 'vectorize', dom)
toc
chebPrice =
   chebfun3 object 
   cols: Inf x 5  chebfun
   rows: Inf x 10 chebfun
  tubes: Inf x 7  chebfun
   core: 5 x 10 x 7
 length: 10, 42, 13
 domain: [0.3, 2] x [0.3, 2] x [-0.9, 0.9]
 vertical scale = 0.19
Elapsed time is 38.815439 seconds.

Let's plot the function at three slices $T = 2$, $K = 0.3$, and $\rho = -0.9$:

slice(chebPrice, 2, 0.3, -0.9)
campos([-10 10 10])
xlabel('T'), ylabel('K'), zlabel('\rho')

Checking the error

Finally, we evaluate both the function handle and the corresponding chebfun3 object. We observe a considerable gain in efficiency: Evaluation of the chebfun3 object is hundreds times faster than evaluation of the function handle.

M = 15;
T = linspace(dom(1), dom(2), M);
K = linspace(dom(3), dom(4), M);
rho = linspace(dom(5), dom(6), M);
V = zeros(M, M, M);
tic
for i =1:M
    for j=1:M
        for k=1:M
            V(i,j,k)=price(T(i), K(j), rho(k));
        end
    end
end
format short
time_price = toc

[xx,yy,zz] = ndgrid(T, K, rho);
tic
VCheb = chebPrice(xx,yy,zz);
time_ChebPrice = toc

err = V - VCheb;
err = max(abs(err(:)))
time_price =
    2.0692
time_ChebPrice =
    0.0076
err =
   2.5650e-04

References

  1. F. Black and M. Scholes, "The pricing of options and corporate liabilities", Journal of Political Economy 81 (1973), 637--654.

  2. R. Carmona and V. Durrleman, "Pricing and Hedging Spread Options", SIAM Review 45 (2003), 627--685.

  3. M. Gass, K. Glau, M. Mahlstedt and M. Mair, "Chebyshev interpolation for parametric option pricing", working paper, 2016, http://arxiv.org/abs/1505.04648.