%% Pricing European and American Spread Options
% This example shows how to price and calculate sensitivities for European and American spread options using various techniques. First, the price and sensitivities for a European spread option is calculated using closed form solutions. Then, price and sensitivities for an American spread option is calculated using finite difference and Monte Carlo simulations. Finally, further analysis is conducted on spread options with a different range of inputs.
% Spread options are options on the difference of two underlying asset prices. For example, a call option on the spread between two assets has the following payoff at maturity:
% max(X_1 - X_2 - K, 0)
% where X_1 is the price of the first underlying asset, X_2 is the price of the second underlying asset, and K is the strike price. At maturity, if the spread X_1 - X_2 is greater than the strike price K, the option holder exercises the option and gains the difference between the spread and the strike price. If the spread is less than 0, the option holder does not exercise the option, and the payoff is 0. Spread options are frequently traded in the energy market. Two examples are:
% Crack spreads: Options on the spread between refined petroleum products and crude oil. The spread represents the refinement margin made by the oil refinery by "cracking" the crude oil into a refined petroleum product.
% Spark spreads: Options on the spread between electricity and some type of fuel. The spread represents the margin of the power plant, which takes fuel to run its generator to produce electricity.
%% Overview of the Pricing Methods
% There are several methods to price spread options, as discussed in [1]. This example uses the closed form, finite difference, and Monte Carlo simulations to price spread options. The advantages and disadvantages of each method are discussed below:
% Closed form solutions and approximations of partial differential equations (PDE) are advantageous because they are very fast, and extend well to computing sensitivities (Greeks). However, closed form solutions are not always available, for example for American spread options.
% The finite difference method is a numerical procedure to solve PDEs by discretizing the price and time variables into a grid. A detailed analysis of this method can be found in [2]. It can handle cases where closed form solutions are not available. Also, finite difference extends well to calculating sensitivities because it outputs a grid of option prices for a range of underlying prices and times. However, it is slower than the closed form solutions.
% Monte Carlo simulation uses random sampling to simulate movements of the underlying asset prices. It handles cases where closed solutions do not exist. However, it usually takes a long time to run, especially if sensitivities are calculated.
%% Pricing a European Spread Option
% The following example demonstrates the pricing of a crack spread option.
% A refiner is concerned about its upcoming maintenance schedule and needs to protect against decreasing crude oil prices and increasing heating oil prices. During the maintenance the refiner needs to continue providing customers with heating oil to meet their demands. The refiner's strategy is to use spread options to manage its hedge.
% On January 2013, the refiner buys a 1:1 crack spread option by purchasing heating oil futures and selling crude oil futures. CLF14 WTI crude oil futures is at $100 per barrel and HOF14 heating oil futures contract is at $2.6190 per gallon.

clear;

% Price, volatility, and dividend of heating oil
Price1gallon = 2.6190;       % $/gallon 
Price1 = Price1gallon*42;    % $/barrel
Vol1 = 0.10;
Div1 = 0.03;

% Price, volatility, and dividend of WTI crude oil
Price2 = 100;     % $/barrel
Vol2 = 0.15;
Div2 = 0.02;

% Correlation of underlying prices
Corr = 0.3;

% Option type
OptSpec = 'call';

% Strike
Strike = 5;

% Settlement date
Settle = '01-Jan-2013';

% Maturity
Maturity = '01-Jan-2014'; 

% Risk free rate
RiskFreeRate = 0.05;

%
% The pricing functions take an interest-rate term structure and stock structure as inputs. Also, you need to specify which outputs are of interest.

% Define RateSpec
Compounding = -1;
Basis = 1;
RateSpec = intenvset('ValuationDate', Settle, 'StartDates', Settle, ...
    'EndDates', Maturity, 'Rates', RiskFreeRate, 'Compounding', ...
    Compounding, 'Basis', Basis);

