- Email: [email protected]

Mike Ludkovski Department of Statistics and Applied Probability, University of California, Santa Barbara, CA 93106 [email protected], http://www.pstat.ucsb.edu/faculty/ludkovski

We consider the valuation of energy storage facilities within the framework of stochastic control. Our two main examples are natural gas dome storage and hydroelectric pumped storage. Focusing on the timing flexibility aspect of the problem we construct an optimal switching model with inventory. Thus, the manager has a constrained compound American option on the inter-temporal spread of the commodity prices. Extending the methodology from Carmona and Ludkovski (2008), we then construct a robust numerical scheme based on Monte Carlo regressions. Our simulation method can handle a generic Markovian price model and easily incorporates many operational features and constraints. To overcome the main challenge of the path-dependent storage levels two numerical approaches are proposed. The resulting scheme is compared to the traditional quasi-variational framework and illustrated with several concrete examples. We also consider related problems of interest, such as supply guarantees and mines management. Key words : gas storage; optimal switching; least squares Monte Carlo; hydro pumped storage; impulse control, commodity derivatives History : First Version: July 2005; This Version: August 24, 2008

Acknowledgments We thank the participants of the Banff BIRS Workshop 07w-5502 “Mathematics and the Environment” and Zhenwei J. Qin for many useful comments and discussions. We also thank the anonymous referees whose feedback led to much improved presentation.

1.

Introduction.

While classical financial contracts such as stocks and bonds are paper assets, ownership of commodities entails physical storage. As a result, the modern commodities industry incorporates an extensive storage infrastructure, including natural gas salt domes, liquified natural gas (LNG) storage tanks, precious metal repositories and hydroelectric reservoirs. In the last decade, with the ongoing deregulation of these industries, storage facilities have also acquired an important role 1

Carmona and Ludkovski: Optimal Switching for Energy Storage

2

in the commodity financial markets. Storage allows for inter-temporal transfer of the commodity and permits exploitation of the fluctuating market prices. The basic principle is to ‘buy low’ and ‘sell high’, such that the realized profit covers the intermediate storage and operating costs. Since the profitable opportunities are driven by price volatility, the storage facility grants its owner a calendar straddle option. Traditionally, storage facilities have been owned by major players in the respective industries who had the enormous capital typically needed to build and maintain them. However, with the liberalized markets, all participants have nowadays the opportunity to rent a storage facility with an eye towards speculation on prices and aggressive profit maximization. For example, with respect to natural gas de Jong and Walet (2003) write that “natural gas storage is unbundled, ... [and] offered as a distinct, separately charged service. ... Buyers and sellers of natural gas have the possibility to use storage capacity to take advantages of the volatility in prices”. The aforementioned price volatility can be either systemic or speculative. For instance, the natural gas market exhibits strong seasonality since the main consumer group is households that use gas for winter heating. Thus, natural gas demand (and prices) has a systematic spike in the cold season, often on the order of 30%-50%. In contrast, in say silver, the volatility in prices is almost entirely speculative, but nevertheless can sometimes lead to price swings of 100% within one year, cf. the Hunt brothers episode in late 1970s (Pirrong 1996). The seasonality effects present in certain markets lead to intrinsic value of storage that can be locked-in by a static purchase and sale of forward contracts. For instance, a typical July-January forward spread in natural gas is on the order of $1/MMBtu and can be easily realized by a simple one-time transaction. However, the presence of an increasingly liquid short-term market permits further dynamic optimization. Intuitively, the manager holds timing options that allow her to optimally exploit opportunities that appear as market prices evolve. Capturing this optionality is crucial in the present competitive markets, especially with the entry of non-energy players who rent facilities with the sole goal of maximizing profit (as opposed to old-fashioned participants who also have strategic aims). Moreover, with the growing importance of energy commodities, sophisticated valuation of energy storage becomes an integral aspect of functioning financial markets.1 Thus, it is necessary to be able to compute the financial extrinsic value of such flexibility. Namely, how much should one pay to gain control of a storage facility for a period of T years? The simple question above hides the associated modeling difficulties. Indeed, the owner faces a multitude of 1

For instance, the $5 billion loss reported by Amaranth LLC. in Fall 2006 was due to a poorly managed bet on the March-April calendar spread in natural gas, which is in turn driven by actions of storage managers during the winter.

Carmona and Ludkovski: Optimal Switching for Energy Storage

3

optionalities and constraints that interact in a nonlinear fashion. First, the purchases and sales can be done immediately, using spot prices, or forward in time using forward prices. Second, the bought commodity must be put on the inventory. As a result, inventory capacity limits, as well as storage costs, delivery charges and other operational and engineering constraints become crucial. However, the latter are intrinsically path-dependent in terms of the storage strategy adopted. Finally, the manager may be exposed to margin requirements (if the commodity is bought with credit), mechanical break-downs and other external events. To overcome these challenges the current literature on commodity storage has largely proceeded in two different directions. The basic practitioner methods have been based on the traditional option-pricing approach. Thus, one makes (often drastic) simplifications to shoehorn the problem into the option pricing framework. For instance, gas storage can be reduced to a collection of calendar Call options paying out the spread between gas prices today and k∆t years from now, with k = 1, . . . , T /∆t (Eydeland and Wolyniec 2003). After such a reduction the extensive existing machinery of derivative pricing can be imported. One gains intuition and computational speed but ignores key operational constraints (such as dynamic capacity limits), since the calendar Calls are priced independently of each other. Furthermore, the method is ad hoc, requiring heuristic adjustments to correct for model assumptions. Alternatively, various stochastic programming algorithms (Nowak and R¨ omisch 2000, Doege et al. 2006) have been considered, especially for hydrothermal systems. These methods maintain the flexibility of incorporating realistic constraints, but instead discretize the set of future scenarios. While powerful, stochastic programming suffers from nonscalability with respect to number of scenarios and time-steps used. To properly account for the interdependence between the timing optionality of the manager in choosing the purchase and sale times and the inventory constraints, one must consider the full stochastic control framework. This leads to a Bellman dynamic programming equation for the value function. From here one may apply the Hamilton-Jacobi-Bellman theory, translating the problem into a quasi-variational partial differential equation (pde) formulation. Such models were studied in Ahn et al. (2002) and Thompson et al. (2003) and solution is then obtained via standard numerical solvers. However, the path-dependency due to presence of inventory implies that the pde is degenerate (convection-dominated) and therefore extra care is necessary. Moreover, the implementation is necessarily price-model dependent and consequently not robust. One approach that mitigates these challenges was recently proposed in Chen and Forsyth (2007). In this paper, we also adopt the stochastic control formulation. However, in contrast to the pde methods above, we proceed to a probabilistic solution based on optimal multiple stopping

Carmona and Ludkovski: Optimal Switching for Energy Storage

4

problems. This perspective allows us to obtain an efficient simulation-based numerical method for valuing energy storage on a finite horizon. The method is flexible and not tied to a particular class of asset prices; in fact we abstract from asset dynamics and take as exogenous the (multidimensional) Markov price process for the commodity. Thus, use of complex price dynamics, such as jump-diffusions, several factors, etc. has only a marginal impact on the efficiency of the algorithm. Compared to previous approaches, our method has several advantages. First, we maintain rigorous modeling of the operational constraints while considering the entire set of future scenarios of commodity prices. Moreover, in contrast to pde solvers which suffer from the curse of dimensionality, our scheme can easily handle multi-dimensional settings. In terms of performance, our scheme is competitive with the pde solvers in one-dimension and is clearly superior in higher dimensions (which are essential in a realistic model, see Section 7). Thanks to its scalability, the algorithm is easily extendable and therefore suitable for realistic use. Thus, our main contribution is a robust numerical method that remains on firm theoretical grounds of stochastic control while bridging the gap to practitioner needs. To be concrete, from now on we focus on the representative example of controlling a natural gas salt dome facility; other applications are addressed in Sections 6 and 7. The rest of the paper is structured as follows. Section 2 describes the stochastic control model we use and its relation to existing literature. Section 3 summarizes the theoretical solution method which is then implemented in Section 4. After outlining the numerical scheme, we proceed to illustrative examples in Section 5. Sections 6 and 7 discuss hydroelectric pumped storage and several problems in natural resource management and demonstrate that our methodology is applicable to a wide variety of real options encountered among commodity derivatives. Finally, Section 8 concludes and outlines future projects.

2.

Stochastic Model.

Natural gas storage is currently the most widespread class of commodity storage infrastructure in the US2 (FERC 2004). A variety of storage options, including depleted gas fields, aquifers, salt domes and artificial caverns are available. In 2006 over 400 such facilities existed in the US and a substantial portion are contracted out for periods of 6–60 months. In the near future, the industry will expand even more with the rolling out of LNG technology and associated storage in North America (see Geman (2005) for the general trends and organization of the gas universe.). In this article we will specifically focus on the case of salt domes which permit the highest rates of injection and withdrawal and therefore contain the most timing optionality (see Table 3). 2

Throughout we focus on the North American markets and use imperial system units.

Carmona and Ludkovski: Optimal Switching for Energy Storage

5

