Modeling and computation of water management by real options

It becomes increasingly important to manage water and improve the efficiency of irrigation under higher temperatures and irregular precipitation patterns. The choice of investment in water saving technologies and its timing play key roles in improving efficiency of water use. In this paper, we use a real option approach to establish a model to handle future uncertainties about the water price. In addition, to match the practical situation, the expiration of the real option is considered to be finite in our model, such that it is difficult to solve the model. Therefore, we reformulate the problem into a linear parabolic variational inequality (Ⅵ) and develop a power penalty method to solve it numerically. Thus, a nonlinear partial differential equation (PDE) is obtained, which is shown to be uniquely solvable and the solution of the nonlinear PDE converges to that of the Ⅵ at the rate of \begin{document}$O(λ^{-\frac{k}{2}})$\end{document} with \begin{document}$λ$\end{document} being the penalty number. Furthermore, a so-called fitted finite volume method is proposed to solve the nonlinear PDE. Finally, several numerical experiments are performed. It is shown that the subjective discount rate will affect the investment boundary mostly, and the flexibility to suspend operation will enlarge the investment region.


1.
Introduction. According to IPCC [39], the scarcity of water is defined as per capita supplies less than 1700 m 3 /year. Obviously, climate change poses an urgent threat to water security. In fact, changes in precipitation and other climatic variables such as higher temperatures may lead to significant changes in water supply in many regions. As a result, about 1.2 billion people worldwide lack access to water, and a total of 2.8 billion face water scarcity for at least one month of the year. In many countries, the agriculture sector is the predominant consumer of water.
Agriculture uses 70% of the world's accessible freshwater, but nearly 60% of this is wasted due to leaky irrigation systems, inefficient application methods as well as the thirsty cultivated land of crops [40]. This wasteful use of water is drying out rivers, lakes and underground aquifers. Irrigated agriculture, which accounts for only 20 percent of cultivated land but over 40 percent of the global harvest, has significant implications for the future of water availability and food security worldwide. The vast quantities of water consumed in irrigation mean that any gains in efficiency can substantially decrease water usage. The average efficiency of irrigation is only 37% worldwide, only about 5% of the irrigated lands use the most efficient irrigation methods available. If farmers could raise total use efficiency of 90%, it would reduce the water demand by upward of 15% [23]. Hence, innovations in agricultural water management have been important and will be even more indispensable in the process of climate change.
Evidently, effective management of agricultural water in water-scarce areas requires efficient approaches. In addition, the uncertainty about water prices should affect the irrigators' decisions regarding the timing of investment in water saving technologies. We select the real options valuation (ROV) as the capital budgeting technique to deal with the decision making problems. The real options valuation analysis has been widely used to manage the project investment in many fields [12], [15], [2], [13], [28], [6], [19], [10], [29]. However, to our best knowledge, most of the real options analysis references about investment in irrigation technology, such as [22], do not take the expiration into consideration as it reduces the complexity of the problem and simplifies the exercise boundary into a constant. This makes the real option turn into a perpetual one, which is not practical in the real world. In this paper, we make a real options analysis on the water management project and assume that the real option depends on both time and water price. Then, the value of the real option is determined by a linear complementarity problem, and the option to suspend application of irrigation technology is considered if the value of water falls down the annual maintenance and operating. Furthermore, it is examined that the impact of subsidy policies to encourage the uptake of technologies.
From the viewpoint of mathematics, our real option model is characterized by a linear complementarity problem (LCP) which can be transformed to a parabolic variational inequality (VI), and an efficient numerical method should be required, since its solution is not available analytically. In recent years, some researchers have paid their attention to solving the VIs arising from different fields [1], [36], [35], [14], [11], [5]. In the case of parabolic VI, Ito and Kunisch [18] propose an approach based on a Lagrange multiplier treatment, and establish the existence and uniqueness of the strong as well as weak solutions. Wang and Yang [30] present a regularization method to a degenerate variational inequality of parabolic type arising from American option pricing. They first use a smoothing technique with small parameter > 0 for the non-smoothing initial value function, and derive the error estimates for the regularized continuous and discrete problems, respectively. In the thesis [34], the author provides a detailed analysis on the power penalty method for the American option pricing problem using the theory of viscosity solutions and gives the convergence rate of the method.
Additionally, in some other literature such as [3], [25], [24], [4], the semi-smooth Newton method combined with the Moreau-Yosida regularization has been used to deal with the optimal control problem and the relevant VI. The semi-smooth Newton method shows the advantages in solving non-smooth problems and has superlinear convergence under some assumptions. However, similar to the general Newton method, the semi-smooth Newton method is a local method and the computation speed heavily depends upon the choice of initial guess.
The power penalty method is easily implemented and widely used in any dimensional problems on both structured and unstructured meshes [7], [8], [33]. In this work, we propose a power penalty approach to approximate the parabolic VI by a nonlinear parabolic partial differential equation (PDE) with an l k penalty term. The power penalty method mainly possesses two advantages. On the one hand, the resulting penalized PDE is of a simple form that is easy to discretize in any dimensions by using the finite element method, finite volume method or finite difference method. On the other hand, it can be shown that the solution to the nonlinear PDE converges to that of the original complementarity problem at the rate of order O(λ − k 2 ) with λ being the penalty number. This exponential convergence rate implies that a desirable accuracy in the approximate solution can be achieved by a judicious choice of the penalty parameters. In other words, we can achieve a much higher level of accuracy with a small penalty parameter by using the l k (0 < k < 1) penalty method.
Moreover, a so-called fitted finite volume method is presented for the numerical solution of the nonlinear partial differential equation. This method is based on a finite volume formulation coupled with a fitted approximation technique [38]. The finite volume method possesses a special feature of the local conservativity of the numerical fluxes, and is becoming more and more popular. See, for instance, [27] for degenerate parabolic problems, [20] for hyperbolic problems, [21] for elliptic problems, and [37,9] for HJB equations.
The paper is organized as follows. In Section 2, a partial differential equation of the real option is established to value the investment opportunity, or the option to invest. Also, the final and boundary conditions are prescribed. A power penalty approximation is proposed for the original problem, and its convergence analysis is presented in Section 3. In Section 4, a so-called fitted volume method is proposed for the discretization of the penalized PDE. In addition, some numerical experiments are performed to illustrate the efficiency and the usefulness of the numerical method in Section 5, and their meanings in economics and management are explained. Finally, some concluding remarks are given in Section 6.

