Carnegie Mellon University

Research Showcase @ CMU

Tepper School of Business

1-2010

Valuation of the Real Option to Store Liquefied

Natural Gas at a Regasification Terminal
Guoming Lai
University of Texas at Austin

Mulan X. Wang
DTE Energy Trading

Sunder Kekre
Carnegie Mellon University, skekre@cmu.edu

Alan Scheller-Wolf
Carnegie Mellon University, awolf@andrew.cmu.edu

Nicola Secomandi
Carnegie Mellon University, ns7@andrew.cmu.edu

Follow this and additional works at: http://repository.cmu.edu/tepper

Part of the Economic Policy Commons, and the Industrial Organization Commons

This Working Paper is brought to you for free and open access by Research Showcase @ CMU. It has been accepted for inclusion in Tepper School of
Business by an authorized administrator of Research Showcase @ CMU. For more information, please contact research-showcase@andrew.cmu.edu.
 Valuation of the Real Option to Store Liquefied Natural Gas at a
Regasification Terminal
Guoming Lai1 , Mulan X. Wang2 , Sunder Kekre3 , Alan Scheller-Wolf3 , Nicola Secomandi3
1 McCombs School of Business, University of Texas at Austin, 1 University Station, B6000, GSB
3.136, Austin, TX 78712-1178, USA
2 DTE Energy Trading, 414 S. Main Street, Suite 200, Ann Arbor, MI 48104, USA
3 Tepper School of Business, Carnegie Mellon University, 5000 Forbes Avenue, Pittsburgh, PA

15213-3890, USA
guoming.lai@mccombs.utexas.edu, wangx@dteenergy.com, {sk0a, awolf, ns7}@andrew.cmu.edu
Tepper Working Paper 2006-E99
Submitted: July 2009; Revised: January 2010

Abstract

The valuation of the real option to store liquefied natural gas (LNG) at the downstream terminal of an LNG
value chain is an important problem in practice. As the exact valuation of this real option is computationally
intractable, we develop a novel and tractable heuristic model for its strategic valuation that integrates models
of LNG shipping, natural gas price evolution, and inventory control and sale into the wholesale natural gas
market. We incorporate real and estimated data to quantify the value of this real option and its dependence
on the throughput of an LNG chain, the type of price variability, the type of inventory control policy
employed, and the level of stochastic variability in both the shipping model and the natural gas price model
used. In addition, we develop an imperfect information dual upper bound to assess the effectiveness of our
heuristic, and find that our method is highly accurate. Our approach also has potential relevance to value
the real option to store other commodities in facilities located downstream from a commodity production or
transportation stage, such as petroleum and agricultural products, chemicals, and metals, or the real option
to store the input used in the production of a commodity, such as electricity.

1 Introduction

Liquefied natural gas (LNG) is natural gas cooled to liquid state at -260F; liquefaction reduces
the volume of natural gas by a factor of more than 600, making storage and shipping practical
(Greenwald [26], Tusiani and Shearer [51]). Special ocean going vessels load LNG at liquefaction
facilities (for example in Trinidad and Tobago, Australia or Qatar), transport it (for days or weeks),
and unload it at terminals (for example in the U.S., Europe or Japan). At these terminals LNG is
pumped into storage tanks, regasified, and then distributed via pipelines or, sometimes, by trucks.
The Energy Information Administration (EIA) has projected that local production of natural
gas will soon be unable to meet its increasing demand in several industrialized countries, and
expects LNG imports to play an important role in bridging this gap (EIA [20, 21]). This long
term projected increase in the world’s natural gas demand is primarily due to natural gas being a
relatively environmentally clean and abundant fuel, which has helped to make it the fuel of choice
for many new power generation projects. Many of these long term forecasts predate the current

1
economic recession, however, as of March 2009 EIA (EIA [22]) forecasts that in 2018 U.S. LNG
imports will peak at 4.9 times their 2008 levels before declining to 2.7 times these levels in 2030.
We have started to see the unfolding of some of these predicted increases, which Jensen [31]
refers to as the “LNG revolution.” Several liquefaction capacity expansion and greenfield projects
have been announced, a number of new terminals have been proposed in North America, and some
of these have been recently completed. Obtaining access to these terminals is necessary to bring
LNG into the natural gas distribution system, and requires leasing storage space and regasification
capacity from the terminal’s operating companies. Hence, industry players face the challenge of
assessing the value of downstream terminal leasing contracts.
The value of such a contract consists of the delivery value and the storage value. In this paper,
we focus on determining the storage value, which requires valuing the real option to store LNG
at a regasification terminal. This topic has not yet been studied in the literature and may not
be well understood in practice. For instance, Holcomb [29] attributes little storage value to LNG
regasification terminals, but Cheniere Energy (www.cheniere.com) attributes strategic importance
to LNG storage at such terminals, having embarked on the construction of a network of LNG
terminals along the U.S. Gulf Coast that, once completed, will feature the largest availability of
LNG storage and regasification capacity in the U.S.
Our interest in this paper is the valuation of the real option to store LNG at regasification
facilities in the presence of a wholesale market for natural gas, which is the case in the U.S., U.K.,
and some parts of Europe (Tusiani and Shearer [51, p. 26]). Exact valuation of this real option is
computationally intractable, especially using an operational model (because of its fine time scale).
Thus, we take a strategic approach and develop a novel and tractable model for the heuristic
valuation of this real option. Our approach integrates models of shipping, commodity price evolu-
tion, inventory control, and storage valuation based on closed queueing networks (CQNs), lattice
approximations of Ito processes, a Markov decision process (MDP), and Monte Carlo simulation,
respectively. Specifically, we extend the CQN model of LNG shipping of Koenigsberg and Lam [33]
and integrate it with (1) lattice approximations of the commodity price models of Jaillet et al. [30]
and Schwartz and Smith [44], modeling seasonality as in Jaillet et al. [30], and (2) an inventory
control MDP, whose policy is evaluated by simulation. We also develop an imperfect information
dual upper bound (Brown et al. [11]) to assess the effectiveness of our practical heuristic model.
We apply our model in a numerical study to estimate and analyze the storage value. We
consider a realistic LNG chain consisting of liquefaction in North Africa, shipping to Lake Charles,
Louisiana, and regasification and sale into the Louisiana natural gas wholesale market using the

2
spot price at Henry Hub, the delivery location of the New York Mercantile Exchange (NYMEX)
natural gas futures contract. We calibrate the price evolution models to prices of traded NYMEX
natural gas futures and options on futures. This application provides the following findings:
(1) The estimates of the value of the real option to store during 12 years range from $89M to
$726M, when the available storage space and LNG chain nominal throughput (ignoring congestion)
vary from 3 billion cubic feet (BCF) and 0.7536 million tons per annum (MTPA) to 24BCF and
7.5362MTPA, respectively, and the regasification capacity is 2BCF/Day (Table 4 in Online Ap-
pendix A reports relevant units of measurements and conversion factors); this takes between 5 and
53 Cpu minutes, depending on the LNG chain configuration.
(2) Comparisons with our upper bound estimates indicate that our model’s valuations are very
accurate, being on average more than 99% of the upper bound estimates.
(3) The value of the real option to store can be nonmonotonic in the throughput of the LNG
chain, which with homogeneous ships is the number of LNG cargos delivered per unit time. This
implies that discretionary regasification capacity is necessary for LNG storage to be most valu-
able, and suggests that the discussion in Holcomb [29] is mainly relevant to situations when the
regasification capacity is comparable to throughput.
(4) Depending on the system configuration and the price model used, approximately between
50% and 62% of the value of storage can be attributed to natural gas price volatility (stochastic
variability), with the remaining part attributable to price seasonality (deterministic variability).
(5) Analysis of our results shows that usage of the forward looking optimal policy of our MDP
is important when the LNG throughput is low relative to the storage space.
(6) The storage value is fairly insensitive to how one models stochastic variability in the shipping
process, due to the low level of congestion in the system configurations that we consider. Thus,
an essentially deterministic model of throughput is often sufficient for valuation, in which case the
relevant Cpu times decrease to 1.4-1.7 minutes. In contrast, the storage value decreases substantially
when we model the price evolution using one factor, rather than two factors. Hence, how one models
price uncertainty has a much stronger impact on the value of storage than how one models stochastic
variability in the processing times in an LNG chain.
Our model and computational results have managerial relevance: Our model is very accurate
and can be used by LNG players to support the negotiation of contracts for access to LNG re-
gasification terminals, and our computational results significantly enhance the understanding of
the value of LNG storage at such a terminal. Our model also has potential applicability, with
suitable modifications, in other settings in which a storage facility is located downstream from the

3
production or shipping stages of a given commodity, or upstream of such processes when it is used
to store their input. Examples include petroleum and agricultural products, chemicals, metals, and
the production of electricity from nonrenewable and renewable sources, such as coal, natural gas,
water, sunlight, and wind (see, e.g., Bannister and Kaye [4], Geman [25], and Baxter [6]).
We proceed by reviewing the relevant literature in §2. We introduce our operational and strate-
gic valuation models in §3 and §4, respectively. We present our heuristic valuation and upper bound
models in §5 and §6, respectively, and quantity the value of storage in §7. We conclude in §8.

2 Literature Review

Our work is unique with respect to the LNG literature: Kaplan et al. [32] and Koenigsberg and
Lam [33] address the modeling of the shipping stage of an LNG chain, but do not investigate the
quantification of the storage value of an LNG regasification terminal, as we do here. Özelkana et
al. [38] present a deterministic optimization model to analyze the design of LNG terminals. In
contrast, our model captures both price and shipping uncertainty. Grønhaug and Christiansen [27]
propose a deterministic optimization model for the tactical management of LNG inventory and ship
routing, whereas we investigate the strategic valuation of downstream LNG storage accounting for
both price and shipping uncertainty.
Our work is related to the real option literature (Trigeorgis [49]) dealing with applications in
commodity and energy industries (see, e.g., Smith and McCardle [47], Clewlow and Strickland [16],
Eydeland and Wolyniec [23], and Geman [25]). To the best of our knowledge, this literature has
not yet studied the valuation of LNG storage: Geman [25, pp. 246-249] briefly describes the LNG
storage setting and Abadie and Chamorro [1] use Monte Carlo simulation to value natural gas
investments, including an LNG plant, but not LNG storage.
Several authors have studied the related real option valuation of natural gas storage, including
Carmona and Ludkovski [12], Chen and Forsyth [14], Boogert and de Jong [10], Lai et al. [34],
Secomandi [45], and Thompson et al. [48], among others. The main difference between natural
gas and LNG storage is that the inflow of commodity into the storage facility is controllable in the
natural gas case, but not in the LNG case. In addition, our work differs from that of Lai et al. [34]
in the natural gas price model, the dynamic programming approximation, and the upper bound
that we use, as well as in some of the managerial insights that we obtain.
Our model is also related to models of hydropower production and sale, especially those with
price uncertainty (see, e.g., the review by Wallace and Fleten [52]). Some of these models are based
on stochastic (mathematical) programming, whereas our approach uses an MDP within a heuristic