A salt dome is an underground natural cave that can store several billion cubic feet of gas (Bcf). It is connected via pumps to the national pipeline system which allows to inject/withdraw gas at a deliverability rate of 0.1 − 0.4Bcf per day. Taking the point of view of the renter, or manager of such a cave, we now wish to maximize economic value by optimizing the dispatching policy, i.e. dynamically deciding when gas is injected and withdrawn, as time and market conditions evolve. We assume that the manager is rational and risk-neutral and aims to maximize total expected revenue over the finite horizon of her rental. We moreover assume that the respective financial markets are liquid and the manager is a price-taker (the situation of price impact is treated in Section 7). The ingredients of our model can be now listed as: • Time horizon T , with a stipulation for the final state of the facility, see (6). • Market gas prices given by a Markov continuous-time stochastic process (Gt ), Gt ∈ Rd , quoted

in dollars per million of British thermal units (MMBtu), with 1 Bcf ≡ 106 MMBtu. • Level of inventory in storage denoted by Ct . • Finite cave capacity represented by cmin ≤ Ct ≤ cmax . • Constant discount (interest) rate r. • Three possible operating regimes of the storage facility: injection, storage and withdrawal. • Denote by ain (Ct ) the injection rate, quoted in Bcf per day. Injection of ain (Ct ) Bcf of gas,

requires the purchase of bin (Ct ) ≥ ain (Ct ) Bcf on the open market. • Similarly the withdrawal rate is labelled aout (Ct ) and causes a market sale of bout (Ct ) ≤

aout (Ct ) Bcf. • Capacity charges Ki (t, Ct ) in each regime that represent direct storage costs, delivery charges,

various O&M costs and seepage losses. The case bi 6= ai indicates gas loss during injection/withdrawal (typically on the scale of 0.25% − 1% for salt dome storage). The transmission rates ai , bi themselves are fixed by the physical characteristics of the facility; they are a function of Ct and are based on gas pressure laws (Thompson et al. 2003). Remark 1. Typically, Gt would represent the price at time t of the near-month forward contract, which is by far the most liquid contract on the market3 . However, given a variety of quoted gas prices (spot, balance-of-the-month, futures, etc.), we remain agnostic about the precise interpretation of the (Gt ) process. The driving process (Gt ) may also include longer maturity forwards. 3 Recent daily volume on NYMEX has been over 90, 000 contracts, with more than 50% of the trades in the nearmonth.

Carmona and Ludkovski: Optimal Switching for Energy Storage

6

Unfortunately, forward selling is problematic, since the sale price is locked-in in advance, while the inventory only changes at delivery time. We assume for simplicity that any purchase or sale is immediately reflected in the current inventory. Label the three regimes above as i ∈ {−1, 0, 1} and denote by ψi (Gt , Ct ) the payoff rate (in $/year) from running the facility in regime i. Then ψi ’s and the corresponding volumetric changes in inventory are given by Inject: ψ−1 (t, Gt , Ct ) = −Gt · bin − K−1 (Ct ), Store: ψ0 (t, Gt , Ct ) = −K0 (Ct ), Withdraw: ψ1 (t, Gt , Ct ) = +Gt · bout − K1 (Ct ),

dCt = ain (Ct ) dt, dCt = a0 (Ct ) dt, dCt = −aout (Ct ) dt.

(1)

In principle, the facility can also be operated at a sub-maximal transfer rate, however when the monetary reward is linear in the pumping rate as in (1), it is always optimal to inject/withdraw at maximum speed. This is the so-called ‘bang-bang’ property of stochastic control problems (Øksendal and Sulem 2005). Many possibilities exist for the form of (Gt ) and there is much recent debate (see Eydeland and Wolyniec (2003), Cartea and Williams (2008)) about appropriate models for gas prices. A standard choice is an Itˆ o diffusion described by a stochastic differential equation (SDE) dGt = µ(t, Gt ) dt + σ(t, Gt ) · dWt ,

(2)

where Wt is a d-dimensional Brownian motion and σ(t, g) is a d × d non-degenerate volatility matrix. A canonical example (see e.g. Jaillet et al. (2004)) is a one-dimensional exponential OrnsteinUhlenbeck process, namely

or

dGt = Gt κ(θ − log Gt ) dt + σ dWt , σ2 − log Gt dt + σ dWt , d(log Gt ) = κ θ − 2κ

(3) G0 = g.

This models the mean-reversion (to the average level eθ ) in gas prices documented by Eydeland and Wolyniec (2003), while keeping log Gt conditionally Gaussian. Upward jumps in (Gt ) can also be considered and may be used to take into account price spikes. The jury is still out whether such jump-diffusion models are appropriate for natural gas. Other possibilities for (Gt ) could include regime-switching, stochastic mean reversion levels, latent factors, L´evy processes, etc. Our method is independent of the assumed model for (Gt ), and in general we only make the following technical Assumption 1 (A) (Gt ) is a d-dimensional, strong Markov, non-exploding process in Rd .

Carmona and Ludkovski: Optimal Switching for Energy Storage

7

(B) The information filtration F = (Ft ) on the stochastic basis (Ω, F , P) is the natural filtration of (Gt ). (C) The reward rate ψi : [0, T ] × Rd × [cmin , cmax ] → R is a jointly Lipschitz-continuous function of (t, g, c) and satisfies "

E

# sup |ψi (t, Gt , Ct )|2 G0 = g, C0 = c < ∞,

∀g, c, i.

t∈[0,T ]

For notational clarity we suppress from now on the dependency of ψi and the coefficients of (2) on time t. Remark 2. Above we have stated the model in continuous-time. This is to conform to classical financial stochastic control models; since the final implementation is computer-based and consequently performed in discrete-time, one could also work in discrete-time from the beginning. 2.1.

Control Problem.

The flexibility available to the manager is specified via the set U of possible storage policies u. For t ∈ [0, T ], ut ∈ {−1, 0, 1} denotes the (dynamically chosen) operating regime of the facility. It is convenient to write u = (ξ1 , ξ2 , . . . ; τ1 , τ2 , . . .), where the variables ξk ∈ {−1, 0, 1} denote the sequence of operating regimes taken by u and τk ≤ τk+1 ≤ T denote the switching times. Thus, P ut = k ξk 1[τk ,τk+1 ) (t), where by convention τ0 = 0, ξ0 = i0 is the initial facility state. Given the initial inventory C0 = c and the storage strategy u, the future inventory C¯t (u) is completely determined. Namely, C¯t (u) satisfies the ordinary differential equation dC¯s (u) = aus (C¯s (u)) ds, In the sequel we will also use the notation C¯t (c, i) , c +

C¯0 (u) = c. Rt 0

(4)

ai (C¯s (c, i)) ds.

Each change of the facility’s regime incurs switching costs. In particular, moving the facility from regime i to regime j costs Ki,j = K(i, j; t, Gt , Ct ). This represents both the effort —one must dispatch workers, coordinate with the outgoing pipeline, stop/start the decompressors, etc.—and the time needed to change the operating mode. We assume that the switching costs are discrete: Ki,j > for all i 6= j and some > 0, and Ki,i = 0. For actual salt dome facilities the switching costs are economically negligible; however, in other applications, such as hydro pumped storage, switching costs may be significant. Also, strictly positive switching costs are needed for technical reasons in our continuous-time model in order to guarantee existence of optimal finite switching strategies (i.e. to rule out chattering, where the owner would repeatedly change the regimes backand-forth over a very short amount of time). Since the ultimate computations are in discrete time, switching costs can be set to zero on implementation-level.

Carmona and Ludkovski: Optimal Switching for Energy Storage

8

A necessary condition for u to belong to the set U of admissible strategies is to be F-adapted, right-continuous and of P-a.s. finite variation on [0, T ]. F-adaptiveness is a standard condition implying that the agent only has access to the observed price process and cannot use any other information. Finite variation means that the number of switching decisions must be finite almost surely. Thus, P[τk < T ∀k > 0] = 0. Other restrictions on u arise from engineering constraints; for example the finite storage constraint requires that C¯t (u) ∈ [cmin , cmax ] for all t ≤ T . Further possibilities are explored in Section 7.2; in the meantime we assume that U (t, c, i), representing the set of all admissible strategies on the time interval [t, T ] starting in regime i and with initial inventory c, is a closed subset of U . Subject to those costs and the operational constraints, the facility manager then maximizes the net expected profit. Given initial conditions at time t, Gt = g, Ct = c, and initial operating regime i, suppose the manager chooses a particular dispatching policy u ∈ U (t, c, i). If we denote by Vu (t, g, c, i) the corresponding expected profit until final date T , then "Z # T X Vu (t, g, c, i) = E e−r(s−t) ψus (Gs , C¯s (u)) ds − e−rτk Kuτk − ,uτk Gt = g, Ct = c . t

(5)

τk

The first term in (5) counts the total revenues and costs from managing the facility up to the horizon T , and the second term counts the incurred switching costs. It remains to specify the terminal condition at T . Typical contracts specify that the facility should be returned with the same inventory C0 as initially held, and in a certain state, e.g. store. To enforce this stipulation, various buy-back provisions are employed. A common condition is ¯ 1 · (C0 − c)+ − K ¯ 2 · (C0 − c)− − K ¯ 3 1i6=0 , V (T, g, c, i; C0 ) = −K

(6)

making the penalty proportional to the difference with stipulated inventory C0 , with multipliers ¯ 1 and K ¯ 2 used for under-delivery and over-delivery respectively, and adding a second penalty of K ¯ 3 if the final regime is not store. Another common choice is V (T, g, c, i; C0 ) = −K ¯ 1 · g · (C0 − c)+ , K which penalizes for having less gas than originally and makes the penalty proportional to current price of gas. Our formal control problem is to compute M

