Commodity Resource Valuation And Extraction: A Pathwise Programming Approach

0.23 MB·

by  Juri Hinz

Commodity Resource Valuation And Extraction:
6 A Pathwise Programming Approach
1
0
2
Juri Hinz Tanya Tarnopolskaya Jeremy Yee*

l
u
J July 14, 2016

3
1
Abstract

]
C Complexityanduncertaintyassociatedwithcommodityresourcevalu-
ationandextractionrequiresstochasticcontrolmethodssuitableforhigh
O
dimensional states. Recent progress in duality and trajectory-wise tech-
.
h niqueshasintroduced avarietyof fresh ideasto thisfieldwith surprising
t results. This paper presents a first application of this promising devel-
a
opment to commodity extraction problems. We introduce efficient algo-
m
rithms for obtaining approximate solutions along with a diagnostic tech-
[ nique,which providesaquantitativemeasurefor solution performancein

2 terms of the distance between the approximate and the optimal control
v policy.
0
Keywords Duality, Markov decision process, natural resource extraction, optimal
1
switching, real option, value function approximation
2
0
0 1 Introduction
.
2
0 Extraction projects for commodities and their valuation can often be formulated as
6 optimalstochasticcontrolproblemsoftheswitchingtype. Suchproblemsusuallyhave
1 no known closed-form solutions and hence numerical methods often remain the only
:
v practical approach to address important operational and investment questions. This
i paperpresentsa family of novelmethodsfor calculating approximate numerical solu-
X
tions based on primal and dual techniques. While the primal methods are based on
r the so-called convex stochastic switching (CSS) and deliver an approximate solution,
a
thedualmethodsutilizerecentadvancesofpathwise dynamicprogrammingtoprovide
solution improvement and its quality assessment. This paper will demonstrate our
algorithms by considering the optimal decision policy and value of natural resource
assets associated with commodity extraction. Numerical results were obtained using
the authors’ software package implemented in the statistical language R. Similar op-
timal switching problems were considered by [2–4,6,16,17]. This paper will depart
from thestandard literature and consider thepresenceof mean reversion in thecom-
modityprice,operationaleffienciesintheformofwastefulextraction, physicaldelivery
obligations to clients and commodity prices subject to stochastic volatility. These
considerations add realism to our study and sheds practical insights.
In the next section, we introduce the problem setting. Section 3 outlines our
numericalapproachwhileSection5containsthenumericalresults. Section6concludes.
1
2 Commodity Extraction As A Stochastic Switch-
ing Problem
Thispaperwillexaminethemanagementofacommodityresourcewithafiniteamount
of commodity as studied by [4]. In this setting, the decision maker aims to maximize
the total expected profit using controls whose costs are given in Table 2.1. Doing
so, the controller dynamically switches the operational mode of the resource between
Opened, Closed, Abandoned . When opened, commodity is extracted, sold at the
{ }
spot market and a revenue is realized based on the commodity price. While closed,
commodity is preserved but a maintenance cost is incurred. The abandoned mode
incursnocost and providesnorevenue. Switchingtotheabandoned modecausesthe
mode to remain permanently there. There is a switching cost between modes.
Table 2.1: Hypothetical Commodity Resource.
Output rate: 10 million pounds/year Inventory level: 150 million pounds
Costs of production: $0.50/pound Opening/closing costs: $200,000
Maintenance costs: $500,000/year Inflation rate: 8%/year
Convenience yield: 1%/year Price variance: 8%/year
Real estate tax: 2%/year Income tax: 50%
Royalty tax: 0% Interest rate: l0%/year
Theoriginalproblemsettingin[4]wasconsideredinacontinuoustimesettingand
solvedusingnumericalpartialdifferentialequations. However,thispaperwillconsider
theproblemasaMarkov decision process (see[11]))indiscretetimeandafinitetime
horizon. The rates, costs and taxes above will beadjusted accordingly.
Given a finite time horizon 0,1,...,T N, consider a controlled Markovian
process (X )T := (P ,Z )T w{hich consist}s⊂of two parts. The discrete component
t t=0 t t t=0
(P )T describes the evolution of a finite-state controlled Markov chain which takes
t t=0
valuesinafinitesetP. Furtherassumethatatanytimet=0,...,T 1thecontroller
−
takes an action a A from a finite set A of all admissible actions in order to cause
the one-step trans∈ition from the mode p P to the mode p′ P with a probability
∈ ∈
Lαep,tp′u(sa)n,owwhteurren(tαop,tph′e(ae)v)po,lpu′t∈iPonaorfetphree-ostpheecrificeodmpstooncehnatst(iZc )mTatroicfetshfeorstaatlleapr∈oceAss.
t t=0
(X )T . Here, we assume that it follows an uncontrolled evolution in the Euclidean
t t=0
space Rd driven as
Z =W Z , t=0,...,T 1
t+1 t+1 t
−
byindependent disturbance matrices (W )T . That is, thetransition kernels a gov-
erningtheevolutionofourcontrolledMatrkto=v1process(X )T :=(P ,Z )T froKmt time
t t=0 t t t=0
t to t+1 is given for each a A by
∈
Ktv(p,z)= αp,p′(a)E(v(p′,Wt+1z)), p∈P, z ∈Rd, t=0,...,T −1
pX′∈P
which acts on each function v : P Rd R where the above expectations are well-
× →
defined.
In this work, we use the discrete component (P )T to convey information re-
t t=0
garding the remaining level of commodity and the current operational mode. More
2
precisely, set P = 0,1,...,I 1,2 to capture all possible positions with the in-
{ }×{ }
terpretationthatthefirstcomponentp(1) oftheposition (p(1),p(2)) Pdescribesthe
∈
remaining commodity level such that p(1) = 0 represent an exhausted or abandoned
resource. Thesecond componentp(2) indicateswhetherthecurrentoperational mode
is closed (p(2) = 1) or opened (p(2) = 2). To introduce the mode control, we set
A = 0,1,2 where a = 0,1,2 stands for the action to abandon, to close, and to
open t{he ass}et, respectively. The second component (Z )T will be used to capture
t t=0
the market dynamics of the the commodity price. The exact form of the disturbance
matrices dependson theparticular choice of theprice process. If thesystem is in the
state (p,z), the costs of applying action a A at time t=0,...,T 1 are expressed
∈ −
throughr (p,z,a). Havingarrivedattimet=T inthestate(p,z),afinalscrap value
t
r (p,z) is collected. Thereby thereward and scrap functions
T
r :P Rd A R, r :P Rd R
t T
× × → × →
are exogenously given for t=0,...,T 1. In our applications, thescrap valueis
−
r (p,z)=maxr (p,z,a), (p,z) P Rd
T T
a∈A ∈ ×
where the reward functions are specified in the following manner. If the asset is
abandoned,there are neithercosts nor revenue
r ((0,p(2),z),a)=0, a A,z Rd.
t
∈ ∈
Forthe case where p(1) >0, we define
r ((p(1),p(2),z),a)=h (z)1 (a)+m 1 (a)+c 1 (a)p(2) a
t t {2} t {1} t {1,2} | − |
withthefollowinginterpretation: Iftheresourceisopened,arevenueiscollectedwhich
dependsoncontinuousstatecomponentz throughapre-specifiedconvexfunction h .
t
Iftheresourceisclosed,thennon-stochasticmaintenancefeem istobepaid. Finally,
t
adeterministicswitchingfeec isincurredwhenevertheoperational modetransitions
t
from opened to closed or vice versa. At each time t = 0,...,T the decision rule π
t
is given by a mapping π :P Rd A, prescribing at time t an action π (p,z) A
t t
for a given state (p,z) P R×d. A→sequenceπ=(π )T−1 of decision rules is calle∈d a
policy. For each policy∈π=×(π )T−1, theso-called potlict=y0value vπ(p ,z ) is definedas
t t=0 0 0 0
thetotal expected reward
T
vπ(p ,z )=E(p0,z0),π r (P ,Z ,π (P ,Z ))+r (P ,Z ) .
0 0 0 t t t t t t t t t
" #
t=0
X
In this formula E(p0,z0),π stands for the expectation with respect to the probability
distributionof(X )T :=(P ,Z )T inducedbytheMarkovtransitionsfrom(P ,Z )
t t=0 t t t=0 t t
to(P ,Z )inducedbythekernels πt(Pt,Zt) fort=0,...,T 1,given theinitial
t+1 t+1 Kt −
value(P ,Z )=(p ,z ).
0 0 0 0
Now we turn to theoptimization goal. A policy π∗ =(π∗)T−1 is called optimal if
t t=0
it maximizes the total expected reward overall policies π vπ(p,z). To obtain such
7→ 0
policy, one introducesfor t=0,...,T 1 the so-called Bellman operator
−
Ttv(p,z)=ma∈aAxrt(p,z,a)+ αp,p′(a)E[v(p′,Wt+1z)] (2.1)
pX′∈P
 