4
decomposition scheme with an imperfect information dual upper bound (Brown et al. [11]), which
differ from MDP type models available in this literature, such as that of Näsäkkälä and Keppo [36]
(see also related work by Harrison and Taksar [28], Drouin et al. [17], and Lamond et al. [35]).

3 Operational Model

In this section we present an operational model for the management of an LNG chain. This model
is computationally intractable; we use it to motivate the strategic valuation model in §4.
Shipping and Terminal Operations. Faithful representation of the interplay between the
LNG shipping, storage, and regasification activities would require modeling them using a continuous
time framework. For ease of exposition, we use a discrete time approach where a given finite time
interval of length T , which we denote by set T := [0, T ], is subdivided into I small time intervals,
each of length ∆t, in set I := {1, . . . , I} (the length ∆t is discussed further below). We index this
set by i and denote time by t and the time when time period i starts by ti (t1 := 0 and tI+1 := T ).
A fleet of N identical LNG ships, each with cargo size C, perform the following activities:
loading at the upstream port, loaded transit to the downstream port, entering the downstream
port, unloading at this port, leaving this port, and ballast transit to the upstream port. Ships are
dedicated and loop between the liquefaction and regasification facilities, by far the most typical
setting in the LNG industry (Greenwald [26], Tusiani and Shearer [51]). Congestion may occur at
the upstream and downstream ports, that is, ships may queue up at these facilities.
We abstract from the details of natural gas liquefaction assuming ample supply. The fact that
LNG liquefaction facilities are designed to run at full capacity, being served by an appropriate
number of ships to satisfy this capacity (Flower [24, p. 96]), supports this. Loading aggregates the
following activities: entering the port (traversing the entry channel), loading the ship, and leaving
the port (traversing the exit channel or the entry channel again if there is only one channel at
the given port). In contrast, as will become apparent soon, it is useful to separately model these
activities, with loading replaced by unloading, for the downstream port.
The state of the shipping system at time ti is si ; we simplify ti to i when used as a subscript.
(Quantities that are underlined will be simplified in later sections.) Each si is an N dimensional
vector of triples that describe the activity performed by each ship, the elapsed time since the start
of this activity, and the position of each ship in any relevant queue. The set of all possible shipping
states at time ti is S i . There is usually uncertainty associated with shipping operations (Kaplan
et al. [32], Ronen [40]). In this section we do not postulate a specific model of this uncertainty,
but we suppose that ∆t is chosen such that it is reasonable to assume that at most one activity

5
can complete with positive probability during a time interval of length ∆t. This assumption would
give rise to daily or even smaller time intervals in applications.
We denote by xi the inventory available at the downstream facility at time ti . Let 0 and X,
respectively, denote the minimum and maximum levels of inventory that can be held in storage
at this facility (a positive minimum inventory level can be easily accommodated at the expense of
additional notation). Hence, the quantity xi is constrained to be in set X := [0, X].
We assume that at most one ship can unload its cargo at any one time. The LNG unloading
and regasification rates are deterministic. We let zi ∈ Z := [0, C] denote the LNG inventory
remaining onboard the ship that is unloading at time ti . Consistent with the objective of running the
liquefaction facility at full capacity, we restrict attention to control policies that do not intentionally
slow down the shipping system. Thus, the amount of LNG unloaded during time period i is not
discretionary and is a deterministic function of xi and zi that we denote by ui (xi , zi ). We denote by
q i the amount of LNG one chooses to regasify during time period i. The set of feasible values that
this quantity can take on depends on xi , zi , and ui (xi , zi ). We denote this set by Q(xi , zi ), since
ui (xi , zi ) is a function of xi and zi . We let z̃i be the random variable that describes the inventory
onboard the ship that unloads during time period i, which depends on the evolution of the shipping
process (we denote random entities as ˜·); zi is equal to C if a ship starts unloading its cargo at time
ti , zi−1 − ui−1 (xi−1 , zi−1 ) otherwise (here i > 1; z1 is determined by s1 ).
Revenue and Cost Structures. We use the Markovian vector-valued stochastic process
{p̃t , t ∈ T } to describe the continuous time evolution of a price state random vector p̃t , with real-
ization pt ∈ <P (P > 1 is an integer). The natural gas spot price at the downstream facility at
time t is the known function gt (pt ) : <P → <+ . We assume an arbitrage free and complete market
for natural gas futures at this location and denote by E[· | pt ] risk neutral conditional expectation
given pt (see, e.g., Duffie [18]). We assume that E[g(p̃t0 ) | pt ] < ∞, ∀t, t0 ∈ T and t 6 t0 . This
ensures that the value functions of the MDPs discussed below and in §§4-6 are finite.
This representation allows us to capture single (P = 1) and multiple (P > 1) factor models
of the evolution of the natural gas spot price, such as the one- and two-factor models of Jaillet et
al. [30] and Schwartz and Smith [44], the latter modified to use deterministic monthly seasonality
factors as in the one factor model. In the two factor model, the natural gas spot price at time
t is gt (p̃t ) := fm(t) exp(p̃1t + p̃2t ), with fm(t) the time t seasonality factor and m(t) the month
corresponding to time t. Using the same notation and terminology of Schwartz and Smith [44], we
let χt := p1t and ξt := p2t , and refer to χt and ξt as the time t values of the short term deviation
factor and long term (equilibrium) factor, respectively. The risk neutral dynamics of these factors

6
are (Schwartz and Smith [44])

dχt = (−κχ χt − λχ )dt + σχ dzχ∗ , (1)

dξt = (µξ − λξ )dt + σξ dzξ∗ , (2)

dχt dξt = ρχξ dt. (3)

Here κχ and σχ are the speed of mean reversion and the volatility of χt ; µξ and σξ are the drift
and the volatility of ξt ; λχ and λξ are the risk premia associated with the two factors; and dzχ∗ and
dzξ∗ are increments to standard Brownian motions with instantaneous correlation ρχξ .
When the equilibrium factor is constant the two factor model (1)-(3) reduces to the one factor
model of Jaillet et al. [30] (see Schwartz and Smith [44, p. 894]), in which the short term factor
pt ≡ χt mean reverts to the constant risk adjusted equilibrium level ξ ∗ , rather than zero, and the
deseasonalized spot price at time t is exp(χt ) (seasonality is modeled as in the two factor model).
The risk neutral dynamics of this short term factor are

dχt = κ(ξ ∗ − χt )dt + σdz ∗ , (4)

where κ and σ are the speed of mean reversion and the volatility of this factor, and dz ∗ is an
increment to a standard Brownian motion.
We account for regasification sales during time period i by multiplying the released quantity,
net of regasification fuel losses (explained in detail below), by the natural gas price prevailing at
time ti . We also assume that the quantity sold does not affect the market price of natural gas and
that the price state vector and the shipping processes evolve independently. These two assumptions
are realistic for modeling an LNG system whose regasification terminal is located in the southern
part of the U.S., e.g., Louisiana, where the natural gas spot market is fairly liquid.
There are operating costs and fuel requirements associated with the physical flows along the
chain (Flower [24]). We denote by h ∈ <+ the per unit and time period physical inventory holding
cost charged against the inventory xi available at time ti at the downstream terminal. We let φ be
the fuel needed to regasify one unit of LNG, that is, the LNG to natural gas yield is 1 − φ. We
denote by c the cost of unloading (handling) one unit of LNG at the downstream terminal.

Remark 1 (Units of measurement). LNG capacity is typically measured in MTPA, cargos in

cubic meters (CM), and LNG downstream storage space and regasification capacity in BCF and
BCF/day. For simplicity, we assume that all the physical quantities are expressed as functions of
million British thermal units (MMBTU), since the natural gas price in the U.S. is expressed in U.S.
dollars per million British thermal units ($/MMBTU). (See Table 4 in Online Appendix A.)

7
 MDP. We now formulate an optimization model to control the inventory level at the down-
stream facility. We assume that the stochastic evolution of the shipping process is uncorrelated with
the evolution of the price of the market portfolio. Thus, the statistical and risk neutral dynamics
of the shipping process can be taken to be identical. This allows us to formulate our model using
risk neutral valuation (see, e.g., Smith [46]).
The state of the system at time ti is (xi , zi , si , pi ). We denote a control policy by π and let Π be
the set of feasible policies. We define the reward obtained in time period i as ri (xi , ui (xi , zi ), q i , pt ) :=
gi (pi )(1 − φ)q i − hxi − cui (xi , zi ). We denote by WI+1 (xI+1 , zI+1 , sI+1 , pI+1 ) the salvage value of
the LNG that is onboard the ships or is available at the regasification facility at time tI+1 . The
objective is to solve the following optimization model:
" #
X
i−1 π π π π I π π π
max E δ ri (x̃i , ui (x̃i , z̃i ), q̃ i , p̃i ) + δ WI+1 (x̃I+1 , z̃I+1 , s̃I+1 , p̃I+1 ) | x1 , z1 , s1 , p1 , (5)
π∈Π
i∈I

where δ denotes the one period risk free discount factor and E[·|x1 , z1 , s1 , p1 ] denotes expectation
with respect to the joint risk neutral probability distribution of the relevant random variables in-
duced by feasible policy π, a dependence indicated by superscripting π, conditional on (x1 , z1 , s1 , p1 ).
Denoting by Wi (xi , zi , si , pi ) the optimal value function in state (xi , zi , si , pi ) at time ti , model
(5) can be reformulated recursively as follows:

Wi (xi , zi , si , pi ) = max ri (xi , ui (xi , zi ), q i , pi )

q i ∈Q(xi ,zi )
h i
+δE Wi+1 (xi + ui (xi , zi ) − q i , z̃i+1 , s̃i+1 , p̃i+1 ) | xi , zi , si , pi ,

∀i ∈ I, xi ∈ X , zi ∈ Z, si ∈ S i , pi ∈ <P . (6)

The value of the regasification terminal at time t1 is W1 (x1 , z1 , s1 , p1 ); the value of its storage
component at this time is this value minus the regasification terminal value under the policy that
regasifies as much as possible in every time period. However, the curse of dimensionality makes
model (5) computationally intractable, and thus approximations are needed for valuation purposes.

4 Strategic Valuation Model