V (t, g, c, i) =

sup

Vu (t, g, c, i),

(7)

u∈U (t,c,i)

with V (t, g, c, i; u) as defined in (5). Besides the value function V we are also interested in explicitly characterizing an optimal policy u∗ (if one exists) that achieves the supremum in (7). Before

Carmona and Ludkovski: Optimal Switching for Energy Storage

9

proceeding, let us emphasize the path-dependent nature of (7). Observe that an optimal policy u∗t0 at intermediate time t < t0 < T depends on the current inventory Ct0 . However, Ct0 = C¯t0 (u∗ ) is itself a function of past strategy (u∗s : t ≤ s ≤ t0 ). Conversely, current Ct0 affects the feasibility of future strategies {us , s ≥ t0 } through the corresponding constraints on U (t0 , c, i). The standard method of solving control problems is by dynamic programming and would proceed backwards in time, from s = T towards s = t0 . However, in our case to find the optimal action at time t0 we need to know optimal actions before t0 , hence the path-dependency and the resulting challenge. Of course, this was abstractly resolved by making Ct a state variable in (7). Unfortunately, from a numerical analysis point of view this is only a superficial fix as Ct is now neither exogenously stochastic, nor directly controlled, creating a degenerate and numerically unstable variable.

3.

Iterative Optimal Stopping.

Without inventory, (7) belongs to the class of Optimal Switching problems. These have been recently extensively studied, both analytically (Zervos 2003, Pham and Ly Vath 2007, Dayanik and Egami 2005), and numerically (Barrera-Esteve et al. 2006, Porchet et al. 2008). In particular, one can exploit the idea of the authors’ earlier paper (Carmona and Ludkovski 2008) to represent (7) as a sequence of optimal stopping (American option) problems. These sub-problems precisely capture the timing flexibility of the manager. Let St denote the set of all F-stopping times between t and T . Recursively construct the functions V k (t, g, c, i) with k = 0, 1, . . ., 0 ≤ t ≤ T , g ∈ Rd , c ∈ [cmin , cmax ] and i ∈ {−1, 0, 1} via hZ M V 0 (t, g, c, i) = E

T

i e−r(s−t) ψi (s, Gs , C¯s (c, i)) ds Gt = g , t " Z τ M k V (t, g, c, i) = sup E e−r(s−t) ψi (s, Gs , C¯s (c, i)) ds τ ∈St

t

# n o + max e−r(τ −t) −Ki,j + V k−1 (τ, Gτ , C¯τ (c, i), j) Gt = g . j6=i

(8)

Recall that C¯s (c, i) denotes the storage level at time s when policy i is constantly applied and initial inventory was c. The results in Carmona and Ludkovski (2008) show that M

Proposition 1. Let U k (t, c, i) = {u ∈ U (t, c, i) : u = (ξ1 , . . . , ξk ; τ1 , . . . , τk )} be the subset of admissible strategies with at most k switches. Then i. V k is equal to the value function for the storage problem with at most k switches allowed: V k (t, g, c, i) = supu∈U k (t,c,i) Vu (t, g, c, i).

Carmona and Ludkovski: Optimal Switching for Energy Storage

10

ii. An optimal strategy u∗ = u∗,k for V k (0, g, c, i) exists, is Markovian and is explicitly defined by τ0∗ = 0, ξ0∗ = i, and for ` = 1, . . . , k by n ( o M ∗ τ`∗ = inf s ≥ τ`−1 : V ` (s, Gs , Cs (u∗ ), i) = maxj6=i −Ki,j + V `−1 (s, Gs , Cs (u∗ ), j) ∧ T, M ξ`∗ = arg maxj6=i −Ki,j + V `−1 (τ`∗ −, Gτ`∗ − , Cτ`∗ − (u∗ ), i) .

(9)

iii. limk→∞ V k (t, g, c, i) = V (t, g, c, i) pointwise, uniformly on compacts. iv. The limit V (t, g, c, i) is continuous and is the minimal solution of the Bellman equation hZ τ e−r(s−t) ψi (Gs , C¯s (c, i)) ds V (t, g, c, i) = sup E τ ∈St

t

i + e−r(τ −t) · max −Ki,j + V (τ, Gτ , C¯τ −t (c, i), j) Gt = g . (10) j6=i

Item (i) says that V k (t, g, c, i) is the maximum expected profit to be had on the time period [t, T ] conditional on the initial state (g, c, i) and at most k switches remaining. This is useful because according to item (iii), for any > 0, there is a K large enough such that an optimal control of V K as defined in (9), generates an -optimal strategy for V . The key insight behind the proposition is the Bellman optimality principle which implies that solving the problem with at most k + 1 switching decisions allowed is equivalent to finding the first optimal decision time τ which maximizes the initial payoff until τ plus the value function at τ corresponding to optimal switching with k switches. Remark 3. We do not have full results showing the uniqueness or existence of optimal control for the original value function V , which is a delicate impulse control problem. On a practical level this makes no difference since an -optimal control is always available. Theoretically, it would be interesting to find a good set of working assumptions to ensure optimality existence/uniqueness. 3.1.

Quasi-variational formulation.

The presented storage model is a special case of stochastic impulse control problems. Hence one can apply the generic quasi-variational method developed by Bensoussan and Lions (1984). The verification theorem presented below states that a suitable smooth candidate function ϕ, which dominates the switching barrier and solves the Kolmogorov pde in the continuation region is indeed the value function of (7). The proof follows from standard techniques, see e.g. Øksendal and Sulem (2005). Proposition 2. Let LG denote the infinitesimal generator of the Markov process (Gt ). Suppose there exists ϕ(t, g, c, i) such that for M

D=

o [n (t, g, c) : ϕ(t, g, c, i) = max {−Ki,j + ϕ(t, g, c, j)} , i

j6=i

Carmona and Ludkovski: Optimal Switching for Energy Storage

11

ϕ(·, i) belongs to C 1,2,2 ([0, T ] × Rd × [cmin , cmax ]) \D ∩ C 1,1,1 (D) for each i ∈ {−1, 0, 1} and satisfies

the following quasi-variational inequality (QVI): n min ϕ(t, g, c, i) − max −Ki,j + ϕ(t, g, c, j) , j6=i o ∂ ∂ − ϕ(t, g, c, i) − L ϕ(t, g, c, i) + a (c) · ϕ(t, g, c, i) − ψ (g, c) + rϕ(t, g, c, i) = 0, G i i ∂t ∂c ¯ 1 · (C0 − c)+ − K ¯ 2 · (C0 − c)− − K ¯ 3 1{i6=0} . ϕ(T, g, c, i) = −K Then ϕ = V is the optimal value function for the storage problem (7). 2

∂ ∂ + 12 σ 2 (g) ∂g If the process (Gt ) is an Itˆ o diffusion as in (2), then LG = µ(g) ∂g 2 is a second-order

differential operator. The derived parabolic pde system with a free boundary can then be solved using standard tools, see for example (Wilmott et al. 1995, Chapter 7). In the context of gas storage this approach has been explored by Ahn et al. (2002). As the simplest choice, consider the basic finite differencing (FD) algorithm. We set up a uniform space-time grid with steps ∆t, ∆g and ∆c in the respective variables, and on this grid solve 0 = ϕt + µ(g)ϕg + σ 2 (g)/2 · ϕgg − ai (c) · ϕc + ψi (g, c) − r · ϕ, ϕ(t, g, c, i) > max −Ki,j + ϕ(t, g, c, j) , j6=i ¯ 1 · (C0 − c)+ − K ¯ 2 · (C0 − c)− − K ¯ 3 1i6=0 , ϕ(T, g, c, i) = −K

(11)

by replacing derivatives with explicit finite differences in the first equation and directly enforcing the barrier condition at each time-step. Using standard properties of the infinitesimal generator LG , one obtains the convergence ϕ(0, g, c, i) − → V (0, g, c, i) as step sizes vanish ∆t → 0, ∆g → 0, ∆c → 0. The FD method is straightforward to implement but will be slow since even in the easiest case, where (Gt ) is one-dimensional and has smooth dynamics, the pde (11) is two-dimensional in space. Furthermore, the degenerate Ct -dynamics cause numerical instability as the pde is convectiondominated (due to absence of ϕcc term). The algorithm is also not robust: for instance, adding jumps to (2) produces a partial integro-differential equation which is non-local and requires special numerical tools, see Thompson et al. (2003) for details. Similarly, pde solvers suffer from the curse of dimensionality —making (Gt ) two-dimensional is still beyond today’s computational power (in the sense of a business-time system). On the other hand, the error analysis of FD algorithms is well-studied and many improvements are possible, including adaptive solution grids, alternating direction implicit schemes, relaxation methods, etc.

4.

Numerical Method.

The benefit of the recursive formulation in (8)-(10) is its suitability for an efficient and scalable numerical implementation. In this section we describe in detail the resulting algorithms.

Carmona and Ludkovski: Optimal Switching for Energy Storage

12

To begin, we discretize time, setting S ∆ , {m∆t, m = 0, 1, . . . , M }, ∆t =

T M

as our discrete time

grid. Managerial decisions are now allowed only at τk ∈ S ∆ . This restriction is similar to looking at Bermudan options as approximation to American exercise rights. Denote by Z t+∆t M ∆ e−r(s−t) · ψi (Gs , Cs ) ds ψi (t, Gt , Ct ) = t

the total cashflows during one time-step and let t1 = m∆t, t2 = (m + 1)∆t be two generic consecutive time steps. In discrete time, the representation of V (t1 , g, c, i) in (10) reduces to deciding between immediate switch at t1 to some other regime j, which must then be maintained until t2 (i.e. τ = t1 in (10)), versus no switching and therefore maintaining regime i until t2 (τ > t1 ⇔ τ ≥ t2 ). In other words, one chooses the best (in terms of continuation value) regime j at t1 , pays the corresponding switching costs, and then waits until t2 . Thus, using the notation of (4), (8) reduces to V (t1 , Gt1 , Ct1 , i) = max −Ki,j + E ψj∆ (t1 , Gt1 , Ct1 ) + e−r∆t · V (t2 , Gt2 , C¯∆t (Ct1 , j), j)| Ft1 . (12) j

The dynamic programming method can now be applied to recursively evaluate (12) backwards in time to obtain the discretized Snell envelopes of the optimal stopping problem (10). Hence, for a numeric evaluation of V it is sufficient to construct an algorithm for evaluating the conditional expectations appearing in (12). 2 Let (Bj (g; t1 , c, i))∞ j=1 be a given orthonormal basis of L (Ft1 ) (selection of (Bj ) is discussed

below). Recall that (Gt ) is Markov, while (Ct ) is determined by u. Thus, we may view the conditional expectation in (12) as a map h i M g 7→ E(g; t1 , c, i) = E ψi∆ (t1 , Gt1 , c) + e−r∆t · V (t2 , Gt2 , C¯∆t (c, i), i) Gt1 = g .

(13)

Nb

The latter may be approximated with a projection on the truncated basis (Bj )j=11 : E(g; t1 , c, i) =

∞ X j=1

N

ˆ t1 , c, i) = αj Bj (g; t1 , c, i) ' E(g;