3
for (p,z) P Rd, acting on all functions v where the stochastic kernel is defined.
∈ ×
Consider the Bellman recursion, also referred toas backward induction:
v (p,z)=r (p,z), v = v for t=T 1,...,0. (2.2)
T T t t t+1
T −
Having assumed that the reward functions are convex and globally Lipschitz and
the disturbances W are integrable, there exists a recursive solution (v∗)T to the
t+1 t t=0
Bellman recursion ( [7]). These functions (v∗)T are called value functions and they
t t=0
determinean optimal policy π∗ =(π∗)T via
t t=0
πt∗(p,z)=argma∈aAxrt(p,z,a)+ αp,p′(a)E[vt∗+1(p′,Wt+1z)], (2.3)
pX′∈P
 
for t=T 1,...,0.
−
3 Primal Approximate Solution
Thefirst step in obtaining anumerical solution to thebackward induction (2.2) is an
appropriate discretization of theBellman operator (2.1) to
n
Ttnv(p,z)=ma∈aAxrt(p,z,a)+ αp,p′(a) νtn+1(k)v(p′,Wt+1(k)z)
pX′∈P Xk=1
 
where the probability weights (νn (k))n corresponds to the distribution sampling
t+1 k=1
(W (k))n ofeachdisturbanceW . Intheresultingmodifiedbackwardinduction
t+1 k=1 t+1
governed by vn = nvn , the modified functions (vn)T need to be described by
t Tt t+1 t t=0
algorithmically tractable objects. Since the reward and scrap functions are convex
in the continuous variable, these modified value functions are also convex and can
be approximated by piecewise linear and convex functions. For this, introduce the
so-called sub-gradient envelope Gmf of a convex function f : Rd R on a grid
Gm Rd with m points i.e. GmS= g1,...,gm →
⊂ { }
Gmf = g∈Gm(▽gf)
S ∨
which is a maximum of the sub-gradients ▽ f of f on all grid points g Gm. Using
g
∈
thesub-gradient envelope operator, definethedouble-modified Bellman operator as
n
Ttm,nv(p,z)=SGmma∈aAxrt(p,z,a)+ αp,p′(a) νtn+1(k)v(p′,Wt+1(k)z).
pX′∈P Xk=1
 