% Define StockSpec for the two assets
StockSpec1 = stockspec(Vol1, Price1, 'Continuous', Div1);
StockSpec2 = stockspec(Vol2, Price2, 'Continuous', Div2);

% Specify price and sensitivity outputs
OutSpec = {'Price', 'Delta', 'Gamma'};

%
% The Financial Instruments Toolbox™ contains two types of closed form approximations for calculating price and sensitivities of European spread options: the Kirk's approximation (spreadbykirk, spreadsensbykirk) and the Bjerksund and Stensland model (spreadbybjs, spreadsensbybjs) [3].
% The function spreadsensbykirk calculates prices and sensitivities for a European spread option using the Kirk's approximation.

% Kirk's approximation
[PriceKirk, DeltaKirk, GammaKirk] = ...
    spreadsensbykirk(RateSpec, StockSpec1, StockSpec2, Settle, ...
    Maturity, OptSpec, Strike, Corr, 'OutSpec', OutSpec)

%
% The function spreadsensbybjs calculates the prices and sensitivities for a European spread option using the Bjerksund and Stensland model. In [3], Bjerksund and Stensland explains that the Kirk's approximation tends to underprice the spread option when the strike is close to zero, and overprice when the strike is further away from zero. In comparison, the model by Bjerksund and Stensland has higher precision.

% Bjerksund and Stensland model
[PriceBJS, DeltaBJS, GammaBJS] = ...
    spreadsensbybjs(RateSpec, StockSpec1, StockSpec2, Settle, ...
    Maturity, OptSpec, Strike, Corr, 'OutSpec', OutSpec)

%
% A comparison of the calculated prices show that the two closed form models produce similar results for price and sensitivities. In addition to delta and gamma, the functions can also calculate theta, vega, lambda, and rho.

displayComparison('Kirk', 'BJS', PriceKirk, PriceBJS, DeltaKirk, DeltaBJS, GammaKirk, GammaBJS)

%
%% Pricing an American Spread Option
% Although the closed form approximations are fast and well suited for pricing European spread options, they cannot price American spread options. Using the finite difference method and the Monte Carlo method, an American spread option can be priced. In this example, an American spread option is priced with the same attributes as the above crack spread option.
% The finite difference method numerically solves a PDE by discretizing the underlying price and time variables into a grid. The Financial Instrument Toolbox™ contains the functions spreadbyfd and spreadsensbyfd, which calculate prices and sensitivities for European and American spread options using the finite difference method. For the finite difference method, the composition of the grid has a large impact on the quality of the output and the execution time. Generally, a finely discretized grid will result in outputs that are closer to the theoretical value, but it comes at the cost of longer execution times. The composition of the grid is controlled using optional parameters PriceGridSize, TimeGridSize, AssetPriceMin and AssetPriceMax.
% To indicate pricing an American option, add an optional input of AmericanOpt with a value of 1 to the argument of the function.

% Finite difference method for American spread option
[PriceFD, DeltaFD, GammaFD, PriceGrid, AssetPrice1, ...
    AssetPrice2] = ...
    spreadsensbyfd(RateSpec, StockSpec1, StockSpec2, Settle, ...
    Maturity, OptSpec, Strike, Corr, 'OutSpec', OutSpec, ...    
    'PriceGridSize', [500 500], 'TimeGridSize', 100, ...
    'AssetPriceMin', [0 0], 'AssetPriceMax', [2000 2000], ...
    'AmericanOpt', 1);

% Display price and sensitivities
PriceFD
DeltaFD
GammaFD

%
% The function spreadsensbyfd also returns a grid that contains the option prices for a range of underlying prices and times. The grid of option prices at time zero, which is the option prices at the settle date, can be plotted for a range of underlying prices.