2.
The PDE real option model and power penalty approach.
2.1. The stochastic process of water price. Following [22], we use the geometric Brownian motion to represent the stochastic process for water price. Let P t denote the water price at time t. Then the dynamical process of the water price is given by where α is the expected price appreciation or drift rate, σ is the standard deviation of the water price, dt is the change in time, and W is a Brownian motion. In this case, water price cannot be negative and its volatility is constant over time.
2.2. The value of investment.
2.2.1. No option to suspend operation when cost exceeds price. Following [6], the expected present value of a project to invest in a water saving technology, denoted by V (P ), is expressed as: where P is the water price, a is the amount of water saved by using the new technology, and C is the operating cost. Here δ denotes the difference between the discount rate and the rate of growth in the asset value expected over time, that is, δ = ρ − α, where ρ is the risk adjusted discount rate. In equation (2), the stochastic variable P is discounted by the risk adjusted discount rate δ and the deterministic variable C is discounted by the risk-free interest rate r.

2.2.2.
Suspending operation when cost exceeds price. When irrigators use the chemical monolayers to save water, they have the option to suspend using the monolayers and purchase some water in the market as the water price falls below the operating cost. In this case, the value of the project is governed by the following ordinary differential equation (ODE) [12]: where Π(P ) = max{P a − C, 0}. We can solve it analytically to obtain where β 1 , β 2 is the two roots of the quadratic equation 1 2 σ 2 β(β − 1) + (r − δ)β − r = 0, and K 1 and B 2 are given by , respectively.
2.3. The value of the real option. There are two parallel methods, namely the dynamic programming method and the contingent claim method, to value a real option. According to [12], the two methods should give the same results of valuation and optimal investment boundary. However, the contingent claim method requires enough traded assets to make the market complete. So, we exploit the dynamic programming method to value the real option and set the discount rate to equal the risk-free interest rate, which implies that the market price of risk is zero.
During each period of time t, the investor can either continue with waiting or exercise the investment to get the payoff F (P ) = V (P ) − I, where I denotes the investment cost. We assume that a constant discount rate ρ equals r at each time point t. Then, the real option value satisfies the following Bellman equation [12,15]: From (5) we know that if then F (P, t) = F (P ), and Multiplying by 1 + ρdt and rearranging, we obtain which, together with Ito's Lemma, leads to Dividing by dt, we get LF > 0, where the partial differential operator L is defined by otherwise, which means that LF = 0 and F − F > 0. In conclusion, we get the following partial differential complementarity problem: for (P, t) ∈ (0, P max ) × [0, T ) with the boundary conditions and the terminal condition We aim at finding out the optimal exercise boundary P (t), which satisfies the value matching condition F (P , t) = F (P ) and the smooth-pasting condition F (P ) = (F ) (P ) (see, for example, [26]).