The corresponding backward induction
vTm,n(p,z) = SGmma∈aAxrT(p,z,a), p∈P,z∈Rd (3.1)
vm,n(p,z) = m,nvm,n(p,z), t=T 1,...0. (3.2)
t Tt t+1 −
yieldstheso-called double-modifiedvaluefunctions (vn,m)T . Underappropriate as-
t t=0
sumptions(see[7]),thedouble-modifiedvaluefunctionsconvergeuniformlytothetrue
valuefunctionsalmostsurelyoncompactsets. Theseassumptionsincludetheconvex-
ityandglobalLipschitzcontinuityoftherewards,theintegrabilityofalldisturbances
and some restrictions on thedistribution sampling and grid density.
4
Sincethedouble-modifiedvaluefunctions(vm,n)T arepiece-wiselinearandcon-
t t=0
vex,theycanbeexpressed inacompact andappealingform usingmatrixrepresenta-
tions. Notethat anypiecewise convexfunction f can bedescribed byamatrix where
each linear functional is represented by a row in the matrix. To denote this relation,
let us agree on the following notation: Given a function f and a matrix F, we write
f F whenever f(z) = max(Fz) holds for all z Rd. Let us emphasize that the
sub∼-gradient envelope operation Gm on a grid Gm∈is reflected in terms of a matrix
S
representativeby a specific row-rearrangement operator
f F Gmf ΥGm[F]
∼ ⇔ S ∼
wheretherow-rearrangementoperatorΥGm associatedwithGm = g1,...,gm Rd
{ }⊂
acts on matrix F with d columns as follows:
(ΥGmF)i,· =Fargmax(Fgi),· for all i=1,...,m. (3.3)
For piecewise convex functions, the result of maximization, summation, and compo-
sition with linear mapping, followed by sub-gradient envelope can be obtained using
their matrix representatives. More precisely, if f F and f F holds, then it
1 1 2 2
∼ ∼
follows that
Gm(f1+f2) ΥGm(F1)+ΥGm(F2)
S ∼
Gm(f1 f2) ΥGm(F1 F2)
S ∨ ∼ ⊔
Gm(f1(W ) ΥGm(F1W)
S · ∼
where W is an arbitrary d d matrix and the operator stands for binding ma-
× ⊔
trices by rows. Therefore, the backward induction (3.1) and (3.2) can be expressed
in terms of the matrix representatives Vm,n(p) of the value functions v(m,n)(p,z) for
t t
p P,z Rd andt=0,...T. Sincethedouble-modifiedbackwardinductioninvolves
∈ ∈
maximization, summations and compositions with linear mappings applied to piece-
wise linear convex functions, it can be rewritten in terms of matrix operations which
results in thefollowing algorithm.
Algorithm 3.1: Convex Stochastic Switching Algorithm
1 for p P do
2 VTm∈,n(p)←ΥGm ⊔a∈ASGmrT(p,.,a)
3 end
4 for t T 1,...,0 do
∈{ − }
5 for p P do
6 V˜tm∈+1,n(p)← nk=1νtn+1(k)ΥGmVtm+,1n(p)Wt+1(k)
7 end
P
8 for p P do
∈
9 Vtm,n(p)←ΥGm ⊔a∈A SGmrt(p,.,a)+ p′∈Pαp,p′(a)V˜tm+,1n(p′)
10 end (cid:16) (cid:17)
P
11 end
The continuation value V˜m,n(p) given by the sixth line in Algorithm 3.1 can be
t+1
further approximated by recalling that the disturbances W are independently and
t+1
5
identically distributed across time. Through the use of a suitable nearest neighbour
algorithmtoconstructpermutationmatrices,theconditionalexpectationonthesixth
line can be approximated to save a significant amount of computational effort. For
furtherdetails, we refer theinterested reader to [8].
A candidate for a nearly optimal policy is given by
πtm,n(p,z)=argma∈aAxrt(p,z,a)+max αp,p′(a)V˜tm+,1n(p′)z.
pX′∈P
  
Given the pointwise convergence of the value functions for all p P, a A, and
∈ ∈
t=0,...,T 1 one can deducethat
−
lim πm,n(p,z)=π∗(p,z) p P,a A, t=0,...,T 1
(m,n)→∞ t ∈ ∈ −
holds for each point (p,z) P Rd. However, this convergence may be of little help
∈ ×
in determining the quality of a candidate policy in practice. To address this issue, a
pathwise dynamic approach will be described in the next subsection. This area has
recently attracted growing research interest on diagnostics and potential a posteriori
justification of an approximatecontrol policy.
4 Dual Solution And Diagnostics
This section utilizes an approach suggested in [12] and [5] to assess the quality of a
given control policy obtained by primal methods. Let us first illustrate their use for
thecaseofaoptimalstoppingproblem. Givenareal-valuedstochasticprocess(Z )T
t t=0
adaptedtoafiltration,considerthefamily ofall 0,1,...,T -valuedstoppingtimes.
V { }
Obviously,theoptimal stopping value is dominated
V∗ =supE(Z ) E( sup Z )
0 τ ≤ t
τ∈V 0≤t≤T
bythe expectation of thepathwise maximum. This dominance still holds
V∗ =supE(Z M ) E( sup (Z M )).
0 τ − τ ≤ t− t
τ∈V 0≤t≤T
whentheoriginalprocess(Z )T isreplacedby(Z M )T withamartingale(M )T
t t=0 t− t t=0 t t=0
starting at theorigin M =0. Furthermore, it turnsout that thisestimate is tight
0
V∗ = inf E( sup (Z M ))=E( sup (Z M∗)) (4.1)
0 (Mt)Tt=0 0≤t≤T t− t 0≤t≤T t− t
where the infimum is taken over the family of all martingales starting at the origin
andinfact isattained at someoptimal martingale(M∗)T from thisfamily. Usually,
t t=0
an optimal martingale is not available, as its knowledge is equivalent to the solution
of the stopping problem. The equation (4.1) yields a practical way to estimate the
optimal value V∗ via a Monte Carlo procedure by determining an average of the
0
pathwisemaximumonanumberofsimulatedtrajectories of(Z M˜ )T . Obviously,
to obtain a tight upper estimate, the martingale (M˜ )T mustt−resemt bt=le0the optimal
t t=0
one. Similarly, a lower estimate for V∗ can be obtained from Monte-Carlo average of
0
independentrealizations ofZ M˜ ,basedonanarbitrarystoppingtimeτ˜. Againto
τ˜ τ˜
−
obtain a tight lower bound, the stopping time τ˜ must be chosen close to the optimal
6
stopping time. This procedure of bound estimate exhibits an interesting self-tuning
property: The closer (M˜ )T and τ˜ are to their optimal counterparts, the lower the
t t=0
varianceoftheMonte-Carlotrialsandthenarrowerthegapbetweenbothbounds. In
a hypothetic case where τ˜ is the optimal stopping time and (M˜ )T is the optimal
t t=0
martingale, the variance of both Monte-Carlo tails reduces to zero. Moreover, in this
case the upperand thelower boundscoincide with the optimal stopping value.
This approach has been successfully applied [1] and generalized to multiple stop-
ping problems in finance [9]. For our switching problem, the pathwise dynamic pro-
gramming approach works similarly (see [8]). Given an approximate control problem
solution,tworandomboundvariablesareconstructed,whoseexpectationsgivealower
andanupperestimatefortheunknownvaluefunction. Theseboundvariablesarede-
termined througha stochastic recursiveproceduresimilar tothebackward induction.
Their calculation is effected in terms of a Monte-Carlo procedure whose in-sample
empirical variance can be used to quote a confidence interval for the unknown value
function. Inthissetting,theself-tuningpropertystatesthatboth,thevarianceofthe
Monte-Carlo trials and the gap between expectations of bound variables decrease, if
theapproximate control policy approaches theoptimal policy.
Suppose that our numerical scheme returns approximate value functions (v )T
t t=0
along with the approximate expected value functions (vE)T , thus we introduce an
t t=0
approximately-optimal policy π=(π )T−1 by
t t=0
πt(p,z)=argmax(rt(p,z,a)+ αp,p′(a)vtE+1(p′,z)))
pX′∈P
To answer the question how far is the strategy π from an optimal strategy π∗, we
estimatetheperformancegap[vπ(p ,z ),vπ∗(p ,z )]atagivenpointapoint(p ,z )
P Rdbyanexplicitconstructio0nof0ran0dom0 var0iab0lesvπ,ϕ(p ,z ),vϕ(p ,z )sat0isfy0in∈g
× 0 0 0 0 0 0
E(vπ,ϕ(p ,z ))=vπ(p ,z ) vπ∗(p ,z ) E(v¯ϕ(p ,z )).
0 0 0 0 0 0 ≤ 0 0 0 ≤ 0 0 0
UsingMonte-Carlosimulations, oneestimatesbothmeansalongwithempirical confi-
dence intervals to estimate the performance gap, which is narrow if the approximate
solutionisclosetooptimal. Thisusefulpropertyisduetothetheself-tuningproperty,
which ensures that the more optimal solutions (v )T (vE)T return narrower gaps
t t=0 t t=0
and lower variance in the Monte-Carlo scheme. This scheme can be implemented as
follows
1. Givenapproximatesolution(v )T (vE)T withthecorrespondingpolicy(π )T−1,
t t=0 t t=0 t t=0
implement control variables (ϕ )T as
t t=1
I
1
ϕt(p,z,a)=pX′∈Pαp,p′(a) I Xi=1vt(p′,Wt(i)z)−vt(p′,Wtz)!,
for all p P, a A, z Rd, where (W(1),...,W(I),W )T are independent
∈ ∈ ∈ t t t t=1
and (W(1),...,W(I),W ) are identically distributed for each t=1,...,T.
t t t
2. Chose a number K N of Monte-Carlo trials and obtain for k = 1,...,K
independent realizatio∈ns (W (ω ))T of disturbances.
t k t=1
3. Startingatzk :=z Rd,definefork=1,...,Ktrajectories(zk)T recursively
0 0∈ t t=0
zk =W (ω )zk, t=0,...,T 1
t+1 t+1 k t −
7
and determinerealizations
ϕ (p,zk ,a)(ω ), t=1,...,T, k=1,...,K.
t t−1 k
4. For each k=1,...,K initialize therecursion at t=T as
vπ,ϕ(p,zk)(ω )=vπ(p,zk)(ω )=r (p,zk) for all p P
T T k T T k T T ∈
and continuefor t=T 1,...,0 and for all p P by
− ∈
vπ,ϕ(p,zk)(ω )=r (p,zk,π (p,zk))+ϕ (p,zk,π (p,zk))(ω )
t t k t t t t t+1 t t t k
+ αp,p′(πt(p,ztk))vπt+,ϕ1(p′,ztk+1)(ωk)
pX′∈P
vϕ(p,zk)(ω )=max r (p,zk,a)+ϕ (p,zk,a)(ω )
t t k a∈A t t t+1 t k
(cid:2) + αp,p′(a)vϕt+1(p′,ztk+1)(ωk)
pX′∈P (cid:3)
5. Calculatesamplemeans 1 K vπ,ϕ(p ,z )(ω ), 1 K vϕ(p ,z )(ω )toes-
K k=1 0 0 0 k K k=1 0 0 0 k
timatethemeansE(vπ,ϕ(p ,z )),E(vϕ(p ,z ))alongwiththeirconfidenceinter-
0 0P0 0 0 0 P
vals. Thesemeanswillbereferredtoastheprimal anddual values,respectively.
5 Numerical Results
This section demonstrates our algorithms on the commodity extraction problem de-
scribed in Section 2. The values given by Table 2.1 will be adjusted in order fit
our discrete time setting with time step ∆ = 0.25 years between the decision times
0,...,T N representing a time horizon [0,T¯] = [0,30] of thirty years, with four
{ } ⊂
decisions made every year. Let us discretize the commodity levels remaining in the
resourcewith thestepsize corresponding tocommodity amountextracted within one
time step ∆. That is, the first entry of the discrete component covers the range
p(1) 0,1,..., p¯ where p¯=15 stands for the minimum numberof years it takes to
∈ ∆
deplete the resource. Further, define the number of decision epochs T = T¯ +1, the
(cid:8) (cid:9) ∆
maintenance costs m = m ∆e(ρ−r−ζ)t∆ and the switching costs c = c e(ρ−r−ζ)t∆
t 0 t 0
with coefficients r = 0.1, ρ = 0.08, ζ = 0.02, m = 0.5 and c = 0.2 standing for
0 0
the interest rate, rate of inflation, real estate tax, initial maintenance cost and initial
switching cost, respectively.
Geometric Brownian Motion: Following [4], the continuous state component
(Z )T is one-dimensional and follows linear state dynamics:
t t=0
σ2
Z =exp µ ∆+σ√∆N Z (5.1)
t+1 − 2 t+1 t
(cid:18)(cid:18) (cid:19) (cid:19)
with σ2 = 0.08 and independent identically standard normally distributed (N )T .
t t=1
Thecontrolledtransitionprobabilitiesofthediscretecomponentgiven(p(1),p(2)) P
∈
is uniquelydetermined by
α (Open)=1,
(p(1),p(2)),(max{p(1)−1,0},Opened)
α (Close)=1,
(p(1),p(2)),(p(1),Closed)
α (Abandon)=1.
(p(1),p(2)),(0,p(2))
8
The mode is controlled deterministically under the above specfication. Now let us
definethecontrolcostsbyafunction,representingtherevenueoftheopenedresource
h (z)=5∆ze−(r+ζ)t∆ 2.5∆e(ρ−r−ζ)t∆
t
−
where thevariable z represents thecommodity price.
An equally spaced grid ranging from 0 to 20 of size 4,001 was used and a dis-
tribution sampling of size 20,000 constructed using equidistant sampling of standard
normalquantileswasappliedtoobtainthefollowingresults. Table5.1liststheprimal
anddualvalues,generatedusingK =1,000pathsandI =1,000sub-simulations. To
compare our results with [4] listed in the first column of the Table 5.1, we assumed
a convenience yield d = 0.01 and set µ = r d = 0.09 in (5.1). We were able to
−
obtain similar results despite [4] using a continuous time formulation. Further, our
results were able to obtain excellent precision as evidenced by the tight bounds and
low standard errors given in theparanthesis.
Table 5.1: Resource Valuation Under Geometric Brownian Motion
B&S Opened Resource Closed Resource
Z0 Open Closed Primal Dual Primal Dual
0.3 1.25 1.45 1.2127(.0026) 1.2156(.0026) 1.4127(.0026) 1.4156(.0026)
0.4 4.15 4.35 4.1059(.0034) 4.1086(.0034) 4.3059(.0034) 4.3086(.0034)
0.5 7.95 8.11 7.9026(.0039) 7.9053(.0039) 8.0752(.0041) 8.0777(.0041)
0.6 12.52 12.49 12.5129(.0042) 12.5153(.0042) 12.4787(.0045) 12.4813(.0045)
0.7 17.56 17.38 17.5869(.0047) 17.5889(.0047) 17.3869(.0047) 17.3904(.0048)
0.8 22.88 22.68 22.9475(.0052) 22.9489(.0052) 22.7475(.0052) 22.7489(.0052)
0.9 28.38 28.18 28.4940(.0057) 28.4957(.0057) 28.2940(.0057) 28.2957(.0057)
1.0 34.01 33.81 34.1667(.0062) 34.1681(.0062) 33.9667(.0062) 33.9681(.0062)
This paperwill illustrate thecontrolpolicies intheform of thecurvesdepictedin
Figure 5.1. We have observed that for all the situations examined in this paper, the
optimalpoliciessharethisstructue. Namely,itisoptimaltooperateanopenresource
when the commodity price is sufficiently high. Below some price level, the resource
should be shut down, either to the closed or abandoned mode. Typically, we have a
curvebifurcationandtheareabetweenbothcurvebranchesrepresentsthepricerange
whereswitchingtotheclosed modeisoptimal. IntherightplotofFigure5.1,theop-
timalpoliciesforaclosedresourceisshownfordifferentvaluesofµ. Recallthatthere
is a strong element of scarcity in our problem. That is, selling a unit of limited com-
modity now means that thesame unit of commodity cannot be sold at a future date.
Thus,weobservefrom Figure5.1,that thedecision makerismore willing topreserve
commodityforfutureuseifthedrifttermµishigh. Inoursubsequentstudies,wewill
assumedriftµ=r d=0.09(unlessotherwisestated)inordertostayinlinewith[4].
−
Ornstein-Uhlenbeck Process: Commodity prices often mimic the business cycle
in practice(see [10,14]). Prices with significant mean reversion often lead todifferent
optimal behavior compared to those following geometric Brownian motions ( [13–
15]). To examine the impact of mean reversion, the logarithm of the commodity
price (logZ˜ )T is assumed to follow an AR(1) process, an auto-regression of order
t t=0
one. With this choice, define the linear state evolution (Z )T in the required form
t t=0
9
Optimal Policy t = 0 Closed Mine
0.6
Open
0.6
0.4 Close
Commodity Price Ore Price mu0000....00006789
Abandon 0.4
0.2
0.0 0.2
0 20 40 60 0 20 40 60
Commodity Reserves Ore Level
Figure 5.1: The optimal policy at t = 0. The horizontal axis represents
remaining reserves while the vertical represents the commodity price. The
right plot displays policy for a commodity resource that begins closed for
µ=0.09,0.08,0.07,0.06.
Z =W Z as
t+1 t+1 t
1 1 0 1
Zt+1 :=(cid:20)logZ˜t+1(cid:21)=" µ− σ22 ∆+σ√∆Nt+1 φ#(cid:20)logZ˜t(cid:21),
(cid:16) (cid:17)
for t = 0,...,T 1, with initial value logZ˜ R and µ = 0.09, σ2 = 0.08. Again,
0
(N )T are indep−endentstandard normally dis∈tributed random variables and thepa-
t t=1
rameter φ [0,1] determines the speed of mean reversion, where φ = 1 gives a
∈
geometric Brownian motion. Using the same parameters as before, merely adjust the
cash flow function accordingly for t=0,...,T 1 to
−
h (z)=5∆exp z(2) exp( (r+ζ)t∆) 2.5∆exp((ρ r ζ)t∆)
t
− − − −
(cid:16) (cid:17)
since now the variable z(2) captures the logarithmic commodity price. Figure 5.2
shows the nearly-optimal policies at initial time t = 0 for an opened asset, which
are computed using 2000 equally spaced grid points ranging from -5 to 5 for log(Z˜ )
t
and 10,000 disturbances constructed from the equidistant quantiles of the normal
distribution. Diagnostics is performed using K =500 sample pathsand I =500 sub-
simulations. Table 5.2contains resultsforφ=1,0.8,0.6 andshows verytight bounds
and low standard errors.
Figure 5.2 provides an interesting insight. We observe that under mean reversion
φ <1,eachadditionalunitofcommodityreserveexhibitdiminishingmarginalvalue
| |
which directly contrasts the behaviour under geometric Brownian motion. This phe-
nomenonwasalsoobservedfortheclosedcommodityresourceandwillhavesignificant
implications later on.
Wasteful Extraction: For w [0,1] and commodity level p(1), assume that
∈
10