Loading [MathJax]/jax/element/mml/optable/GeneralPunctuation.js
Research article Special Issues

Primal-dual active-set method for solving the unilateral pricing problem of American better-of options on two assets

  • In this paper, an efficient numerical algorithm is proposed for the valuation of unilateral American better-of options with two underlying assets. The pricing model can be described as a backward parabolic variational inequality with variable coefficients on a two-dimensional unbounded domain. It can be transformed into a one-dimensional bounded free boundary problem by some conventional transformations and the far-field truncation technique. With appropriate boundary conditions on the free boundary, a bounded linear complementary problem corresponding to the option pricing is established. Furthermore, the full discretization scheme is obtained by applying the backward Euler method and the finite element method in temporal and spatial directions, respectively. Based on the symmetric positive definite property of the discretized matrix, the value of the option and the free boundary are obtained simultaneously by the primal-dual active-set method. The error estimation is established by the variational theory. Numerical experiments are carried out to verify the efficiency of our method at the end.

    Citation: Yiyuan Qian, Haiming Song, Xiaoshen Wang, Kai Zhang. Primal-dual active-set method for solving the unilateral pricing problem of American better-of options on two assets[J]. Electronic Research Archive, 2022, 30(1): 90-115. doi: 10.3934/era.2022005

    Related Papers:

    [1] Yiyuan Qian, Kai Zhang, Jingzhi Li, Xiaoshen Wang . Adaptive neural network surrogate model for solving the implied volatility of time-dependent American option via Bayesian inference. Electronic Research Archive, 2022, 30(6): 2335-2355. doi: 10.3934/era.2022119
    [2] Yinghui Wang, Jiuwei Wang, Xiaobo Song, Yanpeng Hu . Linear convergence of a primal-dual algorithm for distributed interval optimization. Electronic Research Archive, 2024, 32(2): 857-873. doi: 10.3934/era.2024041
    [3] Yan Dong . Local Hölder continuity of nonnegative weak solutions of inverse variation-inequality problems of non-divergence type. Electronic Research Archive, 2024, 32(1): 473-485. doi: 10.3934/era.2024023
    [4] N. Bazarra, J. R. Fernández, R. Quintanilla . A dual-phase-lag porous-thermoelastic problem with microtemperatures. Electronic Research Archive, 2022, 30(4): 1236-1262. doi: 10.3934/era.2022065
    [5] Kaiyu Zhang . Sobolev estimates and inverse Hölder estimates on a class of non-divergence variation-inequality problem arising in American option pricing. Electronic Research Archive, 2024, 32(11): 5975-5987. doi: 10.3934/era.2024277
    [6] Ying Liu, Yanping Chen, Yunqing Huang, Yang Wang . Two-grid method for semiconductor device problem by mixed finite element method and characteristics finite element method. Electronic Research Archive, 2021, 29(1): 1859-1880. doi: 10.3934/era.2020095
    [7] Haiyan Song, Fei Sun . A numerical method for parabolic complementarity problem. Electronic Research Archive, 2023, 31(2): 1048-1064. doi: 10.3934/era.2023052
    [8] Wenya Qi, Padmanabhan Seshaiyer, Junping Wang . A four-field mixed finite element method for Biot's consolidation problems. Electronic Research Archive, 2021, 29(3): 2517-2532. doi: 10.3934/era.2020127
    [9] Xuefei He, Kun Wang, Liwei Xu . Efficient finite difference methods for the nonlinear Helmholtz equation in Kerr medium. Electronic Research Archive, 2020, 28(4): 1503-1528. doi: 10.3934/era.2020079
    [10] Qian He, Wenxin Du, Feng Shi, Jiaping Yu . A fast method for solving time-dependent nonlinear convection diffusion problems. Electronic Research Archive, 2022, 30(6): 2165-2182. doi: 10.3934/era.2022109
  • In this paper, an efficient numerical algorithm is proposed for the valuation of unilateral American better-of options with two underlying assets. The pricing model can be described as a backward parabolic variational inequality with variable coefficients on a two-dimensional unbounded domain. It can be transformed into a one-dimensional bounded free boundary problem by some conventional transformations and the far-field truncation technique. With appropriate boundary conditions on the free boundary, a bounded linear complementary problem corresponding to the option pricing is established. Furthermore, the full discretization scheme is obtained by applying the backward Euler method and the finite element method in temporal and spatial directions, respectively. Based on the symmetric positive definite property of the discretized matrix, the value of the option and the free boundary are obtained simultaneously by the primal-dual active-set method. The error estimation is established by the variational theory. Numerical experiments are carried out to verify the efficiency of our method at the end.



    Due to the increased awareness of risk aversion, the financial derivatives (e.g., options) have attracted more and more attention. An option is a derivative, which can be exercised on the maturity date (European option) or at any time prior to the maturity date (American option). Based on the geometric Brownian motion, the well-known Black-Scholes equation (BS equation) for pricing options was established in 1973. Whereafter, a variety of closed-form expressions have been derived for most European options [1,2,3]. Up to now, there has been no similar results for American options yet. By virtue of the flexibility in the choice of the exercise time, American options have been developed in the option market rapidly. Originally, traders and researchers focused on single-asset options, which are of simple structures [4]. With the development of the financial market and the increasing demand of investors, plentiful multi-asset options have emerged, and the corresponding American options were studied intensively. Although multi-asset American options can be divided into some categories, and may cover any number of assets, the better-of option on two assets, which belongs to the rainbow option family, is typical and popular because of its simplicity and practicality.

    The better-of option was first named "option on the maximum of two risky assets" by Stulz in 1982 [5], and then was called its present name by Jiang [6]. By exploiting the relationship of the pricing formulas for two single asset options, Stulz presented the closed-form solution for the European better-of option. Jiang stated that the American better-of option satisfies a multi-dimensional BS equation with a special payment function. In most instances, the underlying assets of the options pay dividends or have other cash outflows. As is well known, the standard American option written on a single asset that pays dividends can be exercised optimally before the maturity date [7]. And so are multi-asset dividend paying options [8]. When the underlying assets pay continuous dividend yields, the set of optimal prices of multi-asset options with respect to the time forms a continuous surface, which is called the optimal exercise boundary. In this paper, we consider a special two-asset American better-of option pricing problem, where only one underlying asset pays dividends. The optimal exercise boundary for the concerned model divides the solving domain into two parts: one is the bounded region called the holding domain, and the other is the unbounded one called the exercising domain.

    In practice, various methods have been used to price options accurately and efficiently. They can be classified into two categories: analytical approximations [9,10] and numerical solutions [11]. On account of no analytical formula, the pricing problem for American options has to be solved numerically. Binomial method (BM), Monte Carlo (MC) simulations, finite difference methods (FDM) and finite element methods (FEM) [12] are commonly used numerical methods for pricing American options. BM is the most classic numerical technique for options. In 1979, Cox, Ross and Rubinstein first applied the binomial model to price American options [13]. Five years later, Amin and Khanna proved the convergence of this method [14]. With the increase in the number of the underlying assets, the high computational cost and the slow convergence rate incurred by BM. Based on the pioneering works of Boyle, Bossaerts and Tilley, MC simulations are frequently used to the American option pricing problems [15]. The advantage of MC is that it does not depend on the dimension of the problem and thus does not suffer from the curse of dimensionality. As a generalization of BM, FDM is introduced to solve the pricing problems and is widely used because of its simple form. The FDM is easy for implementation and efficient for problems with regular domain. While for problems with irregular domain, FEM is more flexible. It can handle the irregular domain as well as complex boundary conditions efficiently. Meanwhile, FEM has a solid theoretical framework with robust numerical performance. Therefore, FEM has became more and more popular in classical option pricing problem [16], and has been successfully used to deal with better-of options recently [17,18]. In this paper, we shall adopt FEM to discretize the two-asset American better-of option pricing problem when only one asset pays dividend.

    The concerned two-asset American better-of option in this paper satisfies a two-dimension parabolic linear complementary problem (LCP) on an unbounded domain. To solve this model efficiently, we must think about the following issues: 1) In order to better characterize the singularity of the option price near the maturity date, the design of the spatial mesh must be considered carefully. 2) Due to the solution domain being unbounded, it is difficult to design algorithms directly. The truncation method and the corresponding boundary condition, which can guarantee the accuracy of the solution, are important. 3) For the dependency relationship of the option price and the exercise boundary, the pricing model becomes a highly nonlinear problem. The efficiency of numerical algorithm used to settle the option price and the optimal exercise boundary simultaneously is the key issue.

    As regards the first difficulty, the geometric grid for dealing with the classical single asset American options is adopted to guarantee the accuracy of the better-of option around the singularity. For the second issue, there are two mainstream techniques: the transparent boundary condition method and the far-field boundary condition method. The former one was proposed to solve classical American options by Han and Wu in 2003 [19], and improved by Ehrhardt and Mickens in 2008 [20], which is suitable for use in combination with FDM. The latter method was proposed by Kangro and Nicolaides [21]. They presented a reasonable boundary condition along the truncation boundary, and established the pointwise estimate of this truncation method via the comparison principle. Whereafter, many works for pricing American options using this method have been published [16,18]. We shall follow this approach to deal with our problem. Regarding the third issue, the primal-dual active-set (PDAS) method shall be used to solve it efficiently. The PDAS method is a special case of the generalized Moreau-Yosida approximations [22], and is also proved to be a special case of the Newton-type method under some moderate conditions [23]. It is efficient for quadratic programming problems and LCPs. Moreover, the global convergence of PDAS method has been proved by Bergounioux et al. in [24]. For the applications of PDAS on option pricing problems, we refer to [17,25] and references therein for the rich literature.

    The concerned model in this paper is the same as the problem in [18], but the numerical algorithm is very different. Although both methods can be used to obtain the option price and the optimal exercise boundary simultaneously. Compared with the Newton's method in [18], our proposed algorithm has obvious advantages in both theoretical analysis and numerical performance. Firstly, the convergence of our method can be analysed, while the Newton's method in [18] can not be. Secondly, PDAS method as an active-set strategy could effectively reduce the number of the iterations in each time layer, and leads to rapid convergence, which can be verified in the numerical simulations. Moreover, PDAS method is faster than the projected successive overrelaxation method. Therefore, the proposed algorithm in this paper is an efficient method for pricing two-asset American better-of options when only one asset pays dividends.

    The arrangement of this paper is as follows. In Section 2, the variational inequality (VI) model for American better-of options is introduced firstly. Based on some conventional numeraire transformations and far-field technique, the bounded LCP corresponding to the VI is established. In Section 3, the full discretization model of the LCP is constructed by FDM with uniform partition in temporal direction and FEM with geometric mesh in spatial direction, respectively. Then, the PDAS method is adopted to solve the resulting discretized optimization system. The convergence analysis and numerical experiments are implemented in Section 4. The concluding remarks are given in Section 5.

    In this section, we shall introduce the VI model for American better-of options on two underlying assets. Furthermore, based on some changes of variables and the far-field truncation technique, the associated LCP shall be presented.

    For the simplicity of the expression, we introduce some notations firstly. Si, σi and qi(i=1,2) stand for the price, the volatility, and the dividend of the i-th underlying asset, respectively. The correlation coefficient of these two assets is denoted by ρ. It needs to point out that we only consider the case where only one underlying asset pays dividend, as the case in which both underlying assets pay dividends have been discussed exhaustively [17]. Let r, t, and T be the interest rate, the arbitrary time, and the maturity date, respectively. Assume that there is no transaction fee to pay and no arbitrage on the market. Then, the VI model corresponding to the American better-of option price P=P(S1,S2,t) with q1=0 and q2>0 can be described as

    {min{PtLP,Pf(S1,S2)}=0,(S1,S2,t)Σ,P(S1,S2,T)=f(S1,S2),(S1,S2)R2+, (2.1)

    where the payoff function f(S1,S2)=max(S1,S2), the solution domain

    Σ={(S1,S2,t)(S1,S2)R2+,t[0,T)},

    and the operator

    LP=12(σ21S212P2S21+2ρσ1σ2S1S22PS1S2+σ22S222P2S22)+rS1PS1+(rq2)S2PS2rP.

    here, R2+={(S1,S2)Si(0,+),i=1,2}.

    The solution domain Σ of the pricing model (2.1) can be divided into two parts by the optimal exercise boundary, where the payoff of an option holder reaches the maximum. Since q1=0 and q2>0, the exercising domain denoted by Σ2 is on the left of the optimal exercise boundary. When the prices of the underlying assets belong to Σ2, the price of the option is equal to the payoff value by exercising the option, i.e., P(S1,S2,t)=f(S1,S2). Otherwise, the holders may suffer a loss. On the other hand, the holding domain denoted by Σ1 is on the right of the optimal exercise boundary. On this subdomain, the price of the option is unknown, and satisfies P(S1,S2,t)>f(S1,S2). Therefore, the model (2.1) can be rewritten as

    {Pt+LP=0,P>f(S1,S2),(S1,S2,t)Σ1,Pt+LP0,P=f(S1,S2),(S1,S2,t)Σ2. (2.2)

    It is worth to notice that there is an unknown interface between Σ1 and Σ2. In this paper, we denote it by Γ. Hence, the pricing model (2.1) is equivalent to a two-dimensional free boundary problem.

    By means of traditional numeraire transformations

    s=S1S2,p(s,t)=P(S1,S2,t)S2, (2.3)

    the pricing model (2.1) can be transformed into the following one-dimensional VI model

    {min{ptˆLp,pf(s,1)}=0,(s,t)ˆΣ,p(s,T)=f(s,1),s(0,+), (2.4)

    where the solution domain after the transformations

    ˆΣ={(s,t)s(0,+),t[0,T)},

    and the simplified operator

    ˆLp=12ˆσ2s22ps2+q2spsq2p,ˆσ2=σ212ρσ1σ2+σ22.

    Similarly, ˆΣ1, ˆΣ2, and ˆΓ represent the corresponding holding domain, the exercising domain, and the optimal exercise boundary, respectively. The one-dimensional VI model (2.4) can be reformulated as

    {pt+ˆLp=0,p>f(s,1),(s,t)ˆΣ1,pt+ˆLp0,p=f(s,1),(s,t)ˆΣ2. (2.5)

    Here, ˆΣ1={(s,t)s(γ(t),+),t[0,T)}, ˆΣ2=ˆΣˆΣ1 and ˆΓ=γ(t) is the corresponding unknown free boundary. Based on the above discussions, we shall only solve the following unilateral free boundary problem

    {pt+ˆLp=0,(s,t)ˆΣ1,p(s,T)=f(s,1),s(0,+),p(γ(t),t)=1,t[0,T),ps(γ(t),t)=0,t[0,T),lims+(p(s,t)s)=0,t[0,T). (2.6)

    Furthermore, let

    W(s,t)=p(s,t)s, (2.7)

    then the free boundary problem model (2.6) shall be transformed into a standard American put option pricing problem

    {Wt+ˆLW=0,(s,t)ˆΣ1,W(s,T)=f(1s,0),s(0,+),W(γ(t),t)=1γ(t),t[0,T),Ws(γ(t),t)=1,t[0,T),lims+W(s,t)=0,t[0,T), (2.8)

    where the optimal exercise boundary γ(t) is a bounded and nondecreasing function [26], and satisfies

    βα11α2α11β1α2α2α12=γ0:=γ(0)<γ(t)<γ(T)=1. (2.9)

    here, βi=αi1αi(i=1,2), α1<0, and α2>0 are the solutions of the following equation

    12ˆσ2α(α1)+q2αq2=0.

    In order to solve the free boundary problem model (2.8) numerically, the unbounded domain must be truncated firstly. Following the ideas of Holmes and Yang [16], we shall introduce a far-field boundary method to deal with this problem, by which accurate solutions can be obtained efficiently.

    Lemma 1. For a given positive number ε(0,1), let ψ=12ˆσ2, α0=rqˆσ212, then we have

    W(s,t)ε,s[eˆZ,+),t[0,T),

    where

    ˆZ=2ψTα0+2α20ψ2T2ψTlnε.

    Proof. By applying the variable substitutions

    w(z,η)=eαz+βηW(s,t),η=Tt,z=lns, (2.10)

    the free boundary problem model (2.8) can be transformed into

    {wηψ2wz2+νWz+μw=0,z(b(η),+),η(0,T],w(z,0)=h(z,0),z(,+),w(b(η),η)=h(b(η),η),η(0,T],wz(b(η),η)=gz(b(η),η),η(0,T],limz+eαzβηw(z,η)=0,η(0,T], (2.11)

    where the coefficients ψ=12ˆσ2, ν=2αψq2, and μ=(α+1)q2ψα2β, and the functions

    h(z,η)=eαz+βηmax(1ez,0),b(η)=ln(γ(Tη)).

    From the property of γ(t), we can infer that b(η) is a bounded and nonincreasing function with

    lnγ0<b(η)<b(0)=0,η(0,T]. (2.12)

    Taking α=α0 and β=ψα20+q2, yields ν=μ=0. Therefore, the pricing problem model (2.8) becomes a heat equation. By virtue of the property of b(η) in model (2.12), the first equality of model (2.8) holds when z>0. Consequently, by the Theorem 19.3.2 in [27], we obtain

    w(z,η)=z4ψπη0(ηx)32ez24ψ(ηx)w(0,x)dx,z>0.

    The pricing problem model (2.8) is equivalent to a standard American put option pricing problem with the strike price K=1. Based on the fact that the value of the American put option is always less than or equal to the strike price and the transformations model (2.10), we can get w(0,η)eβη, which yields

    w(z,η)zeβη4ψπη0(ηx)32ez24ψ(ηx)w(0,x)dx2eβηz24ψηπ+0ey2dy=eβηz24ψηeβηz24ψT,z>0.

    Furthermore, by the inverse transformations of model (2.10), we have

    W(s,t)=eαzβηw(z,η)eαzz24ψT.

    Finally, from the definition of ˆZ=2ψTα0+2α20ψ2T2ψTlnε, it is easy to verify

    eαzz24ψTε,zˆZ,

    which implies the conclusion of this Lemma.

    Now, using the error bound of Lemma 1, we truncate the solution domain. For a fixed ε(0,1), let

    L=max{lnγ0,ˆZ}, (2.13)

    then the free boundary problem model (2.8) can be approximated by the following LCP

    {(Wt+ˆLW)(Wf(1s,0))=0,s[eL,eL],t(0,T],W(s,T)=f(1s,0),s[eL,eL],W(eL,t)=f(1eL,0),t(0,T],W(eL,t)=f(1eL,0),t(0,T], (2.14)

    with constraints

    Wt+ˆLW0,Wf(1s,0).

    It needs to point out that the choice of L should ensure eLγ(t) and eLeˆZ, which makes sure that the left boundary condition of model (2.14) is accurate, and the right boundary condition is reasonable.

    The LCP model (2.14) is a backward variable coefficients problem on a bounded domain. To solve it efficiently, we change it to a forward constant coefficient problem by the transformations

    v(x,τ)=eax+bτW(s,t),τ=ˆσ22(Tt),x=lns, (2.15)

    with a=2q2ˆσ22ˆσ2, b=12ˆσ2a2+q2. The specific form of the bounded LCP is as follows

    (BLCP){(vτ2vx2)(vq)=0,x[L,L],τ(0,˜T],v(x,0)=q(x,0),x[L,L],v(L,τ)=q(L,τ),τ(0,˜T],v(L,τ)=q(L,τ),τ(0,˜T], (2.16)

    with constraints

    vτ2vx20,vq.

    Here, q(x,τ)=eax+bτf(1ex,0) and ˜T=ˆσ22T. It is easy to get that the relationship of the optimal exercise boundaries between the linear complementary problems (2.14) and (2.16) is

    B(τ)=ln(γ(T2τˆσ2)).

    So far, we have transformed the original variable coefficients pricing model on a two-dimensional unbounded domain into a BLCP on a one-dimensional bounded domain, on which we can design numerical algorithms directly.

    In this section, the variational problem associated with the BLCP model (2.16) is presented firstly. Then, the resulted problem shall be discretized by the FDM and the FEM in temporal direction and spatial direction, respectively. Based on the properties of the discrete system, the PDAS method is proposed to obtain the option price and the optimal exercise boundary simultaneously. The error estimates of the proposed approach are also given at the end of the section.

    Before presenting the variational problem corresponding to the BLCP model (2.16), we ought to introduce some function spaces and function sets for subsequent applications. L2(I) stands for the space of square integrable functions on I=[L,L], and define

    H1(I)={uL2(I)uxL2(I)},H2(I)={uL2(I)uxL2(I),uxxL2(I)},H1τ(I)={uH1(I)uq(x,τ),u(L)=q(L,τ),u(L)=q(L,τ)},H1(I):thedualspaceofH1(I).

    Lemma 2. (cf. [28]) If vL2([0,˜T];H2(I))and vτL2([0,˜T];L2(I)), thenv is the solution of the BLCP model (2.16) if and only ifv is the solution of the following variational inequality

    (VI) Find v(,τ)H1τ(I), such thatv(x,0)=q(x,0) and

    (vτ,uv)+(vx,uxvx)0,uH1τ(I),a.e.τ(0,˜T]. (3.1)

    In order to apply FDM and FEM to the VI model (3.1), some notations must be introduced in advance. The spatial partition of I=[L,L] and the temporal partition of J=[0,˜T] are defined as follows:

    Ix:L=x0<x1<<xN=L,Jτ:0=τ0<τ1<<τM=˜T.

    Let Ii:=(xi1,xi) denote the spatial element and hi:=xixi1,i=1,,N represent the interval length. The grid size of the spatial partition is denoted by h:=max1iNhi. In a similar way, for each temporal element Jj:=(τj1,τj),τj:=τjτj1,j=1,,M, and τ:=max1jMτj represent the local and overall step size, respectively.

    What we should pay attention to is that options are traded more frequently near the optimal exercise boundary Γ. As regards this issue, the most desirable points for the option price are around the optimal exercise boundary with S1=S2. In other words, we ought to set most nodes near the point s=1 or x=0 by the numeraire transformations models (2.3) and (2.15). Therefore, in the spatial direction, we shall use a geometric grid partition and an even number of intervals N, which ensure the nodes should be symmetrical about xN2=0. And we resort an isometric subdivision in temporal direction. That is to say,

    xi=sign(2iN)(2iNN)2L,i=0,1,,N,τj=jτ,τ=˜TM,j=0,1,,M. (3.2)

    The main goal of this subsection is to present the discretization scheme of the VI model (3.1). First, we define the set of piecewise linear functions as follows

    S1τ(I)={uH1τ(I)u(xi)q(xi,τ),u(xi)IiP1,i=1,2,,N},

    where P1 represents the set of polynomials of degree less than or equal to 1. Thus, the semi-discretized approximation of the VI model (3.1) by FEM can be described as

    (SDA) Find vh(x,τ)S1τ(I), such that vh(x,0)=qI(x,0) and

    (vhτ,uhvh)+((vh)x,(uh)x(vh)x)0,uhS1τ(I),τ(0,˜T]. (3.3)

    Here, qI(x,0) stands for the piecewise linear interpolation of q(x,0) in S10(I). Next, the backward Euler method is applied to the SDA model (3.3) at a fixed τ=τj,j=1,,M, the full-discretized approximation of VI model (3.1) is presented as

    (FDA) Find vjh(x)S1τj(I), such that v0h(x)=qI(x,0) and

    (vj1h,uhvjh)+b(vjh,uhvjh)0,uhS1τj(I), (3.4)

    where vjh(x)=vh(x,τj),j=1,,M, the operator and the bilinear function b(,) are defined as follows

    uj:=uj+1ujτ,b(u,v):=(ux,vx).

    For the convenience of algorithm design, we shall derive the matrix-vector form of the FDA model (3.4). Suppose that the basis function set of S1τ(I) is {ϕ0,ϕ1,,ϕN}, and the specific representations of the basis are

    ϕ0(x)={xx1x0x1,x[x0,x1),0,xI[x0,x1),ϕi(x)={xxi1xixi1,x[xi1,xi),xxi+1xixi+1,x[xi,xi+1),0,xI[xi1,xi+1),ϕN(x)={0,xI[xN1,xN),xxN1xNxN1,x[xN1,xN).

    Then, the finite element solutions can be reformulated as

    vjh(x)=N1i=1vjiϕi(x)+q(L,τj)ϕ0(x)+q(L,τj)ϕN(x),j=0,1,,M.

    By substituting the variables vjh in the FDA model (3.4) by the above formula, we obtain the following inequality

    (UVj)T((τ\bmΨ+\bmΦ)Vj\bmΦVj1+Hj)0,UQj, (3.5)

    where j=1,,M,

    U=(u1,u2,,uN1)T,Vj=(vj1,vj2,,vjN1)T,Qj=(q(x1,τj),,q(xN1,τj))T,Hj=((h16τh1)q(L,τj)h16q(L,τj1),0,,0,(hN6τhN)q(L,τj)hN6q(L,τj1))T,

    and

    \bmΦ=(h1+h23h26000h26h2+h33h36000000hN26hN2+hN13hN16000hN16hN1+hN3),\bmΨ=(1h1+1h21h20001h21h2+1h31h30000001hN21hN2+1hN11hN10001hN11hN1+1hN).

    Let C=τ\bmΨ+\bmΦ and Dj=\bmΦVj1Hj, therefore the inequality model (3.5) can be simplified as

    (MVIP)(UVj)T(CVjDj)0,UQj,j=1,,M. (3.6)

    In the above matrix-vector inequality model (3.6), we can verify that the matrix C is positive definite when h2τ is small enough.

    In this subsection, to get the option price and the optimal exercise boundary simultaneously, we adopt PDAS method to solve the MVIP model (3.6). The detailed algorithm is presented in Algorithm 1:

    Algorithm 1 Primal-Dual Active-Set algorithm.
    Input: Vj1, Qj, Dj;
    Output: Vj, xj.
    1: The given parameters: εV>0,δ>0,\bmλ0.
    2: Vpre:=0,Vcur:=max{Vj1,Qj}.
    3: while , do
    4:    .
    5:    ,
    6:    ,
    7:    .
    8:    ,
    9:    ,
    10:    .
    11: end while
    12: .
    13: ,
    14: .

     | Show Table
    DownLoad: CSV

    For , when the numerical solution is obtained, we can get . By means of the transformations models (2.15), (2.7), and (), we can calculate the approximation of , , and . Likewise, we obtain the optimal exercise boundary at . We present the whole process in Algorithm 2 as:

    Algorithm 2 The whole algorithm.

    1: Initial parameters setting: .
    2: Additional parameters calculating: .
    3: Calculate the parameters about boundary: .
    4: Truncate the solving domain: .
    5: Partition:
    6:    ,
    7:    , , ,
    8:    .
    9: Setting:
    10:    ,
    11:   
    12:   
    13: Calculate matrixes: , , .
    14: for , do
    15:    Calculate: , , ,
    16:    Solve: , by PDAS (Algorithm 1).
    17:    ,
    18:    .
    19: end for
    20: .
    21: Calculate: ,
    22: Calculate: ,
    23: Calculate: .

     | Show Table
    DownLoad: CSV

    The error of our proposed method mainly comes twofold: The error between the original problem model (2.1) and the BLCP model (2.16), which has been presented in Lemma 1. The error of the BLCP model (2.16) and the MVIP model (3.6). In this subsection, we focus on the error estimate of . In order to obtain the result of estimation, we introduce some lemmas firstly.

    Lemma 3. (cf. [29]) If and is thepiecewise linear interpolation of on the grid , then wehave

    For any , let , then the following results hold

    Lemma 4. (cf. [30]) Suppose that, , and .For any , let , then we have

    Let represent the linear interpolation of in the temporal direction, based on the Lemmas 3 and 4, we establish the following result.

    Theorem 1. Under the assumptions of Lemma 4, if in addition, , , thenwe have the error estimate as follows

    Proof. For , we can calculate the following equality firstly

    (3.7)

    Let at in VI model (3.1) and in FDA model (3.4), then we obtain

    (3.8)

    Next, we add the two inequalities in model (3.8) to both sides of model (3.7), we have

    (3.9)

    For , define

    then, we can simplify the inequality model (3.9) as below

    Thus, we can get

    (3.10)

    where are defined in Lemma 4. Therefore, the first part of the above inequality can be estimated by

    (3.11)

    Based on the inequalities (3.10) and (3.11), we have

    By virtue of Poincare inequality and , the following result holds

    Moreover, by exploiting the complex trapezoidal formula, we can conclude that

    which leads to the conclusion directly [31].

    In this section, numerical experiments are performed to illustrate the advantages of our proposed algorithm. In order to show the superiority of our algorithm in computational accuracy and speed, we shall compare it with the BM and the Newton's method (NM) [18].

    Let us consider a one-year () American better-of option with . Without loss of generality, the parameters of the first underlying asset are set to be , , the parameters of the other asset are set to be , , the correlation coefficient of these two assets is set to be . Suppose the truncation accuracy in model (2.14) via the definition model (2.13), the estimate model (2.9), and Lemma 1, we obtain a table about the truncation length in deferent cases.

    For convenience, we choose in all experiments, which can ensure the truncated domain cover all the solution domains corresponding to Table 1. Because there are no closed-form solution for American better-of options, we shall apply the solutions obtained by BM with 100,000 points in temporal direction as the benchmarks. The parameters in PDAS method are set to be , , . All the experiments are implemented by Matlab R2014a and on a computer with Intel Core i5 CPU of 2.2GHz.

    Table 1.  The truncated lengths calculated by model (2.13) with .
    0.3 0.4 0.3 0.4 0.3 0.4 0.3 0.4
    2.7220 2.7812 2.6329 2.7167 2.6038 2.6822 2.4712 2.5776

     | Show Table
    DownLoad: CSV

    In this subsection, we mainly verify the superiority of PDAS method in computing the optimal exercise boundaries. The error is measured by the -norm of , where and are the numerical solutions and the benchmark solutions, respectively. With and , Tables 2 and 3 compare the accuracies achieved by running BM, NM, and PDAS method for about the same amount of time (). Similarly, Tables 4 and 5 list the computation costs of the using above three methods in different cases.

    Table 2.  The error estimates on of BM, NM, and PDAS with .
    Error / Time/s
    BM NM PDAS BM NM PDAS BM NM PDAS BM NM PDAS
    0.02 0.3, 0.4 60 100 300 409 100 350 12.156 34.654 1.631 0.135 0.130 0.131
    0.4, 0.4 60 100 300 338 100 350 15.566 58.256 1.833 0.117 0.118 0.120
    0.03 0.3, 0.4 60 100 300 409 100 350 12.596 18.814 1.459 0.138 0.131 0.131
    0.4, 0.4 60 100 300 338 100 350 15.654 30.677 1.641 0.119 0.118 0.114

     | Show Table
    DownLoad: CSV
    Table 3.  The error estimates on of BM, NM, and PDAS with .
    Error / Time/s
    BM NM PDAS BM NM PDAS BM NM PDAS BM NM PDAS
    0.02 0.3, 0.4 60 150 300 528 100 350 10.315 18.142 1.339 0.159 0.151 0.146
    0.4, 0.4 60 150 300 451 100 350 11.580 26.727 1.494 0.145 0.142 0.142
    0.03 0.3, 0.4 60 150 300 528 100 350 09.950 10.508 1.225 0.159 0.159 0.158
    0.4, 0.4 60 150 300 451 100 350 11.039 14.879 1.411 0.146 0.143 0.144

     | Show Table
    DownLoad: CSV
    Table 4.  The computation costs on of BM, NM, and PDAS with .
    Error / Time/s
    BM NM PDAS BM NM PDAS BM NM PDAS BM NM PDAS
    0.02 0.3, 0.4 3000 2000 300 2893 1500 350 1.496 1.678 1.631 48.069 4.225 0.131
    0.4, 0.4 3000 2000 300 2396 2500 350 1.779 1.786 1.833 40.545 6.665 0.120
    0.03 0.3, 0.4 3000 2000 300 2893 900 350 1.524 1.509 1.459 46.583 3.045 0.131
    0.4, 0.4 3000 2000 300 2396 1300 350 1.813 1.709 1.641 41.936 3.624 0.114

     | Show Table
    DownLoad: CSV
    Table 5.  The computation costs on of BM, NM, and PDAS with .
    Error / Time/s
    BM NM PDAS BM NM PDAS BM NM PDAS BM NM PDAS
    0.02 0.3, 0.4 3000 2000 300 3740 1000 350 1.199 1.270 1.333 58.381 3.233 0.142
    0.4, 0.4 3000 2000 300 3195 1300 350 1.375 1.458 1.518 49.031 3.883 0.142
    0.03 0.3, 0.4 3000 2000 300 3740 600 350 1.215 1.265 1.252 59.042 2.327 0.158
    0.4, 0.4 3000 2000 300 3195 800 350 1.387 1.325 1.429 50.168 2.789 0.144

     | Show Table
    DownLoad: CSV

    From the results listed in Tables 25, we can conclude that with the similar computing time (), the errors computed by PDAS method are one order of magnitude smaller than BM and NM. with the similar computing time around (), the computational efficiency of PDAS is almost and times of BM and NM, respectively. Figure 1 shows the optimal exercise boundaries solved by PDAS and the benchmark BM. It is easy to find out that the numerical solutions approximate the benchmarks very well. All of the results confirm that PDAS method is an efficient algorithm to obtain the optimal exercise boundary. At the end of this subsection, to give the investors some intuitive and valuable guidance, we present the original optimal exercise surfaces of model (2.2) in Figure 2 with different parameters, respectively.

    Figure 1.  The optimal exercise boundaries solved by PDAS method and the benchmark BM with . (top left), (top right), (bottom left), (bottom right). and denote the holding domain and the exercising domain, respectively.
    Figure 2.  The optimal exercise boundaries solved by PDAS method with . (top left), (top right), (bottom left), (bottom right).

    In this subsection, we illustrate the efficiency of the PDAS method for solving option prices. Similar to the analysis in previous subsection, we compare the error estimates and computation costs in Tables 69.

    Table 6.  The error estimates on of BM, NM, and PDAS with .
    Error / Time/s
    BM NM PDAS BM NM PDAS BM NM PDAS BM NM PDAS
    0.02 0.3, 0.4 60 100 300 409 100 350 6.073 409.554 0.823 0.135 0.130 0.131
    0.4, 0.4 60 100 250 338 100 350 6.115 444.089 0.975 0.117 0.118 0.115
    0.03 0.3, 0.4 60 100 300 409 100 350 5.885 514.817 0.651 0.138 0.131 0.131
    0.4, 0.4 60 100 300 338 100 350 6.507 405.355 0.984 0.119 0.118 0.114

     | Show Table
    DownLoad: CSV
    Table 7.  The error estimates on of BM, NM, and PDAS with .
    Error / Time/s
    BM NM PDAS BM NM PDAS BM NM PDAS BM NM PDAS
    0.02 0.3, 0.4 60 150 300 528 100 350 2.987 556.106 0.508 0.159 0.151 0.146
    0.4, 0.4 60 150 300 451 100 350 5.079 496.388 0.700 0.145 0.142 0.142
    0.03 0.3, 0.4 60 150 300 528 100 350 3.146 605.635 0.330 0.159 0.159 0.158
    0.4, 0.4 60 150 300 451 100 350 4.978 560.046 0.535 0.146 0.143 0.144

     | Show Table
    DownLoad: CSV
    Table 8.  The computation costs on of BM, NM, and PDAS with .
    Error / Time/s
    BM NM PDAS BM NM PDAS BM NM PDAS BM NM PDAS
    0.02 0.3, 0.4 3000 2000 2000 2893 6000 800 1.261 1.304 1.336 42.699 19.135 1.174
    0.4, 0.4 3000 2000 2000 2396 6000 800 1.169 1.187 1.817 41.764 18.460 1.189
    0.03 0.3, 0.4 3000 2000 2000 2893 6000 800 1.217 1.416 1.072 43.031 19.568 1.164
    0.4, 0.4 3000 2000 2000 2396 6000 800 1.284 1.254 1.613 41.936 18.877 1.168

     | Show Table
    DownLoad: CSV
    Table 9.  The computation costs on of BM, NM and PDAS with .
    Error / Time/s
    BM NM PDAS BM NM PDAS BM NM PDAS BM NM PDAS
    0.02 0.3, 0.4 1500 2000 1200 2644 6000 800 1.169 1.542 1.302 26.770 21.504 0.758
    0.4, 0.4 2000 2000 1200 2608 6000 800 1.130 1.391 1.797 28.930 20.657 0.735
    0.03 0.3, 0.4 1500 2000 1000 2644 6000 800 1.220 1.705 1.137 26.402 21.231 0.669
    0.4, 0.4 2000 2000 1000 2608 6000 800 1.192 1.522 1.613 29.734 20.606 0.625

     | Show Table
    DownLoad: CSV

    The results about option prices listed in Tables 69 can fully illustrate the advantages of PDAS method in computational accuracy and speed. In order to get a more intuitive sense of option pricing, some images of the option price after the transformation model (2.3) and the original option price are shown in Figures 3 and 4, which are the main concerns for financial institutions. From the results in Figure 3, we confirm that the solutions of our proposed method approximate the benchmarks very well, and are all larger than the payoff , which further illustrate the efficiency and rationality of our proposed PDAS method.

    Figure 3.  The option price after the transformation model (2.3) with . (top left), (top right), (bottom left), (bottom right).
    Figure 4.  The option price in the model (2.1) with . (top left), (top right), (bottom left), (bottom right).

    The practicability of PDAS method has been checked, now we shall verify the error estimates obtained in Theorem 1. Taking the temporal direction as an example, suppose that the convergence order in the temporal direction is , then the logarithm of the error between the true and numerical solutions is linear with respect to the logarithm of , which is the degree of freedom in the temporal direction. Therefore, we can conclude that the slope of the logarithm of the error with respect to is equal to . According to the result in Theorem 1, we need to verify in the temporal and spatial direction, respectively. So as to reveal the relationship between the numbers of the spatial nodes and the spatial error, we choose in temporal direction, which can ensure that the spatial error is dominant. Likewise, to verify the relationship between the degrees of freedom in temporal direction and the temporal error, we choose . We also show blue line in each figure as the standard lines to verify the conclusion of theorem. The convergence results of options in spatial and temporal directions are shown in Figures 5 and 6.

    Figure 5.  The convergence rates in spatial direction with , , , (left), (right).
    Figure 6.  The convergence rates in temporal direction with , , , (left), (right).

    It needs to be noted that we have shifted some lines for clarity, but the slopes remain the same, which are the main concerns in our experiments. From Figures 5 and 6, we can conclude that the convergence rates in spatial and temporal directions coincide with error estimate in Theorem 1 are all of order 1.

    In this paper, we propose an efficient numerical method for the unilateral pricing problem of American better-of options with two underlying assets. By some traditional numeraire transformations, the far-field truncation technique and the known information about the free boundary, the original pricing model is approximated by a one-dimensional LCP on a bounded domain. In order to discretize the resulting LCP, the FDM and the FEM are applied in temporal direction and spatial direction, respectively. Based on the positive definite property of the discretized matrix, the discretized system is solved by the PDAS method. Our algorithm can obtain the price and the optimal exercise boundary of American better-of options simultaneously. Furthermore, we use some numerical experiments to verify the superiority of the proposed algorithm on error estimates and computational costs, respectively. In the future work, we shall apply our methods to different stochastic volatility models, regime switching models, and fractional models [32,33,34] to price American better-of options, as well as some inverse problems, such as the implied volatility of American options [35,36].

    The work of K. Zhang was supported by the NSF of China under the grant No. 11871245, and by the Key Laboratory of Symbolic Computation and Knowledge Engineering of Ministry of Education, Jilin University (93K172018Z01). The work of H. Song was supported by the NSF of China under the grant No.11701210, the education department project of Jilin Province under the grant No. JJKH20211031KJ, the NSF of Jilin Province under the grants No.20190103029JH, 20200201269JC and the fundamental research funds for the Central Universities.

    The authors declare there is no conflicts of interest.



    [1] X. J. He, W. Chen, A closed-form pricing formula for European options under a new stochastic volatility model with a stochastic long-term mean, Math. Financ. Econ., 15 (2021), 381–396. https://doi.org/10.1007/s11579-020-00281-y doi: 10.1007/s11579-020-00281-y
    [2] X. J. He, S. Lin, A fractional Black-Scholes model with stochastic volatility and European option pricing, Exp. Syst. Appl., 178 (2021), 114983. https://doi.org/10.1016/j.eswa.2021.114983 doi: 10.1016/j.eswa.2021.114983
    [3] S. L. Heston, A closed-form solution for options with stochastic volatility with applications to bond and currency options, Rev. Financ. Stud., 6 (1993), 327–343. https://doi.org/10.1093/rfs/6.2.327 doi: 10.1093/rfs/6.2.327
    [4] R. Caldana, G. Fusai, A general closed-form spread option pricing formula, J. Banking Finance, 37 (2013), 4893–4906. https://doi.org/10.1016/j.jbankfin.2013.08.016 doi: 10.1016/j.jbankfin.2013.08.016
    [5] R. Stulz, Options on the minimum or the maximum of two risky assets: analysis and applications, J. Financ. Econ., 10 (1982), 161–185. https://doi.org/10.1016/0304-405X(82)90011-3 doi: 10.1016/0304-405X(82)90011-3
    [6] L. Jiang, Mathematical modeling and methods of option pricing, World Scientific Publishing Company, 2005. https://doi.org/10.1142/5855
    [7] I. Karatzas, On the pricing of American options, Appl. Math. Optim., 17 (1988), 37–60. https://doi.org/10.1007/BF01448358 doi: 10.1007/BF01448358
    [8] K. S. Tan, K. R. Vetzal, Early exercise regions for exotic options, J. Deriv., 1 (1995), 42–56. https://doi.org/ssrn.com/abstract=6972
    [9] W. K. Wong, K. Xu, Refining the quadratic approximation formula for an American option, Int. J. Theor. Appl. Financ., 4 (2001), 773–781. https://doi.org/10.1142/S0219024901001243 doi: 10.1142/S0219024901001243
    [10] Y. Gao, H. Y. Liu, X. C. Wang, K. Zhang, On an artificial neural network for inverse scattering problems, J. Comput. Phys., 448 (2022), 110771. https://doi.org/10.1016/j.jcp.2021.110771 doi: 10.1016/j.jcp.2021.110771
    [11] M. Broadie, J. Detemple, American option valuation: new bounds, approximations and a comparison of existing methods, Rev. Financ. Stud., 9 (1996), 1211–1250. https://doi.org/10.1093/rfs/9.4.1211 doi: 10.1093/rfs/9.4.1211
    [12] H. Song, Q. Zhang, J. Li, H. Liu, Finite element method for valuation of American lookback options, Math. Numer. Sin., 38 (2016), 245–256. https://doi.org/10.12286/jssx.2016.3.245 doi: 10.12286/jssx.2016.3.245
    [13] J. C. Cox, S. A. Ross, M. Rubinstein, Option pricing: a simplified approach, J. Financ. Econ., 7 (1979), 229–263. https://doi.org/10.1016/0304-405X(79)90015-1 doi: 10.1016/0304-405X(79)90015-1
    [14] K. Amin, A. Khanna, Convergence of American option values from discrete-to continuous-time financial models, Math. Financ., 4 (1994), 289–304. https://doi.org/10.1111/j.1467-9965.1994.tb00059.x doi: 10.1111/j.1467-9965.1994.tb00059.x
    [15] J. A. Tilley, Valuing American options in a path simulation model, Transactions of the Society of Actuaries, 1993. https://doi.org/10.1.1.577.4823
    [16] A. D. Homes, H. Yang, A front-fixing finite element method for the valuation of American options, SIAM J. Sci. Comput., 30 (2008), 2158–2180. https://doi.org/10.1137/070694442 doi: 10.1137/070694442
    [17] Y. Gao, H. Song, X. Wang, K. Zhang, Primal-dual active set method for pricing American better-of option on two assets, Commun. Nonlinear Sci. Numer. Simul., 80 (2020), 104976. https://doi.org/10.1016/j.cnsns.2019.104976 doi: 10.1016/j.cnsns.2019.104976
    [18] X. Pang, H. Song, X. Wang, K. Zhang, An efficient numerical method for the valuation of American better-of options based on the front-fixing transform and the far field truncation, Adv. Appl. Math. Mech., 12 (2020), 902–919. https://doi.org/aamm.OA-2019-0107
    [19] H. Han, X. Wu, A fast numerical method for the Black-Scholes equation of American options, SIAM J. Numer. Anal., 41 (2003), 2081–2095. https://doi.org/10.1137/S0036142901390238 doi: 10.1137/S0036142901390238
    [20] M. Ehrhardt, R. E. Mickens, A fast, stable and accurate numerical method for the Black-Scholes equation of American options, Int. J. Theor. Appl. Financ., 11 (2008), 471–501. https://doi.org/10.1142/S0219024908004890 doi: 10.1142/S0219024908004890
    [21] R. Kangro, R. Nicolaides, Far field boundary conditions for Black-Scholes equations, SIAM J. Numer. Anal., 38 (2000), 1357–1368. https://doi.org/10.1137/S0036142999355921 doi: 10.1137/S0036142999355921
    [22] K. Ito, K. Kunisch, Augmented lagrangian methods for nonsmooth convex optimization in Hilbert spaces, Nonlinear Anal.: Theory, Methods Appl., 41 (2000), 591–616. https://doi.org/10.1016/S0362-546X(98)00299-5 doi: 10.1016/S0362-546X(98)00299-5
    [23] M. Hintermller, K. Ito, K. Kunisch, The primal-dual active set strategy as a semi-smooth newton method, SIAM J. Optim., 13 (2002), 865–888. https://doi.org/10.1137/S1052623401383558 doi: 10.1137/S1052623401383558
    [24] M. Bergounioux, K. Ito, K. Kunisch, Primal-dual strategy for constrained optimal control problem, SIAM J. Control Optim., 37 (1999), 1176–1194. https://doi.org/10.1137/S0363012997328609 doi: 10.1137/S0363012997328609
    [25] H. Song, X. Wang, K. Zhang, Q. Zhang, Primal-dual active set method for American lookback put option pricing, East Asian J. Appl. Math., 7 (2017), 603–614. https://doi.org/10.4208/eajam.060317.020617a doi: 10.4208/eajam.060317.020617a
    [26] R. Zhang, H. Song, N. Luan, Weak Galerkin finite element method for valuation of American options, Front. Math. China, 9 (2014), 455–476. https://doi.org/10.1007/s11464-014-0358-6 doi: 10.1007/s11464-014-0358-6
    [27] J. R. Cannon, The one-dimensional heat equation, Cambridge University Press, 1984. https://doi.org/10.1017/CBO9781139086967
    [28] H. Song, Q. Zhang, R. Zhang, A fast numerical method for the valuation of American lookback put options, Commun. Nonlinear Sci. Numer. Simul., 27 (2015), 302–313. https://doi.org/10.1016/j.cnsns.2015.03.010 doi: 10.1016/j.cnsns.2015.03.010
    [29] G. Strang, Approximation in the finite element method, Numer. Math., 19 (1972), 81–98. https://doi.org/10.1007/BF01395933 doi: 10.1007/BF01395933
    [30] C. Johnson, A convergence estimate for an approximation of parabolic variational inequalities, SIAM J. Numer. Anal., 13 (1976), 599–606. https://doi.org/10.1137/0713050 doi: 10.1137/0713050
    [31] A. Antuña, J. Guirao, M. Lópeaz, On the perturbations of maps obeying Shannon-Whittaker-Kotel'nikov's theorem generalization, Adv. Diff. Equations, 1 (2021), 1–12. https://doi.org/10.1186/s13662-021-03535-1 doi: 10.1186/s13662-021-03535-1
    [32] X. J. He, W. Chen, Pricing foreign exchange options under a hybrid Heston-Cox-Ingersoll-Ross model with regime switching, IMA J. Manage. Math., (2021), dpab013. https://doi.org/10.1093/imaman/dpab013
    [33] X. J. He, S. Lin, An analytical approximation formula for barrier option prices under the Heston model, Comput. Econ., (2021). https://doi.org/10.1007/s10614-021-10186-7
    [34] D. Ziane, M. H. Cherif, C. Catteni, K. Belghaba, Yang-Laplace decomposition method for nonlinear system of local fractional partial differential equations, Appl. Math. Nonlinear Sci., 4 (2019), 489–502. https://doi.org/10.2478/AMNS.2019.2.00046 doi: 10.2478/AMNS.2019.2.00046
    [35] Y. Deng, H. Liu, X. Wang, D. Wei, L. Zhu, Simultaneous recovery of surface heat flux and thickness of a solid structure by ultrasonic measurements, Electron. Res. Arch., 29 (2021), 3081–3096. https://doi.org/10.3934/era.2021027. doi: 10.3934/era.2021027
    [36] W. Yin, W. Yang, H. Liu, A neural network scheme for recovering scattering obstacles with limited phaseless far-field data, J. Comput. Phys., 417 (2020), 109594. https://doi.org/10.1016/j.jcp.2020.109594 doi: 10.1016/j.jcp.2020.109594
  • Reader Comments
  • © 2022 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(2339) PDF downloads(123) Cited by(0)

Figures and Tables

Figures(6)  /  Tables(9)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog