Abstract

This paper studies the pricing of Asian options when the volatility of the underlying asset is uncertain. We use the nonlinear Feynman-Kac formula in the G-expectation theory to get the two-dimensional nonlinear PDEs. For the arithmetic average fixed strike Asian options, the nonlinear PDEs can be transferred to linear PDEs. For the arithmetic average floating strike Asian options, we use a dimension reduction technique to transfer the two-dimensional nonlinear PDEs to one-dimensional nonlinear PDEs. Then we introduce the applicable numerical computation methods for these two classes of PDEs and analyze the performance of the numerical algorithms.

1. Introduction

Model uncertainty aversion has shown to have important consequences in price behavior in capital markets [13] and macroeconomics [4, 5]. Denis and Martini [6] prove, in the model uncertainty situation the value of any bounded contingent claim can be expressed as the supremum of its expectation over some set of martingale measures.

In finance, there is an important case of model uncertainty called volatility uncertainty in which the uncertainty comes from the volatility coefficient. A major difficulty here is that the probabilities are mutually singular. This type of uncertainty is initially studied by Avellaneda et al. [7] and Lyons [8] for the superhedging of European options with payoffs depending only on the basic assets' terminal values. In this case the values of the derivatives should satisfy the (Black-Scholes-Barenblatt) BSB equations. The corresponding discrete-time case has been studied in Delbaen [9]. But for the path-dependent options, the difficulty is dramatically increased. In this paper, we are going to study the pricing of Asian options in the uncertain volatility model. As far as we know, there is little literature on this problem.

In the model without volatility uncertainty, there is a vast literature on pricing of Asian options.

In the geometric Brownian motion model, Asian options pricing methods are studied extensively. Geman and Yor [10] develop a dimensional Laplace transform method to price Asian options (see also [11]). Kemna and Vorst [12] use Monte Carlo simulation with a specific variance reduction method to compute the prices of fixed strike average-rate options. Ingersoll [13], and Wilmott et al. [14] prove that the two-dimensional PDE for a floating strike Asian option can be reduced to a one-dimensional PDE. Rogers and Shi [15] formulate a one-dimensional PDE that can model the prices of both floating and fixed strike Asian options and propose a numerical technique to compute tight lower bound on European average-rate options' prices. Barraquand and Pudet [16] describe a forward shooting grid algorithm and prove that it is unconditionally convergent. In Zvan et al. [17] and Zvan et al. [18], a flux limiter is used to retain accuracy while preventing oscillations. Simon et al. [19] propose an easy computable upper bound for the prices of arithmetic Asian options. Kaas et al. [20] and Dhaene et al. [21] compute the sharp lower and upper bounds for the prices of arithmetic Asian options by using the method of “comonotonic approximations.” Večeř [22], and Zhang [23, 24] propose various PDE methods. In Marcozzi [25], the first-order hyperbolic term is discretized by using a first-order upwind type method, resulting in at most first-order accuracy. A related approach based on a combination of the (weighted essentially nonoscillatory) WENO discretization and the grid stretching is used for Asian options in Oosterlee et al. [26]. Linetsky [27] investigates a spectral expansion approach.

In the jump diffusion models, methods for pricing Asian options are still under development. Fusai [28] obtains a simple expression for the double transform of the prices of continuously monitored Asian options by means of Fourier and Laplace transform. D’Halluin et al. [29] propose a semi-Lagrangian method to price the continuously observed fixed strike Asian options. Cai and Kou [30] consider the pricing of Asian options for the double exponential jump diffusion model. Bayraktar and Xing [31] obtain a fast numerical approximation scheme and analyze the performance of the numerical algorithm.

In this paper, we use the stochastic processes driven by -Brownian motion to describe the stocks' price processes. -Brownian motion is introduced in Peng [3234], and it has the properties that its increments are zero-mean, independent, and stationary and can be proved to be -normally distributed. The -normal distribution random variables have no mean uncertainty but can have variance uncertainty. The variance uncertainty of the -normal distribution is the source of the volatility uncertainty of the stocks' price processes.

The paper is organized as follows. In Section 2, we use the nonlinear Feynman-Kac formula in the -expectation theory to get the two-dimensional nonlinear PDEs satisfied by the values of the Asian options. For the arithmetic average fixed strike Asian options, we prove the convexity of the options' values with respect to the stock prices and transfer the nonlinear PDEs to linear PDEs with the uncertain volatilities being replaced by the maximum volatilities. For the arithmetic average floating strike Asian options, we use a dimension reduction technique to transfer the two-dimensional nonlinear PDEs to one-dimensional nonlinear PDEs. In Section 3, we introduce the applicable numerical computation methods for these two classes of PDEs and give the numerical results.

