Some Numerical Methods for Solving Stochastic Impulse Control in Natural Gas Storage Facilities

The valuation of gas storage facilities is characterized as a stochastic impulse control problem with finite horizon resulting in Hamilton-Jacobi-Bellman (HJB) equations for the value function. In this context the two catagories of solving schemes for optimal switching are discussed in a stochastic control framework. We reviewed some numerical methods which include approaches related to partial differential equations (PDEs), Markov chain approximation, nonparametric regression, quantization method and some practitioners’ methods. This paper considers optimal switching problem arising in valuation of gas storage contracts for leasing the storage facilities, and investigates the recent developments as well as their advantages and disadvantages of each scheme based on dynamic programming principle (DPP). | Gas storage facility | Stochastic impulse control problems | Optimal switching | Optimal stopping time | Semi-Lagrangian scheme | HJB equation, Viscosity solution | ® 2012 Ibnu Sina Institute. All rights reserved. http://dx.doi.org/10.11113/mjfas.v8n1.121


INTRODUCTION
In the natural gas industry, the modelling and valuation of leases on natural gas storage have been major concerns in the last decade, especially since deregulation of U.S. and Europe energy markets.In this flourishing market, manager/renter of gas storage facilities has been taking advantages of the fluctuation in market prices by releasing natural gas from storage in seasons with high demand.Recently, several authors have focused on solving this operational flexibility problem to maximize the revenue, by exploiting storage operational strategies.These studies, such as [1 through 20], have led to the consideration of the optimal switching models with timing flexibility.In 1997, Pilipovic [16] in her book, which is the first book on energy risk, studied the natural gas price models and storage valuation problem as a price swing straddle, which is an operational method by using the strips of spark-spread options approach such as being used in 2003 by Bringedal [2], Eydeland and Wolyniec [6] and Thompson, Davison and Hassmussen [18].Also, by using real options approach, for instance gas storage can be reduced to a collection of calendar call options.Option pricing methods are intuitive and have fast computational speed of convergence, but it ignores key operational constraints such as dynamic capacity limits [7,18].Moreover, with reference to studies in [1 through 20], the valuation of gas storage facilities is characterized as *Corresponding author at: E-mail addresses: leili_utm@yahoo.com.my(Leyla Ranjbari) a stochastic control problem resulting in Hamilton-Jacobi-Bellman (HJB) equations.As shown in [21,22], the value of a stochastic control problem is identical to the viscosity solution of a HJB equation.As a result, the value of a gas storage contract can be computed by solving the corresponding HJB equation using partial differential equation (PDE)-based approaches such as typically finite difference (FD) methods.It is important to ensure that a numerical scheme converges to the viscosity solution of the corresponding HJB equation.In 2003, Thompson et al. [18] used a Tsitsiklis and Van Roy (TvR) scheme as an explicit time stepping.The TvR scheme suffers from time step restrictions due to stability conditions.To overcome this difficulty, Chen and Forsyth [4, 5 and 23] introduced semi-Lagrangian time-stepping scheme.This implicit finite difference scheme has the same difficulty as the simulationbased methods [3,11] for HJB equation with no bang-bang control, i.e., it is computationally expensive.A control that only takes from a finite set is called bang-bang control.To construct a suitable PDE solver, Ware and Li [20,24] developed a wavelet method coupled with a semi-Lagrangian approach to solve the gas storage HJB equation.Although the wavelet method resolves this difficulty, the convergence proof is problematic.Due to the presence of inventory, the HJB equation is a degenerate PDE and shows the path-dependency of the methods.Also, the PDE methods depend on price models and suffer from dimensionality.
In Ludkovski thesis [11], he introduced the optimal switching approach which reduced the storage facility | 32 | problem to an iterated optimal stopping time over a finite horizon.He developed a robust numerical scheme based on Longstaff-Schwartz Monte Carlo regression (LS) which outperforms the PDE approaches in aspect of dimensionality problems.Besides the capability of easily handling multi-dimensional problems, this scheme considers the operational constraints in the model.This simulation based methods [3,11] can be used directly to solve the stochastic control problem with bang-bang type control.But, for the control which is not of bang-bang type, these methods are computationally expensive.Furthermore, simulation-based schemes have difficulty in achieving high accuracy [5].
In addition to PDE solver and simulation-based approaches, another prevalent approach is lattice or tree models [12,15] which can be considered as stochastic programming approaches.As an example of these approaches, in [15], Parsons presented a two-factor tree approach for valuing natural gas storage leases.His tree model has been tested successfully against historical data to quantify the optionality in storage leases.Meanwhile, lattice/tree methods such as trinomial and binomial tree approaches [7, 15 and 18], are resulted from Bellman's optimal principle and used in real options due to being flexible, relatively simple and readily implemented in a computer program.Thus it can be widely used to evaluate options.However, trees cannot handle optimal exercise strategy problems arisen in natural gas storage [18].There are two reasons for this.First, such trees are just explicit finite difference methods for solving parabolic PDEs.In the case of natural gas storage the operating characteristics lead to equations of a parabolic and hyperbolic nature.Since there are too many state variables in the problem, including the storage level, spot price, and forward prices, it is very difficult to solve the hyperbolic equations by tree models or by solving a finite difference equation [10].Then to solve the hyperbolic equations, we need another technique more sophisticated than tree methods.Secondly, trees cannot handle price spikes properly [18].