In this section we describe our approximate model for the strategic valuation of a downstream
LNG storage and regasification terminal. Given our objective of strategic valuation, instead of
attempting to approximately solve the operational model (6), we formulate a simplified model:
Specifically, we simplify (1) the coupling between the shipping process and the management of the
inventory at the terminal, brought about by the LNG unloading step, and (2) the evolution of the

8
 2. Loaded Transit

1. Loading Port LNG 3. Unloading Port

Service Shipping System Service

4. Ballast Transit

Figure 1: Representation of an LNG shipping system.

price state vector. Although simpler than model (6), the resulting model remains difficult to solve,
but it paves the way for developing the computationally efficient heuristic that we discuss in §5.
Time Aggregation and Shipping and Terminal Operations. The valuation of the termi-
nal involves time horizons exceeding 10 years. Thus, we increase the time scale of the operational
model by aggregating the time periods in set I into J longer time periods in set J := {1, . . . , J};
the start of time period j is time tj and tJ+1 := T . All the aggregate time periods have the same
length, which is application dependent; we take it to be one month in §7.
This aggregate time scale no longer allows us to track the detailed evolution of the interaction
between LNG shipping, storage, and regasification; we are only able to track this interaction at
the level of the number of ships that request to unload their cargos during an aggregate time
period. That is, we treat all cargoes unloaded in this period identically; they all can be regasified
and sold during the period at the prevailing price at the start of the period. Thus, the quantity
of interest to us now is the distribution of the number of cargos or, equivalently, the amount of
LNG unloaded during an aggregate time period. We denote by ũj , with realization uj , the random
amount of LNG unloaded during aggregate time period j and assume that it can take values in set
Uj (xj ) := {0, C, . . . , NjU (xj )C}, where NjU (xj ) is the maximum number of ships that can unload
their cargos during aggregate time period j when the inventory level at time tj is xj . Given this
model, we now need to specify NjU (xj ) and the probability distribution of ũj .
We do this by extending the CQN model of the LNG shipping system proposed by Koenigsberg
and Lam [33] and illustrated in Figure 1. Different from the shipping system described in §3, this
system combines into a single unloading activity the activities of entering/exiting the downstream
port and unloading (we discuss this CQN in detail in §5). We define NjU (xj ) := b(X + Q − xj )/Cc,

9
where Q denotes the regasification capacity available during an aggregate time period. Consistent
with §3, this definition implicitly assumes any operational policy seeks to maximize the utilization
of the available regasification capacity at the downstream facility. The probability distribution of
ũj is determined by the transition laws of the CQN with the restriction that no more than NjU (xj )
ships are allowed to transition from the unloading to the ballast transit stage during aggregate time
period j. If more than NjU (xj ) ships request to unload during this time period, the excess ships
are blocked at the downstream facility at least until time tj+1 . This contrasts with the CQN of
Koenigsberg and Lam [33] that does not include regasification capacity or the maximum inventory
limit at this facility; that is, our CQN models both congestion at the loading and unloading stages,
arising from uncertainty in the CQN as in Koenigsberg and Lam [33], and blocking at the unloading
stage, arising from constraints on the regasification capacity and the maximum inventory space.
The vector sj denotes the state of the shipping system at time tj , which comprises the number
of ships in the various stages, including those blocked; sj is necessary to obtain the probability
distribution of ũj . The set of these possible shipping states is Sj .
We denote by qj ∈ Q(xj , uj ) the amount of LNG that can be feasibly regasified and sold during
aggregate time period j ∈ J given xj ∈ X and uj ∈ Uj (xj ); here Q(xj , uj ) constrains qj as follows:
given tj , xj , uj , and qj , the inventory level at time tj+1 is xj+1 = xj +uj −qj , so that the restrictions
xj+1 ∈ X and qj ∈ [0, Q] imply that Q(xj , uj ) = [max{0, xj + uj − X}, min{Q, xj + uj }].
Discrete Time and Space Price State Process. The domain of the price state vector in
model (6) is unbounded. From now on, we approximate the risk neutral evolution of the vector pt
using a discrete time and space stochastic process that at each time tj , with j ∈ J ∪ {J + 1}, can
take values in the finite set Pj ⊂ <P ; we discuss the generation of this stochastic process in §5.
MDP. Our strategic valuation model is an MDP with stage set J ∪ {J + 1}. Its state in stage
j is (xj , sj , pj ). We let h ∈ <+ denote the per unit physical inventory holding cost during an
aggregate time period and δ the risk free discount factor for this time period. The reward obtained
during aggregate time period j is defined as rj (xj , uj , qj , pj ) := gj (pj )(1 − φ)qj − hxj − cuj . We
denote by Vj (xj , sj , pj ) the optimal value function of our strategic valuation model in stage j and
state (xj , sj , pj ). This function is defined as follows:

Vj (xj , sj , pj ) = E [vj (xj , sj , ũj , pj ) | xj , sj ] , ∀j ∈ J , xj ∈ X , sj ∈ Sj , pj ∈ Pj , (7)

vj (xj , sj , uj , pj ) := max rj (xj , uj , qj , pj ) + δE[Vj+1 (xj + uj − qj , s̃j+1 , p̃j+1 ) | xj , sj , uj , pj ],

qj ∈Q(xj ,uj )
(8)

with boundary conditions VJ+1 (xJ+1 , sJ+1 , pJ+1 ) := [gJ+1 (pJ+1 )(1 − φ) − h]xJ+1 , for all xJ+1 ∈ X ,

10
sJ+1 ∈ SJ+1 , pJ+1 ∈ PJ+1 ; that is, at the final time, tJ+1 , any remaining inventory is released and
sold (for simplicity, we account for the holding cost using the coefficient h). During each aggregate
time period j ∈ J , an amount uj ∈ Uj (xj ) of LNG is unloaded at the downstream terminal at
cost cuj , and becomes available for regasification and sale during this time period; holding cost
hxj is charged against the initial inventory xj ; an optimal amount of LNG from the total available
inventory xj + uj is regasified, and a fraction 1 − φ of this amount is sold into the wholesale
natural gas spot market at the prevailing price gj (pj ). A stochastic transition to the next aggregate
period accounts for the uncertainty in the natural gas spot price and the shipping process, taking
into account the inventory dynamics (as in §3, the stochastic evolution of the shipping process is
independent of that of the price state vector and does not require any risk adjustment).
Comparing models (6) and (7)-(8) reveals two notable effects of our modeling simplification. (1)
The quantity zi is a state variable in the operational model but not in the strategic model, as at the
aggregate time scale the precise status of the inventory onboard an unloading ship is insignificant.
(2) The strategic model has factored expectations, that is, the maximization on the right hand side
of (8) is taken over Q(xj , uj ), assuming the quantity to be unloaded during aggregate time period
j is known; this is supported by the fact that in actual operations one can fairly reliably schedule
ship unloading during an aggregate time period, e.g., one month. This assumption is required
to remove the dependence of the set of feasible regasification decisions from the evolution of the
shipping system within an aggregate time period, that is, to make our time aggregation function.
Storage Valuation. We denote by V1G (x1 , s1 , p1 ) the value in the initial stage and state of a
greedy inventory release policy that in each stage regasifies and sells as much as possible of each
incoming cargo upon receipt. We define the value of storage in the initial stage and state as

S1 (x1 , s1 , p1 ) := V1 (x1 , s1 , p1 ) − V1G (x1 , s1 , p1 ). (9)

5 Heuristic Strategic Valuation Model

In this section we describe our heuristic strategic valuation model.

Heuristic Model Definition. Fundamentally, what continues to make model (7)-(8) com-
putationally difficult is the curse of dimensionality brought about by the coupling between the
shipping/unloading process and the management of the inventory at the terminal, which may re-
sult in ship blocking. Our heuristic model deals with the coupling of these activities in two steps.
First, it decouples them by assuming away the possibility of ship blocking, which eliminates the
dependence between the unloading and inventory processes. This makes computing a heuristic

11
inventory release policy tractable. Second, it reinstates this coupling by including the possibility
of ship blocking at an aggregate time scale when computing the value of storage when using this
policy. Our heuristic model performs these steps by integrating (1) the shipping model, (2) the
price evolution model, (3) the inventory release model, and (4) the storage valuation model.
Shipping Model. The shipping model derives a stochastic description of the number of LNG
cargos unloaded during an aggregate time period. This model decouples the natural gas liquefaction
and LNG shipping activities from the management of downstream LNG storage. In addition, we
simplify our problem by modeling the unloading process as a sequence of independent and identically
distributed random variables, based only on the overall shipping network configuration, not on the
positions of individual ships. Each such random variable, denoted by ũ, represents the amount of
LNG unloaded at the regasification facility by the ships in the network during an aggregate time
period. The distribution of ũ, the set of nonzero probability realizations U and the probability
Pr{ũ = u} for each u ∈ U, is required in the inventory release model.
We determine this distribution in two steps by extending the original CQN of Koenigsberg and
Lam [33]. This extension is a useful and flexible abstraction of LNG shipping: In its basic form
with exponentially distributed processing times, this model provides a conservative estimate of the
throughput of an LNG system (ignoring ship blocking); in a more advanced form with the shipping
times modeled as multistage Coxian distributions, which can approximate the distribution of service
times arbitrarily closely, our model computes a less conservative estimate of this throughput. Our
shipping model can also accommodate other approaches to the modeling of the throughput of an
LNG system. For example, in §7 we discuss results obtained with an essentially deterministic
shipping model that, ignoring congestion, provides an optimistic estimate of throughput.
As an extension to Koenigsberg and Lam [33], suppose that the loading and unloading blocks
in Figure 1 are first come first served (FCFS) exponential queues, and the transit blocks are ample
server (AS) exponential or multistage Coxian queues; the latter case allows more flexibility in
modeling variability than the exponential distribution (see Osogami and Harchol-Balter [37] for a
discussion of Coxian modeling; Koenigsberg and Lam [33] only use exponential queues). With this
representation, the shipping system is a particular CQN called a BCMP network having a closed
form, product form stationary distribution, as defined and proved by Baskett et al. [5]. We denote
by N the set of all possible states of the system, that is, the set of vectors that represent the number
of ships in each stage. The steady state probability that the random variable state of the system,
ñ, is equal to n ∈ N is γ(n) := Pr{ñ = n}, and the CQN throughput can be computed by standard
methods (see Baskett et al. [5]).

12
 We propose the rolling forward method to compute the distribution of ũ. This method uses the