% Plot option prices
figure;
mesh(AssetPrice1, AssetPrice2, PriceGrid(:, :, 1)');
title('American Spread Option Prices for Range of Underlying Prices');
xlabel('Price of underlying asset 1');
ylabel('Price of underlying asset 2');
zlabel('Price of spread option');

%
% An American style option can be priced by Monte Carlo methods using the least square method of Longstaff and Schwartz [4]. The Financial Instruments Toolbox™ contains the functions spreadbyls and spreadsensbyls, that calculate prices and sensitivities of European and American options using simulations. The Monte Carlo simulation method in spreadsensbyls generates multiple paths of simulations according to a geometric Brownian motion (GBM) for the two underlying asset prices. Similar to the finite difference method where the granularity of the grid determined the quality of the output and the execution time, the quality of output and execution time of the Monte Carlo simulation depends on the number of paths (NumTrials) and the number of time periods per path (NumPeriods). Also, the results obtained by Monte Carlo simulations are not deterministic. Each run will have different results depending on the simulation outcomes.

% To indicate that we are pricing an American option using the Longstaff
% and Schwartz method, add an optional input of |AmericanOpt| with a value
% of |1| to the argument of the function.

% Monte Carlo method for American spread option
[PriceMC, DeltaMC, GammaMC] = ...
    spreadsensbyls(RateSpec, StockSpec1, StockSpec2, Settle, ...
    Maturity, OptSpec, Strike, Corr, 'OutSpec', OutSpec, ...
    'NumTrials', 1000, 'Antithetic', true, 'AmericanOpt', 1)

%
% The results of the two models are compared. The prices and sensitivities calculated by the Longstaff and Schwartz method will vary at each run, depending on the outcome of the simulations. It is important to note that the quality of the results from the finite difference method and the Monte Carlo simulation depend on the optional input parameters. For example, increasing the number of paths (NumTrials) for the spreadsensbyls function will result in more precise results at the cost of longer execution times.

displayComparison('Finite Difference', 'Monte Carlo', PriceFD, PriceMC, DeltaFD, DeltaMC, GammaFD, GammaMC)

%
%% Comparing Results for a Range of Strike Prices
% As discussed earlier, the Kirk's approximation tends to overprice spread options when the strike is further away from zero. To confirm this, a spread option is priced with the same attributes as before, but for a range of strike prices.

% Specify outputs
OutSpec = {'Price', 'Delta'};

% Range of strike prices
Strike = [-25; -15; -5; 0; 5; 15; 25];

%
% The results from the Kirk's approximation and the Bjerksund and Stensland model are compared against the numerical approximation from the finite difference method. Since spreadsensbyfd can only price one option at a time, it is called in a loop for each strike value. The Monte Carlo simulation (spreadsensbyls) with a large number of trial paths can also be used as a benchmark, but the finite difference is used for this example.

% Kirk's approximation
[PriceKirk, DeltaKirk] = ...
    spreadsensbykirk(RateSpec, StockSpec1, StockSpec2, Settle, ...
    Maturity, OptSpec, Strike, Corr, 'OutSpec', OutSpec);

% Bjerksund and Stensland model
[PriceBJS, DeltaBJS] = ...
    spreadsensbybjs(RateSpec, StockSpec1, StockSpec2, Settle, ...
    Maturity, OptSpec, Strike, Corr, 'OutSpec', OutSpec);

% Finite difference
PriceFD = zeros(numel(Strike), 1);
DeltaFD = zeros(numel(Strike), 2);
for i = 1:numel(Strike)
    [PriceFD(i), DeltaFD(i,:)] = ...
    spreadsensbyfd(RateSpec, StockSpec1, StockSpec2, Settle, ...
    Maturity, OptSpec, Strike(i), Corr, 'OutSpec', OutSpec, ...    
    'PriceGridSize', [500 500], 'TimeGridSize', 100, ...
    'AssetPriceMin', [0 0], 'AssetPriceMax', [2000 2000]);
end

displayComparisonPrices(PriceKirk, PriceBJS, PriceFD, Strike)

%
% The difference in prices between the closed form and finite difference method is plotted below. It is clear that as the strike moves further away from 0, the difference between the Kirk's approximation and finite difference (red line) increases, while the difference between the Bjerksund and Stensland model and finite difference (blue line) stays at the same level. As stated in [3], the Kirk's approximation is overpricing the spread option when the strike is far away from 0.

% Plot of difference in price against the benchmark
figure;
plot(PriceKirk-PriceFD, 'Color', 'red');
hold on;
plot(PriceBJS-PriceFD, 'Color', 'blue');
hold off;
title('Difference in Price Against Finite Difference');
legend('Kirk', 'BJS', 'Location', 'EastOutside');
xlabel('Strike');
ax = gca;
ax.XTickLabel = Strike;
ylabel('Difference in Price');

%
% Next, the difference in delta between the closed form models and finite difference is plotted. The top plot shows the difference in delta for the first asset, and the bottom plot shows the difference in delta for the second asset. As seen from the small increments in the y-axis of the order 10e-3, it can be seen that all three models (Kirk, BJS, finite difference) produce similar values for delta.

% Plot of difference in delta of first asset against the benchmark
figure;
subplot(2, 1, 1);
plot(DeltaKirk(:,1)-DeltaFD(:,1), 'Color', 'red');
hold on;
plot(DeltaBJS(:,1)-DeltaFD(:,1), 'Color', 'blue');
hold off;
title('Difference in Delta (Asset 1) Against FD');
legend('Kirk', 'BJS', 'Location', 'EastOutside');
xlabel('Strike');
ax = gca;
ax.XTickLabel = Strike;
ylabel('Difference in Delta');

% Plot of difference in delta of second asset against the benchmark
subplot(2, 1, 2);
plot(DeltaKirk(:,2)-DeltaFD(:,2), 'Color', 'red');
hold on;
plot(DeltaBJS(:,2)-DeltaFD(:,2), 'Color', 'blue');
hold off;
title('Difference in Delta (Asset 2) Against FD');
legend('Kirk', 'BJS', 'Location', 'EastOutside');
xlabel('Strike');
ax = gca;
ax.XTickLabel = Strike;
ylabel('Difference in Delta');

%
%% Analyzing Prices and Vega at Different Levels of Volatility
% To further show the type of analysis that can be conducted using these models, the above spread option is priced at different levels of volatility for the first asset. The price and vega are compared at three levels of volatility for the first asset: 0.1, 0.3, and 0.5. The Bjerksund and Stensland model is used for this analysis.

% Strike
Strike = 5;

% Specify output
OutSpec = {'Price', 'Vega'};

% Different levels of volatility for asset 1
Vol1 = [0.1, 0.3, 0.5];

StockSpec1 = stockspec(Vol1, Price1, 'Continuous', Div1);

% Bjerksund and Stensland model
[PriceBJS, VegaBJS] = ...
    spreadsensbybjs(RateSpec, StockSpec1, StockSpec2, Settle, ...
    Maturity, OptSpec, Strike, Corr, 'OutSpec', OutSpec);

displaySummary(Vol1, PriceBJS, VegaBJS)

%
% The change in the price and vega with respect to the volatility of the first asset is plotted below. You can observe that as the volatility of the first asset increases, the price of the spread option also increases. Also, the changes in vega indicate that the price of the spread option becomes more sensitive to the volatility of the first asset and less sensitive to the volatility of the second asset.

figure;

% Plot price for BJS model
subplot(2, 1, 1);
plot(PriceBJS, 'Color', 'red');
title('Price (BJS)');
legend('Price', 'Location', 'EastOutside');
xlabel('Vol of Asset 1');
ax = gca;
ax.XTick = 1:3;
ax.XTickLabel = Vol1;
ylabel('Price');

% Plot of vega for BJS model
subplot(2, 1, 2);
plot(VegaBJS(:,1), 'Color', 'red');
hold on;
plot(VegaBJS(:,2), 'Color', 'blue');
hold off;
title('Vega (BJS)');
legend('Asset 1', 'Asset 2', 'Location', 'EastOutside');
xlabel('Vol of Asset 1');
ax = gca;
ax.XTick = 1:3;
ax.XTickLabel = Vol1;
ax.YLim = [-1 40];
ylabel('Vega');

%
%% Summary
% In this example, European and American spread options are priced and analyzed using various techniques. The Financial Instruments Toolbox™ provides functions for two types of closed form solutions (Kirk, BJS), the finite difference method, and the Monte Carlo simulation method. The closed form solutions are well suited for pricing and sensitivity calculation of European spread options because they are fast. However, they cannot price American spread options. The finite difference method and Monte Carlo method can price both European and American options. However, they are not as fast in pricing European spread options as compared to closed form solutions.
%% References
% [1] Carmona, Rene, Durrleman, Valdo. "Pricing and Hedging Spread Options." SIAM Review. Vol. 45, No. 4, 2003, pp. 627-685.
% [2] Wilmott, Paul, Dewynne, Jeff, Howison, Sam. Option Pricing. Oxford Financial Press, 1993.
% [3] Bjerksund, Petter, Stensland, Gunnar. "Closed form spread option valuation." Department of Finance, NHH, 2006.
% [4] Longstaff, Francis A, Schwartz, Eduardo S. "Valuing American Options by Simulation: A Simple Least-Squares Approach." Anderson Graduate School of Management, UC Los Angeles, 2001.
%
%% Utility Functions

function displayComparison(model1, model2, price1, price2, delta1, delta2, gamma1, gamma2)
% Pad the model name with additional spaces
additionalSpaces = numel(model1) - numel(model2);
if additionalSpaces > 0
    model2 = [model2 repmat(' ', 1, additionalSpaces)];
else
    model1 = [model1 repmat(' ', 1, abs(additionalSpaces))];
end
    
% Comparison of calculated prices
fprintf('Comparison of prices:\n');
fprintf('\n');
fprintf('%s:   % f\n', model1, price1);
fprintf('%s:   % f\n', model2, price2);
fprintf('\n');

% Comparison of Delta
fprintf('Comparison of delta:\n');
fprintf('\n');
fprintf('%s:   % f     % f\n', model1, delta1(1), delta1(2));
fprintf('%s:   % f     % f\n', model2, delta2(1), delta2(2));
fprintf('\n');

% Comparison of Gamma
fprintf('Comparison of gamma:\n');
fprintf('\n');
fprintf('%s:   % f     % f\n', model1, gamma1(1), gamma1(2));
fprintf('%s:   % f     % f\n', model2, gamma2(1), gamma2(2));
fprintf('\n');
end

function displayComparisonPrices(PriceKirk, PriceBJS, PriceFD, Strike)
% Comparison of calculated prices
fprintf('Prices for range of strikes:\n');
fprintf('\n')
fprintf('Kirk     \tBJS      \tFD    \n');
for i = 1:numel(Strike)
    fprintf('%f\t%f\t%f\n', PriceKirk(i), PriceBJS(i), PriceFD(i));
end
end

function displaySummary(Vol1, PriceBJS, VegaBJS)
% Display price
fprintf('Prices for different vol levels in asset 1:\n');
fprintf('\n');
for i = 1:numel(Vol1)
    fprintf('%f\n', PriceBJS(i));
end
fprintf('\n');

% Display vega for first asset
fprintf('Asset 1 vega for different vol levels in asset 1:\n');
fprintf('\n');
for i = 1:numel(Vol1)
    fprintf('%f\n', VegaBJS(i,1));
end
fprintf('\n');

% Display vega for second asset
fprintf('Asset 2 vega for different vol levels in asset 1:\n');
fprintf('\n');
for i = 1:numel(Vol1)
    fprintf('%f\n', VegaBJS(i,2));
end
end

% Copyright 2012-2016 The MathWorks, Inc.