Introduction
The gas storage facility problem discusses about determining the value of gas storage contracts for leasing the storage facilities.Owners of storage facilities lease out space within, and a leaseholder has the right to inject into or withdraw from the facility only for a pre-specified period of time, and within pre-specified volume constraints, which are described through a ratchet schedule [15].A ratchet schedule is a schedule of all possible inventory levels (or ratchets) and their associated daily maximum injection and withdrawal rates.As the lease-holder injects or withdraws, she tends to maximize the revenue by only exploiting storage operational strategies and cover the dealing costs of operating the facility.Owners/operators of storage facilities are not necessarily the owners of the gas held in storage.Indeed, most working gas held in storage facilities is held under lease with shippers, local distribution companies (LDCs), or end users who own the gas.On the other hand, the type of entity that owns/operates the facility will determine to some extent how that facility's storage capacity is utilized.

2.2
Why we choose the salt cavern storage?
Among the underground storage facilities, i.e., depleted reservoirs in oil and/or gas fields, aquifers, salt cavern formations and artificial cavern, salt cavern is a common concern among researchers due to providing very high withdrawal and injection rates relative to their working gas capacity.Since 1993, the natural gas storage industry has attempted to profit from the changes of the market conditions.The deregulation of energy market allows using the storage facility to inject or withdraw as changes in price levels or arbitrage opportunities presents in the market.Therefore, the facilities with more rapid cycle in their inventories, i.e., completely withdraw and refill working gas or vice versa, are more suitable to the flexible operational needs of today's storage users.On the other hand, the maximum profit obtains from high deliverability storage sites, which include salt cavern storage reservoirs as well as some depleted oil or gas reservoirs due to the largest daily withdrawal capability, [18,19].

Storage facility problem
In this section, the corresponding continuous-time stochastic impulse control is introduced, with the notations following Ludkovski's work [11].This operational flexibility problem is modelled as an optimal switching with timing flexibility to maximize the revenue by exploiting storage operational strategies [1 through 20].Before going into the problem, we suppose some assumption on owner and gas market.

Assumption on manager/renter
1.She is rational.The manager/renter makes choices based on her rational outlook, available information and past experiences.And, government policy does not influence on her decisions.In this idea if the owner of storage facility believes that the price for gas will be higher in the future, it will stop or slow supply until the price rises.Consequently, the demand stays the same and gas prices will increase.In sum, the producer believes that the price will rise in the future, makes a rational decision to slow production and this decision partially affect what happens in the future.In conclusion, these rational expectations of the players will partially affect what happens to the economy in the future.2.She is risk-neutral.Risk-neutral investor is only concerned with an investment's expected return and overlooks risk in her investments.
| 33 | 3.She is concerned to maximize the value function of her lease agreement over a finite time.4.She is a price-taker.Her buying or selling transactions have no effect on the market in a form of monopoly.
The owner can alter the injection and withdrawal rates without significantly affecting the market price of natural gas, [3].

Assumption on gas market
1. Market is liquid.In this market, due to high level of trading, trading does not affect the gas price.And, the gas can be easily sold or bought that reflects the ability to convert the gas to cash quickly or being marketability.Then, it is safer to invest in liquid assets than illiquid ones because it is easier for an investor to get his/her money out of the investment.We can say that the higher liquidity, the lower spread and volatility.Although there is no specific liquidity formula, liquidity ratios is used as a measure of liquidity.
2. Market is fairly liberal, i.e., all participants have the opportunity to rent a storage facility to speculate and take the advantages of the volatility in prices to maximize the profit.
3. Gas market shows a strong seasonality which leads to intrinsic value of storage.To buy in low price, the owner locks-in forward contracts and sell in high price in heating/cooling days.