2. The Pricing of Asian Options in Uncertain Volatility Model

2.1. -Brownian Motion and -Expectation

In this subsection, we will introduce the -Brownian motion and -expectation defined in Peng [3336].

Let be a given set and let be a linear space of real valued functions defined on . We suppose that satisfies for each constant and if . The space can be considered as the space of random variables.

A sublinear expectation on is a functional satisfying the following properties: for all , we have(a)Monotonicity: If then .(b)Constant preserving: , for .(c)Sub-additivity: .(d)Positive homogeneity: .

The triple is called a sublinear expectation space.

Definition 1. A -dimensional process on a sublinear expectation space is called a -Brownian motion if the following properties are satisfied:(i);(ii)For each , the increment is -distributed and is independent from , for each and .

Remark 2. Just as in the classical situation, the increments of -Brownian motion are independent from .

Let such that for , be a sequence of partitions of with The quadratic variation process of 1-dimensional -Brownian motion with is defined as is an increasing process with .

Let be the space of all -valued continuous paths , with , equipped with the distance For each fixed , we set .

From Peng [36], the distribution of the quadratic variation process of the -Brownian motion has the following properties.

Lemma 3. For each fixed is identically distributed with and independent from .

Theorem 4. is -distributed; that is, for each ,

The sublinear expectation in the above theorem is called -expectation defined on via the following procedure. Let with for , and

Firstly, construct a sequence of -dimensional random vectors on a sublinear expectation space such that is -normal distributed and is independent from for each . For each with some and , set

And the related conditional expectation of under is defined by where consistently defines a sublinear expectation on and is a -Brownian motion.

We let , denote the completion of under the norm .

2.2. The PDEs

We assume the security market is defined on a sublinear expectation space which opens continuously and offers a constant riskless interest rate to all borrowers and lenders. And we also assume there are no transaction costs and/or taxes. Suppose that the price process of some risky asset is given by where is the -Brownian motion with and , and and are constants. Hence is a geometric -Brownian motion:

Define Then satisfies the following stochastic differential equation:

Let , so where and are matrixes:

The typical Asian options in this security market are the following four types:the arithmetic average fixed strike call option with final payoff, the corresponding arithmetic average fixed strike put option with payoff, the arithmetic average floating strike call option with final payoff, and the corresponding arithmetic average floating strike put option with payoff,

Let denote the price of one of the options at time . Take the arithmetic average fixed strike call option for example, we know .

Here we still assume the market operates in a risk neutral world [7].

Peng [35] introduces the axiomatic conditions satisfied by the valuation mechanism. In the sublinear expectation space , is a filtration consistent valuation operator. By this valuation mechanism, we have Here we assume the riskless interest rate is .

Since we assume the market is risk neutral, we have . It is easy to verify that the discounted prices of the stock satisfy the valuation mechanism ; that is,

Let us assume with

The following result is a corollary of the nonlinear Feynman-Kac formula of Peng [36].

Corollary 5. Let where is a given Lipschitz function and are given Lipschitz functions; that is, , for each ,  ,   and . And denote Then is a viscosity solution of the following PDE:

Proof. The proof of expression (29) is the same as that of Peng [36]. By the proof of Peng's theorem of the nonlinear Feynman-Kac formula and the boundedness of , we can easily prove that is a Lipschitz function; that is, Let us see the continuity of with respect to . Firstly, we have Replace with its Taylor's expansion in the first term on the right side of (33), we have By the following formulas [36] and (32), we have
For the second term on the right side of (33), we have Therefore is continuous in .
Now we are going to prove that is a viscosity solution of PDE (30). For fixed , let such that and . By (29) and Taylor's expansion, for , It is easy to check that Thus is a viscosity subsolution of (30). Similarly we can prove that is a viscosity supersolution of (30).
For the arithmetic average Asian options, we have Then by Corollary 5, satisfies the following PDE: By (18), the above PDE is where

2.3. The Arithmetic Average Fixed Strike Asian Call Option

Let us see the arithmetic average fixed strike Asian call option. From (14), we have Then by the definition of , we have Since By the convexities of function and the sublinear expectation, we can say that is, is convex about . So in the numerical computation, the second order difference of is positive. Then the numerical solutions of (43) are equivalent to those of PDE:

