###### Asian Options Monte Carlo Pricing Using The Lévy Lognormal Approximation:

ABSTRACT— This paper investigates a number of popular methods for pricing Asian options: the Lévy Lognormal approximation and the Black-Scholes Model within a Monte Carlo implementation of arithmetic and geometric averaging methodologies. A multiple control variates technique is also implemented in the Monte Carlo Engine as a means of variance reduction in the price results. Several comparisons are drawn between the various techniques to benchmark the accuracy of each result and the speed of the corresponding algorithms. The Lévy Lognormal approximation provides a quick yet imprecise analytical framework for the geometric Asian option price, while the Monte Carlo implementations of the Black Scholes model with the appropriate variance reduction yields a remarkably accurate result at higher computation times. Multiple iterations of Monte Carlo engine input parameters are presented, and a time-efficient set of variables is benchmarked against popular case methods such as the Linetsky cases for arithmetic average Asian options.

KEYWORDSArithmetic Average, Asian Options, Brownian Motion, Control Variates, Derivatives Pricing, Geometric Average, Geometric Brownian Motion, Lévy Approximation, Lognormal Approximation, Monte Carlo Simulation, Options Pricing, Python, Statistical Modeling, Variance Reduction

I. INTRODUCTION

Asian Options are an exotic derivative where the payoff is determined by the average price of the underlying asset over its maturity T i.e. path dependent contract. The payoffs of a fixed-strike Asian option at maturity T are defined as, $$C\left(K,\ T\right)=\left(K-{\bar{A}}_T,0\right)^+$$ $$P\left(K,\ T\right)=\left({\bar{A}}_T-K,0\right)^+$$

$$\overline{{\bar{A}}_T=\frac{\int_{0}^{T}{S_tdt}}{T}\ i.e.\ \ average\ price\ of\ the\ underlying\ asset}$$
$$K=Strike\ price$$

In terms of time, the averaging can be continuous, or it can be taken in discrete time to assist with computational problems. In terms of degree, the average can be calculated arithmetically or geometrically, a detail that is explicitly stipulated in the option contract. They are defined as, $$A_N=\frac{1}{N}\sum_{i=1}^{N}S_{t_i}$$ $$G_N=(\prod_{i=1}^{N}{S_{t_i})}^{1/N}$$

$$\overline{S_{t_i}=Underlying\ asset\ price\ at\ times\ t_i\ for\ i=1,2,\ldots,N}$$

Asian options using arithmetic averaging procedure are more common than those using geometric averaging. However, the advantages of a geometric average procedure make it a very appealing choice for quantitative modeling and will be discussed further in the following sections.

Most Asian option contracts in industry are known to be averaging-in style. Averaging-in style contracts will sample at discrete and regular time intervals i.e weekly, monthly, quarterly, or annually averaging from inception to maturity date. Conversely, if an Asian option is averaging-out style the average is computed using specific samples near the maturity date. Generally speaking, averaging-out style options are more risky than averaging-in style options since the uncertainty of future spot prices is higher the further in time from the start of the contract. Averaging-out style Asian options have an increased risk exposure to spot price and volatility which is reflected in their higher price.

Arithmetic Asian option prices can only be estimated using numerical methods such as Monte Carlo, because the arithmetic average of lognormal variables do not follow a lognormal distribution. However, an analytical closed form formula for Geometric Asian Options can be adapted from the Black-Scholes-Merton formula, since geometric averages of log-normal variables also follow a lognormal distribution in a risk neutral world. Kemna and Vorst offer a good approximation for an geometric Asian option using the Black-Scholes formula which will be discussed in the pricing algorithm section [2].

Compared to similar ‘vanilla’ exchange traded derivatives Asian options are attractive as their volatility is low due to the averaging mechanics of its payout. Consequently, Asian options are less expensive than their corresponding vanilla derivatives. Furthermore, the higher number of observations, the lower the price of the Asian option. For this reason Asian options are a favored hedging tool of actors in volatile markets such as commodities, foreign exchange, and energy. Consider an airline that steadily buys crude oil whose supply price is not fixed, but set weekly from a particular benchmark. The airline may hedge themselves against a spike in oil by using a tailored Asian option to reflect the weekly purchases. This would be less expensive and more convenient than buying a basket of European options expiring at weekly intervals. Asian options are also be useful in markets with low volume, as the derivative would protect against price manipulation of the underlying asset.

II. MATERIALS AND METHODS

A. Geometric Brownian Motion
The underlying asset price is assumed to follow a stochastic process of a Geometric Brownian Motion (hereinafter GBM), $$S=rSdt+\sigma\ Sdz$$

$$\overline{S=price\ of\ the\ underlying\ asset} \\$$ $$r=risk-free\ interest\ rate \\$$ $$\sigma=volatility \\$$ $$dz=Wiener\ process \\$$

A Wiener process has the following properties:
1. The change $$\mathrm{\Delta\ z}$$ in a short period of time $$\mathrm{\Delta\ t}$$ is, $$\Delta z=\varepsilon\sqrt{\mathrm{\Delta t}}$$ $$\overline{\varepsilon~N\left(0,1\right)}$$
2. The values of $$\mathrm{\Delta\ z}$$ for any two different short intervals of time $$\mathrm{\Delta\ t}$$ are independent.

$$N\left(x,y\right)$$ denotes a normal distribution with mean x and variance y. In discrete time terms, the change in the stock price from Eq. 5 and 6 becomes,

$$\mathrm{\Delta\ S}=rS\mathrm{\Delta\ t}+\sigma\ S\varepsilon\sqrt{\mathrm{\Delta t}}$$ $$\frac{\mathrm{\Delta S}}{S}=r\mathrm{\Delta\ t}+\sigma\varepsilon\sqrt{\mathrm{\Delta t}}$$ $$\fbox{\frac{\mathrm{\Delta S}}{S}~N(r\mathrm{\Delta\ t},\sigma^2\mathrm{\Delta\ t})}$$
Therefore, the rate of return of the stock price over a time interval $$\Delta t$$ follows a normal distribution with mean $$r\Delta t$$, and standard deviation $$\sigma\sqrt{\mathrm{\Delta t}}$$. Applying It$$\hat{o}$$’s lemma, if a stochastic variable $$X$$ follows the It$$\hat{o}$$ process,

$$dX=a\left(X,t\right)dt+b\left(X,t\right)dz$$ A function $$f\left(X,t\right)$$ follows the process, $$df=\left(\frac{\partial f}{\partial X}a+\frac{\partial f}{\partial t}+\frac{1}{2}\frac{\partial^2f}{\partial X^2}b^2\right)dt+\frac{\partial f}{\partial X}bdz$$
This can in turn be applied to the stock price process in Eq. 5 to define the continuous path,

$$dln\left(S\right)=\left(\frac{1}{2}rS+0+\frac{1}{2}\left(\frac{-1}{S^2}\right)\sigma^2S^2\right)dt+\frac{1}{2}\sigma\ Sdz$$ $$dln\left(S\right)=\left(r-\frac{1}{2}\sigma^2\right)dt+\sigma\ dz$$
The equivalent of Eq. 10 in discrete time is therefore,

$$dln\left(S\right)\mathrm{\Delta}\ln{\left(S\right)}=\left(r-\frac{1}{2}\sigma^2\right)\mathrm{\Delta\ t}+\sigma\varepsilon\sqrt{\mathrm{\Delta t}}$$ $$ln\left(S_{t+\mathrm{\Delta t}}\right)-\ln{\left(S_t\right)}=\left(r-\frac{1}{2}\sigma^2\right)\mathrm{\Delta\ t}+\sigma\varepsilon\sqrt{\mathrm{\Delta t}}$$ $$\fbox{S_{t+\mathrm{\Delta t}}={S_t\ast e}^{\left(r-\frac{1}{2}\sigma^2\right)\mathrm{\Delta t}+\sigma\varepsilon\sqrt{\mathrm{\Delta t}}}}$$

Eq. 9 is the GBM price path generator that is used to construct hypothetical trajectories of the price of the underlying asset. The following pseudocode describes the constructor,