stationary distribution γ(·) as a starting point. It transitions the CQN forward through time from
its stationary distribution, tracking the distribution of unloaded amounts over an aggregate time
period: we uniformize (see, e.g., Asmussen [3]), condition on the number of “events” that occur
over this time period, and calculate the distribution of the number of ships unloaded given the
number of total events conditioning on the initial state of the shipping system, as specified by γ(·).
Hence, the distribution of random variable ũ is a one dimensional table.
We compute this distribution analytically, that is, we do not use Monte Carlo simulation.
Specifically, denote by ã the random number of events that occur during an aggregate time period,
a Poisson random variable with appropriate mean, and by η̃ the random number of unloaded ships
during this time period. Let A be such that Pr{ã > A} 6 ² for arbitrarily small ² ∈ <+ . Denote
by Pr{η|a, n} the probability that η ships are unloaded during an aggregate time period given
that the system is initially in state n and a total of a events occur. We have developed analytical
expressions for Pr{η|a, n} through a forward recursion in η, but, since they are somewhat lengthy,
in the interest of space we do not present them here. We compute the distribution of η̃, and hence
P P
that of ũ ≡ η̃C, as Pr{η̃ = η} = [ n∈N γ(n) A a=0 Pr{ã = a} Pr{η|a, n}]/(1 − ²).

Price Evolution Model. The price evolution model generates a lattice representation of the
stochastic evolution of the price state vector during the given time horizon. We use the approach
described by Tseng and Lin [50, §3] to build two trinomial lattices for the two factor model (1)-(3)
and one trinomial lattice for the one factor model (4), with a time step equal to the length of an
aggregate time period. We also use multiplicative adjustment factors to calibrate these lattices
to the natural gas forward curve observed at the beginning of the time horizon. These calibrated
lattices are the stochastic process {pj , j ∈ J ∪ {J + 1}} used in our inventory release model.
Inventory Release Model. The inventory release model incorporates the output of the
shipping and price evolution models to determine an inventory release policy that can be used to
make LNG regasification and natural gas sale decisions. This is a simplified version of model (7)-(8)
whose state in each stage j ∈ J ∪ {J + 1} is the pair (xj , pj ), for all xj ∈ X and pj ∈ Pj ; that is,
the state of the shipping system is no longer part of the state as the aggregate shipping/unloading
process is approximated through the distribution of ũ.
We make the assumption that in each aggregate time period it is possible to regasify all of the
LNG unloaded from ships arriving in this time period, that is, Q > U := max{u : u ∈ U}. (We
could also penalize “blocked” LNG if Q > U were to be unreasonable.) This assumption ensures
that ships cannot be blocked at the downstream facility and makes the set Q(xj , u), defined in §4,

13
 Hold Sell

!
" $ !# "

Figure 2: The structure of the optimal inventory release policy in stage j for a given value of the
price state vector pj .

nonempty for all j ∈ J , xj ∈ X , and u ∈ U.

We denote by VjH (xj , pj ) the value function of the inventory release model in stage j and state
(xj , pj ). This function is defined as
£ ¤
VjH (xj , pj ) = E vjH (xj , pj , ũ) , ∀j ∈ J , xj ∈ X , pj ∈ Pj , (10)
£ H ¤
vjH (xj , pj , u) := max rj (xj , u, qj , pj ) + δE Vj+1 (xj + u − qj , p̃j+1 ) | pj , (11)
qj ∈Q(xj ,u)

H (x
with boundary conditions VJ+1 J+1 , pJ+1 ) := [gJ+1 (pJ+1 )(1 − φ) − h]xJ+1 , for all xJ+1 ∈ X and

pJ+1 ∈ PJ+1 . This formulation is interpreted in a manner analogous to that of model (7)-(8).
Structural Analysis. We now study the structure of the optimal policy of model (10)-(11).
This analysis greatly facilitates the computation of this policy and its efficient use in the storage
valuation model, which makes our integrated model practical. Proposition 1 characterizes the
optimal value function of model (10)-(11). (Online Appendix B includes the proofs for this section.)

Proposition 1 (Optimal value function). In every stage j ∈ J ∪ {J + 1}, the function VjH (xj , pj )
is concave in xj ∈ X for each given pj ∈ Pj .

Turning to the optimal sale action, we define the quantity qj∗ (xj , pj , u) as the largest element that
optimizes the right hand side of (11). Any feasible sale cannot be smaller than max{0, xj + u − X},
because in aggregate time period j one must execute the forced sale qjF (xj , u) := max{0, xj +u−X}
to avoid a tank overflow due to incoming cargos. We call the difference between the feasible
sale qj (xj , pj , u) and the forced sale qjF (xj , u) the optional sale, and denote it by qjO (xj , pj , u) :=
qj (xj , pj , u) − qjF (xj , u); that is, qj (xj , pj , u) ≡ qjF (xj , u) + qjO (xj , pj , u).

Proposition 2 (Optimal policy structure). In every stage j ∈ J , given pj ∈ Pj , define bj (pj ) ∈ X

H (x
as the smallest element of arg maxxj+1 ∈X δE[Vj+1 j+1 , p̃j+1 ) | pj ] − gj (pj )(1 − φ)xj+1 . In stage
∗
j ∈ J , it is optimal to try to sell down to bj (pj ), that is, qj∗ (xj , pj , u) = qjF (xj , u) + qjO (xj , pj , u),
∗
with qjO (xj , pj , u) := min{max{xj + u − qjF (xj , u) − bj (pj ), 0}, Q − qjF (xj , u)}.

The quantity bj (pj ) can be interpreted as a basestock target for optimal optional sales; it is a
target because it is constrained by the limited regasification capacity Q. Given price state vector

14
pj , bj (pj ) partitions the feasible inventory set into two regions, one in which it is optimal to hold
inventory and one in which it is optimal to sell down to the basestock target, as illustrated in
Figure 2. We point out that when it is optimal to sell, it can be optimal to stop selling rather than
draining the terminal as much as possible; formally bj (pj ) ∈ (0, X) is possible, as discussed in §7.
Under the assumptions of Proposition 3, which is related to Propositions 2 and 3 of Secomandi
[45], computing an optimal inventory release policy can be done efficiently.

Proposition 3 (Optimal policy computation). Suppose that X, Q, and each u ∈ U are integer
multiples of some maximal L ∈ <+ . Then, for every j ∈ J ∪ {J + 1}, VjH (xj , pj ) is piecewise
linear and continuous in xj ∈ X for each given pj ∈ Pj , it changes slope at values that are integer
multiples of L, and bj (pj ) is an integer multiple of L (bJ+1 (pJ+1 ) := 0).

The practical implication of this result is that in each stage one needs to compute the optimal
value function only for a finite number of inventory levels, namely 0, L, 2L, . . . , X, and the search
for an optimal basestock target can be limited to one of these values.
Storage Valuation Model. The storage valuation model determines the following lower
bound estimate on the value of storage S1 (x1 , s1 , p1 ), defined by (9):

Ŝ1B (x1 , s1 , p1 ) := V̂1B (x1 , s1 , p1 ) − V̂1G (x1 , s1 , p1 ). (12)

Here, V̂1B (x1 , s1 , p1 ) is the estimate of the value function in the initial stage and state of the optimal
basestock target policy of model (10)-(11) when used as a heuristic policy for model (7)-(8); that
is, when ship blocking at the downstream facility can occur in the manner modeled in §4. This
estimate is obtained by applying this policy within Monte Carlo simulations of the price state
vector process, from the lattice representation of its evolution, and the shipping process with ship
blocking. The quantity V̂1G (x1 , s1 , p1 ) is the estimated value in the initial stage and state of the
greedy policy, discussed at the end of §4, obtained in conjunction with V̂1B (x1 , s1 , p1 ) by using
common random numbers and multiplicative factors computed to adjust the sampled price state
vectors to be consistent with the observed forward curve (these adjustment factors are not the same
ones used in the calibration of the lattices).
The storage valuation model can be used to also estimate the values of the real option to
store due to price seasonality and volatility. The former is the value of this real option with
deterministic natural gas price evolution equal to the forward curve at time 0, which exhibits
significant seasonality, as illustrated in §7. We estimate this value of storage in a manner analogous
to the estimation of the value of storage, except that we compute an optimal inventory release

15
policy using deterministic price dynamics. The latter is the difference between the value of storage
and the value of storage due to price seasonality, and we estimate it accordingly.

6 Upper Bound

In this section we discuss a model that estimates an upper bound on the value of storage to
complement our lower bound of §5. This model is an imperfect information dual upper bound
model, in the sense of Brown et al. [11]. It uses the same periodic review setting as the strategic
model (7)-(8) and employs Monte Carlo simulation to generate a set of sequences of unloaded cargos
of the form {uj , j ∈ J } by sequentially sampling from the stochastic shipping process used by the
shipping model described in §5; that is, ignoring ship blocking. This model then computes an
optimized inventory release policy for each of these sample sequences by using a variant of model
(10)-(11) that assumes perfect foreknowledge of the ship arrival process, as described below.
A tank overflow at the regasification terminal would occur in stage j with inventory level xj
when xj + uj − Q > X. To prevent this, one would have to avoid unloading an amount max{0, xj +
uj − Q − X} of LNG. To simplify the computation of our upper bound, we preprocess the sampled
sequences of unloaded LNG by subtracting from each relevant uj the quantity max{0, uj −Q}, which
is the amount of LNG that would be blocked in stage j if the terminal were full. Then, the upper
0
bound model values this amount of “blocked” LNG at maxj 0 ∈{j,j+1,...,J+1} δ j −j (1−φ)E[gj 0 (p̃j 0 ) | pj ]:
this LNG is allowed to be “virtually stored” for free and sold at the best possible price in the future
(the cost of receiving uj units of LNG is charged in stage j). Assuming perfect knowledge of all
future shipments and truncating allows us to replace the expectation with respect to ũ in (10)
with a degenerate expectation with respect to the deterministic quantity uj − max{0, uj − Q}. We
average the value functions in stage 1 and state (x1 , s1 , p1 ) over the generated sample paths to
obtain the estimate V̂1U B (x1 , s1 , p1 ) of an upper bound on V1 (x1 , s1 , p1 ). Finally, we estimate an
upper bound on the value of storage S1 (x1 , s1 , p1 ), defined by (9), as follows:

Ŝ1U B (x1 , s1 , p1 ) := V̂1U B (x1 , s1 , p1 ) − V1G (x1 , s1 , p1 ). (13)

As estimating the upper bound does not require simulating the evolution of the price state
vector, we do not estimate V1G (x1 , s1 , p1 ) in (13). Instead we take advantage of the linearity in
price of the revenue of the greedy policy in every stage, and the property that the time 0 futures
prices for the relevant stages (maturities), which are available at time 0, are the risk neutral expected
spot prices in these stages given the information available at time 0. Hence, we use these futures
prices directly to compute the value of the greedy policy in the initial stage and state.

