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
-
F. Black and M. Scholes, "The pricing of options and corporate liabilities", Journal of Political Economy 81 (1973), 637--654.
-
R. Carmona and V. Durrleman, "Pricing and Hedging Spread Options", SIAM Review 45 (2003), 627--685.
-
M. Gass, K. Glau, M. Mahlstedt and M. Mair, "Chebyshev interpolation for parametric option pricing", working paper, 2016, http://arxiv.org/abs/1505.04648.