Pricing Asian Options Using Maple - Maple Programming Help

Home : Support : Online Help : Applications and Example Worksheets : Finance : Finance/Examples/AsianOptions

Pricing Asian Options Using Maple

 Overview This example demonstrates the use of Maple for computing the price of an Asian option, a derivative security that has gained popularity in financial markets in recent years. The payoff of an Asian option is based on the difference between an asset's average price over a given time period, and a fixed price called the strike price. Asian options are popular because they tend to have lower volatility than options whose payoffs are based purely on a single price point.  It is also harder for big traders to manipulate an average price over an extended period than a single price, so Asian options offer further protection against risk.   One disadvantage of Asian options is that their prices are very hard to compute using standard techniques.  Unlike European options, which can be priced using the classic Black-Scholes formula, there is no analytical formula for pricing an Asian option when the underlying asset is assumed to have a lognormal distribution, which is par for the course in financial modeling.   In this application, we use two different approaches for computing these prices numerically: (1) Solving a partial differential equation, and (2) Monte Carlo simulation.  We will see that the two numerical solutions that Maple derives are the same, providing strong validation for both the techniques and Maple's numerics.   Arithmetic average Asian options are securities whose payoff depends on the average of the underlying stock price over a certain period of time. To be more precise, the value of the continuous arithmetic Asian call option at time $t\le T$ is given by   ${C}_{t,T}\left(K\right)={ⅇ}^{-r\left(T-t\right)}E{\left[\frac{1}{T}\cdot {{\int }}_{0}^{T}S\left(u\right)\phantom{\rule[-0.0ex]{0.3em}{0.0ex}}{ⅆ}u-K\right]}^{\mathrm{+}}$   where $S\left(t\right)$ is the stock price at time $t$, $T$ is the expiration date, and $K$ is the strike price.   Since no general closed form solution for the price of the arithmetic average Asian option is known, a variety of numerical methods have been developed. These include formulation as a PDE, Monte Carlo simulation, and numeric inversion of the Laplace transform.   This worksheet demonstrates how these methods can be easily and quickly implemented in Maple.

Solution using PDE

 > $\mathrm{restart}$

It was shown by Vecer ("A new PDE approach for pricing arithmetic average Asian options", Journal of Computational Finance, vol. 4, No. 4, 105-113) that the value of the Asian option at time $0$ is given by

${S}_{0}u\left(0,\frac{{S}_{0}-K}{{S}_{0}}\right)$

where the function $u\left(t,z\right)$ satisfies the following PDE

${u}_{t}+r\left(1-\frac{t}{T}-z\right){u}_{z}+\frac{1{\left(1-\frac{t}{T}-z\right)}^{2}{\mathrm{\sigma }}^{2}{u}_{z,z}}{2}=0$

with the boundary condition

$u\left(T,z\right)={z}^{\mathrm{+}}$