The impact of subsidy on investment decision.
There is uncertainty about the value of the water savings because of the fluctuation on water price. Consequently, the return achieved from the investment in a new irrigation technology, for which the main benefit is water savings, will be also uncertain. So, the private incentive to investment in the on-farm technology may not be sufficient and there should be potential to hasten uptake of the technology through the use of cost sharing arrangements to reduce the upfront cost or improve effectiveness of the technology. It can be expected to speed up the technology adoption by reducing the financial barriers, such as shifting or sharing financial risk. Subsidy is the most common form of financial incentives and is needed to defray upfront investment costs and initial productivity losses. It is very flexible, and usually does not require complex institutional arrangements. For example, in line with the Council of Australian Governments (COAG) Intergovernment Agreement on Murray Darling Basin Reform (July 2008), which seeks to improve the use and management of the basin water resources in a way that optimizes economic, social and environmental outcomes, the irrigators would be obligated to transfer to the Commonwealth at least 50 percent of the water savings in the form of permanent entitlement as the return for the subsidy under an agreed timetable (see [41]).
Similar to [22], the subsidy is expressed as a percentage of the technology's capital cost I. We examine a range of scenarios for a subsidy of 5%I to a subsidy of 95%I to find how the investment boundary changes. In this experiment, the constant a in (2) and (4) should also become 0.5a for each subsidy level due to the subsidy under a program of the design examined. We will show the results in Section 6.
3. The power penalty approach. In this section, we approximate the linear complementarity problem by a nonlinear parabolic PDE with an l k penalty term. We then show that the solution to the nonlinear PDE converges to that of the original complementarity problem at the rate of order O(λ −k/2 ) with λ being the penalty number (See Appendix C for the details). To solve the penalized nonlinear equation, the fitted finite volume method is proposed in the next section.
For the convenience of theoretical analysis, we rewrite (6) as the following conservative form: For the discussion convenience, we transform (7) into an equivalent form satisfying homogeneous Dirichlet boundary conditions. Note that this transformation is needed only for the theoretical discussion, but not necessary in computations.
Let F 0 (P ) be a twice differentiable function satisfying the boundary conditions in (8). We introduce a new function where β = σ 2 . Transforming F in (7) into the new function u, we have where It is easy to see that under the transformation (10), the boundary and terminal conditions in (8)-(9) become, respectively, u(0, t) = 0 = u(P max , t) and u(P, T ) = u (P, T ).
3.1. The formulation of complementarity problem into a variational inequality problem. In this subsection, we will reformulate (11) as a variational inequality problem in an appropriate functional setting. Before proceeding, let us first introduce some standard notations to be used in the paper later.
Let Ω = (0, P max ) and let Γ denote the boundary of Ω. For 1 ≤ p ≤ ∞, let L p (Ω) denote the space of all p-integrable functions on Ω with the norm · L p (Ω) , and let H m,p (Ω) denote the usual Sobolev space with the norm · m,p,Ω . When p = 2, we simply use H m (Ω) and · m,Ω to denote H m,2 (Ω) and · m,2,Ω , respectively. We define the weighted Sobolev space H 1 ω (Ω) as follows: We put (12). It is easy to verify that K is a convex and closed subset of H 1 0,ω (Ω). Moreover, in what follows, we will simply write v(t) when we regard v(·, t) as an element of H 1 0,ω (Ω), and will also suppress the independent time variable t (or τ ), when it causes no confusion in doing so. Finally, for any Hilbert . Now, we are in the position to define the following variational inequality problem.
where B(u,v;t) is a bilinear form defined by For this variational inequality problem, we can have the following theorem.
Theorem 3.1. Problem 1 is the variational form of the complementarity problem (11).
See Appendix A for the details of the proof. In order to establish the unique solvability of Problem 1, we study the properties of the operator B(u, v; t).
Lemma 3.2. There exist positive constants C and M, independent of v and w, such that for any v, w ∈ H 1 0,ω (Ω), Using Lemma 3.2 and the theory of abstract variational inequalities, the unique solvability of Problem 1 is established in the following theorem.