function GBM($$S_0$$,K,T,r,$$\sigma$$,$$n_{paths}$$,$$n_{steps}$$):
Define a timestep $$\Delta t$$ by dividing the maturity (T) input by the desired number of discrete steps (price values at time $$t_i$$ along the path, $$n_{steps}$$
Generate an array $$\epsilon$$ of size $$n_{paths}$$ rows and $$n_{steps}$$ columns, where each column represents one price path.
Fill the array will values drawn from a standard normal distribution i.e. $$\varepsilon\ N(0,1)$$.
Perform the operation described in Eq. 9 as a cumulative product of the exponents of the random sample with,
$$Mean=\left(r-\frac{1}{2}\sigma^2\right)\Delta t$$ $$Variance=\sigma\ast\varepsilon\sqrt{\Delta t}$$
The mean and variance of the price path are the risk-free interest rate and the stock volatility of the underlying asset respectively. Fig. 2 shows the result of $$n_{paths} = 1,000$$ price paths with $$n_{steps} = 10,000$$ timesteps.

B. Lognormal (Lévy) approximation and $$A_{N}$$:
The valuation of Arithmetic average Asian options under standard assumptions poses various difficulties. Central to the valuation issue is, with the exception of trivial cases, the lack of closed-form solutions under the conventional geometric diffusion model for the underlying price process. In essence, with no The Lévy approximation [1] solves the issue and enables an accurate approximation of the value of arithmetic average Asian options. More specifically, the hypothesis that “the distribution of an arithmetic average is well-approximated by the lognormal distribution when the underlying price process follows the conventional assumption of a geometric diffusion” [1], is confirmed with tests that simulate volatility similar to real market conditions. This fundamental assumption forms the basis of a performance benchmark of the pricing algorithm implementation against the standard Linetsky test cases [2].

The lognormal approximation arises from the moment matching procedure described in the Appendix wherein the mean $$\left(M_1\right)$$ and variance $$\left(M_2\right)$$ of $$\int_{0}^{T}{S_tdt}$$ in $$\mathcal{P}^\ast$$ is,

$$\mu_{lognormal}=\frac{\left(e^{rT}-1\right)}{rT}S_{0}$$ $$\sigma_{lognormal}=\frac{2S_0^{2}}{T^{2}\left(r+\sigma^2\right)}\left(\frac{e^{2r+σ^{2}Τ}-1}{2r+σ^{2}}-\frac{e^{rT}-1}{r}\right)-\frac{e^{rT}-1}{rTS_{0}}$$
The Monte Carlo Engine that replicates the BS model’s geometric Brownian motion trajectories in repeated succession draws from a log-normal distribution,

$$\ln{\left(X\right)}\sim\mathcal{N}\left(\mu_{lognormal},\sigma_{lognormal}^2\right)$$

and goes beyond the log-normal assumption. The closed-form solution of the price in terms of $$\sigma_{lognormal}$$ and $$\mu_{lognormal}$$ is,

$$A_C\left(K,T\right)=e^{-rT}\left[FN\left(d_1\right)-KN\left(d_2\right)\right] \\\\$$ $$A_P\left(K,T\right)=e^{-rT}\left[KN\left({-d}_2\right)-FN\left({-d}_1\right)\right] \\\\$$
Where,

$$d_{1,2}=\frac{\left(\ln{\frac{S}{K}}+\left(r\pm\frac{\sigma_{lognormal}^2}{2}\right)T\right)}{\sigma_{lognormal}\sqrt T} \\\\$$ $$F={S_0e}^{\mu_{lognormal}} \\\\$$

However, its accuracy is not sufficient and both methods are implemented and compared to each other to determine their strengths and limitations.

C. The modified Black-Scholes-Merton Model:
In the case of Geometric Average Asian Options, the Black-Scholes-Merton options pricing model (hereinafter referred to as BSM model) is fully compatible and can be solved in closed-form. This arises from the inherent lognormality of the two distributions. Therefore, the modified BSM PDE for geometric average call and put options is defined as shown below:

$$G_C\left(K,T\right)=e^{-rT}\ast[G_0\ast\ N\left(d_1\right)-K\ast\ N(d_2)] \\\\$$ $$G_P\left(K,T\right)=e^{-rT}\ast[K\ast\ N\left({-d}_2\right)-G_0\ast\ N\left({-d}_1\right)] \\\\$$ Where, $$d_{1,2}=\frac{1}{\mathrm{\Sigma}_G\sqrt T}(\log{\frac{G_0}{K}}\pm\frac{1}{2}\mathrm{\Sigma}_G^2\ast\ T) \\\\$$ $$\mathrm{\Sigma}_G=\frac{1}{\sqrt3}\sigma \\\\$$ $$G_0={S_0\ast e}^{\frac{1}{2}\left(r-q\right)T-\frac{1}{12}\sigma^2T} \\\\$$
The modified volatility $$\mathrm{\Sigma}_G$$ allows for the treatment of Asian options as European options with a decreased volatility by a factor of $$\sqrt3$$. Under the Lévy approximation [1], the accuracy of both arithmetic and geometric average Asian options can be measured through direct comparisons with the closed form BSM model. This confirms the shared lognormality of the geometric average and enables additional stress test against arithmetic average contracts, which are more commonly used in industry.

D. Monte Carlo Simulator
Monte Carlo experiments (hereinafter referred to as MC) are a class of computational algorithms that optimize a numerical result of a deterministic problem through the use of repeated random sampling. In summary, the following steps are executed in an MC simulation,
• Define the input domain for the task,
• Generate random inputs from a chosen probability distribution,
• Compute the numerical result for a number of iterations,
• Aggregate the results.
The GBM Generator implementation creates the aforementioned random sample and computes the price paths. The particular definition in this paper takes the GBM path array, averages the prices along each path using both an arithmetic and geometric mean, calculates the price of the call or put option and aggregates the option prices by taking the average value. An additional property of the definition allows for variance reduction through the use of the geometric average payoff as a control variate for the arithmetic average. The difference of the two averaging methods is then added back to the arithmetic average to generate narrower confidence intervals and better approximate the true mean of the option price. The functions are divided in separate definitions to optimize computational performance, but Fig. 4 presents the general definition for arithmetic and geometric average Asian Call and Put options. Each simulation invokes the GBM function for a given number of iterations and generates an average payoff over the number of simulations. The final result then becomes the average MC price estimate of the option.

E. Path-wise Greeks Derivations
Equivalent implementations of the Monte Carlo experiment above were constructed in Python’s Tensorflow framework. Tensorflow enables parallel computation for expensive tasks by sending packages to a Graphics Processing Unit (GPU) or a Tensor Processing Unit (TPU). Executing the algorithm within the Google Colaboratory Graphical User Interface allowed the use of Google’s GPU and TPU units for the calculation of the numerical result. This enabled the computation of path-wise Greeks for both arithmetic and geometric average call and put options. Fig. 5 shows the exact python definition of the algorithm, as a pseudocode interpretation of the Tensorflow definitions was extremely verbose. The algorithm generates average values of Greeks across all simulations and stores them in arrays. Table I Illustrates the Greeks generated in a convenient mnemonic rule.

F. Black-Scholes-Merton Greeks
The definition has been updates to include numerical and graphical representations of Greeks in order to better assess the sensitivities of the options contracts under evaluation. Table I shows the BSM Greeks capabilities in the algorithm (underlined). In addition, it includes Rho $$(\rho)$$, or the sensitivity of an option or an options portfolio to changes in the risk-free rate interest rate, as well as Phi $$(\phi)$$, or the expected change in the option premium due to small changes in the foreign currency interest rate, when applicable. All options calculations have been adjusted for the modified BSM model, and $$\mathrm{\Sigma}_G$$ (Eq. 17) in particular. Fig. 7 and 8 present the Greeks surfaces for a range of underlying price $$(S_{0})$$ and maturity $$(T)$$ values of a call and put options under the BSM model.

 Asset Price (S) Volatility (σ) Expiry (T) Option Price Delta (Δ) Vega (V) Theta (Θ) Delta (Δ) Gamma (Γ) Vanna (DΔDσ) Charm (DΔDt) Vega (V) Vanna (DVDS) Volga (DVDσ) Veta (DVDt) Gamma (Γ) Speed (DΓDS) Zomma (DΓDσ) Color (DΓDt) Vomma (DVDσ) N/A Ultima(DVomDσ) Totto (DVomDt)

III. RESULTS AND DISCUSSION

Asian options are rarely traded on an underlying stock of a publicly traded company within NYSE or NASDAQ making information on these options impossible to obtain from Yahoo Finance. However, these options are commonly used on commodities, most commonly on oil futures. This information is commonly available on CME (Chicago Mercantile Exchange). Due to the lack of a personal of corporate/institutional membership to their data service, using CME options data as a benchmark for the options proved to be difficult. However, within the scope of this research, the analytical results are computed and minimized through a series of cross-referencing and variance reduction techniques in order to estimate sufficiently accurate option prices without the need for data extraction with cumbersome data wrangling procedures. Fig. 17 shows the results of the lognormal MC engine. We see that as

A. Linetsky vs. Arithmetic MC
Linetsky [2] offers a set of test cases that has become the leading standard in the research literature with regards to calibration and accuracy test hypotheses. Using the eigenfunction expansion approach, the approximation yields results deemed to be the most accurate arithmetic average computation. The comparison of the Monte Carlo Arithmetic Average Options algorithm agrees with the test results within a range of absolute errors, $$[0.0147\%, 0.7345\%]$$ The error was minimized first, by increasing the number of price paths generated by the GBM function; second, by incrementally increasing the number of Monte Carlo simulations; third, by determining the optimal number of timesteps given the available timeframe. The numerical combination to minimize the errors that stem from insufficient sampling are shown in Fig. 6 and were identified through trial-and-error. The combination was used in all experiments, unless shown otherwise.

B. BSM Model vs. Geometric MC
In the case of the Geometric Average pricing algorithm, the results are compared with the BSM model numerical outputs and graphs of Greeks and graphs of the call option values across time to expiry (T), incremental values of the risk-free interest rate, volatility cases and others. Fig. 10 illustrates the BSM model results and Fig. 12 illustrates the equivalent cases from the MC Call Options under the combination of MC model parameters in Fig. 6. The results for three cases of call and put options with strikes $$K=95, 100, 105$$ are summarized in Table II. In accordance with Ruttiens [3] and Vorst [4], employing the solution to the geometric average problem leads to underpricing in put and/or overpricing in call options. However, the errors are within the range, $$[0.2047\%,0.5438\%]$$ which was deemed sufficient for the purposes of this benchmark.

 Gtype(K, T) BSM Monte Carlo % error (ε) Gc(95,1.0) 12.50853848101788 12.46260400000000 0.367228 Gc(100,1.0) 9.61215966961383 9.638977665500253 0.279001 Gc(105,1.0) 7.185865725921304 7.171154551742235 0.204724 Gp(95,1.0) 2.194652455717919 2.198018800000000 0.153388 Gp(100,1.0) 3.6018135264391495 3.610188200000000 0.232514 Gp(105,1.0) 5.479059464871914 5.489564400000000 0.191729

C. Variance Reduction
A variance reduction technique is employed to define the confidence intervals of the mean option price and to reduce the potential errors of averaging across each simulation and a The geometric average option payoff is used as a control variate to calibrate the accuracy of the arithmetic average component of the algorithm. This can be done by subtracting the geometric average option value from the equivalent arithmetically average counterpart, and then adding the geometrically averaged result back into the final option value. A second experiment was conducted using the difference between the option values of each averaging method and then tracking the confidence intervals at each step to obtain a descriptive set of the deviations between them. The differences become smaller along each timestep, further confirming the efficacy of the variance reduction technique in approximating the true mean of the option value. Using the same visualization script as the BSM results in III.B, the result of the control variates case of the MC method is evaluated and plotted against a range of market variables (Fig. 12).

D. Other MC method error estimates
If the probability density of estimator $$y$$ is $$f\left(y\right)$$, its expected value is, $$\mu_y\ \equiv\ E\left(y\right)=\int_{-\infty}^{\infty}{y\ f\left(y\right)dy}$$ For $$n_{sim}$$ trials, we can use the sample mean as an estimator for $$\mu_y$$, $$\bar{y}\equiv\frac{1}{N}\sum_{i=1}^{N}y_i$$ The moments of the error are, \begin{eqnarray} =E\left(\bar{y}-\mu_y\right)=E\left(\bar{y}\right)-\mu_y \\\\ =E\left(\frac{1}{N}\sum_{i=1}^{N}y_i\right)-\mu_y \\\\ =\frac{1}{N}\sum_{i=1}^{N}E\left(y_i\right)-\mu_y \\\\ \end{eqnarray} Due to the random sampling aspects of Monte Carlo methods, $$E\left(\bar{y}\right)=\mu_y$$. Therefore,

$$Var\left(\bar{y}-\mu_y\right) \\\\$$ $$Var\left(\bar{y}\right)-Var\left(\mu_y\right) \\\\$$ $$Var\left(\frac{\sum_{i=1}^{N}y_i}{N}\right) \\\\$$ $$\frac{1}{N^2}Var\left(\sum_{i=1}^{N}y_i\right) \\\\$$ $$\frac{1}{N^2}\ \sum_{i=1}^{N}Var\left(y_i\right) \\\\$$ $$\frac{1}{N^2}N\sigma_{y}^{2} \\\\$$ $$\fbox{\sigma_{\bar{y}}^2=E\left[\left(\bar{y}-\mu_y\right)^2\right]=\frac{\sigma_y^2}{N}} \\\\$$
Eq. 22 illustrates the interesting relationship between the standard error $$\sigma_{\bar{y}}$$ of the estimator and the sample size N i.e. the standard error decreases with $$\sqrt N$$. In terms of the Monte Carlo Engine, this finding aids in reducing unnecessary computation and setting bounds on the minimum number of paths needed to achieve optimality. To reduce the standard error in the experiment, the sample size was increased by a factor of 1,000 from the original sample. For the variance reduction implementation, the sample size was increased by a factor of 100 for satisfactory results, and to meet computational constraints.

E. Hedging Analysis
The analysis involves a portfolio of 5 Asian call options with stock underlying prices $$S_i$$ generated by the GBM price path generator. The following assumptions are made,
• Each path is refined to approximately 11 updates/second. The calculation assumes that price movements occur evenly throughout a 24-hour cycle.
• After-hours trading does not move the closing price of the security each day i.e. the price of the security at 9:30am (open market) is the same as the closing price the previous day.
• The security pays no dividend.
Each option represents a contract for one of 5 business days. Using the Greek letter derivatives property of the pricing algorithm, the static call option delta, $$\Delta\left(t_i\right)$$ is obtained for each day and is used to create a second portfolio with option values, $$C\left(t_{i}\right)-\Delta\left(t_{i}\right)S_{i}$$ The portfolio construction based on Eq. 23 is called a delta-hedged portfolio and exhibits the unique property of protecting the bundle against directional risk from volatility changes in the price of the underlying asset. The Asian call option delta is shown in Figure 3; it represents the rate of change of the option premium when the underlying asset price changes, and the number of shares required to maintain the overall traders’ position delta neutral.

The $$C\left(t_i\right)$$ plot of the two portfolios across time illustrates the additional protection of a delta-hedged portfolio (Fig. 13 and 14) as seen in the greatly reduced variability of the hedged portfolio throughout the week, as opposed to the naked portfolio. Asian options are easier to hedge than regular options as the payoff from an Asian option becomes more certain with the passage of time. As a result, the amount of uncertainty that needs to be hedged decreases with the passage of time. Delta-hedging is a very common options strategy and is very often paired with gamma hedging in a delta-gamma hedging strategy, which combines both delta and gamma hedges to protect against the risk of price changes in the underlying asset and the delta itself.

IV. NEXT STEPS
Time permitting the team would have liked to incorporate additional properties to the algorithm in order to more accurately reflect on the price dynamics of an Asian stock. One potential path forward is the addition of a volatility surface feature, which would generate more realistic results than those of a Lévy-approximation based model (where the volatility is static across the entire surface i.e. the plane, $$z = \sigma\sqrt{3}$$. Finally, we soon hope to test the pricing algorithm against real data and to take advantage of the unprecedented circumstances in the oil futures markets by incorporating data and features relevant to valuation in similar black-swan events.

Other Asian Options pricing algorithms developed over the years include, but are by no means limited to:

• Path Integral Approach - Effective Classical Potential
• Rogers and Shi’s PDE
• Variance Gamma Model
It would be a fascinating endeavor to incorporate more functionality in the pricing algorithm and be able to price options more accurately and efficiently. The hope of this paper is to generate additional questions and research on the topic and provide more efficient tools and methodologies to industry professionals.

 Case r σ S0 Linetsky Monte Carlo % ε Gc(K,T)1 0.0200 0.10 2.0 0.0559860415 0.0559778385 0.0147 Gc(K,T)2 0.1800 0.30 2.0 0.2183875466 0.2170470505 0.6138 Gc(K,T)3 0.0125 0.25 2.0 0.1722687410 0.1726129894 0.1998 Gc(K,T)4 0.0500 0.50 1.9 0.1931737903 0.1917548655 0.7345 Gc(K,T)5 0.0500 0.50 2.0 0.2464156905 0.2472542514 0.3403 Gc(K,T)6 0.0500 0.50 2.1 0.3062203648 0.3066361640 0.1357 Gc(K,T)7 0.0500 0.50 2.0 0.3500952190 0.3475691276 0.7215

ACKNOWLEDGMENTS
Professor Dan Pirjol for his invaluable guidance throughout this project; the Stevens Financial Engineering Department for providing remote Bloomberg Terminal access; and the Google Colaboratory and Google Compute service providers for enabling large-scale options pricing models at minimal costs for researchers.

AUTHORS

• Theo Dimitrasopoulos– MSc. Financial Engineering, Stevens Institute of Technology 2021; Bachelor of Science in Engineering, Structural Engineering, Princeton University 2017; tdimitr1(at)stevens(dot)edu
• William Kraemer– MSc. Financial Engineering, Stevens Institute of Technology 2020; BEng. Engineering Management, Stevens Institute of Technology 2019; wkraemer(at)stevens(dot)edu
• Snehal Rajguru– MEng. Engineering Management, Stevens Institute of Technology 2021; BEng. Electronics and Telecommunication Engineering, Pune University 2018; srajguru(at)stevens(dot)edu

REFERENCES

• [1] E. Lévy, (1992), “Pricing European average rate currency options”, Journal of International Money md Finance, 1992, vol. 11 pp. 474-491.
• [2] Kemna, A.G.Z., Vorst, A.C.F. (1990), "A pricing method for options based on average asset values". Journal of Banking and Finance, Elsevier, March 1990, vol. 14(1), pp 113-129.
• [3] V. Linetsky (2002), Spectral Decomposition of the Option Value. Published as "Exotic Spectra", RISK, April 2002, pp.85-89.
• [4] Laura Ballotta, Russell Gerrard, Ioannis Kyriakou (2017). Hedging of Asian options under exponential Lévy models: computation and performance, The European Journal of Finance, 23:4, 297-323.
• [5] Yor, M. (2001). Exponential Functionals of Brownian Motion and Related Processes. Springer-Verlag, New York.
• [6] Pirjol, D., Zhu, L. (2017). “Asymptotics for the discrete-time average of the geometric Brownian motion and Asian options”. Advances in Applied Probability, vol. 49(2), pp. 446-480.
• [7] Fu, Michael, Madan, Dilip, Wang, Tong. (1999). Pricing Continuous Asian Options: A Comparison of Monte Carlo and Laplace Transform Inversion Methods. J. Comp. Finance. 2. 10.21314/JCF.1998.024.

APPENDIX

SOURCE CODE


# -*- coding: utf-8 -*-
"""Pricer_v26_0_Release_Build

# **Valuation of Asian Options using the Lévy Approximation**

Course: FE 620: Pricing and Hedging | Stevens Institute of Technology

Group: Theo Dimitrasopoulos, Will Kraemer, Snehal Rajguru

*Version: v26.0*

**Instructions:**

* Click the execution button on the top left corner of each code cell to execute it. To run all cells in descending order, go to the menu bar at the top and click Runtime - Run all (or use the **Ctrl-F9** or **⌘-F9** hotkey for Windows and MacOSX respectively);
* If the code is running slowly, go to Runtime - Change runtime type, and change the Runtime Shape to High-RAM from the dropdown menu;
* The **"!pip install"** lines under the Python packages section (i.e. lines 3-5) only need to be executed the first time you run the notebook. If you receive the message **"Requirement already satisfied:"**, wrap them in treble quotes (add the quotes in lines 2 and 6)*.

## **Python Packages**
"""

#Installations#
'''
!pip install -q quandl
!pip install -q yfinance
!pip install quantsbin
!pip install quantumRandom
'''

# Main modules
import math
import numpy as np
import pandas as pd
'''
import quandl
import yfinance
import quantsbin
'''
# Scipy
from scipy import stats
from scipy.stats import norm
from scipy import special

# Plotters
import matplotlib as mpl
import mpl_toolkits.mplot3d.axes3d as p3
import matplotlib.pyplot as plt
from matplotlib import cm
from matplotlib.animation import FuncAnimation

# Helper modules
import timeit
import datetime
import random
import os
import requests
from datetime import timedelta
#import quantumrandom

# High-performance computing
'''
import tensorflow.compat.v1 as tf
tf.disable_v2_behavior()
'''
import jax
import jax.numpy as jnp
from jax.config import config
from jax.ops import index, index_add, index_update
import numba as nb
from numba import jit, njit, vectorize, cuda

# GPU Pointers for Numba
#os.environ['NUMBAPRO_LIBDEVICE'] = "/usr/local/cuda-10.1/nvvm/libdevice"
#os.environ['NUMBAPRO_NVVM'] = "/usr/local/cuda-10.1/nvvm/lib64/libnvvm.so"

# TPU Pointers for JAX
if 'TPU_DRIVER_MODE' not in globals():
url = 'http://' + os.environ['COLAB_TPU_ADDR'].split(':')[0] + ':8475/requestversion/tpu_driver_nightly'
resp = requests.post(url)
TPU_DRIVER_MODE = 1
config.FLAGS.jax_xla_backend = "tpu_driver"
print(config.FLAGS.jax_backend_target)

"""## **Variables**

### Model Variables
"""

# Initial Underlying Price
S0 = 100                                                #@param {type:"number" }
# Risk-free rate
r = 0.15                                                #@param {type:"number" }
# Dividend Yield Rate
q = 0.0                                                 #@param {type:"number" }
# Valuation Date
t = 0.0                                                 #@param {type:"number" }
# Maturity
T = 1.0                                                 #@param {type:"number" }
# Strike
K = 100                                                 #@param {type:"number" }
# Volatility
sigma = 0.3                                             #@param {type:"number" }
# Number of Price paths for each Monte Carlo simulation
n_paths = 100                                           #@param {type:"integer"}
# Time Steps
n_steps = 504                                           #@param {type:"integer"}
# Number of Monte Carlo simulations
n_sims = 1000                                           #@param {type:"integer"}
# Initial Seed
seed_0 = 42                                             #@param {type:"integer"}
# Initial Seed
seed = 100                                              #@param {type:"integer"}

"""### Plot Variables"""

# Universal Plot width
width = 25 #@param {type:"integer"}

# Universal Plot height
height =  14#@param {type:"integer"}

# Universal xtick size
xtick_size = 8 #@param {type:"integer"}

# Universal ytick size
ytick_size =  8#@param {type:"integer"}

# Universal title font size
title_size = 15 #@param {type:"integer"}

# Universal xlabel font size
xlabel_size = 12 #@param {type:"integer"}

# Universal xlabel font size
ylabel_size = 12 #@param {type:"integer"}

# Universal zlabel font size
zlabel_size = 12 #@param {type:"integer"}

# Universal legend font size
legend_size = 10 #@param {type:"integer"}

# Universal plot font color
color_plots = 'black'

plt.rcParams['figure.figsize'] = (width,height)
params = {'text.color' : 'black',
'xtick.color' : color_plots,
'ytick.color' : color_plots,
'xtick.labelsize' : xtick_size,
'ytick.labelsize' : ytick_size,
'legend.loc' : 'upper left',
}
plt.rcParams.update(params)

"""## **Definitions**

### Brownian Path Generator
"""

def bm_paths(n_paths,n_sims,seed):
if seed != 0:
np.random.seed(seed)
dt = T / n_steps
bt = np.random.randn(int(n_paths), int(n_sims))
W = np.cumprod(bt,axis=1)
return W

"""### Geometric Brownian Path Generator

#### Conventional Implementation:
"""

def gbm_paths_original(S0,K,T,t,r,q,sigma,seed,n_paths,n_steps):
np.random.seed(seed)
dt = T/n_steps
bt = np.random.randn(int(n_paths), int(n_steps))
S = S0 * np.cumprod((np.exp(sigma * np.sqrt(dt) * bt + (r - q - 0.5 * (sigma**2)) * dt)), axis = 1)
for i in range(0,len(S)):
S[i][0] = S0
return S

"""#### JAX Implementation"""

def gbm_paths_jax(S0, r, q, sigma, n_steps,n_paths):
dt = T / n_steps
seed = np.random.randint(0,1000000)
S = jnp.exp((r - q - sigma ** 2 / 2) * dt + np.sqrt(dt) * sigma * jax.random.normal(jax.random.PRNGKey(seed), (n_steps, n_paths)))
S = jnp.vstack([np.ones(n_paths), S])
S = S0 * jnp.cumprod(S, axis=0)
return S

def gbm_paths_jax_iter(S0, r, q, sigma, n_steps,n_paths):
seed = np.random.randint(0,1000000)
key = jax.random.PRNGKey(seed)
dt = T / n_steps
S = np.zeros((n_steps + 1, n_paths))
S[0] = S0
rn = jax.random.normal(key, shape=S.shape)
for t in range(1, n_steps + 1):
S[t] = S[t-1] * np.exp((r - q - sigma ** 2 / 2) * dt + sigma * math.sqrt(dt) * rn[t])
print(S[t])
return S

def gbm_logpaths_jax(S0, r, q, sigma, n_steps,n_paths):
dt = T / n_steps
seed = np.random.randint(0,1000000)
S = jnp.exp((r - q - sigma ** 2 / 2) * dt + np.sqrt(dt) * sigma * np.random.lognormal(mu_lognormal, var_lognormal, (n_steps, n_paths)))
S = jnp.vstack([np.ones(n_paths), S])
S = S0 * jnp.cumprod(S, axis=0)
return S

def gbm_logpaths_jax_iter(S0, r, q, sigma, n_steps,n_paths):
seed = np.random.randint(0,1000000)
key = jax.random.PRNGKey(seed)
dt = T / n_steps
S = np.zeros((n_steps + 1, n_paths))
S[0] = S0
rn = np.random.normal(mu_lognormal, var_lognormal, shape=S.shape)
for t in range(1, n_steps + 1):
S[t] = S[t-1] * np.exp((r - q - sigma ** 2 / 2) * dt + sigma * math.sqrt(dt) * rn[t])
print(S[t])
return S

"""#### Tensorflow Implementation:"""

def tf_graph_gbm_paths():
S0 = tf.placeholder(tf.float32)
K = tf.placeholder(tf.float32)
dt = tf.placeholder(tf.float32)
T = tf.placeholder(tf.float32)
sigma = tf.placeholder(tf.float32)
r = tf.placeholder(tf.float32)
dw = tf.placeholder(tf.float32)
S_i = S0 * tf.cumprod(tf.exp((r-sigma**2/2)*dt+sigma*tf.sqrt(dt)*dw), axis=1)
return (S0, K, dt, T, sigma, r, dw, S_i)

def tf_gbm_paths():
(S0,K, dt, T, sigma, r, dw, S_i) = tf_graph_gbm_paths()
def paths(S_zero, strk, maturity, riskfrate, volatility, seed, n_paths, n_steps):
if seed != 0:
np.random.seed(seed)
stdnorm_random_variates = np.random.randn(n_paths, n_steps)
with tf.Session() as sess:
delta_t = maturity / n_steps
res = sess.run(S_i,
{
S0: S_zero,
K : strk,
r : riskfrate,
sigma: volatility,
dt : delta_t,
T: maturity,
dw : stdnorm_random_variates})
return res
return paths

"""### **Geometric Black-Scholes-Merton**
Call and Put Asian option prices with geometric averaging. The alternative implementation in Tensorlfow offers efficiency and includes the calculations of 1st 2nd and 3rd order Greeks.

#### Conventional Implementation
"""

# Call Options:
def bsm_call(S0, K, T, t, r, q, sigma):
G0 = S0 * np.exp(0.5*(r - q - (sigma**2)/6) * (T - t)
Sigma_G = sigma/np.sqrt(3)
d1 = (1/(Sigma_G * np.sqrt(T - t))) * (np.log(G0/K) + 0.5 * (Sigma_G**2) * (T - t))
d2 = (1/(Sigma_G * np.sqrt(T - t))) * (np.log(G0/K) - 0.5 * (Sigma_G**2) * (T - t))
c = np.exp(-r * (T - t)) * (G0 * N(d1) - K * N(d2))
return c

# Put Options:
def bsm_put(S0, K, T, t, r, q, sigma):
G0 = S0 * np.exp(0.5*(r - q - (sigma**2)/6) * (T - t)
Sigma_G = sigma/np.sqrt(3)
d1 = (1/(Sigma_G * np.sqrt(T - t))) * (np.log(G0/K) + 0.5 * (Sigma_G**2) * (T - t))
d2 = (1/(Sigma_G * np.sqrt(T - t))) * (np.log(G0/K) - 0.5 * (Sigma_G**2) * (T - t))
p = np.exp(-r * (T - t)) * (K * N(-d2) - G0 * N(-d1))
return p

"""#### Tensorflow Implementation"""

# Sample Output
# [NET PRESENT VALUE, [DELTA, VEGA, THETA], [GAMMA, VANNA, CHARM], [VANNA, VOLGA, VETA], [SPEED, ZOMMA, COLOR], [N/A, ULTIMA, TOTTO]]

# Call Options:
def bsm_call_tf(enable_greeks = True):
S0 = tf.placeholder(tf.float32)
K = tf.placeholder(tf.float32)
dt = tf.placeholder(tf.float32)
sigma = tf.placeholder(tf.float32)
r = tf.placeholder(tf.float32)
q = tf.placeholder(tf.float32)
G0 = S0 * tf.exp(0.5 * (r * dt) - ((tf.np.square(sigma)) * dt)/12)
Sigma_G = sigma/tf.sqrt(3.0)
Phi = tf.distributions.Normal(0.,1.).cdf
d_1 = (1/(Sigma_G * tf.sqrt(dt))) * (tf.log(G0/K) + 0.5 * (tf.math.square(Sigma_G)) * dt)
d_2 = (1/(Sigma_G * tf.sqrt(dt))) * (tf.log(G0/K) - 0.5 * (tf.math.square(Sigma_G)) * dt)
npv =  tf.exp(-r * dt) * (G0 * Phi(d_1) - K * Phi(d_2))                # GREEKS TABLE:
target_calc = [npv]                                                    # (e.g. Option Price with respect to Asset Price (S) is delta)
if enable_greeks:                                                      #
grads_greeks = tf.ones([n_paths,1])                                  #                Asset Price (S)   Volatility    Time to Expiry
greeks = tf.gradients(npv, [S0, sigma, dt])                          # Option Price |     delta            vega           theta
dS_2nd = tf.gradients(greeks[0], [S0, sigma, dt])                    # Delta        |     gamma            vanna          charm
dsigma_2nd = tf.gradients(greeks[1], [S0, sigma, dt])                # Vega         |     vanna         vomma/volga       veta
dT_2nd = tf.gradients(dS_2nd[0], [S0, sigma, dt])                    # Gamma        |     speed            zomma          color
dsigma_3rd = tf.gradients(dsigma_2nd[1], [S0, sigma, dt])            # Vomma        |      N/A             ultima         totto
target_calc += [greeks, dS_2nd, dsigma_2nd, dT_2nd, dsigma_3rd]
def execute_graph(S_zero, strk, maturity, riskfrate, volatility):
with tf.Session() as sess:
res = sess.run(target_calc,
{
S0: S_zero,
K : strk,
r : riskfrate,
sigma: volatility,
dt: maturity})
return res
return execute_graph

# Put Options:
def bsm_put_tf(enable_greeks = True):
S0 = tf.placeholder(tf.float32)
K = tf.placeholder(tf.float32)
dt = tf.placeholder(tf.float32)
sigma = tf.placeholder(tf.float32)
r = tf.placeholder(tf.float32)
q = tf.placeholder(tf.float32)
G0 = S0 * tf.exp(0.5 * (r * dt) - ((tf.math.square(sigma)) * dt)/12)
Sigma_G = sigma/tf.sqrt(3.0)
Phi = tf.distributions.Normal(0.,1.).cdf
d_1 = (1/(Sigma_G * tf.sqrt(dt))) * (tf.log(G0/K) + 0.5 * (tf.math.square(Sigma_G)) * dt)
d_2 = (1/(Sigma_G * tf.sqrt(dt))) * (tf.log(G0/K) - 0.5 * (tf.math.square(Sigma_G)) * dt)
npv =  tf.exp(-r * dt) * (K * Phi(-d_2) - G0 * Phi(-d_1))              # GREEKS TABLE:
target_calc = [npv]                                                    # (e.g. Option Price with respect to Asset Price (S) is delta)
if enable_greeks:                                                      #
grads_greeks = tf.ones([n_paths,1])                                  #                Asset Price (S)   Volatility    Time to Expiry
greeks = tf.gradients(npv, [S0, sigma, dt])                          # Option Price |     delta            vega           theta
dS_2nd = tf.gradients(greeks[0], [S0, sigma, dt])                    # Delta        |     gamma            vanna          charm
dsigma_2nd = tf.gradients(greeks[1], [S0, sigma, dt])                # Vega         |     vanna         vomma/volga       veta
dT_2nd = tf.gradients(dS_2nd[0], [S0, sigma, dt])                    # Gamma        |     speed            zomma          color
dsigma_3rd = tf.gradients(dsigma_2nd[1], [S0, sigma, dt])            # Vomma        |      N/A             ultima         totto
target_calc += [greeks, dS_2nd, dsigma_2nd, dT_2nd, dsigma_3rd]
def execute_graph(S_zero, strk, maturity, riskfrate, volatility):
with tf.Session() as sess:
res = sess.run(target_calc,
{
S0: S_zero,
K : strk,
r : riskfrate,
sigma: volatility,
dt: maturity})
return res
return execute_graph

"""### **Monte Carlo Simulator with Arithmetic Average**
Call and Put Asian option prices with arithmetic averaging. The alternative implementation in Tensorlfow offers efficiency and includes the calculations of 1st 2nd and 3rd order Greeks.

#### Conventional Interpretation
"""

# Call Options:
def mc_call_arithm_np(S0, K, T, t, r, q, sigma, seed, n_paths, n_steps):
mc_call_arithm_payoffs = []
for i in range(1,n_paths):
seed = np.random.randint(1,2**31)
S = gbm_paths(S0,K,T,t,r,q,sigma,seed,n_paths,n_steps)
S_arithm_mu = np.mean(S)
arithm_payoff_call = np.exp(-r * T) * max(S_arithm_mu - K, 0)
mc_call_arithm_payoffs.append(arithm_payoff_call)
c = np.mean(mc_call_arithm_payoffs)
return c

# Put Options:
def mc_put_arithm_np(S0, K, T, t, r, q, sigma, seed, n_paths, n_steps):
mc_put_arithm_payoffs = []
for i in range(1,n_paths):
seed = np.random.randint(1,2**31)
S = gbm_paths(S0,K,T,t,r,q,sigma,seed,n_paths,n_steps)
S_arithm_mu = np.mean(S)
arithm_payoff_put = np.exp(-r * T) * max(K - S_arithm_mu, 0)
print(arithm_payoff_put)
mc_put_arithm_payoffs.append(arithm_payoff_put)
p = np.mean(mc_put_arithm_payoffs)
return p

"""#### JAX Implementation"""

def mc_call_arithm(S0, K, T, t, r, q, sigma, n_paths, n_steps, n_sims):
c = []
K = np.full(n_paths,K)
payoff_0 = np.zeros(n_paths)
for i in range(1,n_sims):
S = gbm_paths_jax(S0, r, q, sigma, n_steps, n_paths)
S_arithm = jnp.mean(S, axis=0)
c.append(jnp.mean(jnp.exp(-r*T) * jnp.maximum(S_arithm - K, payoff_0)))
call = np.mean(c)
return call

def mc_put_arithm(S0, K, T, t, r, q, sigma, n_paths, n_steps, n_sims):
p = []
K = np.full(n_paths,K)
payoff_0 = np.zeros(n_paths)
for i in range(1,n_sims):
S = gbm_paths_jax(S0, r, q, sigma, n_steps, n_paths)
S_arithm = jnp.mean(S, axis=0)
p.append(jnp.mean(jnp.exp(-r*T) * jnp.maximum(K - S_arithm, payoff_0)))
put = np.mean(p)
return put

def mc_call_lognormal(S0, K, T, t, r, q, sigma, n_paths, n_steps, n_sims):
c = []
K = np.full(n_paths,K)
payoff_0 = np.zeros(n_paths)
for i in range(1,n_sims):
S = gbm_paths_jax(S0, r, q, sigma, n_steps, n_paths)
S_arithm = jnp.mean(S, axis=0)
c.append(jnp.mean(jnp.exp(-r*T) * jnp.maximum(S_arithm - K, payoff_0)))
call = np.mean(c)
return call

def mc_put_lognormal(S0, K, T, t, r, q, sigma, n_paths, n_steps, n_sims):
p = []
K = np.full(n_paths,K)
payoff_0 = np.zeros(n_paths)
for i in range(1,n_sims):
S = gbm_paths_jax(S0, r, q, sigma, n_steps, n_paths)
S_arithm = jnp.mean(S, axis=0)
p.append(jnp.mean(jnp.exp(-r*T) * jnp.maximum(K - S_arithm, payoff_0)))
put = np.mean(p)
return put

"""#### Tensorflow Implementation:"""

# Sample Output
# [NET PRESENT VALUE, [DELTA, VEGA, THETA], [GAMMA, VANNA, CHARM], [VANNA, VOLGA, VETA], [SPEED, ZOMMA, COLOR], [N/A, ULTIMA, TOTTO]]

# Call Options:
def mc_call_arithm_tf(enable_greeks=False):
(S0, K, dt, T, sigma, r, dw, S_i) = tf_graph_gbm_paths()
A = tf.reduce_sum(S_i, axis=1)/(T/dt)
payout = tf.maximum(A - K, 0)
npv = tf.exp(-r * T) * tf.reduce_mean(payout)
target_calc = [npv]
if enable_greeks:
greeks = tf.gradients(npv, [S0, sigma, dt]) # delta, vega, theta
dS_2nd = tf.gradients(greeks[0], [S0, sigma, dt]) # gamma, vanna, charm
dsigma_2nd = tf.gradients(greeks[1], [S0, sigma, dt]) # vanna, vomma/volga, veta
dT_2nd = tf.gradients(dS_2nd[0], [S0, sigma, dt]) # speed, zomma, color
dsigma_3rd = tf.gradients(dsigma_2nd[1], [S0, sigma, dt]) # N/A, ultima, totto
target_calc += [greeks, dS_2nd, dsigma_2nd, dT_2nd, dsigma_3rd]
def pricer(S_zero, strk, maturity, volatility, riskfrate, seed, n_paths, n_steps):
if seed != 0:
np.random.seed(seed)
stdnorm_random_variates = np.random.randn(n_paths, n_steps)
with tf.Session() as sess:
delta_t = maturity / n_steps
res = sess.run(target_calc,
{
S0: S_zero,
K: strk,
r: riskfrate,
sigma: volatility,
dt: delta_t,
T: maturity,
dw: stdnorm_random_variates})
return res
return pricer

# Put Options:
def mc_put_arithm_tf(enable_greeks=False):
(S0, K, dt, T, sigma, r, dw, S_i) = tf_graph_gbm_paths()
A = tf.reduce_sum(S_i, axis=1)/(T/dt)
payout = tf.maximum(K - A, 0)
npv = tf.exp(-r * T) * tf.reduce_mean(payout)
target_calc = [npv]
if enable_greeks:
greeks = tf.gradients(npv, [S0, sigma, dt]) # delta, vega, theta
dS_2nd = tf.gradients(greeks[0], [S0, sigma, dt]) # gamma, vanna, charm
dsigma_2nd = tf.gradients(greeks[1], [S0, sigma, dt]) # vanna, vomma/volga, veta
dT_2nd = tf.gradients(dS_2nd[0], [S0, sigma, dt]) # speed, zomma, color
dsigma_3rd = tf.gradients(dsigma_2nd[1], [S0, sigma, dt]) # N/A, ultima, totto
target_calc += [greeks, dS_2nd, dsigma_2nd, dT_2nd, dsigma_3rd]
def pricer(S_zero, strk, maturity, volatility, riskfrate, seed, n_paths, n_steps):
if seed != 0:
np.random.seed(seed)
stdnorm_random_variates = np.random.randn(n_paths, n_steps)
with tf.Session() as sess:
delta_t = maturity / n_steps
res = sess.run(target_calc,
{
S0: S_zero,
K: strk,
r: riskfrate,
sigma: volatility,
dt: delta_t,
T: maturity,
dw: stdnorm_random_variates})
return res
return pricer

"""### **Monte Carlo Simulator with Geometric Average**

#### Conventional Implementation
"""

# Call Options:
def mc_call_geom_np(S0, K, T, t, r, q, sigma,n_paths,n_steps,n_sims):
mc_call_geom_payoffs = []
for i in range(1,n_sims):
S = gbm_paths_original(S0,K,T,t,r,q,sigma,seed,n_paths,n_steps)
S_geom_mu = np.exp(np.mean(np.log(S)))
mc_call_geom_payoffs.append(np.exp(-r * T) * max(S_geom_mu - K, 0))
c = np.mean(mc_call_geom_payoffs)
return c

# Put Options:
def mc_put_geom_np(S0, K, T, t, r, q, sigma,n_paths,n_steps,n_sims):
mc_put_geom_payoffs = []
for i in range(1,n_sims):
S = gbm_paths_original(S0,K,T,t,r,q,sigma,seed,n_paths,n_steps)
S_geom = np.exp(np.mean(np.log(S)))
mc_put_geom_payoffs.append(np.exp(-r * T) * max(K - S_geom, 0))
p = np.mean(mc_put_geom_payoffs)
return p

"""#### JAX Implementation"""

def mc_call_geom(S0, K, T, t, r, q, sigma, n_paths, n_steps, n_sims):
c = []
K = np.full(n_paths,K)
payoff_0 = np.zeros(n_paths)
for i in range(1,n_sims):
S = gbm_paths_jax(S0, r, q, sigma, n_steps, n_paths)
S_geom = jnp.exp(jnp.mean(jnp.log(S), axis=0))
c.append(jnp.mean(jnp.exp(-r*T) * jnp.maximum(S_geom - K, payoff_0)))
call = np.mean(c)
return call

def mc_put_geom(S0, K, T, t, r, q, sigma, n_paths, n_steps, n_sims):
p = []
K = np.full(n_paths,K)
payoff_0 = np.zeros(n_paths)
for i in range(1,n_sims):
S = gbm_paths_jax(S0, r, q, sigma, n_steps, n_paths)
S_geom = jnp.exp(jnp.mean(jnp.log(S), axis=0))
p.append(jnp.mean(jnp.exp(-r*T) * jnp.maximum(K - S_geom, payoff_0)))
put = np.mean(p)
return put

"""#### Tensorflow Implementation"""

# Sample Output
# [NET PRESENT VALUE, [DELTA, VEGA, THETA], [GAMMA, VANNA, CHARM], [VANNA, VOLGA, VETA], [SPEED, ZOMMA, COLOR], [N/A, ULTIMA, TOTTO]]

# Call Options:
def mc_call_geom_tf(enable_greeks=True):
(S0, K, dt, T, sigma, r, dw, S_i) = tf_graph_gbm_paths()
A = tf.pow(tf.reduce_prod(S_i, axis=1), dt / T)
payout = tf.maximum(A - K, 0)
npv = tf.exp(-r * T) * tf.reduce_mean(payout)                          # GREEKS TABLE:
target_calc = [npv]                                                    # (e.g. Option Price with respect to Asset Price (S) is delta)
if enable_greeks:                                                      #
grads_greeks = tf.ones([n_paths,1])                                  #                Asset Price (S)   Volatility    Time to Expiry
greeks = tf.gradients(npv, [S0, sigma, dt])                          # Option Price |     delta            vega           theta
dS_2nd = tf.gradients(greeks[0], [S0, sigma, dt])                    # Delta        |     gamma            vanna          charm
dsigma_2nd = tf.gradients(greeks[1], [S0, sigma, dt])                # Vega         |     vanna         vomma/volga       veta
dT_2nd = tf.gradients(dS_2nd[0], [S0, sigma, dt])                    # Gamma        |     speed            zomma          color
dsigma_3rd = tf.gradients(dsigma_2nd[1], [S0, sigma, dt])            # Vomma        |      N/A             ultima         totto
target_calc += [greeks, dS_2nd, dsigma_2nd, dT_2nd, dsigma_3rd]
def pricer(S_zero, strk, maturity, riskfrate, volatility, seed, n_paths, n_steps):
if seed != 0:
np.random.seed(seed)
stdnorm_random_variates = np.random.randn(n_paths, n_steps)
with tf.Session() as sess:
delta_t = maturity / n_steps
res = sess.run(target_calc,
{
S0: S_zero,
K: strk,
r: riskfrate,
sigma: volatility,
dt: delta_t,
T: maturity,
dw: stdnorm_random_variates})
return res
return pricer

# Put Options:
def mc_put_geom_tf(enable_greeks=True):
(S0, K, dt, T, sigma, r, dw, S_i) = tf_graph_gbm_paths()
A = tf.pow(tf.reduce_prod(S_i, axis=1), dt / T)
payout = tf.maximum(K - A, 0)
npv = tf.exp(-r * T) * tf.reduce_mean(payout)                          # GREEKS TABLE:
target_calc = [npv]                                                    # (e.g. Option Price with respect to Asset Price (S) is delta)
if enable_greeks:                                                      #
grads_greeks = tf.ones([n_paths,1])                                  #                Asset Price (S)   Volatility    Time to Expiry
greeks = tf.gradients(npv, [S0, sigma, dt])                          # Option Price |     delta            vega           theta
dS_2nd = tf.gradients(greeks[0], [S0, sigma, dt])                    # Delta        |     gamma            vanna          charm
dsigma_2nd = tf.gradients(greeks[1], [S0, sigma, dt])                # Vega         |     vanna         vomma/volga       veta
dT_2nd = tf.gradients(dS_2nd[0], [S0, sigma, dt])                    # Gamma        |     speed            zomma          color
dsigma_3rd = tf.gradients(dsigma_2nd[1], [S0, sigma, dt])            # Vomma        |      N/A             ultima         totto
target_calc += [greeks, dS_2nd, dsigma_2nd, dT_2nd, dsigma_3rd]
def pricer(S_zero, strk, maturity, riskfrate, volatility, seed, n_paths, n_steps):
if seed != 0:
np.random.seed(seed)
stdnorm_random_variates = np.random.randn(n_paths, n_steps)
with tf.Session() as sess:
delta_t = maturity / n_steps
res = sess.run(target_calc,
{
S0: S_zero,
K: strk,
r: riskfrate,
sigma: volatility,
dt: delta_t,
T: maturity,
dw: stdnorm_random_variates})
return res
return pricer

"""### Lognormal Monte Carlo

#### Conventional Implementation
"""

def Monte_Carlo_lognormal(iterations = 1000000, K = 110, T = 1, r = 0.1):
Price = np.exp(-r*T) * np.maximum(lognorm - K, 0)
lognorm = np.random.lognormal(mu_lognormal, sigma_lognormal, (n_paths)
std = np.std(Price)
mean = np.mean(Price)
return mean+(1.96*std/np.sqrt(iterations)), mean, mean-(1.96*std/np.sqrt(iterations))

"""#### JAX Implementation"""

def mc_call_log(S0, K, T, t, r, q, sigma, n_paths, n_steps, n_sims):
c = []
K = np.full(n_paths,K)
payoff_0 = np.zeros(n_paths)
mean = []
std = []
for i in range(1,n_sims):
S = gbm_paths_jax(S0, r, q, sigma, n_steps, n_paths)
S_geom = jnp.exp(jnp.mean(jnp.log(S), axis=0))
c.append(jnp.mean(jnp.exp(-r*T) * jnp.maximum(S_geom - K, payoff_0)))
std.append(np.std(c))
mean.append(np.mean(c))
call = np.mean(c)
return call, mean+(1.96*std/np.sqrt(n_sims)), mean, mean-(1.96*std/np.sqrt(n_sims))

"""### PDF, CDF"""

# Probability density function of standard normal
def dN(x):
return np.exp(-0.5 * x ** 2) / np.sqrt(2 * np.pi)

# Cumulative density function of standard normal
def N(u):
q = special.erf(u / np.sqrt(2.0))
return (1.0 + q) / 2.0

"""### **Greeks**

#### Greeks Derivations
"""

# 1st Order Greeks:

# Delta
def bsm_delta(S0, K, T, t, r, q, sigma, optiontype):
G0 = S0 * np.exp(0.5 * (T - t) * (r - q - (sigma**2)/6))
Sigma_G = sigma/np.sqrt(3)
d1 = (1/(Sigma_G * np.sqrt(T - t))) * (np.log(G0/K) + 0.5 * (Sigma_G**2) * (T - t))
if(optiontype == "Call"):
delta = np.exp(-(T - t)) * N(d1)
elif(optiontype == "Put"):
delta = -np.exp(-(T - t)) * N(-d1)
return delta

# Vega
def bsm_vega(S0, K, T, t, r, q, sigma, optiontype):
G0 = S0 * np.exp(0.5 * (T - t) * (r - q - (sigma**2)/6))
Sigma_G = sigma/np.sqrt(3)
d1 = (1/(Sigma_G * np.sqrt(T - t))) * (np.log(G0/K) + 0.5 * (Sigma_G**2) * (T - t))
vega = S0 * dN(d1) * np.sqrt(T - t)
return vega

# Rho
def bsm_rho(S0, K, T, t, r, q, sigma, optiontype):
G0 = S0 * np.exp(0.5 * (T - t) * (r - q - (sigma**2)/6))
Sigma_G = sigma/np.sqrt(3)
d1 = (1/(Sigma_G * np.sqrt(T - t))) * (np.log(G0/K) + 0.5 * (Sigma_G**2) * (T - t))
d2 = (1/(Sigma_G * np.sqrt(T - t))) * (np.log(G0/K) - 0.5 * (Sigma_G**2) * (T - t))
if(optiontype == "Call"):
rho = K * (T - t) * np.exp(-r * (T - t)) * N(d2)
if(optiontype == "Put"):
rho = -K * (T - t) * np.exp(-r * (T - t)) * N(-d2)
return rho

# Theta
def bsm_theta(S0, K, T, t, r, q, sigma, optiontype):
G0 = S0 * np.exp(0.5 * (T - t) * (r - q - (sigma**2)/6))
Sigma_G = sigma/np.sqrt(3)
d1 = (1/(Sigma_G * np.sqrt(T - t))) * (np.log(G0/K) + 0.5 * (Sigma_G**2) * (T - t))
d2 = (1/(Sigma_G * np.sqrt(T - t))) * (np.log(G0/K) - 0.5 * (Sigma_G**2) * (T - t))
if(optiontype == "Call"):
theta = -(S0 * dN(d1) * Sigma_G / (2 * np.sqrt(T - t)) - r * K * np.exp(-r * (T - t)) * N(d2))
if(optiontype == "Put"):
theta = -(S0 * dN(d1) * Sigma_G / (2 * np.sqrt(T - t)) + r * K * np.exp(-r * (T - t)) * N(-d2))
return theta

# 2nd Order Greeks:

# Gamma
def bsm_gamma(S0, K, T, t, r, q, sigma, optiontype):
G0 = S0 * np.exp(0.5 * (T - t) * (r - q - (sigma**2)/6))
Sigma_G = sigma/np.sqrt(3)
d1 = (1/(Sigma_G * np.sqrt(T - t))) * (np.log(G0/K) + 0.5 * (Sigma_G**2) * (T - t))
gamma = dN(d1) / (S0 * Sigma_G * math.sqrt(T - t))
return gamma

# Charm
def bsm_charm(S0, K, T, t, r, q, sigma, optiontype):
G0 = S0 * np.exp(0.5 * (T - t) * (r - q - (sigma**2)/6))
Sigma_G = sigma/np.sqrt(3)
d1 = (1/(Sigma_G * np.sqrt(T - t))) * (np.log(G0/K) + 0.5 * (Sigma_G**2) * (T - t))
d2 = (1/(Sigma_G * np.sqrt(T - t))) * (np.log(G0/K) - 0.5 * (Sigma_G**2) * (T - t))
charm = -dN(d1) * (2 * r * (T-t) - d2 * Sigma_G * np.sqrt(T-t)) / (2 * (T - t) * Sigma_G * np.sqrt(T-t))
return charm

# Phi
def bsm_phi(S0, K, T, t, r, q, sigma, optiontype):
G0 = S0 * np.exp(0.5 * (T - t) * (r - q - (sigma**2)/6))
Sigma_G = sigma/np.sqrt(3)
d1 = (1/(Sigma_G * np.sqrt(T - t))) * (np.log(G0/K) + 0.5 * (Sigma_G**2) * (T - t))
d2 = (1/(Sigma_G * np.sqrt(T - t))) * (np.log(G0/K) - 0.5 * (Sigma_G**2) * (T - t))
phi = 0.01 * T * S0 * np.exp(-q * T) * N(d1)
return phi

# Vanna
def bsm_vanna(S0, K, T, t, r, q, sigma, optiontype):
G0 = S0 * np.exp(0.5 * (T - t) * (r - q - (sigma**2)/6))
Sigma_G = sigma/np.sqrt(3)
d1 = (1/(Sigma_G * np.sqrt(T - t))) * (np.log(G0/K) + 0.5 * (Sigma_G**2) * (T - t))
d2 = (1/(Sigma_G * np.sqrt(T - t))) * (np.log(G0/K) - 0.5 * (Sigma_G**2) * (T - t))
vanna = S0 * np.exp(-q * T) * d2 / Sigma_G * dN(d1)
return vanna

# Vomma
def bsm_vomma(S0, K, T, t, r, q, sigma, optiontype):
G0 = S0 * np.exp(0.5 * (T - t) * (r - q - (sigma**2)/6))
Sigma_G = sigma/np.sqrt(3)
d1 = (1/(Sigma_G * np.sqrt(T - t))) * (np.log(G0/K) + 0.5 * (Sigma_G**2) * (T - t))
d2 = (1/(Sigma_G * np.sqrt(T - t))) * (np.log(G0/K) - 0.5 * (Sigma_G**2) * (T - t))
vomma = -0.01 * np.exp(-q * T) * d2 / Sigma_G * dN(d1)
return vomma

"""#### Greeks: Sensitivity Plotter"""

def greeks_plot_tool(greek_function, x_var_name, y_var_name,S0, K, T, t, r, q, sigma, x, y, optiontype, greek, plot):
# Initialise vector to store our option values and then iterate over
# Assumption that we're using a constant sized vector length for each variable
# Need to change the variables being iterated over here for each update (possibly a better way to do this)
V = np.zeros((len(S0), len(S0)), dtype=np.float)
for i in range(len(S0)):
for j in range(len(S0)):
V[i, j] = greek_function(S0[i], K[i], T[j], t[i], r[i], q[i], sigma[i], optiontype)

# Initiliase plotting canvas
surf = plot.plot_surface(x, y, V, rstride=1, cstride=1, alpha=0.7, cmap=cm.RdYlGn, antialiased = True) #cmap = cm.RdYlBu, for surface plotting
plot.set_xlabel(x_var_name, color = color_plots, fontsize = xlabel_size)
plot.set_ylabel(y_var_name, color = color_plots, fontsize = ylabel_size)
plot.set_zlabel("%s(K, T)" % greek, color = color_plots, fontsize = zlabel_size)
plot.set_title("%s %s" % (optiontype, greek), color = color_plots, fontsize=title_size, fontweight='bold')

# Calculate colour levels based on our meshgrid
Vlevels = np.linspace(V.min(), V.max(), num=8, endpoint=True)
xlevels = np.linspace(x.min(), x.max(), num=8, endpoint=True)
ylevels = np.linspace(y.min(), y.max(), num=8, endpoint=True)

cset = plot.contourf(x, y, V, Vlevels, zdir='z',offset=V.min(), cmap=cm.RdYlGn, alpha = 0.5, antialiased = True)
cset = plot.contourf(x, y, V, xlevels, zdir='x',offset=x.min(), cmap=cm.RdYlGn, alpha = 0.5, antialiased = True)
cset = plot.contourf(x, y, V, ylevels, zdir='y',offset=y.max(), cmap=cm.RdYlGn, alpha = 0.5, antialiased = True)

# Set our viewing constraints
plt.clabel(cset,fontsize=10, inline=1, color = color_plots, antialiased = True)
plot.set_xlim(x.min(), x.max())
plot.set_ylim(y.min(), y.max())
plot.set_zlim(V.min(), V.max())

# Colorbar
colbar = plt.colorbar(surf, shrink=1.0, extend='both', aspect = 10)
l,b,w,h = plt.gca().get_position().bounds
ll,bb,ww,hh = colbar.ax.get_position().bounds
colbar.ax.set_position([ll, b+0.1*h, ww, h*0.8])

"""### Option Sensitivity

#### BSM Option Sensitivity
"""

def bsm_plot_values(function,S0, K, T, t, r, q, sigma,optiontype):
fig = plt.figure(figsize=(30,90))
points = 100

# Option(K,T) vs. Strike
klist = np.linspace(K-30, K+30, points)
vlist = [function(S0, K, T, t, r, q, sigma) for K in klist]
fig1.plot(klist, vlist)
fig1.grid()
fig1.set_title('BSM %s Option Value vs. Strike' % optiontype, color = color_plots, fontsize = title_size)
fig1.set_xlabel('Strike $K$', color = color_plots, fontsize = xlabel_size)
fig1.set_ylabel('%s Option Value' % optiontype, color = color_plots, fontsize = ylabel_size)

# Option(K,T) vs. Strike vs. Underlying Price
klist = np.linspace(K-30, K+30, points)
s0list = np.linspace(S0 - 20, S0 + 20, points)
V = np.zeros((len(s0list), len(klist)), dtype=np.float)
for j in range(len(klist)):
for i in range(len(s0list)):
V[i, j] = function(s0list[i], klist[j], T, t, r, q, sigma)

x, y = np.meshgrid(klist, s0list)
fig2.patch.set_alpha(0.0)
fig2.plot_wireframe(x, y, V, linewidth=1.0, color = color_plots) #cmap = cm.RdYlGn, for surface plotting
fig2.set_title('BSM %s Option Value vs. Strike vs. Underlying Price' % optiontype, color = color_plots, fontsize = title_size)
fig2.set_xlabel('Strike $K$', color = color_plots, fontsize = xlabel_size)
fig2.set_ylabel('Stock/Underlying Price ($)', color = color_plots, fontsize = ylabel_size) fig2.set_zlabel('%s Option Value' % optiontype, color = color_plots, fontsize = zlabel_size) # Option(K,T) vs. Time fig3 = fig.add_subplot(822) tlist = np.linspace(0.0001, T, points) vlist = [function(S0, K, T, t, r, q, sigma) for T in tlist] fig3.plot(tlist, vlist) fig3.grid() fig3.set_title('BSM %s Option Value vs. Time' % optiontype, color = color_plots, fontsize = title_size) fig3.set_xlabel('Maturity$T$', color = color_plots, fontsize = xlabel_size) fig3.set_ylabel('%s Option Value' % optiontype, color = color_plots, fontsize = ylabel_size) # Option(K,T) vs. Time vs. Underlying Price tlist = np.linspace(0.0001, T, points) s0list = np.linspace(S0 - 20, S0 + 10, points) V = np.zeros((len(s0list), len(tlist)), dtype=np.float) for j in range(len(tlist)): for i in range(len(s0list)): V[i, j] = function(s0list[i], K, tlist[j], t, r, q, sigma) fig4 = fig.add_subplot(824, projection="3d") x, y = np.meshgrid(tlist, s0list) fig4.patch.set_alpha(0.0) fig4.plot_wireframe(x, y, V, linewidth=1.0, color = color_plots) #cmap = cm.RdYlGn, for surface plotting fig4.set_title('BSM %s Option Value vs. Time vs. Underlying Price' % optiontype, color = color_plots, fontsize = title_size) fig4.set_xlabel('Maturity$T$', color = color_plots, fontsize = xlabel_size) fig4.set_ylabel('Stock/Underlying Price ($)', color = color_plots, fontsize = ylabel_size)
fig4.set_zlabel('%s Option Value' % optiontype, color = color_plots, fontsize = zlabel_size)

# Option(K,T) vs. r
rlist = np.linspace(0, r, points)
vlist = [function(S0, K, T, t, r, q, sigma) for r in rlist]
fig5.plot(rlist, vlist)
fig5.grid()
fig5.set_title('BSM %s Option Value vs. r' % optiontype, color = color_plots, fontsize = title_size)
fig5.set_xlabel('Risk-free rate $r$', color = color_plots, fontsize = xlabel_size)
fig5.set_ylabel('%s Option Value' % optiontype, color = color_plots, fontsize = ylabel_size)

# Option(K,T) vs. r vs. Underlying Price
rlist = np.linspace(0, r, points)
s0list = np.linspace(S0 - 20, S0 + 20, points)
V = np.zeros((len(s0list), len(rlist)), dtype=np.float)
for j in range(len(rlist)):
for i in range(len(s0list)):
V[i, j] = function(s0list[i], K, T, t, rlist[j], q, sigma)

x, y = np.meshgrid(rlist, s0list)
fig6.patch.set_alpha(0.0)
fig6.plot_wireframe(x, y, V, linewidth=1.0, color = color_plots) #cmap = cm.RdYlGn, for surface plotting
fig6.set_title('BSM %s Option Value vs. r vs. Underlying Price' % optiontype, color = color_plots, fontsize = title_size)
fig6.set_xlabel('Risk-free rate $r$', color = color_plots, fontsize = xlabel_size)
fig6.set_ylabel('Stock/Underlying Price ($)', color = color_plots, fontsize = ylabel_size) fig6.set_zlabel('%s Option Value' % optiontype, color = color_plots, fontsize = zlabel_size) # Option(K,T) vs. Implied Vol. fig7 = fig.add_subplot(826) slist = np.linspace(0.01, sigma, points) vlist = [function(S0, K, T, t, r, q, sigma) for sigma in slist] fig7.plot(slist, vlist) fig7.grid() fig7.set_title('BSM %s Option Value vs. Volatility' % optiontype, color = color_plots, fontsize = title_size) fig7.set_xlabel('Volatility$\sigma$', color = color_plots, fontsize = xlabel_size) fig7.set_ylabel('%s Option Value' % optiontype, color = color_plots, fontsize = ylabel_size) # Option(K,T) vs. Volatility vs. Underlying Price slist = np.linspace(0.01, sigma, points) s0list = np.linspace(S0 - 20, S0 + 20, points) V = np.zeros((len(s0list), len(slist)), dtype=np.float) for j in range(len(slist)): for i in range(len(s0list)): V[i, j] = function(s0list[i], K, T, t, r, q, slist[j]) fig8 = fig.add_subplot(828, projection="3d") x, y = np.meshgrid(slist, s0list) fig8.patch.set_alpha(0.0) fig8.plot_wireframe(x, y, V, linewidth=1.0, color = color_plots) #cmap = cm.RdYlGn, for surface plotting fig8.set_title('BSM %s Option Value vs. Volatility vs. Underlying Price' % optiontype, color = color_plots, fontsize = title_size) fig8.set_xlabel('Volatility$\sigma$', color = color_plots, fontsize = xlabel_size) fig8.set_ylabel('Stock/Underlying Price ($)', color = color_plots, fontsize = ylabel_size)
fig8.set_zlabel('%s Option Value' % optiontype, color = color_plots, fontsize = zlabel_size)

"""#### Monte Carlo Option Sensitivity"""

def mc_plot_values(function,S0, K, T, t, r, q, sigma, n_paths, n_steps, n_sims, optiontype):
fig = plt.figure(figsize=(30,90))
points = 100

# Option(K,T) vs. Strike
klist = np.linspace(K-30, K+30, points)
vlist = [function(S0, K, T, t, r, q, sigma,n_paths,n_steps, n_sims) for K in klist]
fig1.plot(klist, vlist)
fig1.grid()
fig1.set_title('Monte Carlo %s Option Value vs. Strike' % optiontype, color = color_plots, fontsize = title_size)
fig1.set_xlabel('Strike $K$', color = color_plots, fontsize = xlabel_size)
fig1.set_ylabel('%s Option Value' % optiontype, color = color_plots, fontsize = ylabel_size)

# Option(K,T) vs. Strike vs. Underlying Price
klist = np.linspace(K-20, K+20, points)
s0list = np.linspace(S0 - 20, S0 + 20, points)
V = np.zeros((len(s0list), len(klist)), dtype=np.float)
for j in range(len(klist)):
for i in range(len(s0list)):
V[i, j] = function(s0list[i], klist[j], T, t, r, q, sigma,n_paths,n_steps, n_sims)

x, y = np.meshgrid(klist, s0list)
fig2.patch.set_alpha(0.0)
fig2.plot_wireframe(x, y, V, linewidth=1.0, color = color_plots) #cmap = cm.RdYlGn, for surface plotting
fig2.set_title('Monte Carlo %s Option Value vs. Strike vs. Underlying Price' % optiontype, color = color_plots, fontsize = title_size)
fig2.set_xlabel('Strike $K$', color = color_plots, fontsize = xlabel_size)
fig2.set_ylabel('Stock/Underlying Price ($)', color = color_plots, fontsize = ylabel_size) fig2.set_zlabel('%s Option Value' % optiontype, color = color_plots, fontsize = zlabel_size) # Option(K,T) vs. Time fig3 = fig.add_subplot(822) tlist = np.linspace(0.0001, T, points) vlist = [function(S0, K, T, t, r, q, sigma,n_paths,n_steps, n_sims) for T in tlist] fig3.plot(tlist, vlist) fig3.grid() fig3.set_title('Monte Carlo %s Option Value vs. Time' % optiontype, color = color_plots, fontsize = title_size) fig3.set_xlabel('Maturity$T$', color = color_plots, fontsize = xlabel_size) fig3.set_ylabel('%s Option Value' % optiontype, color = color_plots, fontsize = ylabel_size) # Option(K,T) vs. Time vs. Underlying Price tlist = np.linspace(0.0001, T, points) s0list = np.linspace(S0 - 20, S0 + 10, points) V = np.zeros((len(s0list), len(tlist)), dtype=np.float) for j in range(len(tlist)): for i in range(len(s0list)): V[i, j] = function(s0list[i], K, tlist[j], t, r, q, sigma,n_paths,n_steps, n_sims) fig4 = fig.add_subplot(824, projection="3d") x, y = np.meshgrid(tlist, s0list) fig4.patch.set_alpha(0.0) fig4.plot_wireframe(x, y, V, linewidth=1.0, color = color_plots) #cmap = cm.RdYlGn, for surface plotting fig4.set_title('Monte Carlo %s Option Value vs. Time vs. Underlying Price' % optiontype, color = color_plots, fontsize = title_size) fig4.set_xlabel('Maturity$T$', color = color_plots, fontsize = xlabel_size) fig4.set_ylabel('Stock/Underlying Price ($)', color = color_plots, fontsize = ylabel_size)
fig4.set_zlabel('%s Option Value' % optiontype, color = color_plots, fontsize = zlabel_size)

# Option(K,T) vs. r
rlist = np.linspace(0, r, points)
vlist = [function(S0, K, T, t, r, q, sigma,n_paths,n_steps, n_sims) for r in rlist]
fig5.plot(rlist, vlist)
fig5.grid()
fig5.set_title('Monte Carlo %s Option Value vs. r' % optiontype, color = color_plots, fontsize = title_size)
fig5.set_xlabel('Risk-free rate $r$', color = color_plots, fontsize = xlabel_size)
fig5.set_ylabel('%s Option Value' % optiontype, color = color_plots, fontsize = ylabel_size)

# Option(K,T) vs. r vs. Underlying Price
rlist = np.linspace(0, r, points)
s0list = np.linspace(S0 - 20, S0 + 20, points)
V = np.zeros((len(s0list), len(rlist)), dtype=np.float)
for j in range(len(rlist)):
for i in range(len(s0list)):
V[i, j] = function(s0list[i], K, T, t, rlist[j], q, sigma,n_paths,n_steps, n_sims)

x, y = np.meshgrid(rlist, s0list)
fig6.patch.set_alpha(0.0)
fig6.plot_wireframe(x, y, V, linewidth=1.0, color = color_plots) #cmap = cm.RdYlGn, for surface plotting
fig6.set_title('Monte Carlo %s Option Value vs. r vs. Underlying Price' % optiontype, color = color_plots, fontsize = title_size)
fig6.set_xlabel('Risk-free rate $r$', color = color_plots, fontsize = xlabel_size)
fig6.set_ylabel('Stock/Underlying Price ($)', color = color_plots, fontsize = ylabel_size) fig6.set_zlabel('%s Option Value' % optiontype, color = color_plots, fontsize = zlabel_size) # Option(K,T) vs. Implied Vol. fig7 = fig.add_subplot(826) slist = np.linspace(0.01, sigma, points) vlist = [function(S0, K, T, t, r, q, sigma,n_paths,n_steps, n_sims) for sigma in slist] fig7.plot(slist, vlist) fig7.grid() fig7.set_title('Monte Carlo %s Option Value vs. Volatility' % optiontype, color = color_plots, fontsize = title_size) fig7.set_xlabel('Volatility$\sigma$', color = color_plots, fontsize = xlabel_size) fig7.set_ylabel('%s Option Value' % optiontype, color = color_plots, fontsize = ylabel_size) # Option(K,T) vs. Volatility vs. Underlying Price slist = np.linspace(0.01, sigma, points) s0list = np.linspace(S0 - 20, S0 + 20, points) V = np.zeros((len(s0list), len(slist)), dtype=np.float) for j in range(len(slist)): for i in range(len(s0list)): V[i, j] = function(s0list[i], K, T, t, r, q, slist[j],n_paths,n_steps, n_sims) fig8 = fig.add_subplot(828, projection="3d") x, y = np.meshgrid(slist, s0list) fig8.patch.set_alpha(0.0) fig8.plot_wireframe(x, y, V, linewidth=1.0, color = color_plots) #cmap = cm.RdYlGn, for surface plotting fig8.set_title('Monte Carlo %s Option Value vs. Volatility vs. Underlying Price' % optiontype, color = color_plots, fontsize = title_size) fig8.set_xlabel('Volatility$\sigma$', color = color_plots, fontsize = xlabel_size) fig8.set_ylabel('Stock/Underlying Price ($)', color = color_plots, fontsize = ylabel_size)
fig8.set_zlabel('%s Option Value' % optiontype, color = color_plots, fontsize = zlabel_size)

"""### Lévy Approximation"""

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
import time
r = 0.1
sigma = 0.3
T = 1
K = 90
S0 = 100

############### Exact Formulas ###############

def ExactPricing(r,T,K,S0,sigma):
d1 = (np.log(S0/K) + (r + (sigma**2)/2)*T)/(sigma*np.sqrt(T))
d2 = d1 - sigma*np.sqrt(T)
return S0*norm.cdf(d1) - K * np.exp(-r*T)* norm.cdf(d2)

def geometricExactPricing(r,T,K,S0,sigma):
d1 = (np.log(S0/K) + (r+(sigma**2)/6)*T/2)/(sigma*np.sqrt(T/3))
d2 = (np.log(S0/K) + (r-(sigma**2)/2)*T/2)/(sigma*np.sqrt(T/3))
return S0*np.exp(-(r + (sigma**2)/6)*T/2) * norm.cdf(d1) - K * np.exp(-r*T)* norm.cdf(d2)

############## Calcul de l'espérance pour la lognormal

def esperance(r = 0.15, sigma = 0.3, T = 1, S0 = 100):
return ((np.exp(r*T)-1)*S0)/(r*T)

############## Calcul de la variance pour la lognormal

def variance(r = 0.15, sigma = 0.3, T = 1, S0 = 100):
temp1 = (2*S0**2/((T**2)*(r+sigma**2)))
temp2 = (np.exp((2*r+sigma**2)*T)-1)/(2*r+sigma**2)
temp3 = (np.exp(r*T)-1)/r
temp = temp1 * (temp2 - temp3)
return temp - (((np.exp(r*T)-1)*S0)/(r*T))**2

def ecarttype(r = 0.15, sigma = 0.3, T = 1, S0 = 100):
return np.sqrt(variance(r, sigma, T, S0))

esp = esperance()
var = variance()
EC = ecarttype()

mulog = np.log(esp**2/np.sqrt(var + esp**2))
sigmalog = np.sqrt(np.log(1 + (var/(esp**2))))

################# Monte Carlo Simulation ##########

def Monte_Carlo_lognormal(iterations = 1000000, K = 95, T = 1, r = 0.15):
lognorm = np.random.lognormal(mulog, sigmalog, iterations)
Price = np.exp(-r*T) * np.maximum(lognorm - K, 0)
std = np.std(Price)
mean = np.mean(Price)
return mean+(1.96*std/np.sqrt(iterations)), mean, mean-(1.96*std/np.sqrt(iterations))

def Monte_Carlo_lognormal_it(iterations = 1000000, K = 95, T = 1, r = 0.15):
lognorm = np.random.lognormal(mulog, sigmalog, iterations)
Price = np.exp(-r*T) * np.maximum(lognorm - K, 0)
lower_list = []
upper_list = []
mean_list = []
for i in range(1, iterations, 100):
std = np.std(Price[:i])
mean = np.mean(Price[:i])
mean_list.append(mean)
upper_list.append(mean+(1.96*std/np.sqrt(i)))
lower_list.append(mean-(1.96*std/np.sqrt(i)))
return upper_list, mean_list, lower_list

#dataplot = Monte_Carlo_lognormal_it(iterations = 1000000, K = 95, T = 1, r = 0.15)
plt.plot(dataplot[0])
plt.plot(dataplot[1])
plt.plot(dataplot[2])
plt.xlabel('$n_(simulations)$',fontsize = xlabel_size, color=color_plots)
plt.ylabel('C(K=95,T=1)',fontsize = ylabel_size, color=color_plots)
plt.title('Expected Call Option Price using the Lévy Approximation ($S_0=100, K=95, T=1$')
data = Monte_Carlo_lognormal(iterations = 1000000, K = 95, T = 1, r = 0.15)
print(data)

"""### **Control Variates**

#### Conventional Implementation
"""

# Call Options:
def mc_call_control_variates(S0, K, T, t, r, q, sigma,iterations,timesteps, simulations):
c_cv_temp = []
c_cv = []
c_upper_cv = []
c_lower_cv = []
for i in range(1,simulations):
seed = random.randint(1,100000)
S = gbm_paths(S0, K, T, t, r, q, sigma, seed, iterations, timesteps)
S_arithm = np.mean(S)
S_geom = np.exp(np.mean(np.log(S)))
payoff_arithm = np.exp(-r * T) * max(S_arithm - K, 0)
payoff_geom = np.exp(-r * T) * max(S_geom - K, 0)
c_cv_temp.append(payoff_arithm - payoff_geom)
c_cv.append(np.mean(c_cv_temp))
c_upper_cv.append(np.mean(c_cv) + 1.96 * np.std(c_cv)/np.sqrt(i))
c_lower_cv.append(np.mean(c_cv) - 1.96 * np.std(c_cv)/np.sqrt(i))
return c_cv, c_upper_cv, c_lower_cv

# Put Options:
def mc_put_control_variates(S0, K, T, t, r, q, sigma,iterations,timesteps,simulations):
p_cv_temp = []
p_cv = []
p_upper_cv = []
p_lower_cv = []
for i in range(1,simulations):
seed = random.randint(1,100000)
S = gbm_paths(S0, K, T, r, q, sigma, seed, iterations, timesteps)
S_arithm = np.mean(S)
S_geom = np.exp(np.mean(np.log(S)))
payoff_arithm = np.exp(-r * T) * max(K - S_arithm, 0)
payoff_geom = np.exp(-r * T) * max(K - S_geom, 0)
p_cv.append(payoff_geom - payoff_arithm)
p_cv.append(np.mean(p_cv_temp))
p_upper_cv.append(np.mean(p_cv) + 1.96 * np.std(p_cv)/np.sqrt(i))
p_lower_cv.append(np.mean(p_cv) - 1.96 * np.std(p_cv)/np.sqrt(i))
return p_cv, p_upper_cv, p_lower_cv

"""#### JAX Implementation"""

# Call Options:
def mc_call_cv(S0, K, T, t, r, q, sigma, n_paths, n_steps, n_sims):
c_cv = []
c_cv_temp = []
c_upper_cv = []
c_lower_cv = []
K = np.full(n_paths,K)
payoff_0 = np.zeros(n_paths)
for i in range(1,n_sims):
S = gbm_paths_jax(S0, r, q, sigma, n_steps, n_paths)
S_arithm = jnp.mean(S, axis=0)
S_geom = jnp.exp(jnp.mean(jnp.log(S), axis=0))
payoff_arithm = jnp.mean(jnp.exp(-r*T) * jnp.maximum(S_arithm - K, payoff_0))
payoff_geom = jnp.mean(jnp.exp(-r*T) * jnp.maximum(S_geom - K, payoff_0))
c_cv_temp.append(jnp.mean(payoff_arithm - payoff_geom))
c_cv.append(c_cv_temp)
c_upper_cv.append(np.mean(c_cv) + 1.96 * np.std(c_cv)/np.sqrt(i))
c_lower_cv.append(np.mean(c_cv) - 1.96 * np.std(c_cv)/np.sqrt(i))
return c_cv, c_upper_cv, c_lower_cv

# Put Options:
def mc_put_cv(S0, K, T, t, r, q, sigma, n_paths, n_steps, n_sims):
p_cv = []
p_cv_temp = []
p_upper_cv = []
p_lower_cv = []
K = np.full(n_paths,K)
payoff_0 = np.zeros(n_paths)
for i in range(1,n_sims):
S = gbm_paths_jax(S0, r, q, sigma, n_steps, n_paths)
S_arithm = jnp.mean(S, axis=0)
S_geom = jnp.exp(jnp.mean(jnp.log(S), axis=0))
payoff_arithm = jnp.mean(jnp.exp(-r*T) * jnp.maximum(K - S_arithm, payoff_0))
payoff_geom = jnp.mean(jnp.exp(-r*T) * jnp.maximum(K - S_geom, payoff_0))
p_cv_temp.append(jnp.mean(payoff_arithm - payoff_geom))
p_cv.append(p_cv_temp)
p_upper_cv.append(np.mean(p_cv_temp) + 1.96 * np.std(p_cv_temp)/np.sqrt(i))
p_lower_cv.append(np.mean(p_cv_temp) - 1.96 * np.std(p_cv_temp)/np.sqrt(i))
return p_cv, p_upper_cv, p_lower_cv

"""## **Performance Tests**

### Black-Scholes vs. Geometric Monte Carlo

#### Test Cases:
"""

# Call Options
print('- - - - - - - - - - - - - - - - -')
n_paths = 100000  # Number of paths per simulation
n_steps = 252     # Number of steps N (1 monitoring point per business day)
n_sims = 100      # Number of simulations

for i in range(95,110,5):
bs = bsm_call(S0, i, T, t, r, q, sigma)
mc = mc_call_geom(S0, i, T, t, r, q, sigma, n_paths, n_steps, n_sims)
error = np.abs((bs-mc)/bs)*100
print('G_c(K=%d,T=%d):'%(i,T))
print('Black-Scholes:',bs)
print('Monte Carlo:', mc)
print('Error: %f%%' % error)
print('- - - - - - - - - - - - - - - - -')

# Put Options
print('- - - - - - - - - - - - - - - - -')
n_paths = 100000  # Number of paths per simulation
n_steps = 504     # Number of steps N (2 monitoring points per day)
n_sims = 100      # Number of simulations

for i in range(95,110,5):
bs = bsm_put(S0, i, T, t, r, q, sigma)
mc = mc_put_geom(S0, i, T, t, r, q, sigma, n_paths, n_steps, n_sims)
error = np.abs((bs-mc)/bs)*100
print('G_p(K=%d,T=%d):'%(i,T))
print('Black-Scholes:',repr(bs))
print('Monte Carlo:',repr(mc))
print('Error:',repr(error),'%')
print('- - - - - - - - - - - - - - - - -')

"""#### Final Results:

**Call Options:**
- - - - - - - - - - - - - - - - -
G_c(K=95,T=1):
Black-Scholes: 12.508538481017892
Monte Carlo: 12.447389441700023
Error: 0.488858%
- - - - - - - - - - - - - - - - -
G_c(K=100,T=1):
Black-Scholes: 9.61215966961383
Monte Carlo: 9.638977665500253
Error: 0.279001%
- - - - - - - - - - - - - - - - -
G_c(K=105,T=1):
Black-Scholes: 7.185865725921304
Monte Carlo: 7.171154551742235
Error: 0.204724%
- - - - - - - - - - - - - - - - -
**Put Options:**
- - - - - - - - - - - - - - - - -
G_p(K=95,T=1):
Black-Scholes: 2.194652455717922
Monte Carlo: 2.1827188484716005
Error: 0.543758407634428 %
- - - - - - - - - - - - - - - - -
G_p(K=100,T=1):
Black-Scholes: 3.6018135264391438
Monte Carlo: 3.5923289672198133
Error: 0.26332732524071434 %
- - - - - - - - - - - - - - - - -
G_p(K=105,T=1):
Black-Scholes: 5.479059464871914
Monte Carlo: 5.455180840115154
Error: 0.4358161270169417 %
- - - - - - - - - - - - - - - - -

### Asian call options with Arithmetic Average: Linetsky Test Cases

#### Introduction:

The arithmetic average option pricer is benchmarked against the test cases in Table B of the Linetsky paper 1 *Exotic spectra, Risk magazine, April 2002. V. Linetsky* (reproduced below):

**%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%**

**C   r       σ     T   S0   EE       MC           %Err**

1   0.0200   0.10   1   2.0   0.05602   0.0559860415   0.017

2   0.1800   0.30   1   2.0   0.21850   0.2183875466   0.059

3   0.0125   0.25   2   2.0   0.17250   0.1722687410   0.063

4   0.0500   0.50   1   1.9   0.19330   0.1931737903   0.084

5   0.0500   0.50   1   2.0   0.24650   0.2464156905   0.095

6   0.0500   0.50   1   2.1   0.30640   0.3062203648   0.106

7   0.0500   0.50   2   2.0   0.35030   0.3500952190   0.146

**%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%**

**Notes:**
* *(EE = Eigenfunction Expansion i.e. the Black-Scholes analytic result in
this algorithm; MC = Monte-Carlo estimate);*
* *All test cases have a strike K = 2.0 and a dividend yield q = 0.0;*

#### Test Cases:
"""

print('- - - - - - - - - - - - - - - - -')
# Test Case 1 Parameters:
# r = 0.02
# sigma = 0.1
# T = 1
# S0 = 2.0
n_paths = 100000  # Number of paths per simulation
n_steps = 252     # Number of steps M (fixed for T = 1 year)
n_sims = 1000     # Number of simulations
lntk1 = 0.0559860415

mc = mc_call_arithm(2.0, 2.0, 1, 0.0, 0.02, 0.0, 0.1,n_paths,n_steps,n_sims)
error = np.abs((lntk1-mc)/lntk1)*100
print('Case 1 (r=0.02, sigma=0.10, T=1, S0=2.0):')
print('G_c(K=%d,T=%d):'%(2.0,1))
print('Linetsky:',repr(lntk1))
print('Monte Carlo:',repr(mc))
print('Error:',repr(error),'%')

print('- - - - - - - - - - - - - - - - -')

print('- - - - - - - - - - - - - - - - -')
# Test Case 2 Parameters:
# r = 0.18
# sigma = 0.30
# T = 1.0
# S0 = 2.0
n_paths = 100000  # Number of paths per simulation
n_steps = 252     # Number of steps M (fixed for T = 1 year)
n_sims = 1000     # Number of simulations
lntk2 = 0.2183875466

mc = mc_call_arithm(2.0, 2.0, 1.0, 0.0, 0.18, 0.0, 0.3,n_paths,n_steps,n_sims)
error = np.abs((lntk2-mc)/lntk2)*100
print('Case 2 (r = 0.18, sigma = 0.30, T = 1, S0 = 2.0):')
print('G_c(K=%d,T=%d):'%(2.0,1))
print('Linetsky:',repr(lntk2))
print('Monte Carlo:',repr(mc))
print('Error:',repr(error),'%')

print('- - - - - - - - - - - - - - - - -')

print('- - - - - - - - - - - - - - - - -')
# Test Case 3 Parameters:
# r = 0.0125
# sigma = 0.25
# T = 2.0
# S0 = 2.0
n_paths = 100000  # Number of paths per simulation
n_steps = 504     # Number of steps M (fixed for T = 2 years)
n_sims = 1000     # Number of simulations
lntk3 = 0.1722687410

mc = mc_call_arithm(2.0, 2.0, 2.0, 0.0, 0.0125, 0.0, 0.25,n_paths,n_steps,n_sims)
error = np.abs((lntk3-mc)/lntk3)*100
print('Case 3 (r = 0.0125, sigma = 0.25, T = 2, S0 = 2.0):')
print('G_c(K=%d,T=%d):'%(2.0,2))
print('Linetsky:',repr(lntk3))
print('Monte Carlo:',repr(mc))
print('Error:',repr(error),'%')

print('- - - - - - - - - - - - - - - - -')

print('- - - - - - - - - - - - - - - - -')
# Test Case 4 Parameters:
# r = 0.05
# sigma = 0.50
# T = 1.0
# S0 = 1.9
n_paths = 100000  # Number of paths per simulation
n_steps = 252     # Number of steps M (fixed for T = 1 year)
n_sims = 1000     # Number of simulations
lntk4 = 0.1931737903

mc = mc_call_arithm(1.9, 2.0, 1.0, 0.0, 0.05, 0.0, 0.50,n_paths,n_steps,n_sims)
error = np.abs((lntk4-mc)/lntk4)*100
print('Case 4 (r = 0.05, sigma = 0.50, T = 1, S0 = 1.9):')
print('G_c(K=%d,T=%d):'%(2.0,1))
print('Linetsky:',repr(lntk4))
print('Monte Carlo:',repr(mc))
print('Error:',repr(error),'%')

print('- - - - - - - - - - - - - - - - -')

print('- - - - - - - - - - - - - - - - -')
# Test Case 5 Parameters:
# r = 0.05
# sigma = 0.50
# T = 1.0
# S0 = 2.0
n_paths = 100000  # Number of paths per simulation
n_steps = 252     # Number of steps M (fixed for T = 1 year)
n_sims = 1000     # Number of simulations
lntk5 = 0.2464156905

mc = mc_call_arithm(2.0, 2.0, 1.0, 0.0, 0.05, q, 0.50,n_paths,n_steps,n_sims)
error = np.abs((lntk5-mc)/lntk5)*100
print('Case 5 (r = 0.05, sigma = 0.50, T = 1, S0 = 2.0):')
print('G_c(K=%d,T=%d):'%(2.0,1))
print('Linetsky:',repr(lntk5))
print('Monte Carlo:',repr(mc))
print('Error:',repr(error),'%')

print('- - - - - - - - - - - - - - - - -')

print('- - - - - - - - - - - - - - - - -')
# Test Case 6 Parameters:
# r = 0.05
# sigma = 0.50
# T = 1.0
# S0 = 2.1
n_paths = 100000  # Number of paths per simulation
n_steps = 252     # Number of steps M (fixed for T = 1 year)
n_sims = 1000     # Number of simulations
lntk6 = 0.3062203648

mc = mc_call_arithm(2.1, 2.0, 1.0, 0.0, 0.05, 0.0, 0.50,n_paths,n_steps,n_sims)
error = np.abs((lntk6-mc)/lntk6)*100
print('Case 6 (r = 0.05, sigma = 0.50, T = 1, S0 = 2.1):')
print('G_c(K=%d,T=%d):'%(2.0,1))
print('Linetsky:',repr(lntk6))
print('Monte Carlo:',repr(mc))
print('Error:',repr(error),'%')

print('- - - - - - - - - - - - - - - - -')

print('- - - - - - - - - - - - - - - - -')
# Test Case 7 Parameters:
# r = 0.05
# sigma = 0.50
# T = 2.0
# S0 = 2.0
n_paths = 100000  # Number of paths per simulation
n_steps = 504     # Number of steps M (fixed for T = 1 year)
n_sims = 1000     # Number of simulations
lntk7 = 0.3500952190

mc = mc_call_arithm(2.0, 2.0, 2.0, 0.0, 0.05, q, 0.50,n_paths,n_steps,n_sims)
error = np.abs((lntk7-mc)/lntk7)*100
print('Case 7 (r = 0.05, sigma = 0.50, T = 2.0, S0 = 2.0):')
print('G_c(K=%d,T=%d):'%(2.0,2))
print('Linetsky:',repr(lntk7))
print('Monte Carlo:',repr(mc))
print('Error:',repr(error),'%')

print('- - - - - - - - - - - - - - - - -')

"""#### Final Results:

- - - - - - - - - - - - - - - - -
Case 1 (r=0.02, sigma=0.10, T=1, S0=2.0):
G_c(K=2,T=1):
Linetsky: 0.0559860415
Monte Carlo: 0.0559778385271486
Error: 0.014651817902507759 %
- - - - - - - - - - - - - - - - -
Case 2 (r = 0.18, sigma = 0.30, T = 1, S0 = 2.0):
G_c(K=2,T=1):
Linetsky: 0.2183875466
Monte Carlo: 0.2170470505500451
Error: 0.6138152430505384 %
- - - - - - - - - - - - - - - - -
Case 3 (r = 0.0125, sigma = 0.25, T = 2, S0 = 2.0):
G_c(K=2,T=2):
Linetsky: 0.172268741
Monte Carlo: 0.1726129894916586
Error: 0.19983224446888723 %
- - - - - - - - - - - - - - - - -
Case 4 (r = 0.05, sigma = 0.50, T = 1, S0 = 1.9):
G_c(K=2,T=1):
Linetsky: 0.1931737903
Monte Carlo: 0.1917548655071399
Error: 0.7345327700287397 %
- - - - - - - - - - - - - - - - -
Case 5 (r = 0.05, sigma = 0.50, T = 1, S0 = 2.0):
G_c(K=2,T=1):
Linetsky: 0.2464156905
Monte Carlo: 0.2472542514106814
Error: 0.3403033747485359 %
- - - - - - - - - - - - - - - - -
Case 6 (r = 0.05, sigma = 0.50, T = 1, S0 = 2.1):
G_c(K=2,T=1):
Linetsky: 0.3062203648
Monte Carlo: 0.3066361640492379
Error: 0.1357843230019895 %
- - - - - - - - - - - - - - - - -
Case 7 (r = 0.05, sigma = 0.50, T = 2.0, S0 = 2.0):
G_c(K=2,T=2):
Linetsky: 0.350095219
Monte Carlo: 0.3475691276119277
Error: 0.7215440974280498 %
- - - - - - - - - - - - - - - - -

### Monte Carlo vs. Number of Trials

The results of 1,50000,100000,150000, ... ,1000000 paths are compared to the Black-Scholes theoretical price. The errors are then log-plotted against the number of paths in order to determine the sensitivity of the pricing algorithm to the number of trials.
"""

# Iterations Cases:
it_paths = [1,100,1000,10000,100000,1000000]

# Call Options:
bsm_c = bsm_call(S0, K, T, t, r, q, sigma)
print('The Black-Scholes Theoretical Call Price is:',bsm_c)
errors_c = []

for i in range(1,len(it_cases)):
mc_c = mc_call_geom(S0, K, T, t, r, q, sigma, it_paths[i], n_steps, n_sims)
print('For %d trials, the Monte Carlo Call Price estimate is:' %it_paths[i],mc_c)
error_c = 100 * (np.abs(mc_c[0]-bsm_c)/bsm_c)
errors_c.append(error_c)

# Errors v. Simulations Log-Plot:
errors_v_runs = plt.figure()
ax.set_title("Errors vs. Simulations (Call)")
ax.plot(np.log(it_cases), np.log(errors_c), "o-")

# Put Options:
bsm_p = bsm_put(S0, K, T, t, r, q, sigma)
print('The Black-Scholes Theoretical Call Price is:',bsm_p)
errors_p = []

for i in range(0,len(it_cases)):
mc_p = mc_put_geom(S0, K, T, t, r, q, sigma, it_cases[i], n_steps, n_sims)
print('For %d iterations, the Monte Carlo Call Price estimate is:' %it_cases[i],mc_p)
error_p = 100 * (np.abs(mc_p - bsm_p)/bsm_p)
errors_p.append(error_p)

# Errors v. Simulations Log-Plot:
errors_v_runs = plt.figure()
ax.set_title("Errors vs. Simulations (Call)")
ax.plot(np.log(it_cases), np.log(errors_p), "o-")

"""### Multiple Control Variates"""

# Base parameters:
S0 = 100          # Spot price
K = 95            # Strike price
T = 1.0           # Maturity in years
r = 0.15          # Risk-free rate
q = 0.0           # Option dividend yield
sigma = 0.3       # Volatility
n_paths = 100000       # Number of paths per simulation
n_steps = 252   # Number of steps M
n_sims = 100     # Number of simulations

# Call Options:
mc_call_cv_0 = mc_call_cv(S0, K, T, t, r, q, sigma, n_paths, n_steps, n_sims)
c_cv_values = mc_call_cv_0[0]
c_cv_upper = mc_call_cv_0[1]
c_cv_lower = mc_call_cv_0[2]
data_call = [c_cv_upper, c_cv_lower, c_cv_values]
for i in range(0,len(data_call)):
plt.plot(data_call[i])
plt.title('$Arithmetic - Geometric Call Option Value Difference',color=color_plots,fontsize=title_size) plt.xlabel('$\Delta$t', color=color_plots, fontsize=xlabel_size) plt.ylabel('Difference', color=color_plots, fontsize=ylabel_size) plt.show() # Put Options: mc_put_cv = mc_put_control_variates(S0, K, T, t, r, q, sigma, n_paths, n_steps, n_sims) p_cv_values = mc_put_cv[0] p_cv_upper = mc_put_cv[1] p_cvv_lower = mc_put_cv[2] data_put = [p_cv_upper, p_cvv_lower, p_cv_values] for i in range(0,len(data_put)): plt.plot(data_put[i]) plt.title('$Arithmetic - Geometric Put Option Value Difference',color=color_plots,fontsize=title_size)
plt.xlabel('$\Delta$t', color=color_plots, fontsize=xlabel_size)
plt.ylabel('Difference', color=color_plots, fontsize=ylabel_size)
plt.show()

"""## **Hedging Test**

### Daily Data Generation

The timesteps factor was set to timesteps = 1,000,000.

Assumptions:
*  Each path is refined to approximately 11 updates/second. The calculation assumes that price movements occur evenly throughout a 24-hour cycle.
* After-hours trading does not move the closing price of the security each day i.e. the price of the security at 9:30am (open market) is the same as the closing price the previous day.
* The security pays no dividend.
"""

# Define parameters:
n_steps = 1000000 # Granularity of ~11 movements/sec
n_paths = 5 # each path represents 1 business day, and a full business week of trading is generated for a hypothetical security.
S_geometric = []
call = []
call_deltas = []
call_gammas = []
time_vector = np.arange(1,n_paths+1,1) # For plotting

# Initialize values:
S0 = 100
K = 95
t = 0.0
r = 0.15
q = 0.0
sigma = 0.3
n_paths = 100000  # Number of paths per simulation
n_steps = 252     # Number of steps N (1 monitoring point per business day)
n_sims = 100      # Number of simulations

# Modified Geometric Brownian Path Generator; The last element of each price path becomes the first element of the next:
def gbm_paths_hedging(S0,T,r,q,sigma,seed,n_paths,n_steps):
np.random.seed(seed)
dt = T/n_steps
bt = np.random.randn(int(n_paths), int(n_steps))
S = S0 * np.cumprod((np.exp(sigma * np.sqrt(dt) * bt + (r - q - 0.5 * (sigma**2)) * dt)), axis = 1)
S[0][0] = S0
for i in range(0,len(S)-1):
S[i+1][0] = S[i][-1]
return S

# Plot Price Paths:
for i in range(0,5):
plt.plot(np.transpose(daily_data[i]), label = ('Day %s' % str(i+1)))
plt.title('Price Path, Day %s' % str(i+1), color = color_plots, fontsize = title_size)
plt.xlabel('Timesteps', color = color_plots, fontsize = xlabel_size)
plt.ylabel("Underlying Price ($)", color = color_plots, fontsize = ylabel_size) plt.legend(fontsize = legend_size) plt.show() """### Portfolio Sensitivity #### Options on Daily Stock Data """ # Calculate option prices: n_paths = 100000 n_steps = 252 for i in range(1,5): S = gbm_paths_jax(100,0.15,0.0,0.3,n_steps, n_paths) S_geom = np.exp(np.mean(np.log(S))) S_geometric.append(S_geom) call.append(np.exp(-r * (T-t) * max(S_geom - K, 0))) call_deltas.append(bsm_delta(S0,K,T,t,r,q,sigma,'Call')) t = i*T/n_steps print(S_geometric) """#### BSM Values (Control) **1st Order Greeks: Sensitivity to Strike and Stock/Underlying Price** """ fig, ax = plt.subplots(nrows=3, ncols=2, sharex=False, sharey=False, figsize=(30, 30)) fig.suptitle('1st Order Greeks: Sensitivity to Strike and Stock/Underlying Price.$S_0 = %d$' % S0, color = color_plots, fontsize=title_size) #Variables: klist = [85, 95, 100, 105, 115] S0list = np.arange(50,150) plt.subplot(321) for i in klist: c = [bsm_delta(S0, i, T, t, r, q, sigma, "Call") for S0 in S0list] p = [bsm_delta(S0, i, T, t, r, q, sigma, "Put") for S0 in S0list] plt.plot(S0list, c, label = ("Delta Call K=%i" % i )) plt.plot(S0list, p, label = ("Delta Put K=%i" % i )) plt.xlabel('Stock/Underlying Price ($)', color = color_plots, fontsize = xlabel_size)
plt.ylabel("Delta", color = color_plots, fontsize = ylabel_size)
plt.legend(fontsize = legend_size)

plt.subplot(322)
for i in klist:
c = [bsm_gamma(S0, i, T, t, r, q, sigma, "Call") for S0 in S0list]
p = [bsm_gamma(S0, i, T, t, r, q, sigma, "Put") for S0 in S0list]
plt.plot(S0list, c, label = ("Gamma Call K=%i" % i ))
plt.plot(S0list, p, label = ("Gamma Put K=%i" % i ))

plt.xlabel('Stock/Underlying Price ($)', color = color_plots, fontsize = xlabel_size) plt.ylabel("Gamma", color = color_plots, fontsize = ylabel_size) plt.legend(fontsize = legend_size) plt.subplot(323) for i in klist: c = [bsm_vega(S0, i, T, t, r, q, sigma, "Call") for S0 in S0list] p = [bsm_vega(S0, i, T, t, r, q, sigma, "Put") for S0 in S0list] plt.plot(S0list, c, label = ("Vega Call K=%i" % i )) plt.plot(S0list, p, label = ("Vega Put K=%i" % i )) plt.xlabel('Stock/Underlying Price ($)', color = color_plots, fontsize = xlabel_size)
plt.ylabel("Vega", color = color_plots, fontsize = ylabel_size)
plt.legend(fontsize = legend_size)

plt.subplot(324)
for i in klist:
c = [bsm_rho(S0, i, T, t, r, q, sigma, "Call") for S0 in S0list]
p = [bsm_rho(S0, i, T, t, r, q, sigma, "Put") for S0 in S0list]
plt.plot(S0list, c, label = ("Rho Call K=%i" % i ))
plt.plot(S0list, p, label = ("Rho Put K=%i" % i ))

plt.xlabel('Stock/Underlying Price ($)', color = color_plots, fontsize = xlabel_size) plt.ylabel("Rho", color = color_plots, fontsize = ylabel_size) plt.legend(fontsize = legend_size) plt.subplot(325) for i in klist: c = [bsm_theta(S0, i, T, t, r, q, sigma, "Call") for S0 in S0list] p = [bsm_theta(S0, i, T, t, r, q, sigma, "Put") for S0 in S0list] plt.plot(S0list, c, label = ("Theta Call K=%i" % i )) plt.plot(S0list, p, label = ("Theta Put K=%i" % i )) plt.xlabel('Stock/Underlying Price ($)', color = color_plots, fontsize = xlabel_size)
plt.ylabel("Theta", color = color_plots, fontsize = ylabel_size)
plt.legend(fontsize = legend_size)

plt.subplot(326)
for i in klist:
c = [bsm_charm(S0, i, T, t, r, q, sigma, "Call") for S0 in S0list]
p = [bsm_charm(S0, i, T, t, r, q, sigma, "Put") for S0 in S0list]
plt.plot(S0list, c, label = ("Charm Call K=%i" % i ))
plt.plot(S0list, p, label = ("Charm Put K=%i" % i ))

plt.xlabel('Stock/Underlying Price ($)', color = color_plots, fontsize = xlabel_size) plt.ylabel("Charm", color = color_plots, fontsize = ylabel_size) plt.legend(fontsize = legend_size) plt.show() """**1st Order Greeks: Sensitivity to Risk-Free Rate + Stock/Underlying Price**""" fig, ax = plt.subplots(nrows=3, ncols=2, sharex=False, sharey=False, figsize=(30, 30)) fig.suptitle('1st Order Greeks: Sensitivity to Risk-Free Rate + Stock/Underlying Price', color = color_plots, fontsize=title_size, fontweight='bold') # Variables: rlist = [0.0,0.01,0.1] S0list = np.arange(50,150) K = 95 r = 0.15 sigma = 0.3 T = 1.0 t = 0.0 q = 0.0 plt.subplot(321) for i in rlist: c = [bsm_delta(S0, K, T, t, i, q, sigma, "Call") for S0 in S0list] p = [bsm_delta(S0, K, T, t, i, q, sigma, "Put") for S0 in S0list] plt.plot(c, label = ("Delta Call r=%.2f" % i )) plt.plot(p, label = ("Delta Put r=%.2f" % i )) plt.xlabel('Stock/Underlying Price ($)', color = color_plots, fontsize = xlabel_size)
plt.ylabel("Delta", color = color_plots, fontsize = ylabel_size)
plt.legend(fontsize = legend_size)

plt.subplot(322)
for i in rlist:
c = [bsm_gamma(S0, K, T, t, i, q, sigma, "Call") for S0 in S0list]
p = [bsm_gamma(S0, K, T, t, i, q, sigma, "Put") for S0 in S0list]
plt.plot(c, label = ("Gamma Call r=%.2f" % i ))
plt.plot(p, label = ("Gamma Put r=%.2f" % i ))

plt.xlabel('Stock/Underlying Price ($)', color = color_plots, fontsize = xlabel_size) plt.ylabel("Gamma", color = color_plots, fontsize = ylabel_size) plt.legend(fontsize = legend_size) plt.subplot(323) for i in rlist: c = [bsm_vega(S0, K, T, t, i, q, sigma, "Call") for S0 in S0list] p = [bsm_vega(S0, K, T, t, i, q, sigma, "Put") for S0 in S0list] plt.plot(c, label = ("Vega Call r=%.2f" % i )) plt.plot(p, label = ("Vega Put r=%.2f" % i )) plt.xlabel('Stock/Underlying Price ($)', color = color_plots, fontsize = xlabel_size)
plt.ylabel("Vega", color = color_plots, fontsize = ylabel_size)
plt.legend(fontsize = legend_size)

plt.subplot(324)
for i in rlist:
c = [bsm_rho(S0, K, T, t, i, q, sigma, "Call") for S0 in S0list]
p = [bsm_rho(S0, K, T, t, i, q, sigma, "Put") for S0 in S0list]
plt.plot(c, label = ("Rho Call r=%.2f" % i ))
plt.plot(p, label = ("Rho Put r=%.2f" % i ))

plt.xlabel('Stock/Underlying Price ($)', color = color_plots, fontsize = xlabel_size) plt.ylabel("Rho", color = color_plots, fontsize = ylabel_size) plt.legend(fontsize = legend_size) plt.subplot(325) for i in rlist: c = [bsm_theta(S0, K, T, t, i, q, sigma, "Call") for S0 in S0list] p = [bsm_theta(S0, K, T, t, i, q, sigma, "Put") for S0 in S0list] plt.plot(c, label = ("Theta Call r=%.2f" % i )) plt.plot(p, label = ("Theta Put r=%.2f" % i )) plt.xlabel('Stock/Underlying Price ($)', color = color_plots, fontsize = xlabel_size)
plt.ylabel("Theta", color = color_plots, fontsize = ylabel_size)
plt.legend(fontsize = legend_size)

plt.subplot(326)
for i in rlist:
c = [bsm_charm(S0, K, T, t, i, q, sigma, "Call") for S0 in S0list]
p = [bsm_charm(S0, K, T, t, i, q, sigma, "Put") for S0 in S0list]
plt.plot(c, label = ("Charm Call r=%.2f" % i ))
plt.plot(p, label = ("Charm Put r=%.2f" % i ))

plt.xlabel('Stock/Underlying Price ($)', color = color_plots, fontsize = xlabel_size) plt.ylabel("Charm", color = color_plots, fontsize = ylabel_size) plt.legend(fontsize = legend_size) plt.show() """**BSM Option Sensitivity Definitions Tests**""" # Variables: S0 = 100 K = 95 T = 1.0 t = 0.0 r = 0.01 q = 0.0 sigma = 0.1 n_paths = 1000# Number of paths per simulation n_steps = 252 # Number of steps M n_sims = 7 # Number of simulations mc_call_cv(S0,K,T,t,r,q,sigma,n_paths,n_steps,n_sims) """**Monte Carlo Option Sensitivity Definitions Tests**""" # Greeks Plot tool for all Greeks: S0 = np.linspace(50, 150, 100) K = np.linspace(95.0,95.0, 100) T = np.linspace(0.1, 1.0, 100) t = np.linspace(0.0, 0.0, 100) r = np.linspace(0.15, 0.15, 100) q = np.linspace(0.0,0.0,100) sigma = np.linspace(0.30, 0.30, 100) x, y = np.meshgrid(S0, T) fig = plt.figure(figsize=(30,90)) greeks = [bsm_delta,bsm_gamma,bsm_vega,bsm_theta,bsm_rho,bsm_vanna,bsm_vomma, bsm_phi, bsm_charm] greeks_names = ['Delta','Gamma','Vega','Theta','Rho','Vanna','Volga','Phi','Charm'] for i in range(len(greeks)): ax = fig.add_subplot(len(greeks), 2, i+1, projection='3d') ax.patch.set_alpha(0.0) greeks_plot_tool(greeks[i],"Stock/Underlying Price ($)", "Expiry (T)", S0, K, T, t, r, q, sigma, x, y, "Call", greeks_names[i], ax)
plt.show()

"""### Hedged-Unhedged Portfolio Variability"""

# Shift values:
def convert(lst):
return [ i - 50 for i in lst ]

# Calculate Deltas

for j in range(0,len(daily_data)):
for i in range(0,n_steps):
call_delta = bsm_delta(daily_data[j][i],K,T,,r,q,sigma,'Call')
call_delta = call_delta/n_steps
call_deltas.append(call_delta)
print(call_deltas)

for j in range(0,len(daily_data)):
for i in range(0,n_steps):
call_gamma = bsm_gamma(daily_data[j][i],K,T,,r,q,sigma,'Call')
call_gamma = call_delta/n_steps
call_gammas.append(call_gamma)
print(call_gamams)
# Results
#print('Call Deltas:')
#print(call_deltas)

# Hedged Portfolio
hedged = []
for i in range(0,len(daily_data)):
hedged.append(call[i]-call_deltas[i]*S_geometric[i])

#hedged = convert(hedged)
std_unhedged = np.std(call)
std_hedged = np.std(hedged)
print('The standard deviation of the unhedged options portfolio is:',std_unhedged)
print('The standard deviation of the hedged options portfolio is:',std_hedged)

# Plots
plt.plot(time_vector,call, label = 'Unhedged Portfolio')
plt.plot(time_vector,hedged, label = 'Hedged Portfolio')
plt.title('Variability in a Hedged and Unhedged Options Portfolio', color = color_plots, fontsize = title_size)
plt.xlabel('Time (days)', color = color_plots, fontsize = xlabel_size)
plt.ylabel('Option Price', color = color_plots, fontsize = ylabel_size)
plt.legend(fontsize = legend_size)

"""## **Final Project Deliverables**"""

# Test MC Engine:

S0 = 100          # Spot price
K = 95            # Strike price
T = 1.0           # Maturity in years
r = 0.15          # Risk-free rate
q = 0.0           # Option dividend yield
sigma = 0.3       # Volatility
n_paths = 100000       # Number of paths per simulation
n_steps = 252   # Number of steps M
n_sims = 7     # Number of simulations

# mc_plot_values(mc_call_geom,S0, K, T, t, r, q, sigma, n_paths, n_steps, n_sims,'Call')
#for i in range(1,10):
#   print(mc_call_geom(S0, K, T, t, r, q, sigma, n_paths, n_steps, n_sims))
#   print(bsm_call(S0,K,T,t,r,q,sigma))

# Test MC Engine:

S0 = 100          # Spot price
K = 95            # Strike price
T = 1.0           # Maturity in years
r = 0.15          # Risk-free rate
q = 0.0           # Option dividend yield
sigma = 0.3       # Volatility
n_paths = 100000       # Number of paths per simulation
n_steps = 252   # Number of steps M
n_sims = 100     # Number of simulations

plt.plot(mc_call_cv(S0,K,T,t,r,q,sigma,n_paths,n_steps,n_sims))


DOCUMENTS

All requests for copies of the research, please forward to my university email address listed on my homepage.