b1 X

αj Bj (g; t1 , c, i),

(14)

j=1

where αj are the R-valued projection coefficients. The right hand side of (14) is a finite-dimensional projection of the continuation values onto the basis functions and can be replaced with an empirical regression based on a Monte Carlo simulation. This then gives a method for implementing (12) on a computer. n We begin by generating N sample paths (gm∆t ), n = 1, . . . , N of the discretized (Gt ) process with

a fixed initial condition G0 = g = g0n . As mentioned before, the inventory depends on the policy choice, so it cannot be directly simulated. To overcome this problem, we shall construct a grid in the C-variable and compute V (t, g, c, i) only for c ∈ {c0 = cmin , c1 , . . . , cNC = cmax }.

Carmona and Ludkovski: Optimal Switching for Energy Storage

13

We will approximate the value function by the empirical average of the pathwise quasi-values PN (from now on simply values) V (0, g, c, i) ' N1 n=1 v(0, g0n , c, i). The values v(t, gtn , c, i) are computed ¯1 · recursively in a backward fashion, starting with the terminal condition of (6): v(T, gTn , c, i) = −K ¯ 2 · (C0 − c)− − K ¯ 3 1i6=0 . Consider again two consecutive time steps t1 , t2 and suppose (C0 − c)+ − K inductively that we know v(t2 , gtn2 , c, i) along the paths (gtn2 )N n=1 and for c ∈ {c` , ` = 1, . . . , NC }. Our n ˆ m∆t goal is to compute v(t1 , gtn1 , c, i). To obtain the prediction E(g ; t1 , c, i) of the continuation value,

one first computes v(t2 , gtn2 , C¯∆t (c, i), i). Note that in general C¯∆t (c, i) does not belong to the grid {c` }, and interpolation is needed. Then one regresses v(t2 , gtn2 , C¯∆t (c, i), i) against the basis functions Nb

(Bj (gtn1 ; t1 , c, i))j=11 to find the corresponding αj ≡ αj (t1 , c, i) and applies (14). By analogue of (12), the estimate for v(t1 , gtn1 , c, i) is then ˆ tn ; t1 , c, i) ∨ max −Ki,j + E(g ˆ tn ; t1 , c, j) . v(t1 , gtn1 , c, i) = E(g 1 1 j6=i

(15)

ˆ The Observe that (15) performs pathwise computations, while using across-the-paths projection E. scheme (15) first appeared in Tsitsiklis and van Roy (2001) in the context of American option pricing. In our setting we call it a mixed interpolation-Tsitsiklis-van Roy scheme (MITvR). It is also useful to think in terms of the optimal storage strategy. Let ˆn (t1 ; i) ∈ {−1, 0, 1} represent the optimal decision on the n-th path at time t = t1 and current regime i. The analogue of (12) implies that (recall Ki,i ≡ 0) ˆ tn ; t1 , c, j) . ˆn (t1 ; i) = arg max −Ki,j + E(g 1 j

(16)

Thus, the set of paths on which it is optimal to switch at time t = m∆t is given by {n : ˆn (m∆t; i) 6= i}. Equation (16) can be used to construct the switching boundaries which partition [0, ∞] × [cmin , cmax ] into regions where injection, withdrawal, etc., is optimal at date t. The efficiency of (15) is enhanced by using the same set of paths to compute all the conditional expectations. Nevertheless, because of the capacity variable C the above approach is still timeintensive. Indeed, at every time step m∆t and regime i, we must run a separate regression for each inventory grid point c` . Hence, in terms of computational complexity, the above method is equivalent to solving N c optimal switching problems. The choice of appropriate basis functions (Bj (·; t, c, i)) in (14) is user-defined. A detailed analysis of different orthogonal families is available in Stentoft (2004). Empirically, basis choice has only a mild effect on numerical precision, but strongly affects the variance of the algorithm. Thus, customization is desirable and it helps to use basis functions that resemble the expected shape of the value function. In practice, Nb1 as small as five or six normally suffices, and having more bases

Carmona and Ludkovski: Optimal Switching for Energy Storage

14

can often lead to worse numerical results due to overfitting. Let us mention that the requirement of an orthonormal basis is purely theoretical and any set of linearly independent functions will suffice. Some of our favorite choices are exponential functions eαg and the polynomials g m ; this choice is essentially heuristic. Also, we typically select the bases independent of parameters (t, c, i), though the latter offer a wide scope for additional fine-tuning. 4.1.

Quasi-Simulation of Inventory Levels.

To maintain numerical efficiency it is desirable to avoid the fixed discretization in the C-variable that resembles the slow lattice schemes. Accordingly, we propose the following alternative that uses pathwise and regime-dependent inventory levels (cnm∆t (i)). The idea is to perform a bivariate regression in (14) of tomorrow’s value against the (price, inventory) pair. The paths (cnm∆t (i))M m=1 are generated backwards during the dynamic programming procedure by combining randomization and guesses of today’s optimal strategy. Besides added efficiency, we are also guided by considerations of accuracy. Quasi-simulation of inventory allows us to use the Longstaff and Schwartz (2001) scheme of computing pathwise value functions of optimal stopping problems. From simpler problems of American option pricing and plain optimal switching we know that the LSM scheme typically has less bias (though more variance) than the TvR scheme (15) (Ludkovski 2005). We inductively assume again that we are given the 3N values v(t2 , gtn2 , cnt2 (i), i), i ∈ {−1, 0, 1}, n = Nb

¯j (g, c; t1 , i))j=12 . For a given path n, regime i, and a 1, . . . , N , as well as bivariate basis functions (B given inventory cnt1 (i) (see below about obtaining cnt1 (i)) we make the optimal switching decision as follows (compare with (15)): 1. For each k ∈ {−1, 0, 1}, regress {e−r∆t · v(t2 , gtn2 , cnt2 (k), k)}N n=1 against the basis functions ¯j (gtn , cnt (k); t1 , k), j = 1, . . . , Nb ). This gives a prediction (B 2 1 2 N

˜ : (g, c, k) 7→ E

b2 X

i h ∆ −r∆t ¯ ¯ α ¯ j Bj (g, c; t1 , k) ' E ψk (t1 , g, C−∆t (c, k)) + e · v(t2 , Gt2 , c, k) Gm∆t = g (17)

j=1

of the value tomorrow given today’s prices and tomorrow’s inventory. 2. Compute C¯∆t (cnt1 (i), j), the inventory tomorrow given today’s inventory cnt1 (i) and the decision to switch to j. 3. The optimal decision is the regime ˆn (t1 , i) maximizing the approximate continuation value, cf. (12), o n ˜ tn , C¯∆t cnt (i), j , j) − Ki,j . ˆn (t1 , i) = arg max E(g 1 1 j

(18)

Carmona and Ludkovski: Optimal Switching for Energy Storage

15

4. If C¯∆t (cnt1 (i), ˆn ) = cnt2 (ˆ n ), then the Longstaff-Schwartz update is used: v(t, gtn1 , cnt1 (i), i)

Z

t2

= t1

e−r(t−t1 ) ψˆn (Gt , C¯t−t1 (cnt1 (i), ˆn )) dt + e−r∆t · v(t2 , gtn2 , cnt2 (ˆn ), ˆn ) − Ki,ˆ . (LSM)

Else, one updates via ˜ tn , C¯∆t (cnt (i), ˆ), ˆ) − Ki,ˆ . v(t1 , gtn1 , cnt1 (i), i) = E(g 1 1

(TvR)

The first case (LSM) stands for Least Squares Monte Carlo. In that version the across-the-paths regression is used primarily to make the optimal switching decision, but is not necessarily fed into the pathwise values. This helps to eliminate potential biases from the regressions by preventing error accumulation across time-steps. To preserve this beneficial look-ahead property of the Longstaff and Schwartz (2001) algorithm, we attempt to speculatively pick cnm∆t (i) such that the first case (LSM) occurs as much as possible. In other words, as we move back in time, we try to select inventory levels that form an optimal (price, inventory)-path on the remaining time interval. When this is not possible (due to capacity or other constraints, or if our guess of ˆ is incorrect), we fall back onto the basic (TvR) scheme. Accurate guessing of ˆ means that we correctly select the optimal strategy (up to the errors resulting from the projection). In such a case, we have RT P v(t1 , gtn1 , cnt1 (i), i) = t1 e−r(s−t) ψu∗s (Gt , C¯s (u∗ )) ds − τ ∗

k

k

path. Observe also that the method can be used iteratively over several simulation runs, improving the guesses of ˆ over time. The terminal inventory levels cnT are randomized and obtained by independent and uniform samples from [cmin , cmax ]. At each step t = m∆t, some randomization in (cnt1 (i)) is also desirable in order to avoid clustering and allow for good fit during the regression step. Of course, randomization reduces the number of paths satisfying (LSM), and balancing the two objectives is a detail that we leave to implementation. We christen this scheme Bivariate Least Squares Monte Carlo (BLSM). 4.2.

Algorithm Summary

¯j ) and algorithm 1. Select a set of univariate basis functions (Bj ), bivariate basis functions (B parameters ∆t, M, N, Nb1 , Nb2 . n 2. Generate N paths of the price process: {gm∆t , m = 0, 1, . . . , M, n = 1, 2, . . . , N } with fixed

initial condition g0n = g0 . Generate a random terminal inventory level cnT (i) for each path and each regime i. 3. Initialize the pathwise values v(T, gTn , cnT (i), i) from (6).

Carmona and Ludkovski: Optimal Switching for Energy Storage

16

4. Moving backward in time with t = m∆t, m = M, . . . , 0 repeat the Loop, where the computations are based on (10): i) Guess Current Inventory: generate (cnm∆t (i)) by guessing the optimal decision ˆn (m∆t, i) and solving cnm∆t (i) = C¯−∆t (cn(m+1)∆t (ˆ n (m∆t, i), ˆn (m∆t, i)). ii) Regression Step: do the bivariate regression of (17). iii) Optimal Decision Step: find the optimal decision using (18). n iv) Update Step: compute v(m∆t, gm∆t , cnm∆t (i), i) via (LSM) and (TvR).