3.2.
The power penalty approach. To derive the power penalty approach, we first consider the following nonlinear variational inequality problem: where and [z] + = max{0, z} for any z.
The unique solvability of this problem is guaranteed by the coerciveness and continuity of the bilinear operator B and the lower semi-continuity of j. It is easy to show that by virtue of Lemma 3.2 and (16), all the required conditions are satisfied by the bilinear form B. Hence, (15) is uniquely solvable.
From (16) we can see that j(v) is differentiable. Thus, (15) is equivalent to the following problem. where We remark that (17)-(18) is a penalized variational equation corresponding to (13). The strong form of (17)- (18), which defines the penalized equation approximating (11), is given by with the given boundary and final conditions u λ (P, t)| Γ = 0 and u λ (P, T ) = u (P, T ).
If k = 1 2 , this penalty approach corresponds to the quadratic penalty approach. While k = 1, the typical l 1 penalty approach is obtained. When k > 1, it is the socalled lower order penalty approach. Instead of solving the original problem, we will present a so-called fitted finite volume method for (19), which is an approximation of the original problem, in the next section. In Appendix C, we will investigate the convergence rate of u λ to u as λ → ∞.

4.
The fitted finite volume method. In this section, we will present the socalled fitted finite volume method for (19). The idea of this method is based on a finite volume formulation coupled with a fitted approximation technique. This fitting technique is to approximate the flux of a given function locally by a constant, yielding a locally nonlinear approximation to the function [16,17,31,32]. In what follows, we will develop the method for our nonlinear partial differential equation with penalty term. To avoid an overloading of symbols, we omit the subscript λ of F λ used in the above, but bear in mind that this F is the solution of (19).