Stochastic impulse control
We study the following optimal switching problem from Carmona and Ludkovski's work [3] as the storage facility model, with no changes or supplements.In advance, we label the three above operational strategies in the way that injection (-1), store or doing nothing (0) and withdrawal (1), i.e., we have three operational regimes  ∈ {−1,0,1} by the payoff rate   (,   ,   ) in $  from running the facility in regime  .
Given the initial inventory  0 =  and the storage strategy u, the future inventory   � () is completely determined.Namely,   � () satisfies the ordinary differential equation In addition, we also use the following notation. 0 Each change of the facility's regime incurs switching costs.In particular, moving the facility from regime i to regime j costs  , = (, ; ,   ,   ).We assume that the switching costs are discrete:  , ≥  for all  ≠  and some  > 0 , and  ,=0 .In this research, among the storage facilities we only focus on salt dome facilities the switching costs are economically negligible; however, in other facilities, switching costs may be significant.Also, switching costs are assumed strictly positive so that it prohibits the owner would repeatedly change the regimes back-and-forth over a very short time.Since the ultimate computations are in discrete time, switching costs can be set to zero on implementation-level.
| 35 | resulted from applying Bellman's optimal principle.In fact, a smooth solution of the quasi-variational Hamilton-Jacobi-Bellman (HJB) inequality is the value function of the impulse control problem.To prove the uniqueness of solutions, the notion of viscosity solutions is used as weak solutions.Indeed, the value function is always the unique viscosity solution of the system in verification theorem [11, 22 and 26].The PDE solver relies on the quasi-variational formulation of verification theorem.Generally speaking, the PDE methods transform the stochastic control problem into a parabolic partial differential equation with a free boundary.To solve this parabolic PDE, there are a variety of tools such as the basic finite difference scheme (FD).By setting up a space-time grid, we solve the PDE by replacing derivatives with finite differences in the HJB equation and directly enforcing the barrier condition at each step.In spite of being the simplest setting, the approaches pose a multitude of challenges that prevent rigorous solutions.An FD method is easy to implement but suffers several major drawbacks.First, FD often needs a large number of time steps to be of numerical stability.Second, the method suffers from dimensionality problem: the size of the space grid is exponential in number of dimensions of price process d and generally speaking d > 3 is computationally infeasible.Finally, the switching boundary will inevitably be jagged due to the presence of a grid in the x-space.The last point can be alleviated with the use of an adaptive grid.The method's accuracy depends on the order of the approximation used for the derivatives of u, which should be at least (∆ + ∆ 2 ).

Nonparametric regressions
This approach is related to the third strategy as we mentioned before to compute the conditional expectations.To compute the conditional expectation, there are several methods, including Malliavin calculus, Monte Carlo simulations and regression.In regression schemes we have a choice between regression against basis functions and fully non-parametric regression.Choosing a non-parametric regression relieves us of the concerns regarding selecting appropriate basis functions and may produce smoother conditional distributions [11].The simplest version of this approach is k-nearest neighbours multivariate kernel regression.Thus, this regression is replaced by a local linear combination of the other paths' values with the weights proportional to the distance.The use of nearest neighbours reduces the dimensionality difficulty.It is attractive due to its robustness and involving no additional error.However, the two major difficulties with the kernel method are selecting an appropriate bandwidth and the computations around the edges.The bandwidth h controls the peak of the weights around (x, y).Thus, as the bandwidth increases, more distant points carry more weight and the estimate becomes more 'global'.However, choosing h is heuristic and may require a lot of trial-and-error.The other difficulty is when the regressed value is extreme; the nearest neighbours have bias for the extreme values in contradiction to regression algorithm.

Quantization method
A powerful non-Markovian version of the second strategy, i.e., replacing the continuous-space process by a discrete approximation, with an adaptive approximating grid is the so-called Quantization Scheme (QS).The main motivation of quantization is to find a small and efficient approximating grid in exchange for giving up the Markov property and closed-form formulae for transition probability matrix ℙ  =   (, ).This method instead of computing ℙ  via PDE solver such as FD, uses Monte Carlo simulation.Whiles this is likely to be slow, it can be done just once off-line and stored for later calculations.The gain will increase the robustness and results in much better dimensional scaling.The transition matrix is required to approximate closely the dynamics of between cells of the adjacent quantization grids.The quantization scheme can be summarized in these following stages: 1. Construct the quantization grid by Voronoi tessellations (cells) as a partition of state space E.