Theoretically, the above PDE is defined on domain . But generally, closed form solutions cannot be found; we have to find numerical solutions. We have no interest in finding the solution in the entire infinite domain. Therefore we consider the domain .

Now we determine the boundary conditions of (48). If , that is, , for the positiveness of , we have . Therefore the final payoff on the option is certain to be positive. At time , this payoff can be expressed as This payoff can also be obtained by using the following self-financing duplicating portfolio strategy. Assume that an investor commits into riskless bonds in order to secure the return promised by the first part of the right of (49), that is, at time . If the investor is to be certain of accruing the return promised by the second part of the right of (49), he is obliged to transfer a certain portion of stock (equal to ) to a riskless bond in every time interval . Overall this strategy requires a sum equal to plus the portion of a stock which can be expressed as follows: Then the price of the option when must be equal to

The above discussion is similar to that in the model without volatility uncertainty. Equation (51) obviously satisfies (48). So we only need to find the solution of (48) in the domain .

On the boundary , the boundary condition from (51) can be written as

We note that if , (13) tells us that is also equal to zero for and is thus equal to . By (46), we can write the following boundary condition:

We need another boundary condition at . Hugger [37] gives the boundary condition at in the model without volatility uncertainty. By the same discussion, we can get the same boundary condition in uncertain volatility model:

So, the values of the arithmetic average fixed strike Asian call option satisfy

For the arithmetic average fixed strike Asian put option with terminal payoff , we will value this option by the following put-call parity.

It is easily to see that we have So

is the time value of time claim , where the stock has certain volatility . To get claim at , let us see the replicating strategy: in any time interval , transfer a certain portion of stock to a riskless bond. This strategy ensures the final payoff of . The stocks that should be transferred are Then the time value of equals . So the put-call parity can be written as

2.4. The Arithmetic Average Floating Strike Asian Call Option

For the arithmetic average floating strike Asian call (put) option, the PDE (43) still applies with terminal condition

In Ingersoll [13] where there is no volatility uncertainty, due to the linear homogeneity of the PDE and its terminal condition, their PDE can be reduced to a one-dimensional PDE. In our uncertain volatility model, the PDE (43) is not linear homogeneous anymore. But we can still reduce PDE (43) to a one-dimensional PDE by the following argument.

For the arithmetic average floating strike Asian call option, we know

Denote Then by (14), we have and . So with .

Then the function satisfies the following PDE: which is a one-dimensional PDE.

For the European style Asian option, we must impose boundary conditions both at and as . The boundary condition as is simple. The only way that can tend to infinity is that tends to zero. In this case the option will not be exercised, and so For computational purpose, we truncate the asset region into , where is sufficiently large to ensure the accuracy of the solution. If , the PDE degenerates into the first-order hyperbolic PDE

So satisfies

By the same argument, for the arithmetic average floating strike Asian put option, its prices could be calculated by with satisfying the following PDE:

Generally, the put-call parity will not hold any more in the uncertain volatility model. Let us see the following argument. The mature payoff of a portfolio with one European arithmetic average floating strike Asian call holding long and one put holding short is The value of this portfolio is equal to one consisting of one stock and a financial product whose payoff is . The values of the financial product with final payoff satisfy the following PDE: with terminal condition . We seek a solution of (71) of the form with and . Substituting (72) into (71) and satisfying the boundary conditions, we find that For by the valuation mechanism in space , we take the discounted -conditional expectation on both sides and get That is where is the time price of the arithmetic average floating strike Asian call option and is the time price of the arithmetic average floating strike Asian put option.

3. Numerical Computation

In the model without volatility uncertainty, Asian options are much more difficult to value than regular stock options. Standard techniques tend to be impractical, inaccurate, or slow. For example, traditional binomial lattice methods require such enormous amounts of computer memory for the necessity of keeping track of every possible path throughout the tree that they are effectively unusable. Monte Carlo simulation works well for European style options (see [12]) but is relatively slow. Sun [38] proposed to use the weighted upwind method [39] to price Asian options, but he did not give numerical results. The Asian option pricing in the uncertain volatility model will be even more difficult than in the model without volatility uncertainty.

For the arithmetic average fixed strike Asian call options, due to the convexity of the value with respect to variable , (43) can be written as (48) which are the equations satisfied by the arithmetic average fixed strike Asian call options with certain volatility .