16
7 Quantification of the Value of the Real Option to Store

In this section we report the results of a numerical study of the real option value of downstream
LNG storage. After introducing the setting of the study and the estimation of the parameters of
the natural gas price models used, we discuss our valuation results and related managerial insights.
Operational Parameters and Operating Costs. Table 1 summarizes the operational pa-
rameters and operating costs of the liquefaction, shipping, and regasification stages used in our
numerical experiments. We consider a cargo size of 145, 000CM, a common size in the LNG indus-
try (Flower [24]). We set the distance between the liquefaction and regasification facilities equal to
7,000 nautical miles (NM), which is roughly equal to the distance between Egypt, an LNG exporting
country, and Lake Charles, Louisiana, which hosts an LNG terminal operated by Trunkline LNG
(http://infopost.panhandleenergy.com/InfoPost/jsp/frameSet.jsp?pipe=tlng). We assume that the
speed of each ship is 19 knots, a realistic value (Flower [24, p. 100], Cho et al. [15]), which makes
a one way trip approximately 15 days long. The mean service times at the liquefaction and regasi-
fication facilities are one day each, which are representative of typical operations (EIA [19]). The
unloading (handling) charge and the regasification fuel loss are the “Currently Effective Rates” of
Trunkline LNG for firm terminal service as of 5/29/2009 for the unloading charge and 6/26/2006
for the regasification fuel loss (the 5/29/2009 rate sheet does not display a numerical figure for this
quantity, so we employ the analogous figure from the 6/26/2006 rate sheet).
According to EIA [19], two cargos is the industry rule of thumb for the size of the receiving
tanks. Since the Lake Charles terminal storage space amounts to roughly three cargos and some of
the newly developed terminals in the U.S. have even larger sizes, we also consider larger values for
this parameter. Moreover, for completeness we include one cargo storage size in our analysis. The
send out capacity is 2BCF/day, which is consistent with the 2.1BCF/day peak capacity of the Lake
Charles terminal (the capacity of the Sabine Pass, Texas, terminal, operated by Cheniere Energy,
is 2.6BCF/day; Cheniere Energy has also proposed two other terminals with capacity equal to
2.6BCF/day and 3.3BCF/day, respectively). Apparently, Trunkline LNG and the other companies
that manage the active regasification terminals in the U.S. do not charge a holding cost, so we set
this to zero in our experiments.
We consider fleet sizes ranging from 1 to 10 ships in unitary increments. Their throughput
levels vary from 0.7536MTPA to 7.5362MTPA, assuming that the ships are operated 365 days
per year with deterministic processing (loading, unloading, and transit) times; with exponentially
distributed processing times, ignoring ship blocking, these throughput figures are no more than

17
 Table 1: Operational parameters and costs.
Liquefaction Average Loading Time
1 Day
Shipping
Average One Way Transit Time Distance Speed Ship Size
15 Days 7,000NM 19 Knots 145,000CM
Regasification
Average Unloading Time Capacity Fuel Loss Unloading Cost
1 Day 2BCF/Day 1.69% $0.0017/MMBTU

NYMEX Natural Gas Futures Prices

9
8.5

7.5
7
$/MMBTU

6.5
6

5.5

5
4.5
4

3.5
0 20 40 60 80 100 120 140 160
Maturity

Figure 3: NYMEX natural gas futures prices traded on 5/29/2009.

2.22% lower, indicating a very low level of congestion in the systems that we analyze.
Estimation of the Parameters of the Natural Gas Price Models. We need to estimate
the following parameters pertaining to the natural gas price models (1)-(3) and (4): κχ , λχ , σχ ,
µ∗ξ := µξ − λξ , σξ , ρχξ , χ0 , and ξ0 , where the latter two are initial values for the two factors, and
f1 , . . . , f12 . We assume that regasified LNG is sold into the Louisiana wholesale natural gas spot
market at the Henry Hub price. Thus, we use NYMEX data for estimation purposes and employ
a dataset that includes natural gas futures prices and prices of call and put options on natural gas
futures from 5/29/2009. We now describe our estimation approach.
Consider model (1)-(3). Denote by F (t, t0 ) the time t price of a futures contract for delivery at
time t0 > t. Under this model, ln F (t, t0 ) can be expressed as

0 0 λχ σ̌ 2 (t, t0 )
ln F (t, t0 ) = ln fm(t0 ) + e−κχ (t −t) χt + ξt + µ∗ξ (t0 − t) − [1 − e−κχ (t −t) ] + , (14)
κχ 2

18
 Black Implied Volatilities of NYMEX Options on Natural Gas Futures Prices
0.7

0.65

0.6

0.55

0.5

0.45

0.4

0.35

0.3

0.25
0 5 10 15 20 25 30 35 40 45
Maturity

Figure 4: Black implied volatilities of NYMEX options on natural gas futures prices traded on
5/29/2009.

0 σχ2 0 ρχξ σχ σξ
σ̌ 2 (t, t0 ) := [1 − e−κχ (t −t) ] + σξ2 (t0 − t) + 2[1 − e−2κχ (t −t) ] . (15)
2κχ κχ
Closed form expressions for the time 0 prices of European call and put options on futures price
F (t, t0 ) under model (1)-(3) depend on (15). Given the time 0 market price of a traded (call or put)
European option on a futures and the futures price, one can use the well known Black [9] formula
for the option price to compute, by means of standard techniques (Eydeland and Wolyniec [23, pp.
147-150]), the so called implied volatility parameter σ̂B , where subscript B stands for Black. We
numerically compute a Black volatility for each option in our dataset. Under model (1)-(3) the
price of a European call/put option on a futures price would match the market price of the traded
2 = σ̌ 2 (t, t0 )/t.
option if σ̂B
Figures 3 and 4 illustrate the futures prices and Black implied volatilities in our dataset. The
marked seasonality in futures prices is notable. The overall decline of the volatilities with increasing
futures price maturity is typical; it is known as the Samuelson [42] effect. Following Clewlow and
Strickland [16, p. 160], we estimate the price model parameters by minimizing the sum of the
squared percent deviations of observed futures prices and implied volatilities subject to appropriate
P
nonnegativity constraints and the constraint 12m=1 ln fm = 0 (this is a normalization condition that

uses a yearly cycle with constant seasonality factors within each month as in Jaillet et al. [30]).
Table 2 displays the relevant estimates. The root mean squared errors (RMSEs) between the
observed and estimated futures prices and implied volatilities are 0.0584 and 0.0280, respectively.

19
Table 2: Estimates of the parameters of the natural gas price models using NYMEX data from
5/29/2009.
Two factor model One factor model
Initial short term level (χ0 ) 0.7417 Initial level (χ0 ) 1.2721
Speed of mean reversion (κχ ) 1.5245 Speed of mean reversion (κ) 1.0547
Short term volatility (σχ ) 0.7388 Volatility (σ) 0.6696
Short term risk premium (λχ ) −2.1219 Risk adjusted long term level (ξ ∗ ) −2.0421
Initial long term level (ξ0 ) 0.4724
Risk adjusted long term drift (µ∗ξ ) 0.0038
Long term volatility (σξ ) 0.1300
Correlation (ρχξ ) −0.0886
January factor (f1 ) 1.0826 January factor (f1 ) 1.0809
February factor (f2 ) 1.0773 February factor (f2 ) 1.0761
March factor (f3 ) 1.0441 March factor (f3 ) 1.0433
April factor (f4 ) 0.9537 April factor (f4 ) 0.9532
May factor (f5 ) 0.9467 May factor (f5 ) 0.9464
June factor (f6 ) 0.9547 June factor (f6 ) 0.9547
July factor (f7 ) 0.9693 July factor (f7 ) 0.9674
August factor (f8 ) 0.9717 August factor (f8 ) 0.9711
September factor (f9 ) 0.9690 September factor (f9 ) 0.9697
October factor (f10 ) 0.9750 October factor (f10 ) 0.9767
November factor (f11 ) 1.0133 November factor (f11 ) 1.0156
December factor (f12 ) 1.0565 December factor (f12 ) 1.0590

We estimate the relevant parameters of model (4) in a manner analogous to the estimation of those
of model (1)-(3). Table 2 reports also these estimates. The RMSEs of the futures prices and implied
volatilities for this estimated model are 0.2266 and 0.0332, respectively. Thus, using two factors,
instead of one, yields a more accurate fit of the data.
Valuation Results. We employ a valuation period of twelve years divided into monthly time
periods. Although a time horizon of twenty years would be more in line with industry practices
regarding the valuation of LNG projects (Flower [24]), our choice stems on the fact that the NYMEX
natural gas forward curve spans approximately twelve years with monthly maturities. However, if
one were willing to extend the forward curve beyond this time horizon our model could be applied
to this longer time horizon at the expense of additional run time. Thus, we set J := 143. We let
the time 0 price of natural gas be equal to the Henry Hub spot price traded on 5/29/2009, which
is $3.92/MMBTU. We use an annual risk free rate equal to 0.47%, the one year U.S. treasury rate
on 5/29/2009. We build the lattices described in §5 using the parameter estimates shown in Table
2 and the futures prices displayed in Figure 3. At time t1 no inventory is available in storage and
all the ships are in the ballast stage; we set x1 and s1 accordingly.
The software implementation of our model features inventory and action sets expressed in num-

20
 Value of Storage
800
X = 1C
X = 2C
X = 3C
700 X = 4C
X = 5C
X = 6C
X = 7C
600
X = 8C

500
$M

400

300

200

100

0
1 2 3 4 5 6 7 8 9 10
Fleet Size

Figure 5: The estimated value of the option to store, Ŝ1B (x1 , s1 , p1 ), using the two factor price
model and exponentially distributed processing times ($M).

ber of cargos. This is justified by Proposition 3, whose assumptions are satisfied by the values of
the relevant parameters used in this section. The Cpu times reported below pertain to computa-
tions performed on a 64 bits Monarch Empro 4-Way Tower Server with four AMD Opteron 852
2.6GHz processors, each with eight DDR-400 SDRAM of 2 GB and running Linux Fedora 11 (all
the reported results were obtained using only one processor). Our model was coded in C++ and
compiled using the compiler g++ version 4.3.0 20080428 (Red Hat 4.3.0-8).
We first discuss the valuation results obtained using the two factor price model and a shipping
model with exponentially distributed processing times (Table 5 in Online Appendix C reports the
probability mass functions of the number of unloaded cargos per aggregate time period computed
by the rolling forward method of §5 for different fleet sizes; computing this table takes less than 1
Cpu second). We compare these with results obtained with an essentially deterministic shipping
system when discussing the shipping model effect, and with those obtained with the one factor
model in our discussion of the price model effect.
The value of storage. Figure 5 displays the estimated value of the option to store for different
fleet sizes and levels of storage space. The system configurations that we consider satisfy our
assumption made in §5 that Pr{ũ > U } = 0 holds. The estimates of the value of storage, obtained
using 500,000 price and unloaded cargo sample paths, vary from $89M (1 ship and 1 cargo of

21
 Value of Storage
800
X = 1C
X = 2C
X = 3C
700 X = 4C
X = 5C
X = 6C
X = 7C
600
X = 8C

500
$M

400

300

200

100

0
1 2 3 4 5 6 7 8 9 10
Fleet Size

Figure 6: The estimated value of the option to store, Ŝ1B (x1 , s1 , p1 ), using the two factor price
model and the essentially deterministic shipping system when Q = 1BCF/day ($M).

storage space) to $726M (10 ships and 8 cargos of storage space), and their relative standard errors
are 0.06% (we refer to a standard error expressed as a fraction of the estimate it pertains to as a
relative standard error). The Cpu times required to compute these values vary from approximately
5 to 53 minutes; these times also include the estimation of the value of storage due to seasonality,
discussed below, using 500,000 samples.
With the two factor price model, we use 1,000 unloaded cargo sample paths to estimate the
upper bound, which takes between 61 and 144 minutes of Cpu time depending on the instance.
We reject the hypothesis that S1U B (x1 , s1 , p1 ) is equal to S1B (x1 , s1 , p1 ) in favor of the alternative
hypothesis that S1U B (x1 , s1 , , p1 ) is greater than S1B (x1 , s1 , p1 ) at the 5% significance level in 50 out
of 80 combinations of fleet size and storage space. Moreover, the estimates of the value of storage
are never lower than 99.11% of their upper bound estimates; on average this figure is 99.73%. These
results suggest that our estimates of the value of storage are of very high quality, and reflect the
very low ship blocking probability and low dependence between unloaded amounts in successive
aggregate time periods in the systems that we consider.
Throughput and space effects. Figure 5 illustrates that the value of storage increases in the
available space. This is intuitive, as more available space allows more effective exploitation of high
natural gas prices. In addition, as one would expect, the marginal benefit of additional storage is

22
Table 3: Percent improvements of the basestock target policy relative to the myopic policy, using
the two factor price model and exponentially distributed processing times.
Storage Size (# of Cargos)
# of Ships 1 2 3 4 5 6 7 8
1 0.36 2.27 5.16 8.60 12.42 16.60 21.14 26.05
2 0.03 0.26 0.95 2.00 3.31 4.83 6.51 8.29
3 0.01 0.04 0.21 0.60 1.18 1.89 2.73 3.67
4 0.00 0.00 0.05 0.18 0.45 0.84 1.30 1.84
5 0.00 0.00 0.01 0.05 0.17 0.38 0.66 0.99
6 0.00 0.00 0.00 0.01 0.06 0.16 0.33 0.56
7 0.00 0.00 0.00 0.00 0.02 0.06 0.16 0.31
8 0.00 0.00 0.00 0.00 0.01 0.02 0.06 0.15
9 0.00 0.00 0.00 0.00 0.00 0.01 0.02 0.06
10 0.00 0.00 0.00 0.00 0.00 0.00 0.01 0.03

decreasing, as each additional unit of space is used less often. Although not immediately evident
from Figure 5, it is more interesting that the value of storage can be nonmonotonic in the number of
ships, or, equivalently, throughput: with 1 cargo and 2 cargos of storage space the value of storage
decreases from $91.3052M to $91.3044M and from $182.607M to $182.601M when the number of
ships increases from 9 to 10. These drops in storage value are minimal, but Figure 6 more forcefully
illustrates this possibility (we obtain the values displayed in this figure by using the essentially
deterministic shipping model described below when discussing the shipping model effect). This
suggests that, due to regasification capacity constraints, after a critical level of throughput the
value of storage decreases. Intuitively, as throughput grows large, forced sales approach the system
capacity. As optional sales shrink, so does the ability to store, and consequently the value of storage
(see Online Appendix C for a more detailed justification of this statement). The fundamental insight
here is that discretionary regasification capacity must be available for an LNG terminal to have
storage value. Otherwise, as throughput approaches the regasification capacity, an LNG terminal
becomes a delivery mechanism with little or no storage value. This suggests that the discussion in
Holcomb [29] (see §1) is pertinent to situations that satisfy this condition.
The storage values due to price seasonality and volatility. There is almost an equal split between
the storage values due to price seasonality and volatility, with the former varying between 46% and
51% of the total option value (the relative standard error of the value due to seasonality is no more
than 0.0186%). Considering the pronounced seasonality displayed by the natural gas forward curve
(see Figure 3), it is remarkable that more value is not attributable to seasonality, or, put another
way, that there is such significant value in adapting the inventory release policy to price volatility.
Forward looking optimization effect. As mentioned in §5, an optimal basestock target policy

23
can be nontrivial, that is, the basestock targets may take values different from 0 and X. This is
a consequence of the forward looking optimization nature of this policy. To assess the relevance
of this feature of this policy, we compare the percent improvements of this policy relative to a
myopic policy that in every stage sells as much as possible if the difference between the discounted
price of the futures with maturity in the next stage and the spot price is positive, and holds as
much inventory as possible otherwise (except for the forced sales); this is a basestock target policy
whose targets are set naı̈vely either equal to 0 or X in every stage. Computing the myopic policy
does not require solving an MDP. Table 3 reports these figures using the two factor price model
and the exponential shipping system. The basestock target policy outperforms the myopic policy
by no more than 26.05%, 8.29%, 3.67%, 1.84%, and 0.99% with fleets growing from 1 to 5 ships,
respectively, and by less than 1% with 6 or more ships. For each fleet size, these improvements
increase (weakly) in the available storage space. These results suggest that the simple myopic
policy seems adequate for storage valuation when the throughput is sufficiently high relative to
the storage space. Otherwise, there appears to be significant value from optimizing the inventory
release policy by taking into account the entire future consequences of a current action.
Shipping model effect. As pointed out by Koenigsberg and Lam [33], exponential processing
times, with a coefficient of variation equal to 1, are unrealistic, especially for shipping times whose
averages are of the order of two or more weeks. Since our CQN model is a BCMP network, we could
reduce the variability of the transit times by increasing the number of Erlang stages in the two
shipping blocks (exponential times correspond to the case of a single Erlang stage). For simplicity,
we consider the extreme case in which all the processing times are deterministic and equal to their
means. In this case the throughput of the LNG chain, expressed in number of cargos per day, is
equal to the fleet size divided by the sum of the average processing times.
In computing the value of storage in this case, we modify our deterministic assumption slightly
to satisfy the conditions of Proposition 3 in §5: we use a two point unloading random variable
that takes values equal to the floor and the ceiling of the throughput with probabilities determined
to make the mean of this random variable equal to this throughput. We refer to this case as the
essentially deterministic system. For consistency, we continue to use 500,000 samples to estimate
the value of storage and that due to price seasonality; this takes between 1.4 and 1.7 Cpu minutes.
The estimate of the value of storage obtained under the exponential assumption is at least 98%
of the value of storage computed in the essentially deterministic case (the relative standard errors of
the value of storage in the latter case are roughly 0.06%); on average this figure is 99.37%. That the
value of storage is lower in the former case is intuitive, since an essentially deterministic shipping

24
system allows for easier planning of released inventory. That this value is not dramatically lower
stems from the low ship blocking probabilities and congestion in the exponential shipping system.
Similar to the exponential case, the estimated values due to price seasonality approximately range
from 47% to 51% of the option value; their relative standard errors are no more than 0.0031%.
Thus, a simple version of the shipping model may be adequate to value storage for the range of
parameters that we consider.
Price model effect. We reexamine all the previous comparisons, but now for the one factor
price model. In this case the estimates of the value of storage are between 83% and 84% of those
obtained using the two factor price model, both with the exponential and essentially deterministic
shipping systems (the relative standard errors with 500,000 samples are 0.05% for both systems).
Hence, the value of storage due to price seasonality, which is the same with the one and two factor
models, is relatively more important with a single factor, varying between 55% and 62% of the
value of storage (the relative standard errors of the values due to price seasonality are no more
than 0.0189% in the exponential case and 0.0031% in the essentially deterministic case).
Estimating the value of storage is less computationally demanding with a one factor model, and
so is estimating our upper bound. Thus, we also estimate the upper bound using 500,000 samples
of unloaded cargo sequences, rather than using 1,000 samples as in the two factor model case. With
the exponential shipping system it takes from about 2 to 38 Cpu minutes to run our valuation model
and from about 7 to 28 Cpu minutes to compute the upper bound; in the essentially deterministic
shipping system it takes about 26 Cpu seconds to run our valuation model. For the exponential
shipping system, we reject the hypothesis that S1U B (x1 , s1 , p1 ) is equal to S1B (x1 , s1 , p1 ) in favor
of the alternative hypothesis that S1U B (x1 , s1 , , p1 ) is greater than S1B (x1 , s1 , p1 ) at the 5% level of
significance in 53 out of 80 cases; the estimated value of storage is at least 99% of its upper bound
and on average it is 99.71% of this value. The value of storage so computed is at least 97.84% of
that obtained with the essentially deterministic system and on average it is 99.33% of this value.
These results suggest that a quicker but appreciably lower estimate of the value of storage can
be obtained by using a one factor model, in which case our integrated model continues to deliver
very high quality storage valuations relative to our upper bound estimate.

8 Conclusions

Motivated by current developments in the LNG industry, we develop a real option model for the
strategic valuation of downstream LNG storage. Unique to our model is the integration of models
of natural gas liquefaction and LNG shipping, natural gas price evolution, and LNG inventory

25
regasification and sale into the wholesale spot market. This provides a heuristic strategic valuation
of the real option to store LNG at a regasification terminal. We apply our model to real and
estimated data and find these valuations to be highly accurate. We also investigate how these
values depend on the level of stochastic variability in the shipping model, the type of natural gas
price model used, the LNG throughput, and the type of inventory control policy employed.
Our model has the potential to be applied in practice as it is both computationally manageable
and highly accurate. It could be used by LNG players to assess the value of leasing contracts
on regasification facilities, or as an economic valuation model by different parties involved in the
development of LNG projects (see the discussion in Flower [24, p. 120]). Moreover, while our
focus has been on LNG, our model and analysis have potential applicability in other commodity
industries that exhibit uncertainty in the commodity production or shipping processes, for example
those characterized by random yield and/or spot price fluctuations.
For further research, one could assess the dependence of our valuation results on the type of
multifactor model used to represent the evolution of the price of natural gas. In particular, it would
be interesting to study the case when this evolution is captured using an equilibrium model, such
as that of Routledge et al. [41], or reduced form models with more than two factors, such as those
presented by Schwartz [43], Casassus and Collin-Dufresne [13], and Lai et al. [34].
Another additional research area is the development of approximate dynamic programming
algorithms (Bertsekas [7, Chapter 6], Adelman [2], Powell [39]) for solving the operational model
discussed in §3, perhaps extended to include the status of the liquefaction facility and the inventory
available at this location. This would allow one to obtain a policy for tactical and operational
control. The work of Besbes and Savin [8] is pertinent here.

Acknowledgments

This research was conducted under CART grants from the Carnegie Mellon Tepper School of
Business. We thank Christine Parlour for insightful discussions and the Area Editor, Bert Zwart,
and the anonymous review team for constructive feedback that led to significant improvements to
this paper.

References

[1] L. M. Abadie and J. M. Chamorro. Monte Carlo valuation of natural gas investments. Review
of Financial Economics, 18:10–22, 2009.

26
 [2] D. Adelman. Math programming approaches to approximate dynamic programming. Tutorial,
INFORMS Annual Meeting, Pittsburgh, PA, USA, 2006.

[3] S. Asmussen. Applied Probability and Queues. Springer, New York, NY, USA, 2003.

[4] C. H. Bannister and R. J. Kaye. A rapid method for optimization of linear systems with
storage. Operations Research, 39:220–232, 1991.

[5] F. Baskett, K. M. Chandy, R. R. Muntz, and F. G. Palacios. Open, closed, and mixed networks
of queues with different classes of customers. Journal of the ACM, 22:248–260, 1975.

[6] R. Baxter. Energy Storage: A Nontechnical Guide. PennWell Corporation, Tulsa, OK, USA,
2006.

[7] D. P. Bertsekas. Dynamic Programming and Optimal Control – Volume I. Athena Scientific,
Belmont, MA, USA, third edition, 2000.

[8] O. Besbes and S. Savin. Going bunkers: The joint route selection and refueling problem.
Manufacturing & Service Operations Management, Articles in Advance, 2009.

[9] F. Black. The pricing of commodity contracts. Journal of Financial Economics, 3:167–179,
1976.

[10] A. Boogert and C. de Jong. Gas storage valuation using a Monte Carlo method. Journal of
Derivatives, 15:81–98, 2008.

[11] D. B. Brown, J. E. Smith, and P. Sung. Information relaxation and duality in stochastic
dynamic programs. Operations Research, Forthcoming, 2009.

[12] R. Carmona and M. Ludkovski. Valuation of energy storage: An optimal switching approach.
Working Paper, Princeton University, Princeton, NJ, USA, 2007.

[13] J. Casassus and P. Collin-Dufresne. Stochastic convenience yield implied from commodity
futures and interest rates. Journal of Finance, 60:2283–2331, 2005.

[14] Z. Chen and P. A. Forsyth. A semi-Lagrangian approach for natural gas storage valuation and
optimal operation. SIAM Journal on Scientific Computing, 30:339–368, 2007.

[15] J. H. Cho, H. Kotzot, F. de la Vega, and C. Durr. Large LNG carrier poses economic advan-
tages, technical challenges. LNG Observer, 2, October 2005. On-line supplement to Oil & Gas
Journal.

27
[16] L. Clewlow and C. Strickland. Energy Derivatives: Pricing and Risk Management. Lacima
Publications, London, UK, 2000.

[17] N. Drouin, A. Gautier, B. F. Lamond, and P. Lang. Piecewise affine approximations for the
control of a one-reservoir hydroelectric system. European Journal of Operational Research,
89:53–69, 1996.

[18] D. Duffie. Dynamic Asset Pricing Theory. Princeton University Press, Princeton, NJ, USA,
2001.

[19] Energy Information Administration. U.S. Natural Gas Markets: Mid-Term Prospects for Nat-
ural Gas Supply. U.S. Department of Energy, Washington, DC, USA, 2001.

[20] Energy Information Administration. The Global Liquefied Natural Gas Market: Status &
Outlook. U.S. Department of Energy, Washington, DC, USA, 2003.

[21] Energy Information Administration. Natural Gas, in International Energy Outlook 2006. U.S.
Department of Energy, Washington, DC, USA, 2006.

[22] Energy Information Administration. Annual Energy Outlook 2009 with Projections to 2030.
Report #: DOE/EIA-0383(2009), U.S. Department of Energy, Washington, DC, USA, March
2009.

[23] A. Eydeland and K. Wolyniec. Energy and Power Risk Management. John Wiley & Sons,
Inc., Hoboken, NJ, USA, 2003.

[24] A. R. Flower. LNG project feasibility. In G. B. Greenwald, editor, Liquefied Natural Gas:
Developing and Financing International Projects, pages 73–124. Kluwer Law International,
London, UK, 1998.

[25] H. Geman. Commodities and Commodity Derivatives: Modeling and Pricing for Agriculturals,
Metals and Energy. John Wiley & Sons, Chichester, UK, 2005.

[26] B. Greenwald, G., editor. Liquefied Natural Gas: Developing and Financing International
Projects. Kluwer Law International, London, UK, 1998.

[27] R. Grønhaug and M. Christiansen. Supply chain optimization for the liquefied natural gas
business. In L. Bertazzi, M. G. Speranza, and J. A. E. E. van Nunen, editors, Innovation in
distribution logistics, volume 619 of Lecture Notes in Economics and Mathematical Systems,
pages 195–218. Springer, Berlin, Germany, 2009.

28
[28] J. M. Harrison and M. I. Taksar. Instantaneous control of Brownian motion. Mathematics of
Operations Research, 8:439–453, 1983.

[29] J. H. Holcomb. Storage downstream of regasification optimizes LNGs value. LNG Observer,
3, July 2006. On-line supplement to Oil & Gas Journal.

[30] P. Jaillet, E. I. Ronn, and S. Tompaidis. Valuation of commodity-based swing options. Man-
agement Science, 50:909–921, 2004.

[31] J. T. Jensen. The LNG revolution. The Energy Journal, 24:1–45, 2003.

[32] M. Kaplan, R. C. Wentworth, and R. J. Hischier. Simulation and optimization of LNG shipping
systems. In ASME Transaction Petroleum Mechanical Engineering and Pressure Vessels and
Piping Conference, pages 1–11, New Orleans, LA, USA, September 17-21 1972.

[33] E. Koenigsberg and R. C. Lam. Cyclic queue models of fleet operations. Operations Research,
24:516–529, 1976.

[34] G. Lai, F. Margot, and N. Secomandi. An approximate dynamic programming approach to

benchmark practice-based heuristics for natural gas storage valuation. Operations Research,
Forthcoming, 2009.

[35] B. F. Lamond, S. L. Monroe, and M. J. Sobel. A reservoir hydroelectric system: Exactly and
approximately optimal policies. European Journal of Operational Research, 81:535–542, 1995.

[36] E. Näsäkkälä and J. Keppo. Hydropower with financial information. Applied Mathematical
Finance, 15:503–529, 2008.

[37] M. H.-B. Osogami, T. A closed-form solution for mapping general distributions to minimal
PH distributions. Performance Evaluation, 63:524–552, 2006.

[38] E. C. Özelkana, A. DAmbrosio, and S. G. Teng. Optimizing liquefied natural gas terminal
design for effective supply-chain operations. International Journal of Production Economics,
111:529–542, 2009.

[39] W. B. Powell. Approximate Dynamic Programming: Solving the Curses of Dimensionality.

Wiley-Interscience, Hoboken, NJ, USA, 2007.

[40] D. Ronen. Maritime transportation of oil and gas: The role of quantitative methods and
information systems in improving profitability. In A. Poirier and G. Zaccour, editors, Maritime

29
 and Pipeline Transportation of Oil and Gas, pages 236–246. Éditions Technip, Paris, France,
1991.

[41] B. Routledge, D. J. Seppi, and C. Spatt. Equilibrium forward curves for commodities. Journal
of Finance, 55:1297–1338, 2000.

[42] P. Samuelson. Proof that properly anticipated prices fluctuate randomly. Industrial Manage-
ment Review, 6:41–49, 1965.

[43] E. S. Schwartz. The stochastic behavior of commodity prices: Implications for valuation and
hedging. Journal of Finance, 52:923–973, 1997.

[44] E. S. Schwartz and J. E. Smith. Short-term variations and long-term dynamics in commodity
prices. Management Science, 46:893–911, 2000.

[45] N. Secomandi. Optimal commodity trading with a capacitated storage asset. Management
Science, Forthcoming, 2009.

[46] J. E. Smith. Alternative approaches for solving real-options problems. Decisions Analysis,
2:89–102, 2005.

[47] J. E. Smith and K. F. McCardle. Options in the real world: Lessons learned in evaluating oil
and gas investments. Operations Research, 47:1–15, 1999.

[48] M. Thompson, M. Davison, and H. Rasmussen. Natural gas storage valuation and optimization:
A real options application. Naval Research Logistics, 56:226–238, 2009.

[49] L. Trigeorgis. Real Options: Managerial Flexibility and Strategy in Resource Allocation. The
MIT Press, Cambridge, MA, USA, 1996.

[50] C. L. Tseng and K. Y. Lin. A framework using two-factor price lattices for generation asset
valuation. Operations Research, 55:234–251, 2007.

[51] M. D. Tusiani and G. Shearer. LNG: A Nontechnical Guide. PennWell Corporation, Tulsa,
OK, USA, 2007.

[52] S. W. Wallace and S.-E. Fleten. Stochastic programming models in energy. In A. Ruszczyn-
ski and A. Shapiro, editors, Stochastic programming, volume 10 of Handbooks in Operations
Research and Management Science, pages 637–677. Elsevier, Amsterdam, The Netherlands,
2003.

30
 Online Appendix
A Units of Measurement and Conversion Factors

Table 4 reports relevant units of measurements and conversion factors.

B Proofs for §5

Proof of Proposition 1 (Optimal value function). By induction. The claimed property

clearly holds in stage J + 1. Make the induction hypothesis that this property holds also in stages
H (x
j + 1, . . . , J. Consider stage j. The function Vj+1 j+1 , pj+1 ) is concave in xj+1 ∈ X for each
H (x
given pj+1 ∈ Pj+1 . This implies that E[Vj+1 j+1 , p̃j+1 ) | pj ] is concave in xj+1 ∈ X for each given

pj ∈ Pj . Fix pj ∈ Pj , u ∈ U, x1j , x2j ∈ X with x1j 6= x2j , and qj1 , qj2 ∈ Q(xj , u) with qj1 6= qj2 . Define
x1j+1 := x1j + u − qj1 , x2j+1 := x2j + u − qj2 , and xθj+1 := θx1j+1 + (1 − θ)x2j+1 for some θ ∈ [0, 1]. Since
x1j+1 , x2j+1 ∈ X , the convexity of X implies that xθj+1 ∈ X . The concavity of E[Vj+1
H (x
j+1 , p̃j+1 ) | pj ]

in xj+1 for given pj and the linearity of expectation imply that

H
E[Vj+1 (xθj+1 , p̃j+1 ) | pj ] > θE[Vj+1
H
(x1j+1 , p̃j+1 ) | pj ] + (1 − θ)E[Vj+1
H
(x2j+1 , p̃j+1 ) | pj ].

Thus, the definitions of x1j+1 and x2j+1 imply that E[Vj+1

H (x + u − q , p̃
j j j+1 ) | pj ] is jointly concave in
H (x + u − q , p̃
xj and qj for given pj and u, and so is gj (pj )(1 − φ)qj − hxj − cu + δE[Vj+1 j j j+1 ) | pj ].

This property, the convexity of set A(u) := {(x, q) : x ∈ X , q ∈ Q(x, u)}, and Proposition B-4
in Heyman and Sobel [53, p. 525] imply that vjH (xj , pj , u) is concave in xj for given pj and u.
The concavity of VjH (xj , pj ) in xj for given pj follows since VjH (xj , pj ) = E[vjH (xj , pj , ũ)], and the
property holds in all stages by the principle of mathematical induction. ¤
Proof of Proposition 2 (Optimal policy structure). Consider arbitrary stage j ∈ J . Since
finding an optimal action is equivalent to finding an optimal optional sale after having performed

Table 4: Units of measurement and conversion factors.

T Metric Tons (LNG)
MTPA Million Metric Tons per Annum
CM Cubic Meters (LNG)
NM Nautical Miles
MMBTU Million British Thermal Units
BCF Billion Cubic Feet
1T = 51.98237MMBTU
1Knot = 1NM per Hour
1CM = 23.6863MMBTU
1BCF = 1,100,000MMBTU

OA-1
the forced sale, we denote by yj (xj , u) the post forced sale inventory level given xj and u:
½
X xj + u ∈ (X, X + U ],
yj (xj , u) :=
xj + u xj + u ∈ [0, X].

This inventory level can only take values in set X . The costs hxj and cuj do not affect the choice
of an optimal optional sale in stage j, and we can restrict our attention for this purpose to state
(yj , pj ) ∈ X ×Pj . To find an optimal optional sale in this state, we first consider the relaxed problem
of finding an optimal optional sale by ignoring the capacity restriction. Thus, the feasibility set is
simply equal to [0, yj ] for each yj ∈ X , and the problem to be solved is maxqO ∈[0,yj ] νj (yj , qjO , pj ),
j

where νj (yj , qjO , pj ) := gj (pj )(1 − φ)qjO + H (y

δE[Vj+1 j − qjO , p̃j+1 ) | pj ].
Define q̌jO (yj , pj ) as the largest quantity in set arg maxqO ∈[0,yj ] νj (yj , qjO , pj ). Following Porteus
j

[54, p. 67], by letting xj+1 = yj − qjO , it clearly holds that

max νj (yj , qjO , pj ) = max νj (yj , yj − xj+1 , pj )

qjO ∈[0,yj ] xj+1 ∈[0,yj ]

= gj (pj )(1 − φ)yj

© H
ª
+ max δE[Vj+1 (xj+1 , p̃j+1 ) | pj ] − gj (pj )(1 − φ)xj+1 .
xj+1 ∈[0,yj ]

Notice that any element of arg maxxj+1 ∈X νj (yj , yj − xj+1 , pj ), and so its smallest one denoted
by x̌j+1 (pj ), does not depend on yj . In particular, it holds that
½
O 0 yj ∈ [0, x̌j+1 (pj )],
q̌j (yj , pj ) =
yj − x̌j+1 (pj ) yj ∈ (x̌j+1 (pj ), X].

This implies that we can define bj (pj ) := x̌j+1 (pj ). Imposing the capacity constraint on optimal
∗
optional sale q̌jO (yj , pj ) yields the stated expression for qjO (xj , pj , uj ). ¤
Proof of Proposition 3 (Optimal policy computation). By induction. The claimed
properties clearly hold in stage J + 1. Make the induction hypothesis that they hold also in stages
j + 1, . . . , J. Consider stage j. Fix pj ∈ Pj and u ∈ U. Recall that set Pj+1 is finite. This
H (x
and the induction hypothesis imply that, given pj , δE[Vj+1 j+1 , p̃j+1 ) | pj ] is piecewise linear

and continuous in xj+1 ∈ X and changes slope in xj+1 only at integer multiples of L. It is
H (X + u − q , p̃
easy to show that, given pj and u, gj (pj )(1 − φ)qj − hX − cu + δE[Vj+1 j j+1 ) | pj ]

changes slope in qj only at integer multiples of L, which implies that bj (pj ) is an integer multiple
of L. It is now shown that the function VjH (xj , pj ) is piecewise linear and continuous in xj ∈ X
and changes slope in xj only at integer multiples of L. By Proposition 2, the following cases
need to be considered: (H) xj + u ∈ [0, bj (pj )) and (S) xj + u ∈ [bj (pj ), X + U ]. Case (H):
vjH (xj , pj , u) = −cu − hxj + δE[Vj+1
H (x + u, p̃
j j+1 ) | pj ]. Case (S): if bj (pj ) can be reached from

OA-2
Table 5: Probability mass functions of the number of unloaded cargos in one aggregate time
period (1 month) computed by the rolling forward method for the shipping model with exponential
processing times.
Unloaded # of Ships in the Fleet
Cargos 1 2 3 4 5 6 7 8 9 10
0 0.2809 0.0791 0.0223 0.0063 0.0018 0.0005 0.0002 0.0000 0.0000 0.0000
1 0.5216 0.2935 0.1241 0.0468 0.0166 0.0057 0.0019 0.0006 0.0003 0.0001
2 0.1776 0.3726 0.2728 0.1462 0.0672 0.0282 0.0111 0.0042 0.0016 0.0006
3 0.0189 0.1956 0.3037 0.2514 0.1565 0.0829 0.0395 0.0176 0.0074 0.0030
4 0.0010 0.0512 0.1877 0.2627 0.2332 0.1609 0.0947 0.0501 0.0245 0.0113
5 0.0000 0.0074 0.0699 0.1762 0.2346 0.2181 0.1621 0.1037 0.0595 0.0316
6 0.0000 0.0006 0.0166 0.0792 0.1649 0.2136 0.2054 0.1618 0.1105 0.0680
7 0.0000 0.0000 0.0026 0.0248 0.0834 0.1548 0.1971 0.1946 0.1604 0.1158
8 0.0000 0.0000 0.0003 0.0055 0.0310 0.0846 0.1458 0.1836 0.1852 0.1586
9 0.0000 0.0000 0.0000 0.0008 0.0087 0.0355 0.0842 0.1376 0.1722 0.1769
10 0.0000 0.0000 0.0000 0.0001 0.0018 0.0116 0.0385 0.0828 0.1302 0.1624
11 0.0000 0.0000 0.0000 0.0000 0.0003 0.0030 0.0141 0.0403 0.0807 0.1234
12 0.0000 0.0000 0.0000 0.0000 0.0000 0.0005 0.0041 0.0160 0.0413 0.0782
13 0.0000 0.0000 0.0000 0.0000 0.0000 0.0001 0.0010 0.0052 0.0175 0.0415
14 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0002 0.0014 0.0062 0.0186
15 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0001 0.0004 0.0019 0.0070
16 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0001 0.0005 0.0023
17 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0001 0.0006
18 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0001

xj + u then vjH (xj , pj , u) = gj (pj )(1 − φ)[xj + u − bj (p)] − cu − hxj + δE[Vj+1

H (b (p ), p̃
j j j+1 ) | pj ],

otherwise vjH (xj , pj , u) = gj (pj )(1 − φ)Q − cu − hxj + δE[Vj+1

H (x + u − Q, p̃
j j+1 ) | pj ]. It is an easy

task to verify that, for given pj and u, vjH (xj , pj , u) is piecewise linear and continuous in xj and
changes slope in xj only at integer multiples of L. Hence, VjH (xj , pj ) satisfies this property because
it is a convex combination of a finite number of functions that also satisfy this property. By the
principle of mathematical induction, the claimed properties hold in every stage. ¤

C Supporting Material for §7

Table 5. Table 5 displays the probability mass functions of the number of unloaded cargos per
aggregate time period (1 month) computed by the rolling forward method when the loading, un-
loading, and shipping times are exponentially distributed according to the parameters displayed in
Table 1.
Explanation of the Decreasing Value of Storage for High Throughput. Consider a
degenerate unloading random variable equal to u. Momentarily, impose the constraint that the
inventory level at time tJ+1 ≡ T be zero, i.e., xJ+1 = 0. Suppose that u = Q, so that the LNG rate

OA-3
into the terminal is equal to the maximum rate out of it. In this case, the value of storage must
be zero because no amount of LNG can be stored in any aggregate time period. Thus, when u is
sufficiently high, as it approaches Q from below, the value of storage decreases to zero.
Now, remove the constraint xJ+1 = 0, so that the only conditions imposed on xJ+1 are 0 6
xJ+1 6 X. Finally, make the realistic assumption that QJ > X, that is, a full terminal at time
t1 ≡ 0 can be emptied by time T . If u = Q, any amount of LNG not released in some aggregate
time period j ∈ J must be stored until time T , and the maximum amount of stored LNG during
the entire planning horizon is min{uJ, X} = min{QJ, X} = X. Thus, the value of storage for
any level of throughput that allows one to store at least an amount X of LNG during the entire
planning horizon (that is, for any u 6 Q such that uJ > X) must be at least the one obtainable
when u = Q. In other words, u = Q is the level of throughput that minimizes the value of storage
among all those that satisfy uJ > X. Therefore, as u > X/J approaches Q, obviously from below,
the value of storage decreases in a neighborhood of Q.
When the unloading random variable is not degenerate, explaining the decreasing value of
storage after some level of throughput is more involved, but the main intuition provided here
remains relevant.

References

[53] D. P. Heyman and M. J. Sobel. Stochastic Models in Operations Research: Volume II, Stochas-
tic Optimization. Dover Publications, Inc., Mineola, NY, USA, 2004.

[54] E. L. Porteus. Foundations of Stochastic Inventory Theory. Stanford Business Books, Stanford,
CA, USA, 2002.

OA-4