MATLAB Examples

Intrinsic Storage Optimization

This example demonstrates optimizing a storage facility and valuing a storage contract using intrinsic valuation. The optimization involves finding the optimal positions in a set of forward natural gas contracts that maximize the storage value while meeting all the injection, withdrawal and volume constraints.

Contents

Import Futures Quotes

The dataset command is used to import NYMEX Natural Gas futures market quotes from an Excel spreadsheet.

data = dataset('xlsfile','NGFutQuote.xlsx');

% Display first 24 quotes
display(data(1:24,:))
Warning: Variable names were modified to make them valid MATLAB identifiers. 

ans = 

    Symbol         Month              Last     Chg       Bid      Ask      DayHigh    DayLow    Vol      OpInt 
    'NGQ10'        '8/1/2010'         4.657    -0.197    4.656    4.657    4.871      4.613     75166    127454
    'NGU10'        '9/1/2010'         4.694    -0.186    4.687    4.689    4.893      4.648     15262    143591
    'NGV10'        '10/1/2010'        4.777    -0.162    4.765    4.767    4.954      4.729     15972     82795
    'NGX10'        '11/1/2010'            5    -0.158    4.995        5    5.166       4.96      6209     41194
    'NGZ10'        '12/1/2010'        5.259    -0.142    5.255    5.261      5.4       5.22      4344     38750
    'NGF11'        '1/1/2011'         5.428    -0.137     5.42    5.427    5.557      5.393      7859     57108
    'NGG11'        '2/1/2011'         5.393    -0.132    5.392    5.402    5.513      5.377      1109     18747
    'NGH11'        '3/1/2011'         5.303    -0.112    5.294    5.305    5.412      5.273      2373     53074
    'NGJ11'        '4/1/2011'         5.072    -0.083     5.06    5.073    5.139      5.038      2083     40850
    'NGK11'        '5/1/2011'          5.07    -0.098    5.073    5.094    5.137      5.053       491     19745
    'NGM11'        '6/1/2011'         5.125    -0.087    5.124    5.141    5.202      5.099       168      8344
    'NGN11'        '7/1/2011'         5.193    -0.079     5.18    5.202     5.21      5.162        99      5742
    'NGQ11'        '8/1/2011'           5.2    -0.116    5.218    5.248    5.295        5.2        49      5307
    'NGU11'        '9/1/2011'         5.265    -0.081     5.25    5.278    5.265       5.23       184      5292
    'NGV11'        '10/1/2011'         5.35    -0.081    5.342    5.366     5.38       5.32      1154     18513
    'NGX11'        '11/1/2011'        5.555    -0.096    5.551      5.6      5.6      5.555       141      4856
    'NGZ11'        '12/1/2011'        5.805      -0.1      5.8    5.851     5.85       5.79        19     12430
    'NGF12'        '1/1/2012'         6.005    -0.066    5.962     6.01     6.01       5.99        33     10575
    'NGG12'        '2/1/2012'         5.964    -0.044     5.88     5.93    5.964      5.964         2      3895
    'NGH12'        '3/1/2012'         5.845         0    5.721     5.79        0          0         0     10795
    'NGJ12'        '4/1/2012'         5.285      -0.1    5.261     5.34     5.29      5.285         5     11288
    'NGK12'        '5/1/2012'           5.3    -0.093     5.26    5.362    5.301        5.3         5      3399
    'NGM12'        '6/1/2012'          5.34    -0.091    5.289    5.407     5.35       5.34        58      1686
    'NGN12'        '7/1/2012'          5.39    -0.101      5.2        0     5.39       5.39         4      1556

Storage parameters

These are parameters that define the storage facility, its constraints and costs.

facility.Vmax = 1e6; % MMBtu (Maximum gas volume for facility)
facility.Vmin = 0; % MMBtu (Minimum gas volume for facility)
facility.Imax = 8197; % MMBtu/day (max injection rate)
facility.Wmax = 16393; % MMBtu/day (max withdrawal rate)
facility.cI = .01; % $/MMBtu (cost of injection)
facility.cW = .01; % $/MMBtu (cost of withdrawal)

display(facility)
facility = 

    Vmax: 1000000
    Vmin: 0
    Imax: 8197
    Wmax: 16393
      cI: 0.010000000000000
      cW: 0.010000000000000

Valuation Parameters

These are parameters that define the storage valuation

valDate = '1 Jul 2010'; % Valuation date
startDate = '1 Aug 2010'; % Contract start time
endDate = '1 Aug 2011'; % Contract end time
startVol = 0; % Start volume
endVol = 0; % End volume
r = 3.5/100; % Risk free interest rate
liquiditySpread = .02; % If either the bid and ask are 0, assume this is the market spread