Zvan et al. [17] point out that the Crank-Nicolson scheme in conjunction with the van Leer flux limiter method can be used to numerically value Asian options with certain volatilities. The van Leer limiter has the property that it is total variation diminishing and thus produces oscillation free second order accurate solutions. But the van Leer limiter is nonlinear, so the solutions should be obtained by using full Newton iteration at each time level. In this paper, we will use the alternating direction implicit (ADI) methods (see Duffy [40]).

The arithmetic average fixed strike Asian put options can be valued by the put-call parity.

For the arithmetic average floating strike Asian call (put) options, we need to solve (67) (see (69)). These PDEs are one-dimensional in space. And they are not linear equations but HJB equations in some sense. For such PDEs, we will use the fully implicit positive coefficient finite different method [41] to solve them. We will see that this method can ensure the numerical solutions converge to the viscosity solutions, and is accurate, efficient, and quick to be implemented.

All the calculations were performed by Matlab 2009, on a computer with 1.83 GHz hard disk and 512 MB memory.

3.1. Numerical Computation of Fixed Strike Asian Options

An equivalent formulation of (48) in terms of the average is given in Barraquand and Pudet [16] in certain volatility model with volatility :

Equation (77) is convection dominated in the direction because there is no diffusion effect in this dimension. Barraquand and Pudet [16] find that an explicit centrally weighted scheme for (77) is unstable. In particular, the convective term in the dimension becomes very large as . Barraquand and Pudet also note that implicit centrally weighted schemes will generally produce unsatisfactory results because of the numerical diffusion introduced by this first order accurate in time scheme.

Zvan et al. [17] point out that under the condition of convection dominated, solutions generated by using a centrally weighted scheme on the convective term for (77) cannot be ensured to be free of oscillations. And solutions generated by the first order upstream weighting scheme on the convective term are no longer oscillatory, but their profiles are too diffuse. To produce oscillation free solutions without the excessive diffusion of first order upstream weighting, Zvan et al. [17] suggest the nonlinear van Leer flux limiter in conjunction with the Crank-Nicolson scheme. Since the flux limiter is nonlinear, the solutions should be obtained by using full Newton iteration at each time level.

The ADI method, pioneered in the United States by Douglas, Rachford, Peaceman, Gunn, and others, has a number of advantages. First, explicit difference methods are rarely used to solve initial boundary value problems because of their poor stability. Implicit methods have superior stability properties but unfortunately they are difficult to solve in two and more dimensions. Consequently, ADI methods become an alternative because they can be programmed by solving a simple trigonal system of equations. For the convection-dominated problems, we could use the exponentially fitted schemes in each leg of the approximate scheme to ensure the stability of the results [40]. Since Barraquand and Pudet [16] have used the forward shooting grid method and Zvan et al. [17] have used the van Leer flux limiter, we will test the ADI method in this paper by comparing our results with theirs.

For convenience, let us first convert (55) to be the following forward equation in time by substituting with :

We let ,  , be a set of portion points in satisfying and use even partition of time; that is, ,  . Define the partitions of space variables as and .

With ADI, we march from time level to time level and then from time level to time level . In this case we use exponential fitting in all space variables and implicit Euler in time. The first leg is given by the scheme with Let Then the scheme (79) can be written as We let Then we can write (82) as the following equivalent matrix form

The second leg of the discretized PDE in (78) is given by the scheme Let Then the scheme (85) can be written as

We let The matrix form of (87) is The matrix equations (84) and (89) can be solved by using LU decomposition.

The goal of the computation is to find the fair price of the Asian option at emission, that is, . We need to make sense of any value and stock price at time . There are no problems with the stock prices, but the integral at time zero can only be 0 for continuous averages; that is, .

Figure 1 is the results we obtained by using the above discretization scheme with ,  ,  ,  , and . The point marked by is the price of the Asian option. We choose to ensure the desirable accuracy.

Table 1 shows the convergence of the results with the refinement of the grid. From (78), the values of variable lie in the interval , where is the boundary and 0 is the point where the continuous averages Asian options are valued, so there is no point that is less important on the interval . Hence we use uniform grids in direction .

In our uncertain volatility model, the pricing PDEs (48) are similar to the PDEs in certain volatility model with volatility . Barraquand and Pudet [16] and Zvan et al. [17] give numerical results of Asian options with certain volatilities. In the uncertain volatility model, when ,  , the Asian option PDE is the same as that of the Asian option with certain volatility . When the parameters ,  ,  , and grid ,  , the price of the Asian call option we calculate is 13.82. The results with the same parameters of Barraquand and Pudet [16] and Zvan et al. [17] are 13.825 and 13.721, respectively.