Simulate
� � ,  � � with, for instance, Competitive Learning Vector Quantizer (CLVQ) algorithm.3. Solve the optimal switching problem by computing the pseudo-Snell envelops of the non-Markov �  � �.The advantageous of the quantization scheme can be mentioned as: 1.Since the entire algorithm for constructing�  � � is Monte Carlo based, it is considered to be robust.2. Simulation can be parallelized.
3. Once the grids and transition are computed, any control problem is solved because the approximation of �  � � is independent from the optimization steps in contrast to regression algorithm.However, there are some drawbacks to the method: 1. Estimation of the transition probability  � is problematic.2. To achieve an acceptable accuracy, the method needs to run in large number of simulations.
Ludkovski [3,11] replaced the optimal switching problem with finite horizon by a recursive optimal stopping problems in finite time.Then, the exponential maturity randomization was used to construct the iterative sequence of infinite horizon stopping times.In stead of directly solving the obtained iterative problems, he used LS Monte Carlo regression scheme to compute the conditional expectation by exploiting the concept of Snell envelops for revenue functions.

PRACTITIONER METHODS
Besides these methods, there are some practitioner methods to solve the operational flexibility problem of natural gas storage facility which are mentioned in the following:

Classical net present value (NPV) theory
In Dixit [27], he used discounted cash flow analysis to estimate the value of the asset based on projections of future prices and proper weighing and discounting of possible cases.Uncertainty was essentially eliminated, as static scenarios were used to forecast the future and select pre-determined optimal behaviour.The opportunity of dynamically responding to prices was ignored and as a result the contracts were consistently underpriced.

Markov decision process (MDP)
This is a stochastic programming approach implemented in Yushkevich [28].In MDP, we have the tree-based versions of the stochastic control formulation, and there exist sequential decision problems under uncertainty.The problem is discretized in time and the path space of (  ) is broken into a finite number of actions.Then the transition probabilities (, , ∆) are estimated for each current outcome x and possible transition action y.Finally, the problem is solved via dynamic programming that corresponds to a lattice discretization of the Quasi-Variational Inequality (QVI).Indeed, MDP or stochastic programming approach is a dynamic construction of the decision tree.The goals of using MDP are mentioned in below: (A) In finite horizon problem, the aim is to maximize the expected reward for the next n steps.(B) In infinite horizon, maximizing the expected discounted reward is considered.
The MDP tree is simple and intuitive to understand.However, if one must solve numerically then the computational complexity explodes for long horizons with much optionality.

Strips of spark-spread options approach
This approach is being carried out in Bringedal [2], Eydeland and Wolyniec [6] and Thompson, Davidson and Rassmussen [18].The motivation is to reduce the problem to pricing standard financial options whose valuation is well understood.Accordingly, the payoff from the power plant is represented as a collection of European options that pay the maximum value to be obtained during each period.This method is intuitive and have fast computational speed of convergence, but it ignores key operational constraints such as dynamic capacity limits.

CONCLUSION
As a comparison for the numerical methods considered in this paper, one can mention that in small dimensions (d < 3) the best algorithm is the PDE approache such as semi-Lagrangian and TvR schemes.Otherwise, for d > 2, the LS and Quantization schemes should be used due to scalability to high dimensional problems.An industrial-strength implementation should be very fast and produce provable accurate results.And, to choose between LS and quantization schemes, if one is looking for a quick tool to solve optimal switching for a variety of (  ), the LS scheme is the best.On the other hand, if (  ) is fixed, the quantization may be better.
• Level of inventory in storage denoted by   .• Finite cave capacity represented by   ≤   ≤   .• Constant discount (interest) rate r.• Denote the injection rate by   (  ), quoted in Bcf per day.Injection of   (  ) Bcf of gas requires the purchase of   (  ) ≥   (  ) Bcf on the open market.• Similarly the withdrawal rate is labelled   (  ) and causes a market sale of   (  ) ≤  Bcf.• Capacity charges   (,   ) in each regime.• The case   ≠   indicates gas loss during injection/withdrawal (typically on the scale of 0.25%−1% for salt dome storage).The transmission rates   ,   are a function of   and are based on gas pressure laws.