First set up the equation:

 > $q:=t→1-\frac{t}{T}$
 ${q}{≔}{t}{↦}{1}{-}\frac{{t}}{{T}}$ (2.1)
 > $\mathrm{pde}:=\frac{\partial }{\partial t}u\left(t,z\right)+r\left(1-\frac{t}{T}-z\right)\left(\frac{\partial }{\partial z}u\left(t,z\right)\right)+\frac{1{\left(1-\frac{t}{T}-z\right)}^{2}{\mathrm{σ}}^{2}\left(\frac{{\partial }^{2}}{\partial {z}^{2}}u\left(t,z\right)\right)}{2}=0$
 ${\mathrm{pde}}{≔}\frac{{\partial }}{{\partial }{t}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{u}\left({t}{,}{z}\right){+}{r}{}\left({1}{-}\frac{{t}}{{T}}{-}{z}\right){}\left(\frac{{\partial }}{{\partial }{z}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{u}\left({t}{,}{z}\right)\right){+}\frac{{\left({1}{-}\frac{{t}}{{T}}{-}{z}\right)}^{{2}}{}{{\mathrm{\sigma }}}^{{2}}{}\left(\frac{{{\partial }}^{{2}}}{{\partial }{{z}}^{{2}}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{u}\left({t}{,}{z}\right)\right)}{{2}}{=}{0}$ (2.2)

State the boundary conditions:

 > ${\mathrm{ibc}}_{1}:=u\left(T,z\right)=\mathrm{max}\left(z,0\right);$
 ${{\mathrm{ibc}}}_{{1}}{≔}{u}\left({T}{,}{z}\right){=}{\mathrm{max}}\left({0}{,}{z}\right)$ (2.3)
 > ${\mathrm{ibc}}_{2}:=u\left(t,-1\right)=0;$
 ${{\mathrm{ibc}}}_{{2}}{≔}{u}\left({t}{,}{-1}\right){=}{0}$ (2.4)
 > ${\mathrm{ibc}}_{3}≔\mathrm{D}\left[2\right]\left(u\right)\left(t,1\right)=1$
 ${{\mathrm{ibc}}}_{{3}}{≔}{{\mathrm{D}}}_{{2}}\left({u}\right)\left({t}{,}{1}\right){=}{1}$ (2.5)

Assume a fixed time period of 1 year, a starting price of $100, a strike price of$100, and .

 > $T:=1;$
 ${T}{≔}{1}$ (2.6)
 > ${S}_{0}≔100.;$
 ${{S}}_{{0}}{≔}{100.}$ (2.7)
 > $K:=100.;$
 ${K}{≔}{100.}$ (2.8)
 > $\mathrm{σ}:=0.2;$
 ${\mathrm{\sigma }}{≔}{0.2}$ (2.9)
 > $r:=0.15;$
 ${r}{≔}{0.15}$ (2.10)

Now solve the PDE numerically using Maple's built-in solver.

 > $\mathrm{pds}:=\mathrm{pdsolve}\left(\mathrm{pde},\left\{{\mathrm{ibc}}_{1},{\mathrm{ibc}}_{2},{\mathrm{ibc}}_{3}\right\},\mathrm{numeric},\mathrm{spacestep}=0.001,\mathrm{timestep}=0.001\right);$
 ${\mathrm{pds}}{≔}{\mathbf{module}}\left({}\right)\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{...}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathbf{end module}}$ (2.11)
 > $\mathrm{val}:=\mathrm{pds}:-\mathrm{value}\left(t=0,\mathrm{output}=\mathrm{listprocedure}\right):$
 > $v:=\mathrm{eval}\left(u\left(t,z\right),\mathrm{val}\right):$

Finally, evaluate ${S}_{0}u\left(0,\frac{{S}_{0}-K}{{S}_{0}}\right)$ at Maple's numerical solution of $u\left(z,t\right)$.

 >
 ${8.40878796725886}$ (2.12)

Monte Carlo Simulation

Alternatively, you can use Monte Carlo methods to price the option.

 > $\mathrm{restart}:$
 > $\mathrm{with}\left(\mathrm{Finance}\right):$

Use the same parameters as before.

 > $T:=1;$
 ${T}{≔}{1}$ (3.1)
 > ${S}_{0}≔100.;$
 ${{S}}_{{0}}{≔}{100.}$ (3.2)
 > $K:=100.;$
 ${K}{≔}{100.}$ (3.3)
 > $\mathrm{σ}:=0.2;$
 ${\mathrm{\sigma }}{≔}{0.2}$ (3.4)
 > $r:=0.15;$
 ${r}{≔}{0.15}$ (3.5)

In the classical Black-Scholes setting the underlying process follows a geometric Brownian motion with drift $r$ and volatility $\mathrm{σ}.$

 > $X:=\mathrm{GeometricBrownianMotion}\left({S}_{0},r,\mathrm{σ}\right):$

Option payoff is given by the following formula:

 > $P:=t→\mathrm{max}\left(\frac{1{{∫}}_{0}^{t}X\left(u\right)\phantom{\rule[-0.0ex]{0.3em}{0.0ex}}{ⅆ}u}{t}-K,0\right)$
 ${P}{≔}{t}{↦}{\mathrm{max}}\left(\frac{{{\int }}_{{0}}^{{t}}{X}\left({u}\right)\phantom{\rule[-0.0ex]{0.3em}{0.0ex}}{ⅆ}{u}}{{t}}{-}{K}{,}{0}\right)$ (3.6)

You can compute the risk-neutral price of the option as expected payoff discounted using the risk-free rate.

 > $\mathrm{DiscountFactor}\left(r,T\right)\mathrm{ExpectedValue}\left(P\left(T\right),\mathrm{timesteps}=100,\mathrm{replications}={10}^{4},\mathrm{output}=\mathrm{value}\right)$
 ${8.342020802}$ (3.7)

You can also compute the various statistics related to this simulation.

 >
 value  =  9.569798629        standarderror  =  .1022532908    standarddeviation  =  10.22532908             skewness  =  1.172246041             kurtosis  =  1.209793258              minimum  =  0.              maximum  =  66.07060386        standarderror  =  .1022532908

Alternatively, you can generate sample data for the option payoff and perform statistical analysis of this data.

 > $A:=\mathrm{SampleValues}\left(P\left(T\right),'\mathrm{timesteps}'=100,'\mathrm{replications}'={10}^{4}\right);$
 ${{\mathrm{_rtable}}}_{{18446744074799623334}}$ (3.8)
 > eval(NextAfter);
 ${\mathbf{proc}}\left({}\right)\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathbf{option}}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathrm{builtin}}{=}{\mathrm{NextAfter}}{;}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathbf{end proc}}$ (3.9)
 >
 ${{\mathrm{_rtable}}}_{{18446744074799664774}}$ (3.10)
 >