Table 2 is the comparison of our results with those of Zvan et al. [17]. Z, F, and V in Table 2 refer to the result of Zvan et al. [17]. Our results are denoted by Asian fixed strike, which calculated with ,  , and for maturities of one quarter, half a year, and one year, respectively. Since , the partitions in direction are related with the maturity . The execution time is much less than theirs. The grid spacing was chosen to achieve an accuracy of at least 0.10% of . Our results are close to that of Zvan et al. [17] and the difference is less than 0.10% of .

Our calculation procedure is different from Zvan et al. [17], as well as Barraquand and Pudet [16], in the following four aspects. First, our pricing PDE (78) is in terms of the running sum . Zvan et al. [17] and Barraquand and Pudet [16] use the pricing PDE where is defined as . But PDE (78) and (90) are equivalent for pricing fixed strike European Asian options. Second, since the PDEs are different, PDE (78) and (90) must have different boundary conditions. We do not know the boundary conditions used by Zvan et al. [17]. The forward shooting grid method used by Barraquand and Pudet [16] proceeds without any boundary conditions. Third, our PDE (78) is defined on the domain , and is the boundary and is the point where the options are valued. So there is no point that is less important on the interval . Hence our numerical schemes use uniform grid spacing, while Zvan et al. [17] use nonuniform grid spacing of on domain . To achieve an accuracy of at least 0.1% of , our grids are much finer than Zvan et al. [17].

3.2. Numerical Computation of Floating Strike Asian Options

To value the arithmetic average floating strike Asian call (put) options, we need to solve (67) (see (69)) for , and the values of the options will be with . Actually, what we really care about is the option values at time , while, from the definition of , it is obvious that and therefore , so our goal of the computation is to find the prices of the Asian options at ; that is, .

To solve (67), we use the fully implicit positive coefficient discretization, which will ensure convergence to the viscosity solution.

Let , and then the PDE in (67) can be changed into the following forward equation in time:

Now, let us consider the discretization of (91). Define a grid and a set of timesteps . Let the discretization parameters be given by we use even discretization in time; that is, , .

Let be the approximate values of the solution at ,  ; that is, . If , we use forward difference to term . Otherwise, we use backward difference to term .

The general form of the discretized PDE (91) is with where

Of course, we have and . So the discretization equation (93) satisfies the positive coefficient condition of Forsyth and Labahn [41].

Let Then (93) can be written as the following matrix form

By the boundary condition we replace in with . Then (97) can be written as with

It is clear that is -matrix. From the property of -matrix, the solution of (99) exists and is unique.

It is important to ensure that we generate a numerical solution which converges to the viscosity solution. It has been shown that (67) satisfies the strong comparison property [4244]. Then, from Barles and Souganidis [45] and Barles [46], a numerical scheme converges to the viscosity solution if the method is consistent, stable, and monotone. Thus, we will show that our numerical scheme satisfies these conditions.

Lemma 6 (Stability). The discretization (93) is stable.

Proof. It follows from (93) that since .
Denote If , then we have Otherwise, or . Therefore, with ,   being the given boundary conditions.

Lemma 7 (Monotonicity). The discretization (93) is monotonic.

Proof. The lemma is trivially true on the boundary, so let us just see the cases and . Denote
For any , since and .
Similarly,

Lemma 8 (Consistency). The discretization (93) is consistent.

Proof. Suppose is a smooth test function with bounded derivatives of all orders with respect to . Let Then where if ; otherwise Using Taylor series expansions to (109), we have that is, (93) is consistent.

Therefore, we have the following.

Theorem 9. The solution of the discretization scheme (93) converges uniformly to the unique viscosity solution of PDE (91).

Table 3 shows the convergence of the values of the arithmetic average floating strike Asian call and put options with the refinement of the grid. We use uniform grids. To ensure the desirable accuracy, we choose . The execution times are less than 40 seconds.

Figure 2 is the values of function for the arithmetic average floating strike Asian call option with the given parameters specified under the figure. The corresponding option values are .

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

This work is supported by National Natural Science Foundation of China (11301011) and Beijing Natural Science Foundation (1112009). The authors would like to thank the anonymous referees for the comments.