v) Switching Sets: the points n Cm∆t (i, j) , {(gm∆t , cnm∆t ) : n is such that ˆn (m∆t, i) = i}

(19)

define the empirical region in the (G, C)-space where switching from regime i to regime j is optimal. This defines the optimal strategy at t = m∆t. 5. end Loop 6. Interpolate V (0, g0 , c, i) from the N values v(0, g0n , cn0 (i), i) for the desired inventory level c (using splines, kernel regression, etc.). Remark 4. As mentioned before, in the discrete-time version we allow switching costs to be zero, Ki,j ≡ 0. In that case V (t, g, c, i) does not depend on the current regime i and so one can save on the corresponding computations. 4.3.

Algorithm Complexity.

The BLSM algorithm requires O (N · M · ((Nb1 )3 + (Nb2 )3 )) operations. The most computationally intensive operation is the regression step where we face matrices of size N × Nb1 and N × Nb2 , and which make the algorithm linear in the larger dimension N and cubic in the smaller dimensions Nb1 , Nb2 . In contrast, the algorithm complexity of the MITvR scheme is O (N · M · Nc · (Nb1 )3 ) where Nc is the grid size in the C-variable. However, these expressions hide the relationship between Nb1 , Nb2 and N , because more basis functions require more paths for accurate evaluation of the regression step. In fact, according to Glasserman and Yu (2004), N must be asymptotically exponential in the number of basis functions. On the other hand, to perform the bivariate regression (17), it is likely that a large number of basis functions Nb2 is needed, about 12 − 15 in our experience. Hence in the BLSM algorithm N must be taken larger than in the the MITvR case. Precise comparison is hard because the BLSM scheme inherently generates more variance and we have no hard benchmark to go by. For our examples we find that N = 16, 000 for the MITvR scheme is reasonable, while N = 40, 000 is needed for BLSM. Practically speaking this implies that BLSM is about twice as fast as MITvR, see Section 5. The memory requirements of both schemes are n O(N · M ) corresponding to the need to store the entire sample paths (gm∆t )N n=1 in memory.

Carmona and Ludkovski: Optimal Switching for Energy Storage

4.4.

17

Convergence.

The presented algorithms have several layers of approximations. Three major types of errors can be identified: error due to time discretization and the corresponding restriction of strategies to U ∆ , projection error and Monte Carlo sampling error. Detailed error analysis has been performed