You can consider a more complicated process for modeling the underlying asset. For example, you can incorporate jumps or regime switching into the model.

 > $\mathrm{X1}:=\mathrm{GeometricBrownianMotion}\left({S}_{0},r,\mathrm{σ}\right);$
 ${\mathrm{X1}}{≔}{\mathrm{_X4}}$ (3.11)
 > $\mathrm{X2}:=\mathrm{GeometricBrownianMotion}\left({S}_{0},0.1,0.5\right);$
 ${\mathrm{X2}}{≔}{\mathrm{_X5}}$ (3.12)
 > $P:=⟨⟨0.5,0.5⟩|⟨0.5,0.5⟩⟩;$
 $\left[\begin{array}{cc}0.5& 0.5\\ 0.5& 0.5\end{array}\right]$ (3.13)
 > $\mathrm{X3}:=\mathrm{RegimeSwitchingProcess}\left(P,⟨\mathrm{X1},\mathrm{X2}⟩,1,12\right);$
 ${\mathrm{X3}}{≔}{\mathrm{_X6}}$ (3.14)
 > $P:=\left(X,t\right)→\mathrm{max}\left(\frac{{{∫}}_{0}^{t}X\left(u\right)\phantom{\rule[-0.0ex]{0.3em}{0.0ex}}{ⅆ}u}{t}-K,0\right);$
 ${P}{≔}\left({X}{,}{t}\right){↦}{\mathrm{max}}\left(\frac{{{\int }}_{{0}}^{{t}}{X}\left({u}\right)\phantom{\rule[-0.0ex]{0.3em}{0.0ex}}{ⅆ}{u}}{{t}}{-}{K}{,}{0}\right)$ (3.15)
 > $\mathrm{DiscountFactor}\left(r,T\right)\mathrm{ExpectedValue}\left(P\left(\mathrm{X1},T\right),\mathrm{timesteps}=120,\mathrm{replications}={10}^{4},\mathrm{output}=\mathrm{value}\right);$
 ${8.497423429}$ (3.16)
 > $\mathrm{DiscountFactor}\left(r,T\right)\mathrm{ExpectedValue}\left(P\left(\mathrm{X2},T\right),\mathrm{timesteps}=120,\mathrm{replications}={10}^{4},\mathrm{output}=\mathrm{value}\right);$
 ${12.45400511}$ (3.17)
 > $\mathrm{DiscountFactor}\left(r,T\right)\mathrm{ExpectedValue}\left(P\left(\mathrm{X3},T\right),\mathrm{timesteps}=120,\mathrm{replications}={10}^{4},\mathrm{output}=\mathrm{value}\right);$
 ${10.85502698}$ (3.18)

Here is another example. You can consider another asset, which also follows a geometric Brownian motion but with different parameters. Consider the following contract: if the average price of $X$ is higher than the average price of $Y$, we get the difference, otherwise we get nothing.

 > $Y:=\mathrm{GeometricBrownianMotion}\left({S}_{0},0.1,0.5\right):$
 > $P:=t→\mathrm{max}\left(\frac{1{{∫}}_{0}^{t}X\left(u\right)\phantom{\rule[-0.0ex]{0.3em}{0.0ex}}{ⅆ}u}{t}-\frac{1{{∫}}_{0}^{t}Y\left(u\right)\phantom{\rule[-0.0ex]{0.3em}{0.0ex}}{ⅆ}u}{t},0\right):$
 > $\mathrm{ExpectedValue}\left(P\left(T\right),\mathrm{timesteps}=100,\mathrm{replications}={10}^{4}\right);$
 $\left[{\mathrm{value}}{=}{14.85155399}{,}{\mathrm{standarderror}}{=}{0.1795756877}\right]$ (3.19)

Again, in order to get the price of the contract, you have to discount this expected value using the risk-free rate.

 > $\mathrm{DiscountFactor}\left(r,T\right)\mathrm{ExpectedValue}\left(P\left(T\right),\mathrm{timesteps}=100,\mathrm{replications}={10}^{4},\mathrm{output}=\mathrm{value}\right);$
 ${12.83559553}$ (3.20)
 >