The Pricing of Asian Options in Uncertain Volatility Model

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.


Introduction
Model uncertainty aversion has shown to have important consequences in price behavior in capital markets [1][2][3] 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 firstorder 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.

Mathematical Problems in Engineering
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 [32][33][34], 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 twodimensional 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.

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 [33][34][35][36].
Let Ω be a given set and let H be a linear space of real valued functions defined on Ω.We suppose that H satisfies  ∈ H for each constant  and || ∈ H if  ∈ H.The space H can be considered as the space of random variables.
A sublinear expectation E on H is a functional E : H → R satisfying the following properties: for all ,  ∈ H, we have The triple (Ω, H, E) is called a sublinear expectation space.
From Peng [36], the distribution of the quadratic variation process of the -Brownian motion has the following properties.
The sublinear expectation Ê in the above theorem is called -expectation defined on   (Ω) via the following procedure.Let with   =   ,  ∈ [0, ∞) for  ∈ Ω, and Firstly, construct a sequence of -dimensional random vectors () ∞ =1 on a sublinear expectation space ( Ω, H, Ẽ) such that   is -normal distributed and  +1 is independent from ( 1 , . . .,   ) for each  = 1, 2, . ... For each  ∈   (Ω) with some  ∈  , ( × ) and And the related conditional expectation of under Ω   is defined by where under the norm 2.2.The PDEs.We assume the security market is defined on a sublinear expectation space (Ω, H, Ê) 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 (  ) ≥0 of some risky asset is given by where  is the -Brownian motion with  2 = − Ê[− 2  1 ] and  2 = Ê[ 2  1 ], and  and σ are constants.Hence ( ] ) ]≥0 is a geometric -Brownian motion: Define Then  ] satisfies the following stochastic differential equation: Let  ] = ( ] ,  ] )  , so where   and σ are 2 × 2 matrixes: The Let (,   ,   ) denote the price of one of the options at time .Take the arithmetic average fixed strike call option for example, we know (,   ,   ) = max(  / − , 0).
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 (Ω, H, Ê), Ê [⋅] 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 { Ê [⋅]} ≥0 ; that is, Let us assume with The following result is a corollary of the nonlinear Feynman-Kac formula of Peng [36].
where Φ : R 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 (, ) ∈ (0, ) × R 2 , let  ∈  2,3   ((0, ) × R 2 ) such that  ≥  and (, ) = (, ).By (29) and Taylor's expansion, for  ∈ (0,  − ), 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 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  = {(, , ) ∈ [0, ] × [0, ∞) × [0, ∞)}.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 (1/) −(−]) Δ]) 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  * = {(, , ) ∈ [0, ] × [0,  * ) × [0, )}.
On the boundary  = , the boundary condition from (51) can be written as We note that if   = 0, (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  =  max = S * .Hugger [37] gives the boundary condition at  =  max in the model without volatility uncertainty.By the same discussion, we can get the same boundary condition in uncertain volatility model: For the arithmetic average fixed strike Asian put option with terminal payoff (,   ,   ) = max( −   /, 0), we will value this option by the following put-call parity.
It is easily to see that we have max So ) .(60) 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 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  = 0 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 lim For computational purpose, we truncate the asset region (0, +∞) into (0,  max ), where  max is sufficiently large to ensure the accuracy of the solution.If  = 0, 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   = 0 and   = −1/.Substituting (72) into (71) and satisfying the boundary conditions, we find that For by the valuation mechanism in space (Ω, H, Ê), 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.

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 onedimensional 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.

Numerical Computation of Fixed Strike Asian Options.
An equivalent formulation of (48) in terms of the average   = (1/) ∫  0    is given in Barraquand and Pudet [16]  (77) 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  → 0. 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 convectiondominated 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.
With ADI, we march from time level  to time level  + 1/2 and then from time level  + 1/2 to time level  + 1.
In this case we use exponential fitting in all space variables Mathematical Problems in Engineering 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, ( 0 ,  0 , 0).We need to make sense of any value  ∈ [0,  max ] and stock price  ∈ [0,  max ] at time  = 0.There are no problems with the stock prices, but the integral at time zero can only be 0 for continuous averages; that is,  0 = 0.
Figure 1 is the results we obtained by using the above discretization scheme with σ = 0.8,  = 0.5,  = 0.10,  = 1, and  = 100.The point marked by * is the price (0, , 0) of the Asian option.We choose  max = 3 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 [0, ], 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 [0, ].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 σ = 0.8,  = 0.5, the Asian option PDE is the same as that of the Asian option with certain volatility  = 0.4.When the parameters  = 0.10,  = 1,  = 95, and grid Δ = 0.11875, ΔS = 2.85, 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 Δ = 0.01, 0.02, 0.04, Δ = 0.059375, 0.11875, 0.2375, and Δ = 2.85 for maturities of one quarter, half a year, and one year, respectively.Since  max = , 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]

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  = 0, while, from the definition of   , it is obvious that  0 = 0 and therefore  0 = 0, so our goal of the computation is to find the prices of the Asian options at  = 0; that is, (0, , 0) = (, 0).
To solve (67), we use the fully implicit positive coefficient discretization, which will ensure convergence to the viscosity solution.
Let    be the approximate values of the solution at   ,   ; that is,    ≃ (  ,   ).If 1 −   ≥ 0, we use forward difference to term   (, ).Otherwise, we use backward difference to term   (, ).
The general form of the discretized PDE (91) is with where
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 [42][43][44].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.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  max = 2.The execution times are less than 40 seconds.

Lemma 6 (Stability
Figure 2 is the values of function (0, ) for the arithmetic average floating strike Asian call option with the given parameters specified under the figure.The corresponding option values are (0, 0).

Table 2 :
[17]d strike Asian call values when  = 0.10 and  = 95.The values are the option prices at  = 100.Our results are denoted by Asian fixed strike, which calculated with Δ = 0.01, 0.02, 0.04, Δ = 0.059375, 0.11875, 0.2375, and Δ = 2.85 for maturities of one quarter, half a year, and one year, respectively.Z, F, and V refer to the result of Zvan et al.[17].Execution times (in seconds) are in parentheses.