Process Forward Curve Dataset

Extract only the required contracts for the valuation and compute the adjusted bid and ask prices that take into account injection/withdrawal costs and market liquidity.

sDate = datenum(startDate);
eDate = datenum(endDate);
monthN = datenum(data.Month);

% Slice data set to obtain only the contracts for the required months
data = data(monthN >= sDate & monthN < eDate, {'Month', 'Bid', 'Ask', 'Last'});

% Add a spread to illiquid contracts, if necessary
data.Bid(data.Bid == 0) = data.Last(data.Bid == 0) - liquiditySpread/2;
data.Ask(data.Ask == 0) = data.Last(data.Ask == 0) + liquiditySpread/2;

% Incorporate injection and withdrawal costs
adjBid = data.Bid - facility.cW;
adjAsk = data.Ask + facility.cI;

% Display the adjusted forward curve bid and ask prices.
plot(datenum(data.Month), [adjBid adjAsk] , 'o-');
datetick('x','mmmyy'); legend('Adj. Bid', 'Adj. Ask'); grid on
xlabel('Expiration'); ylabel('Adjusted Forward Price $/MMBtu');
title('Contract costs adjusted for market spread and cost of injection/withdrawal');

Set up the Optimization

Compute the objective function and constraints for the linear program that will be used to optimize the storage facility and compute an intrinsic value. The variables in the optimization are N long and N short positions in the futures contracts

% Objective function
% This means we are maximizing the sum of prices we pay for injections and
% those we receive for withdrawals. -f' * x, where x is the number of
% long and short contracts for each forward.
f = [adjAsk ; -adjBid];

% Parameters
N = length(adjBid);
daysInMonth = monthDays(datenum(data.Month));

% Constraints
% These are inequality, equality and bound constraints. These ensure the
% volume of the facility is always between Vmin and Vmax and that the
% volume at the end of the period is endVol.
[A, b, Aeq, beq, lb, ub] = createConstraints(facility, N, daysInMonth, startVol, endVol);

Run optimization

Feed objective function and constraints to linear program solver

optimPos = linprog(f, A, b, Aeq, beq, lb, ub);
Optimization terminated.

Display solution

Display storage value, positions, volumes and profits by contract

optimPos = round(optimPos);

soln = data(:,{'Month','Bid','Ask'});
soln.Injection = optimPos(1:N);
soln.Withdrawal = optimPos(N+1:end);
soln.Volume = startVol + cumsum(soln.Injection - soln.Withdrawal);
soln.Capital = round(cumsum(-f(1:N) .* soln.Injection - f(N+1:end) .* soln.Withdrawal));

% Calculate optimized intrinsic storage value
storageValue     =  -f(1:N)' * soln.Injection - f(N+1:end)' * soln.Withdrawal;
discStorageValue = (-f(1:N) .* soln.Injection - f(N+1:end) .* soln.Withdrawal)'...
                   * exp(-r*(datenum(data.Month) - datenum(valDate))/365);

format long
disp('--------- Intrinsic Valuation ------------');
fprintf('Total Storage Contract Value: $%0.2f\n', storageValue);
fprintf('Discounted Storage Contract Value: $%0.2f\n\n', discStorageValue);
disp('Schedule:');
fprintf('Start Volume: %d\n', startVol);
disp(soln);
--------- Intrinsic Valuation ------------
Total Storage Contract Value: $644177.48
Discounted Storage Contract Value: $557807.90

Schedule:
Start Volume: 0

    Month              Bid      Ask      Injection    Withdrawal    Volume     Capital 
    '8/1/2010'         4.656    4.657    254107            0         254107    -1185917
    '9/1/2010'         4.687    4.689    245910            0         500017    -2341448
    '10/1/2010'        4.765    4.767    254107            0         754124    -3555318
    '11/1/2010'        4.995        5    245876            0        1000000    -4787156
    '12/1/2010'        5.255    5.261         0            0        1000000    -4787156
    '1/1/2011'          5.42    5.427         0       508183         491817    -2037886
    '2/1/2011'         5.392    5.402         0       459004          32813      432473
    '3/1/2011'         5.294    5.305         0        32813              0      605857
    '4/1/2011'          5.06    5.073    245910            0         245910     -644103
    '5/1/2011'         5.073    5.094    254107            0         500017    -1941066
    '6/1/2011'         5.124    5.141      8166            0         508183    -1983129
    '7/1/2011'          5.18    5.202         0       508183              0      644177

Visualize results

Display the results visually in a custom plot.

plotResults(soln, facility, startVol, endVol)