in Carmona and Ludkovski (2008) for the case of the TvR scheme with no inventory. Taking Ct to be a ‘dummy’ variable determined by the dynamics of Gt and policy u, the results carry over without change. Analysis of the BLSM scheme is too involved, however see partial results in this direction in Egloff (2005). This lack of provable convergence results is typical for Monte Carlo optimal stopping methods, largely due to the nonlinearity introduced by the stopping boundaries. Nevertheless, extensive empirical experiments (Stentoft 2004) have strongly supported the general TvR/LSM methodology. The error from discretizing (Gt ) and simultaneously restricting the switching times to occur √ only at the discrete time grid points is known to be O( ∆t). The error from approximating ˆ k)) when the conditional expectations with a projection (14) is on the order of O(∆t−k · (kE − E computing V k (Carmona and Ludkovski 2008). This suggests that the projection errors multiply in the number of decisions taken k. However, empirically the dependence on ∆t is much better, especially under the BLSM scheme, so this upper bound is probably not tight. In any case for practical examples, the typical number of switches is in single digits. Finally, the third source of error is due to approximating the projections with an empirical regression using N realizations of n , n = 1, . . . , N ). This error is difficult to analyze due to interactions between the the paths (gm∆t

path-by-path maximum taken in (16) and the across-the-paths regression. No convergence behavior is known; however numerical experiments suggest that it is close to O((∆t · N )−1/2 ), which is the expected rate for Monte Carlo methods. Table 1 illustrates this conjecture on Example 2 below. We run the BLSM algorithm using 8000 − 40000 Monte Carlo paths and tabulate the resulting standard errors. At least in this case we see in fact a faster than O(N −1/2 ) convergence. Some residual variance in V remains even for large N due to the sensitivity of the algorithm n to outliers in generated paths of gm∆t for small m. Since all paths begin from the same

initial value g0n , the regression step is sensitive in the first few time steps.

5.

Numerical Results.

In this section we present several examples to show the structure of the storage problem and the scalability of our algorithm.

Carmona and Ludkovski: Optimal Switching for Energy Storage

18 Table 1

Convergence of Monte Carlo error for Example 1 under the BLSM scheme. The initial gas price is G0 = 3 $/MMBtu, initial inventory is C0 = 4 Bcf, initial regime is i0 = store and the horizon is T = 1 year with ∆t = 0.005. The table shows the mean and standard deviation of V (0, G0 , C0 , i0 ) (in MM$/MMBtu) obtained by running the algorithm 50 times.

Paths N 8000 16000 24000 32000 40000 Table 2

Mean 9.63 9.41 9.37 9.38 9.35

Std. Dev 0.4179 0.1362 0.0961 0.0663 0.0647

Comparison of numerical results for Example 1. The initial gas price is G0 = 3 $/MMBtu, initial inventory is C0 = 4 Bcf, initial regime is i0 = store and the horizon is T = 1 year. The table shows V (0, G0 , C0 , i0 ) (in MM$/MMBtu) computed using four different methods. Standard deviations for the Monte Carlo methods were obtained by running the respective algorithms 50 times.

Method Coarse FD Fine FD MITvR BLSM

V (0, G0 , C0 , i0 ) 9.32 9.44 9.86 9.35

Std. Dev – – 0.021 0.067

Time (min) 24 65 47 32

Example 1. As a first illustration of our approach, consider a facility with a total capacity of cmax = 8 Bcf rented out for one year, T = 1. The price process is taken from the data of de Jong and Walet (2003), d log Gt = 17.1 · (log 3 − log Gt ) dt + 1.33 dWt .

(20)

Observe the very fast mean-reversion κ = 17.1 of the prices, with a half-life of 15 days. The initial inventory is 4Bcf and the terminal condition is V (T, g, c, i) = −2 · g · max(4 − c, 0). Thus, the manager is penalized at double the market price for final inventory being less than 4 Bcf and receives no compensation for any excess. The other parameters (in yearly units) in (1) are ain (c) = 0.06 · 365, Ki (c) ≡ 0.1c, Ki,j ≡ 0.25 for i 6= j, aout (c) = 0.25 · 365, r = 0.06, bi (c) ≡ ai (c). Thus, it takes about 8/0.06 = 133 days to fill the facility and 8/0.25 = 32 days to empty it. In this simple example we have taken the injection/withdrawal rates to be independent of inventory levels and assumed lossless pumping. We solve this storage problem using three different solvers: an explicit finite-difference pde solver discretizing (11), the MITvR scheme of (15) and the BLSM scheme of (LSM). The results are

Carmona and Ludkovski: Optimal Switching for Energy Storage

Figure 1

19

Value function surface for Example 1 showing V (0.5, g, c, store; T = 1) as a function of current gas price Gt = g and current inventory Ct = c.

summarized in Table 2. As an extra check we used two different grid sizes for the pde solver: a coarse 250 × 250 (G, C)-grid with 10000 time steps and a finer 500 × 500 (G, C)-grid with 20000 time steps. The MITvR scheme used 200 time steps, 10000 paths with six basis functions and 80 grid points in the C-variable. The quasi-simulation BLSM scheme used 200 time steps and 40,000 paths with fifteen basis functions. Taking the fine pde solver as the benchmark value, we see that the simulation methods are within 5% of the optimal value and seem to have an upper bias. In particular, the BLSM algorithm with 40,000 paths produces a storage value that is within 1% of the fine FD value, giving an acceptable level of accuracy. The computational challenges involved are indicated by the long running times of the algorithms.4 . In this light, the 45% time savings obtained by the joint (G, C)-regression are significant from a practical point of view. Figure 1 shows the value function V (t, g, c, i) as a function of current price and inventory for an intermediate time t = 0.5 and i = store regime. Not surprisingly, higher inventory increases the value function since one has the opportunity to simply sell the extra gas on the market. In the Gt -variable we observe a parabolic shape with a minimum around the long-term mean 3$/MMBtu. Thus, deviations of Gt from its mean imply higher future profits, confirming our intuition about storage acting as a financial straddle. 4

The simulation methods were run in Matlab on a 1.6GHz desktop. The pde solver was written in C++ and run on the same machine.

Carmona and Ludkovski: Optimal Switching for Energy Storage

20

We next investigate the sensitivities of the value function. Table 3 shows the effect of storage flexibility on V . Higher transmission rates increase the extrinsic value of storage, since the manager can move more gas in and out of the facility under “favorable” circumstances. In the example considered, the smaller injection rate acts as a bottleneck on the manager’s flexibility, so the derived extrinsic value is more sensitive to ainj than to aout . Table 4 studies the effect of other parameters on the extrinsic value. We find that switching costs Ki,j have a major impact on V . High Ki,j makes the manager risk-averse and unwilling to change the pumping regime until a very good opportunity comes along (since each switch has a high upfront fixed cost). We also find that because of the limited transmission rates and the mean-reverting nature of the prices, there are dis-economies of scale with respect to facility size. Thus, cutting the facility size to 6Bcf reduces value by nearly 16%, but an increase from 10Bcf to 12Bcf produces a benefit of just 3%. The last panel of Table 4 studies the impact of price mean-reversion on value of storage. Intuition suggests that stronger mean-reversion would reduce price fluctuations, and therefore also the opportunities to capitalize on temporal price spreads. This is illustrated in Table 4 which shows that increasing the mean-reversion strength κ in (20) decreases storage value. Table 3

Effect of Storage Flexibility on the Value Function in Example 1. Extrinsic value corresponds to V (0, 3, 4, 0). Results obtained using the BLSM algorithm with N = 40000 paths and ∆t = 0.005.

ain Bcf/day 0.06 0.03 0.12 0.18 0.12

aout Bcf/day

Extrinsic Value

0.25 0.125 0.50 0.75 0.25

9.35 4.75 14.50 17.28 12.33

Example 2. Our second example is based on the situation presented in Thompson et al. (2003). A mean-reverting model is taken for gas prices, with a seasonally-adjusted mean-reverting level. The gas prices satisfy dGt = 4 · (6 + sin(4πt) − Gt ) dt + 0.5Gt dWt . Thus, the mean-reversion level has a seasonal component representing the summer/winter price increases in the North American markets. This seasonality implies an approximate trough-to-peak calendar spread of $1 for each half year.

Carmona and Ludkovski: Optimal Switching for Energy Storage

Table 4

21

Effect of Engineering Characteristics on the Value Function in Example 1. Results obtained using the BLSM algorithm with 40,000 paths.

Effect of Facility Size Capacity (Bcf) V (0, 3, 4, 0) 6 7.78 8 9.35 10 10.26 12 10.58 Effect of Switching Costs Ki,j , i 6= j V (0, 3, 4, 0) 0.01 13.25 0.1 11.40 0.25 9.35 0.5 6.73 Effect of Storage Costs Ki V (0, 3, 4, 0) 0 9.77 0.05 · c 9.56 0.1 · c 9.35 Effect of Mean Reversion Strength κ V (0, 3, 4, 0) 7.1 10.47 12.1 10.32 17.1 9.35 22.1 8.73 Secondly, the injection and withdrawal rates are ratcheted in terms of current inventory. Thus, injection rate decreases as the amount of gas in the facility grows, and conversely withdrawal rate decreases as less gas is on inventory. More precisely, the facility capacity is cmax = 2 Bcf and q √ 1 1 the yearly rates are aout (c) = 0.177 · 365 c, ain (c) = 0.0632 · 365 c+0.5 − 2.5 . These rate functions are related to the ideal gas law. Namely, gas transmission rate is proportional to pressure in the reservoir, which in turn is inversely quadratically related to gas volume. Thus, maximum injectivity at cmin = 0 is 0.08Bcf/day, and maximum withdrawal at cmax = 2 is 0.25Bcf/day. The monetary reward functions are given by ψ−1 (g, c) = −(ain (c) + 0.6205)g, dCt = ain (Ct ) dt, Inject: Store: ψ0 = 0, dCt = 0 dt, Withdraw: ψ1 (g, c) = aout (c)g, dCt = −aout (Ct ) dt.

(21)

Note that there is gas loss during injection represented by the term 0.6205g. We take a horizon of one year T = 1 with no terminal penalty, V (T, g, c, i) = 0. There are no switching costs and r = 0.1. Figure 2 presents the optimal control for Example 2 at different times to maturity. The three shades indicate the regions of injection (on the left), storage (dark region in the middle) and

Carmona and Ludkovski: Optimal Switching for Energy Storage

22

withdrawal (on the right) respectively. We observe that the optimal action regions Cm∆t (i, j) of (19) are strongly dependent on time to maturity t. Namely, as the terminal horizon approaches, the no-action regions widen, as the fixed switching costs begin to dominate the potential benefit from changing regimes. Also, the terminal condition starts to influence the manager’s optimal decision. Figure 2 also shows that Cm∆t (i, j) depends on i, since the positive switching costs induce inertia and make the manager reluctant to change policies. Finally, the ratcheting transmission rates cause the curvatures in Cm∆t (i, j) as a function of current inventory c.

Switching regions C0.25 (0, ·) from regime store at 3 months to maturity.

Figure 2

Switching regions C0.5 (1, ·) from regime withdraw at 6 months to maturity.

Optimal Controls for Example 2 using the BLSM algorithm with 32,000 paths. Each point corresponds to a simulated (gtn , cn ˆ of (18) (red for inject, black for t (i)) pair, and the color indicates the optimal store, blue for withdraw).

Example 3. Finally, in our third example we illustrate the flexibility of the simulation approach with regards to more complex price processes. It is well known that a one-factor diffusive model does not provide a good fit to gas prices. Accordingly, let us consider a two-factor model with jumps; namely a log-mean-reverting diffusive factor and a second mean-reverting pure jump factor. The second factor captures spikes in natural gas prices without making the mean-reversion rate unnecessarily high (Kluge (2004)): ( 1 dGt = 4(log 6 − log G1t )G1t dt + 0.5G1t dWt , dG2t = 26(0 − G2t ) dt + ξt dNt − λµJ dt,

(22)

where (Nt ) is an independent Poisson process with intensity λ, and the jump size ξt ∼ N (µJ , σJ ) has normal distribution. The total gas price is the product Gt , G1t · exp(G2t ), and G2 can be interpreted as the multiplicative jump factor that causes price spikes on the scale of µJ %. The possibility of price spikes makes storage much more valuable since it increases the volatility of

Carmona and Ludkovski: Optimal Switching for Energy Storage

23

inter-temporal spreads. We pick λ = 12, µJ = 0.02, σJ = 0.1 for the jump component, as well as T = 2, r = 0.06, Ki,j = Ki ≡ 0, V (T, ·) = 0. Implementing Example 3 requires only minor modifications to the implementation of Example 2, which essentially reduce to simulation of the bivariate price process (G1t , G2t ) and selection of ¯j (g 1 , g 2 ). These changes only take a few minutes to make, and the bivariate basis functions B resulting algorithm will take only a little longer to run (depending on how many basis functions are added to deal with (G2t )). In contrast, with a pde approach, the new state dimension would require an extensive rewrite of the code and would slow the performance by an order of magnitude. Since the value function V (t, g1 , g2 , c, i) now has three space variables, in Figure 3 we visualize the dependence of V just on the two price factors (g1 , g2 ) for different inventory levels c. Since each factor is mean-reverting and the total price is a product of the two, the value function will exhibit a parabolic straddle shape in each factor. Thus, in Figure 3, when g2 < 0 one can expect prices to rise back to their “normal” level, and so this presents an opportunity for injection, at least when inventory is low. Conversely, when g2 > 0 and inventory is high, we are in an upward spike with prices expected to fall and an attractive withdrawal opportunity. The dependence of V on g1 is similar to that of Figure 1. Overall, as a function of the price and the spike factor, the value function exhibits an asymmetric “bowl” shape, which in turn is dependent on the current level of inventory.

c = 0.4Bcf .

Figure 3

6.

c = 1Bcf .

c = 1.6Bcf .

Value function V (1, (g 1 , g 2 ), c, −1) for Example 3 for different inventory levels c.

Hydroelectric Pumped Storage.

Another important practical application of our model is hydroelectric pumped storage. Pumped storage consists of a large reservoir of water held by a hydroelectric dam at a higher elevation.

24

Carmona and Ludkovski: Optimal Switching for Energy Storage

When desired, the dam can be opened which activates the turbines and moves the water to another, lower reservoir. The generated electricity is sold to the power grid. As the water flows, the upper reservoir is depleted. Conversely, in times of low electricity demand (weekends, etc.), the water can be pumped back into the upper reservoir using special electricity-operated pumps (with the required energy purchased from the grid). The efficiency of the system is about 80%, and the capacity of such pumped storage facilities is typically on the order of several hundred megawatthours (MWh). Currently pumped storage is the dominant type of electricity storage with more than a hundred facilities around the world (ASCE 1996). The major use of pumped storage is to capture the daily off-peak/on-peak spread in electricity prices. However, by adjusting pump rates, the reservoir manager can also take advantage of seasonal price/water level fluctuations like in the gas storage problem, see Bl¨ ochlinger et al. (2004). Beyond direct losses from upstream pumping, stored water is subject to evaporation. At the same time, precipitation and/or upper river run-off provide reservoir replenishment. Realistic modeling is complicated by the need to compute the potential energies of the reservoirs which depend on the relative levels of the water and in turn modify generation rates ai (Ct ). We abstract from these concerns and treat the problem in our framework of commodity storage (5), with an addition of another variable modeling weather. Let Lt represent the Markovian weather state at time t (e.g. Lt can be a humidity index or river flow rate). (Lt ) controls reservoir gains/losses, so that inventory depletes at rate d(Lt , Ct ) irrespective of the storage regime. The inventory Ct represents water level ¯ > 1, bin = Ka ¯ in , in the upper reservoir5 . The pumping inefficiency is represented by a multiplier K bout = aout that affects the cost of pumping. The overall model is thus: ¯ · g · ain − K−1 (c), ψ−1 (g, c) = −K dCt = [ain − d(Lt , Ct )] dt, Pump: Store: ψ0 (g, c) = −K0 (c), dCt = [a0 − d(Lt , Ct )] dt, Generate: ψ1 (g, c) = +g · aout − K1 (c), dCt = [−aout − d(Lt , Ct )] dt.

(23)

Once a suitable model is chosen for (Lt ) (see e.g. Cao et al. (2004)), implementation would be similar to Example 2 above, and would require minimal changes to the simulation algorithm. Note that one could mix-and-match different model types for different variables, for instance a jump-diffusion model for gas prices, and a seasonal AR(1) model for river run-off. Such flexibility would be hard to achieve outside of simulation paradigm (compare to the stochastic programming approach of Nowak and R¨ omisch (2000), Doege et al. (2006)). 5

A full model should also specify the lower reservoir inventory since the latter also depletes over time.

Carmona and Ludkovski: Optimal Switching for Energy Storage

7.

25

Extensions.

In this section we discuss various extensions and modifications that can be made to our model. First, let us remark that many other resource management problems can be recast in our framework. Such problems all feature fluctuating commodity prices, finite inventory constraints and a small number of operating regimes describing the facility state. Below we elaborate on some of the possibilities. 7.1.

Other Applications.

7.1.1.

Mine management A producer extracts metal from a mine with initial capacity C0 .

As the resource is mined, the inventory Ct declines. In the meantime, the producer can control the mine operating regime to time the extraction with high commodity prices represented by (Gt ). In this situation, the remaining resource amount Ct is non-decreasing, since only extraction is possible. Exhaustion implies that no profit is available when Ct = 0: V (t, g, 0, i) = 0. Armed with our methodology we can e.g. redo in a more efficient manner (see Ludkovski (2005) for the computation) the copper mine example analyzed in Brennan and Schwartz (1985). Moreover, we can easily add further constraints to their model. A related application is production of oil from oilfields. In the latter context, Ct can be increased by further exploration; at the same time fixed extraction costs increase as the field gets depleted, so that Ki ≡ Ki (Ct ). One may also add a termination option that allows total shutdown and avoids the ongoing O&M costs Ki . 7.1.2.

Hydroelectric Power Generation This setting is similar to the pumped storage

problem; however no pumping is available and the dammed reservoir is replenished solely with river run-off. The latter is modelled by a stochastic process (Lt ). When the turbines are running the produced electricity is sold at the spot power price Gt . As before, the inventory Ct is the current amount of water in the dammed reservoir. Reservoir management has been already studied by ancient Egyptians and Mesopotamians; related stochastic control models have recently been considered by Keppo (2002) and McNickle et al. (2004). Note that on a practical level a major challenge is the long-memory hydrological features of (Lt ). 7.1.3.

Power Supply Guarantees Yet another possibility similar to the pumped storage

above is the case of power supply guarantees. The latter involve a hybrid energy storage/power generation setting. By law, the North American Load-Serving Entities (LSE), i.e. the local power utilities, are obligated to provide power irrespective of demand. The latter is stochastic so that the LSE faces uncertain demand (volume risk), coupled with uncertain fuel prices (price risk). To insulate against risk, the LSE might operate an energy storage facility (e.g. a natural gas aquifer)

Carmona and Ludkovski: Optimal Switching for Energy Storage

26

that acts as a buffer between risky supply costs and risky demand needs. Letting (Dt ) represent the demand at time t, one obtains a model similar to (23) where the reservoir is depleted at rate Dt due to the requirement of producing fuel. Note that typically (Dt ) is highly correlated with the fuel price (Gt ) as high demand drives up the spot prices. Thus, the marginal cost of storage is high precisely when prices are high. See Deng et al. (2006) for further details. 7.1.4.

Emissions Trading Another application area is emissions trading, cf. Insley (2003).

A firm running a factory is subject to emission laws and must account for its pollution by buying publicly traded emission permits with current price Gt . The non-increasing inventory Ct in this case corresponds to the total number of remaining factory orders that must filled in the current quarter. Hence, the management must satisfy all the orders while minimizing emission costs. The firm opportunistically runs its production given emission price Gt and current shipment backlog. This setup is similar to supply guarantees, with an additional constraint of Ct ≤ C0 − Ot where Ot is the (deterministic) shipping timetable supplied by the customers. Violation of this constraint causes a severe penalty as the firm misses its shipment. Many other situations can be imagined—forest management, oilfield development, pipeline shipping, etc. From the above descriptions it should be evident that our numerical algorithm would carry over easily to the new settings. Summarizing, optimal switching with inventory is a widespread financial setting with many practical applications. 7.2.

Incorporating Other Features.

From a practical standpoint the presented models are gross simplifications. However, as already advertised, the simulation framework permits great flexibility. To illustrate the possibilities, we briefly discuss additional features that a practitioner is likely to implement. First of all, one is likely to use a more general price model for (Gt ) than (2). As already mentioned, all that is absolutely necessary is to satisfy Assumption 1, thus extra features such as regime-switching or latent factors are easily implementable. As already shown in Examples 2 and 3, seasonality effects, models with jumps and multi-factor models can be implemented straightforwardly. The dynamics of (Gt ) might also be affected by the choice of strategy. Indeed, since the manager tends to buy when prices are low and sell when prices are high, her influence is to smooth out the price fluctuations of (Gt ). This effect can be quite pronounced in segmented markets based on supply and demand (e.g. gas markets with little outside connectivity). As long as the effect is limited to the coefficients of (2), µ = µ(ut , ·), σ = σ(ut , ·), such market impact can be treated by independently simulating price paths under each separate regime and then modifying the LSM

Carmona and Ludkovski: Optimal Switching for Energy Storage

27

algorithm like in Carmona and Ludkovski (2008). If the transmission losses and engineering costs Ki are nonlinear in the pumping rates, it may become optimal to inject/withdraw gas at submaximal rates. In such a case, the optimal control u∗ will take on a continuum of values. Our method relies heavily on u∗ belonging to a (small) finite number of regimes; however a reasonable first-order correction would be to add a few more regimes (i.e. to discretize the range of u∗ ) to the original three considered here. Another challenge is proper modeling of borrowing constraints faced by the manager. In the North American gas industry, the facility typically borrows money in the summer to inject gas and then repays its loans in the winter as gas is withdrawn. In the meantime, the creditors often impose margin requirements regarding the value of stored gas versus the original loan. Thus, if prices drop, the manager might receive a margin call that would require him to sell off some of the inventory (at a loss) to raise capital. To account for this, one can let Bt be the cumulative borrowed capital taken out for inventory purchase. The margin constraints are then imposed on the net equity Bt − Ct · Gt ; alternatively some absolute borrowing limits Bt > −K could be required. The option of forward sales is another crucial feature. Forward sales allow the manager to lock-in future profits while reducing earnings volatility and form a bread-and-butter business in the gas storage industry. From the modeling perspective, a forward contract is challenging due to its nonMarkovian nature, which necessitates complex account-keeping for gas already sold but still sitting in the facility (or gas already bought but still not on inventory). Indeed, imagine that at time t = t0 , the manager forward-sells some quantity C0 of gas at future date t = t1 . This now affects her possible future strategies: the manager must have at least C0 in inventory at t = t1 , and must also start to withdraw C0 after t1 . Such constraints are modelled easily enough in our framework; however because they affect the future, they are not easily implementable in the dynamic programming method —to find the value of the forward sale at t = t0 we must recompute the optimal strategy ¯ (t, c, i), which is computationally challenging (essentially under the new admissibility restrictions U requiring as much effort as the original computation). Hence, modeling of forward sales would lead to a “tree” of simulations with a separate branch for each possible forward sale or purchase. This might still be practical to do if the forward sales occur infrequently. Finally, an interesting research direction would be incorporation of realistic risk objectives for the manager. In this paper we have assumed that the manager maximizes total (discounted) earnings from the asset. In practice this would lead to overly aggressive strategies and high earnings volatility. Thus, it is desirable to impose risk constraints that lead to more conservative speculation. One method for doing so can be achieved by replacing the linear conditional expectations in (12) with

Carmona and Ludkovski: Optimal Switching for Energy Storage

28

non-linear expectations that take into account risk preferences Musiela and Zariphopoulou (2004). Alternatively, one may penalize the variance in cumulative gains/losses which were denoted by BT above.

8.

Conclusion.

This paper presented a simple model for energy storage that emphasizes the intertemporal optionality of the asset. Assuming that the commodity is bought and sold on the spot market, we have maximized the expected profit given operational constraints, in particular inventory limits and switching costs. While the model sidesteps the possibilities of forward trades, it properly accounts for the dynamic nature of the problem, which is a crucial aspect of revenue maximization. Our approach is scalable and robust and we provide a detailed description of implementation. As our numerical examples attest, the model is computationally efficient and we believe better than any other proposed in the literature. We hope it can fill in the gap between current practitioner needs and academic models. Moreover, our strategy is applicable to many related problems, such as hydroelectric pumped storage, power supply guarantees, natural resource management and emissions trading. As the next step in improving our model, one should study more advanced risk objectives. In our model we have assumed linear risk-preferences of the agent who maximizes expected revenue under a given pricing measure. A real-life storage manager is in fact risk-averse, and would therefore take additional steps to hedge her risks (by e.g. reducing the variance of cashflows). However, direct hedging is likely to be of limited use in practice, because the facility buys gas based on local prices that are different from the liquidly traded indices. For instance, for gas storage the liquid instruments are the Henry Hub contracts available on the New York Mercantile Exchange (NYMEX). Such NYMEX trading would expose the agent to basis risk due to the local-index price spread. Consequently, to study hedging one must consider the risk-preferences of the manager, an issue that was alluded to in Section 7.2. On the theoretical level, financial hedging would require analysis of a combined control problem, namely the mixture of optimal switching and portfolio optimization in an incomplete market.

References Ahn, H., A. Danilova, G. Swindle. 2002. Storing arb. Wilmott 1.

Carmona and Ludkovski: Optimal Switching for Energy Storage

29

ASCE. 1996. Hydroelectric pumped storage technology : international experience. New York, NY. Prepared by Task Committee on Pumped Storage of the Committee on Hydropower of the Energy Division of the American Society of Civil Engineers. Barrera-Esteve, C., F. Bergeret, C. Dossal, E. Gobet, A. Meziou, R. Munos, D. Reboul-Salze. 2006. Numerical methods for the pricing of Swing options: a stochastic control approach. Methodology and Computing in Applied Probability 8(4) 517–540. Bensoussan, A., J-L. Lions. 1984. Impulse Control and Quasi-Variational Inequalities. Gauthier-Villars, Paris. Bl¨ochlinger, L., T. Bollinger, J. Maurer, M. Semadeni. 2004. An evaluation model for the marked-to-market value of hydropower plants. Proceedings of the 6th IAEE European Conference on Modelling in Energy Economics and Policy. Brennan, M., E. Schwartz. 1985. Evaluating natural resource investments. Journal of Business 58 135–157. Cao, M., A. Li, J. Wei. 2004. Precipitation modeling and contract valuation: A frontier in weather derivatives. Journal of Alternative Investments 7(2) 92–99. Carmona, R., M. Ludkovski. 2008. Pricing asset scheduling flexibility using optimal switching. Applied Mathematical Finance in press. URL http://www.informaworld.com/10.1080/13504860802170507. Cartea, Alvaro, Thomas Williams. 2008. Uk gas markets: The market price of risk and applications to multiple interruptible supply contracts. Energy Economics 30(3) 829–846. Chen, Z., P. A. Forsyth. 2007. A semi-lagrangian approach for natural gas storage valuation and optimal operations. SIAM Journal on Scientific Computing 30(1) 339–368. Dayanik, S., M. Egami. 2005. A direct solution method for optimal switching problems of one-dimensional diffusions. Working paper. de Jong, C., K. Walet. 2003. To store or not to store. Energy Power Risk Management 8–11. Deng, S., Y. Oum, S.S. Oren. 2006. Hedging quantity risks with standard power options in a competitive wholesale electricity market. Naval Research Logistics 53(7) 697–712. Doege, J., H.-J. L¨ uthi, Ph. Schiltknecht. 2006. Risk management of power portfolios and valuation of flexibility. OR Spectrum 28(2) 267–287. Egloff, D. 2005. Monte Carlo algorithms for optimal stopping and statistical learning. Annals of Applied Probability 15(2) 1396–1432. Eydeland, A., K. Wolyniec. 2003. Energy and Power Risk Management: New Developments in Modeling, Pricing and Hedging. John Wiley& Sons, Hoboken, NJ. FERC. 2004. State of and issues concerning underground natural gas storage. Tech. rep., Federal Energy Regulatory Commission. Packet no. AD04-11-0000.

Carmona and Ludkovski: Optimal Switching for Energy Storage

30

Geman, H. 2005. Commodities and Commodity Derivatives : Modelling and Pricing for Agriculturals, Metals and Energy. John Wiley& Sons, Hoboken, NJ. Glasserman, P., B. Yu. 2004. Number of paths versus number of basis functions in American option pricing. Annals of Applied Probability 14(4) 2090–2119. Insley, M. 2003. On the option to invest in pollution control under a regime of tradable emissions allowances. Canadian Journal of Economics 35(4) 860–883. Jaillet, P., E. Ronn, S. Tompaidis. 2004. Valuation of commodity based Swing options. Management Science 50(7) 909–921. Keppo, J. 2002. Optimality with hydropower system. IEEE Transactions on Power Systems. 17(3) 583–589. Kluge, T. 2004. Pricing options in electricity markets. Preprint. Longstaff, F.A., E.S. Schwartz. 2001. Valuing American options by simulations: a simple least squares approach. Review of Financial Studies 14 113–148. Ludkovski, M. 2005. Optimal switching with application to energy tolling agreements. Ph.D. thesis, Princeton University. McNickle, D., E. Read, J. Tipping. 2004. The incorporation of hydro storage into a spot price model for the New Zealand electricity market. Sixth European Energy Conference: Modelling in Energy Economics and Policy. Zurich. Musiela, M., T. Zariphopoulou. 2004. A valuation algorithm for indifference prices in incomplete markets. Finance & Stochastics 8(3) 399–414. Nowak, M., W. R¨ omisch. 2000. Stochastic Lagrangian relaxation applied to power scheduling in a hydrothermal system under uncertainty. Annals of Operations Research 100(4) 251–272. Øksendal, B., A. Sulem. 2005. Applied stochastic control of jump diffusions. Springer-Verlag, Berlin. Pham, H., V. Ly Vath. 2007. Explicit solution to an optimal switching problem in the two-regime case. SIAM J. Control Optim. 46(2) 395–426. Pirrong, C. 1996. Corners and squeezes: the economics, law, and public policy of financial and commodity market manipulation. Kluwer Academic Publishers. Porchet, A., N. Touzi, X. Warin. 2008. Valuation of a power plant under production constraints and market incompleteness. Mathematical Methods of Operations research (to appear). Stentoft, L. 2004. Assessing the least squares Monte Carlo approach to American option valuation. Review of Derivatives Research 7(3) 129–168. Thompson, M., M. Davison, H. Rasmussen. 2003. Natural gas storage valuation and optimization: A real options approach. Tech. rep., University of Western Ontario. Tsitsiklis, J.N., B. van Roy. 2001. Regression methods for pricing complex American-style options. IEEE Transactions on Neural Networks 12(4) 694–703.

Carmona and Ludkovski: Optimal Switching for Energy Storage

31

Wilmott, P., S. Howison, J. Dewynne. 1995. The mathematics of financial derivatives. Cambridge University Press, Cambridge. A student introduction. Zervos, M. 2003. A problem of sequential entry and exit decisions combined with discretionary stopping. SIAM Journal on Control and Optimization 42(2) 397–421.

Copyright © 2019 PROPERTIBAZAR.COM. All rights reserved.