4.1.
Discretization. It is easy to show that under the inverse of the transformation (10), (19) can be rewritten as where and P i+ 1 2 = Pi+Pi+1 2 for each i = 1, 2, . . . , N − 1. These mid-points form a second partition of (0, P max ) if we define P − 1 2 = P 0 and P N + 1 2 = P N . Integrating both sides of (21) over ( Applying the mid-point quadrature rule to the first, third and the last terms, we obtain from the above that We ignore the complex deduction of the approximation to the flux at the two points P i+ 1 2 and P i− 1 2 , since further details of the discussion can be found in [27] and [32]. The flux has the following analytical expression We define a global piecewise constant approximation to f (F ) by f h (F ) satisfying Substituting (24) into (22), we can get for i = 1, 2, . . ., N − 1, where This is a system of first order linear ordinary differential equations (ODEs). To discretize this system, we let t i (i = 0, 1, . . . , M ) be a set of partition points in [0, T ] satisfying T = t 0 > t 1 > · · · > t M = 0 with subinterval length ∆t m = t m+1 − t m < 0. Then, we apply the two-level implicit time-stepping method with a splitting parameter θ ∈ [0, 1] to (25) to gain the following matrix: (26) for m = 0, 1, . . . , M − 1, where E m and E m+1 denote the values of E(t) at t m and t m+1 , respectively, and . Apparently, the scheme is the Crank-Nicolson scheme when θ = 0.5 and the implicit Euler scheme when θ = 1, both of which are absolutely stable, and they are of second and first order accuracy, respectively.

4.2.
The solution of the nonlinear system. Since the system (26) is nonlinear, we apply the classical Newton method to solve it. Note that when k > 1, we see that d (F m i ) → ∞ as F m i − F i → 0 + . To overcome this difficulty, we smooth the penalty term by the following way: Now, applying the Newton's method to (26), we get for l = 1, 2, . . . , with ω 0 being a given initial guess, where J D (ω) denotes the Jacobian of the column vector D(ω), and γ ∈ (0, 1] denotes a damping parameter. We then choose From [32] we can know that the coefficient matrix of system (28) is an M-matrix, such that (28) has a unique solution.

Numerical experiments.
In this section, we demonstrate the efficiency and the usefulness of our numerical method discussed above by solving the following test problems with time step splitting parameter θ = 0.5, power penalty method parameters λ = 10 and k = 4. Especially, the meanings of these examples in economics and management are also discussed. Test 1: European real option with parameters [22]: P max = $1000, a = 86ML per year, C = $48091, I = $50121, α = 0.02, T = 1, σ = 0.2, r = 0.07, ρ = 0.07, δ = 0.05. The value of the project is chosen as (2). Although a real option is similar to the American option, we test a European real option first here to demonstrate the convergence rate of the discretization method. The valuation PDE of a European option is (21) without the penalty term. To solve the European real option with the parameters in Test 1, we divide (0, P max ) and (0, T ) uniformly into 50 and 50 sub-intervals, respectively. The numerical values of the European option are plotted in Figure 1.
To examine the convergence rate, we present the errors in the discrete L ∞ -norm at the final time step t = 0 in Table 1. The numerical solution on the mesh with N = 2 14 and M = 2 13 is regarded as the "exact" solution F . The ratio is defined as: Moreover, linear regression is used to show that these data obey the following error estimate:   which demonstrates numerically that Test 2: American real option with parameters [22]: P max = $1000, a = 6.1ML per year, C = $315, I = $36750, α = 0.02, T = 1, σ = 0.2, r = 0.07, ρ = 0.07, δ = 0.05. The value of the project is chosen as (2).
To solve the American real option with the parameters in Test 2, we divide (0, P max ) and (0, T ) uniformly into 200 and 50 sub-intervals, respectively. We note that the local approximation to the first order partial derivative ∂F ∂P can be easily obtained from the real option value. This quantity, known as the ∆ of an option, is important in practice, and the optimal exercise boundary can be also obtained from the value of ∆. The option value, ∆, and the optimal exercise boundary are plotted in Figure 2 and Figure 3, respectively.
To examine how the parameters impact the timing of investments, we plot the optimal exercise boundary about different parameters in Figure 4.
Evidently, from Figure 4 we can see that the investment boundary moves towards right with higher σ, ρ, I and C, and with lower a , which means that the non investment region (or continuation region in option's version) becomes bigger and the investment region (or stopping region in option's version) becomes smaller. In addition, the investment boundary is less sensitive to C and a than it is to σ, ρ  and I. Furthermore, the subjective discount rate ρ effects the investment boundary mostly compared with other parameters. Test 3: American real option with parameters [22]: P max = $1000, a = 86ML per year, C = $48091, I = $50121, α = 0.02, T = 1, σ = 0.2, r = 0.07, ρ = 0.07, δ = 0.05. The value of the project is chosen as (4).
Similarly, we divide (0, P max ) and (0, T ) uniformly into 200 and 50 sub-intervals, respectively. Note that in the case of (4), the derivative of V (P ) with respect to P is still a function of water price P . We plot the real option value, ∆ of the option, and the optimal exercise boundary in the following Figures 5 and 6.
Furthermore, the flexibility to suspend operation will affect the investment boundary, and the result is plotted in the following Figure 7. Clearly, from Figure 7 we can see that there is a bigger investment region for the case "option to suspend operation" than the case "no option to suspend operation". The potential to suspend    operation of the new technology helps to offset some of the risk associated with the technology, because the user is not locked into incurring the high operating costs. Test 4: American real option with subsidy. The parameters and the value of project are identical to Test 2. Note that the irrigators would have to relinquish 50% of the water savings from the use of new technology due to accepting the subsidy program, and we examine a range of scenarios for government cost sharing of capital costs: a 5% subsidy through to a 95% subsidy to the capital cost of the technology, so we set parameter a as 0.5a and parameter I in the interval [0.05I, 0.95I], respectively.
While the subsidy reduces the upfront capital cost of the technology, the irrigator must weigh this against the value of the water savings returned to the government. In this case, the irrigator's investment decision changes compared to the base case discussed in Test 2, where a lower capital cost but also a lower level of water savings would be faced after accepting the subsidy.
We plot the investment trigger price under different levels of subsidy at four time points: t = 0, 0.25, 0.5, 0.75 in Figure 8. The real option approach is equivalent to the NPV at maturity and the solid line can be seen as the trigger price at t = 1. Assuming that a temporary water price is $500 per ML, the NPV analysis would suggest a 30% subsidy to the cost to encourage investment. With the real options analysis, the subsidy to fasten the uptake of the technology is about 40%, 45%, 50%, and 55% at t = 0.75, t = 0.5, t = 0.25, and t = 0, respectively. This result potentially explains that the adoption rates sometimes remain lower than expected even when a subsidy to the capital cost of the investment is offered. At time point t = 0, the trigger price without the subsidy program in the base case Test 2 is about $450 per ML. To maintain this trigger price, a subsidy of almost 60% of the upfront cost would need to be offered. More than 60% of the capital cost would be offered in order to encourage the uptake of the new technology. The results are similar at other time points.
The fact that the program increases the trigger price for investment decisions implies that there may not be a strong intention to accept the subsidy that requires a return of permanent entitlement. The incentive to invest in the new technology and to retain all the on-farm water savings would be greater than participating in a program that offers a subsidy to capital costs in return for 50% of the water savings. This is because that the 50% water savings return can be seen as a permanent investment and then the irrigators lose the flexibility to manage water price risk.
6. Concluding remarks. Obviously, climate changes affect agricultural activities mostly and directly, such that it is increasingly important to manage water efficiently under higher temperatures and irregular precipitation patterns. For instance, developing the water saving technologies is one of the efforts. To this purpose, the choice of investment and its timing, which is affected by the uncertainty in water price, is essential in improving efficiency of water use.
In this paper, we utilize a real option approach to establish a linear parabolic variational inequality model to deal with future uncertainties about the water price. Especially, from the viewpoint of practice, instead of infinity, the expiration of the real option is considered to be finite in our model. Thus, solving the built model is quite difficult. Therefore, a power penalty approach is employed to solve it, and a nonlinear partial differential equation is obtained. Furthermore, a so-called fitted finite volume method is proposed to solve the resulting nonlinear PDE. To illustrate the efficiency and the usefulness of the method as well as the impacts of main parameters on the timing of investments, we perform several numerical experiments, which show that the subjective discount rate effects the investment boundary mostly, and the flexibility to suspend operation enlarges the investment region. In addition, the rates of subsidy, which are used to reduce the upfront cost of the new technology, are likely to be higher than those implied by the NPV analysis. Proof. For any w ∈ K, it follows from the definition of K that w − u ≤ 0 a.e. on Θ = Ω × (0, T ).
(30) On the other hand, from the third relationship of (11) we have Therefore, (30) Since θ ∈ [0, 1], we can see that (31) leads to which, together with the definition of the bilinear form B(·, ·; t) in (14), yields This proves the first inequality. The proof of the second inequality is standard and thus is omitted here.
Appendix C. Convergence analysis of the power penalty approach. We now show that as λ → ∞, the solution to Problem 2 converges to that of Problem 1 at the rate of order O(λ −k/2 ) in a proper norm. We start this discussion by the following lemma.
Lemma C.1. Let u λ be the solution to Problem 2. If u λ ∈ L p (Θ), then there exists a positive constant C, independent of u λ and λ, such that where k is the power of the power penalty function and p = 1 + 1/k.
Moreover, it follows from integration from t to T and using the above two inequalities as well as (36) that Thus, substituting the above inequality and (35) into (34), we obtain Now, from (37) and (38) we have from which the first inequality in (32) follows. Furthermore, from the above inequality we can also gain (φ(t), φ(t)) which implies the second inequality in (32). Using Lemma C.1, we are able to show that the solution to Problem 2 converges to that of Problem 1 at the rate of order λ −k/2 , as stated in the next theorem.