\documentclass[11pt]{article}
\usepackage[sort]{natbib}
\usepackage{fancyhdr}
\usepackage{bm,amsmath,tikz,bbm,amsfonts,nicefrac,latexsym,amsmath,amsfonts,amsbsy,amscd,amsxtra,amsgen,amsopn,bbm,amsthm,amssymb, enumerate, appendix, listings,subcaption}
\usepackage{pifont}
\usepackage{graphicx}
%\usepackage[margin=0.5in]{geometry}
\usepackage[nottoc,numbib]{tocbibind}
\usepackage[a4paper,left=2cm,right=2cm, bottom=1cm]{geometry}
%\usepackage[toc,page]{appendix}
\lstset{
frame=single,
breaklines=true,
postbreak=\raisebox{0ex}[0ex][0ex]{\ensuremath{\color{red}\hookrightarrow\space}}
}
%-----------------------------------------------------------------------------------USER INTRODUCED PACKAGES BEGINS
%-----------------------------------------------------------------------------------USER INTRODUCED PACKAGES ENDS
%\oddsidemargin 0.1cm
%\topmargin -2.0cm
%\bottommargin -2.0cm
%\textheight 24.0cm
%\textwidth 15.25cm
%\parindent=0pt
%\parskip 1ex
\renewcommand{\baselinestretch}{1.1}
\pagestyle{fancy}
\pagenumbering{roman}
%-----------------------------------------------------------------------------------USER DEFINED COMMANDS/ENVIRONMENTS BEGINS
%\newtheorem{propos}{Propostion}[section]
%\newtheorem{lemma}{Lemma}[section]
%\newtheorem{theorem}{Theorem}[section]
%\newtheorem{corolloary}{Corollory}
%\newtheorem{defn}{Definition}[section]
%\newtheorem{ex}{Example}
%\newtheorem{recall}{Recall:}
%\newtheorem{nt}{Note:}
%\newcommand{\aut}{\textnormal{Aut}}
%\newcommand{\inn}{\textnormal{Inn}}
%\newcommand{\tick}{\ding{51}}
%\newcommand{\cross}{\ding{55}}
%\newcommand*{\rom}[1]{\expandafter\@slowromancap\romannumeral #1@}
\usepackage{amsthm}
%Normal theorem enviroment is for theorems, corollaries, lemmas, propositions, conjectures, criteria,
\newtheorem{thrm}{Theorem}[section]
\newtheorem{prop}{Proposition}[section]
\newtheorem{lem}{Lemma}[section]
\newtheorem{defn}{Definition}[section]
\theoremstyle{definition} %definitions, conditions, problems, and examples
\newtheorem{ex}{Example}[section]
\theoremstyle{remark} %remarks, notes, notation, claims, summaries, acknowledgments, cases, and conclusions
\newtheorem*{claim*}{Claim:}
%-----------------------------------------------------------------------------------USER DEFINED COMMANDS/ENVIRONMENTS ENDS
%\newcounter{mycounter}
%\newtheorem{thrm}{Theorem}[numberby]
%-----------------------------------------------------------------------------------HEADER AND TITLE INFORMATION BEGINS
\lhead{\normalsize \textrm{Numerical and Analytic Methods in Option Pricing}}
\chead{}
\rhead{\normalsize D. Edwards }
\lfoot{\normalsize \textrm{MA4XA}}
\cfoot{\thepage}
\rfoot{Dr. A. Chernov}
\setlength{\fboxrule}{4pt}\setlength{\fboxsep}{2ex}
\renewcommand{\headrulewidth}{0.4pt}
\renewcommand{\footrulewidth}{0.4pt}
\setlength\headheight{14pt}
\addtolength{\textheight}{-54pt}
\title{Numerical and Analytic Methods in Option Pricing}
\author{D. Edwards \\\phantom{} \\ Student Number: $20006274$ \\ \phantom{} \\ Department of Mathematics and Statistics \\ University of Reading }
%-----------------------------------------------------------------------------------HEADER AND TITLE INFORMATION ENDS
\begin{document}
\clearpage\maketitle
\begin{abstract}
In this paper we discuss how to price American, European and Asian options using a geometric Brownian motion model for stock price. We investigate the analytic solution for Black-Scholes differential equation for European options and consider numerical methods for approximating the price of other types of options. These numerical methods include Monte Carlo, binomial trees, trinomial trees and finite difference methods. We conclude our discussion with an investigation of how these methods perform with respect to the changes in different Greeks. Further analysing how the value of a certain Greeks affect the price of a given option.
\end{abstract}
\thispagestyle{empty}
\newpage
\tableofcontents
\listoffigures
\listoftables
\newpage
\pagenumbering{arabic}
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
\section{Introduction}
%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
\subsection{The History and Origin of Options}
Options, in one form or another, have been around for at least centuries. There are reports that they have been around for millennia, being used as early as the sixth century B.C. in Ancient Greece \citep{ES}. Further there are many other instances where we see options crop up in history, clauses in marine cargo contracts, such as those of the Romans and Phoenicians, would now be considered options.
The reason for their wide spread use is, in part, due to the guarantees that they afford to those who buy them. If we were to be traveling to a foreign country to buy something to trade, it is helpful to know when we get back that we are going to be able to sell our product for a price agreed upon prior to the trip. This allows us insurance on the trip, the knowledge of how much we will make when we return. Alternatively, if someone else was to go on a voyage to acquire something we wanted, and we believed that the price of that commodity was to go up, then making sure we could buy it at the current price would allow us to make profit.
The interest then became how to price these contracts. If you want to sell grain for twice its current price in one months’ time, it is obvious that the price of grain is unlikely to double in one month, thus the price of this contract would likely be very high. Conversely, if we were wanting to buy grain at twice the price in one months’ time, the cost of this contract would be very low. Owing again to the unlikely event that the grain does double in price. Thus someone would happily take on the contract cheaply knowing they can sell their grain to you for what is likely to be much more than it would be worth in a months’ time.
%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
\subsection{Option Basics}
Before we may discuss the pricing of \emph{options} we need to introduce a few terms. The most key of these being the definition of \emph{options} and, we introduce this and some other key ideas here following the structure of \citet[Chap. 8]{JCH}.
\begin{defn}[Option]
An \textbf{option} is a contract that gives the buyer the right, but not the obligation, to buy or sell an underlying asset at a specific price (known as the \textbf{strike price}) on or before a certain date.
\end{defn}
By an \emph{underlying asset} we refer to the financial instrument on which the options price is based. This could be futures, stocks, commodities or currency, noting a change in the price of the underlying asset causes a change in the price of the option.
In the introduction we discussed the two types of options in an imprecise way, the buy side and sell side of an option. These concepts are extended within the next definition, to give rise to the two types of options that are used today, which we define here.
\begin{defn}[Call]
A \textbf{call} option is a contract that gives the buyer the right, but not the obligation, to buy an underlying asset at a specific price on or before a certain date (known as the \textbf{expiry date}).
\end{defn}
\begin{defn}[Put]
A \textbf{put} option is a contract that gives the buyer the right, but not the obligation, to sell an underlying asset at a specific price on or before a certain date (known as the \textbf{expiry date}).
\end{defn}
There are many different variations of puts and calls, the most popular being \emph{European, American} and \emph{Asian} option types. Other types of options are often referred to as \emph{exotic options}, these include; basket options, barrier options, binary options, and down-and-out options. Each of these types has a different structure in the way the pay-out is calculated or when it may be \emph{exercised}. This term is often used, and before we can differentiate between the types of options we must introduce it.
\begin{defn}[Exercise]
We say that we \textbf{exercise} an option when we decide to use that option. That is, to buy or sell the underlying asset the option corresponds to.
\end{defn}
The main difference between European and American options is when you can exercise them. The term \textbf{European} is given to those options that may be exercised \textbf{only} at the date of expiration. Alternatively \textbf{American} options may be exercised \textbf{at any time} before, or on, the date of expiration.
For both European and American options, as call options allow us to buy the stock at a specific price, the payoff would be given by $\max(S-K,0)$ and for a put $\max(K-S,0)$, where the stock price at the time the option is exercised (sometimes termed when it reaches maturity) is $S$ and the strike price is $K$.
Asian options are very different from the European and American options. This is because the pay-out of this type of option is not dependent on a single value on or before the option expires.
\begin{defn}[Asian Options]
An option is said to be \textbf{Asian} or sometimes an \textbf{average option} if the pay-out depends on the average price of the underlying asset over a period of time.
\end{defn}
Insofar, we have discussed many different terms relating to options. We seek to give substance to these ideas with two brief examples, one for a put and one for a call. These are similar example as to the ones found in \citet[P. 193]{JCH} but have been adapted to better suit our discussion.
\begin{ex}
Assume that a stock is trading at $\$50$ per share. We could buy a European call option, the right to buy the stock, with a strike price of $\$55$. This is essentially ``betting the price goes up''; here we are hoping that the price rises to above $\$55$, and if this happens we make profit. Let’s say that this option entitles us to buy $100$ shares of the stock and cost us $\$11$. For a European option the pay-out, that is the amount of money we make on this transaction, is dependent on the price at the date of expiration. So we have three cases,
\begin{enumerate}
\item The stock ends below $\$55$. In this case we have lost money, that is the cost of the option $\$11$. This is because we are able to buy the stock for $\$55$, but it is trading at less than that, so exercising this option is worthless. So we have made no money here, yet only lost the cost of the option.
\item The stock rises to or above $\$55.12$. Here we have made a profit, we are able to buy $100$ shares of the stock at $\$55.12$ or above. Hence, we buy these for our strike price $\$55$ and sell them for the price it is trading at. As the stock is on or above $\$55.12$ per share we make $\$55.12 - \$55 = \$0.12$, so for all $100$ shares we make $\$12$. Finally, subtracting the cost of the option $\$11$ we have made at least a dollar. It is worth noting that as the price of the stock could increase up to any amount, our theoretical maximum profit here is infinite.
\item The stock ends between $\$55$ and $\$55.11$. Here we have made a profit. As with the profitable case, we buy the stock for $\$55$ and sell it for the price it is trading at. We see the maximum profit we make is $(\$55 - \$55.11) \times 100 = \$11$. However we have already paid $\$11$ for the option, so we have lost money, yet we still exercise the option. This is because we make some money on the cost of the option back, i.e. if the stock was at $\$55.07$, then we only lose $\$4$ due to the $\$7$ dollar we made from the option.
\end{enumerate}
\end{ex}
\begin{ex}
As our second example we will consider an American put on a stock. Here we have the right to sell a stock for a specified strike price, let’s assume that this is $\$45$. Much like with the call this is a kind of bet that we are placing, here we are ``betting that the stock will decrease in price''. So let’s say this option cost us $\$15$ then if the stock is trading at $\$50$, as with the previous example we either make nothing, reduce our losses or make profit. However this is an American option and as such can be exercised at any time between purchase and expiry dates. Hence, the only way we are guaranteed to lose the cost of the option is if the stock does not ever drop below $\$44.85$.
If it only drops below this value once and we do not exercise we would loose out. This shows the difficulty when choosing to exercise American options.
\end{ex}
%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
\subsection{Preliminaries}
%probability distributions
%
\label{pel}
Lastly before we begin our full discussion of the methods used to price options we must introduce a few basic mathematical results and notations that will be used throughout. These are mainly either results from stochastic calculus \citep{GSS} or probability \citep{DR}.
\begin{defn}[Stochastic variable]
Given a probability space, with events $x$, we can introduce a \textbf{stochastic (or random) variable} as a function of $x$, denoted $f(x)$. In particular the identity function $I(x) = x$ is one such random variable. Continuing on we will use capital letters to denote random variables and lower case letters to denote their values.
Note that the random variable can be either continuous or discrete.
\end{defn}
\begin{defn}[Stochastic process]
If a random variable $X$ is dependent on time, so that it is defined at different instances of time $t_1, t_2, \dots, t_n$, then $X(t)$ is called a \textbf{stochastic process}.
\end{defn}
\begin{defn}[Gaussian process]
A process is said to be \textbf{Gaussian} if all possible distributions $X_t$ are Gaussian. This means that a Gaussian process is characterised fully by the mean and variance of $X_t$.
\end{defn}
\begin{defn}[Characteristic function]
If $X$ is a stochastic variable taking a continuous range of real numbers, its \textbf{characteristic function}, $G(s)$, is defined as,
\begin{align*}
G(s) = \left< \mathrm{e}^{\mathrm{i} sX}\right> = \int P(x)\mathrm{e}^{\mathrm{i} sx}\mathrm{d}x
\end{align*}
where the integral varies over the range of $x$.
%NEEDS TO BE IMAGINARY
\end{defn}
\begin{ex}[Characteristic function of a Gaussian]
The \textbf{Gaussian} or \textbf{normal distribution} is defined as,
\begin{align*}
P(x) = \frac{1}{\sqrt{2\pi\sigma^2}}\exp\left({\frac{-(x-\mu)^2}{2\sigma^2}}\right) \qquad x \in (-\infty,\infty),
\end{align*}
where $\mu$ is the mean and $\sigma^2$ the variance. Hence, from the definition of a characteristic function we have that,
\begin{align*}
G(s)&= \int P(x)\mathrm{e}^{\mathrm{i}sx}\mathrm{d}x \\
&= \int^\infty_{-\infty} \frac{1}{\sqrt{2\pi\sigma^2}}\exp\left({\frac{-(x-\mu)^2}{2\sigma^2} + \mathrm{i}sx}\right) \\
&= \exp\left({-\frac{\sigma^2s^2}{2} + \mathrm{i}s\mu}\right)
\end{align*}
\end{ex}
\begin{defn}[The It\^o integral]
The \textbf{It\^o integral} between two points $t$ and $t_0$ is defined similarly to the Riemann integral. We first begin by discretizing time, a typical type would be $t_i = t_0 + \frac{i}{n}(t - t_0)$ for $i = 0, 1, \dots, n$ %check this!! may be from 1
. Then the It\^o integral of a function $g(X_t,t)$ will be defined as,
\begin{align*}
\int^t_{t_0}g(X_t)\mathrm{d}W_t = \lim_{n \rightarrow \infty}\sum^n_{i=1}g(X_{t_i},t)(W_{t_i} - W_{t_{i-1}})
\end{align*}
where $W_t$ is the Wiener process which we will introduce later. We will use the notation, $\Delta W_i = W_{t_i} - W_{t_{i-1}}$ and $\Delta t = t_i -t_{i-1}$
\end{defn}
\begin{defn}[Moments and their properties]
The \textbf{moment} of order $m$ is defined as,
\begin{align*}
\mu_m =
\begin{cases}
\int x^mP(x)\mathrm{d}x \quad \textnormal{for continuous variables} \\
\sum_i x^m_iP(x_i) \quad\textnormal{for discreet variables.}
\end{cases}
\end{align*}
The \textbf{average} of any function of a stochastic variable $X$, denoted $\mathbb{E}[f(X)]$, and is defined as,
\begin{align*}
\mathbb{E}[f(X)] =
\begin{cases}
\int f(x)P(x)\mathrm{d}x \quad\textnormal{for continuous variables} \\
\sum_i f(x_i)P(x_i) \quad\textnormal{for discreet variables} %check f(x_i) or f(x)
\end{cases}
\end{align*}
Then we may now write the definition of a moment as,
\begin{align*}
\mu_m = \mathbb{E}[x^m]
\end{align*}
Lastly from the our definition, due to the linearity of sums and integrals as well as the fact the measure of a probability space is one, the following hold, for random variables $X$ and $Y$ and $a \in \mathbb{R}$:
\begin{enumerate}
\item $\mathbb{E}[aX]$ $=$ a$\mathbb{E}[X]$,
\item $\mathbb{E}[X + a$ = $\mathbb{E}[X]+ a$,
\item $\mathbb{E}[X+Y] = \mathbb{E}[X] + \mathbb{E}[Y]$,
\item $\mathbb{E}[XY] = \mathbb{E}[X] \mathbb{E}[Y]$ if $X$ and $Y$ are independent.
\end{enumerate}
The last property is derived from the fact that if they are independent then their \textbf{covariance} is zero, which we define below.
\end{defn}
\begin{defn}[Variance and Covariance]
The \textbf{variance} of a stochastic variable $X$ is given by,
\begin{align*}
\textnormal{Var}[X]&=\mathbb{E}[(X-\mathbb{E}[X])^2] \\
&=\mathbb{E}[X^2] - (\mathbb{E}[X])^2.
\end{align*}
Furthermore we have the \textbf{covariance} of two stochastic variables $X$ and $Y$ is defined by,
\begin{align*}
\textnormal{Covar}&=\mathbb{E}[(X-\mathbb{E}[X])(Y-\mathbb{E}[Y])] \\
&=\mathbb{E}[XY] - \mathbb{E}[X]\mathbb{E}[Y]
\end{align*}
\end{defn}
\begin{defn}[Lognormal Distribution]
\label{lnd}
A variable, $X$, is said to be \textbf{lognormally distributed} if its logarithm is normally distributed. Thus, $Y = \ln(X)$ is normally distributed. Further, if $Y$ is normally distributed with mean $a$, and standard deviation $b$, we have that,
\begin{align*}
\mathbb{E}[X] &= \mathrm{e}^{a + \frac{b^2}{2}}, \\
\textnormal{Var}[X] &= (\mathrm{e}^{\sigma^2}-1)\mathrm{e}^{2a + b^2}.
\end{align*}
\end{defn}
\begin{defn}[Cumulative distribution function]
The \textbf{cumulative distribution function}, $C_X(x)$, of a continuous random variable $X$ is given by,
\begin{align*}
C_X(x) = P(X\leq x),
\end{align*}
Where $P(X\leq x)$ is the probability of $X<x$. The cumulative distribution function of the standard normal distribution is given by,
\begin{align*}
\Phi(x) = \frac{1}{\sqrt{2\pi}}\int^x_{-\infty}\mathrm{e}^{-t^2/2}\mathrm{d}t.
\end{align*}
\end{defn}
\begin{defn}[Probability density function]
The \textbf{probability density function} (P.D.F.) of a continuous distribution is the derivative of the cumulative distribution function. This satisfies that if $X$ is a random variables and has P.D.F. $f$ then the expected value of $X$ is given by,
\begin{align*}
\mathbb{E}[X] = \int^\infty_{-\infty}xf(x)\mathrm{d}x
\end{align*}
\end{defn}
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
\section{The Black-Scholes Model for Stocks}
The trouble with pricing options is without knowing the path the stock is likely to travel, it is difficult to price the option, as clearly the price of an option must be dependent on the stock price. Here we will develop a continuous time stochastic model for stocks and through this derive an equation for the price of options. We may then solve this explicitly for European options to gain the formulae for the pricing of European options; this model was first developed by \citet{BS} and further improved by \citet{RM}. Before we develop this model, there are a number of key ideas and concepts we will need to introduce. We will follow much the same derivation as \citet[Chap.9,12-13]{JCH}.
%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
\subsection{Factors Affecting the Price of Stock Options}
\label{fasp1}
As we expect the price of the option to reflect the price of the underlying asset, it seems logical that a factor that affects the price of the underlying asset would also affect the price of the option \citep[Chap. 9]{JCH}. So for stock options we consider the factors that affect stock prices. There are six major factors. These are;
\begin{enumerate}
\item The current price the stock is trading at, $S_0$,
\item The strike price, $K$,
\item The time until expiration, $T$,
\item The volatility of the price of the stock, $\sigma$,
\item The risk-free interest rate, $r$,
\item The dividends that are expected to be paid, $q$.
\end{enumerate}
Here the risk-free interest rate is the theoretical rate of return on a completely risk free investment. The volatility represents how the price varies over time, which we will give a more rigorous definition of this later. How each of these factors affects the price of European and American puts and calls is given in Table \ref{tab:factors}. We will discuss how each of these factors affects the option. Note that in Table \ref{tab:factors} we use $+$ to mean the factor increases the price of the option, a $-$ for a decrease and a $?$ when the relationship is unknown.
\begin{table}[ht]
\centering
\begin{tabular}{|l||c|c|c|c|}
\hline
Variable & European Call & European Put & American Call & American Put \\ \hline
Current stock price & + & - & + & - \\
Strike Price & - & + & - & + \\
Time to expiration & ? & ? & + & + \\
Volatility & + & + & + & + \\
Risk-free rate & + & - & + & - \\
Amount of future dividends & - & + & - & + \\\hline
\end{tabular}
\caption{\label{tab:factors} The effect different factors have on the price of different options }
\end{table}
Most of the ways each factor affects each type of option is intuitivee. If the price goes up, the price of a call goes up as we are likely to see greater increases in the price of the stock and the price of a put goes down as the decrease will not be as fast as the increase. This is due to the fact that, as we will see later, the path the stock takes is the initial price multiplied by some variables. Hence increases tend to happen faster than decreases.
The risk free interest rate is a more complex idea. Within the economy, as interest rates increase, investors expect more return from the option, however the value of any money earnt in the future decreases, due to these interest rates, termed inflation. This increases the price of stocks slightly resulting in a higher chance of calls paying off and less of puts paying off.
With relation to the dividends, as the ex-dividend date approaches, that is the date at which entitlement to dividends changes, the stock price decreases. Hence, due to the relationships for calls and puts, the price of calls decreases and the price of a put increases.
The most anomalous observations are those of the time to expiration and the volatility. The time to expiration for American options increases for both puts and calls. This is due to that, for two options, if the only difference is that one has a longer time to expiration, then the owner of the option with a longer time has all of the same opportunities to exercise and more. This increases the value of the option. For European options, this is not necessarily the case. If we have two options that straddle the ex-dividend date, the one with the shorter life would be worth more.
Finally we consider volatility, we have not defined volatility until this point, we will discuss definition in later sections. For now we merely remark that volatility is a measure of how uncertain we are of the stocks future price. If volatility increases, the chance the stock will do very well or very poorly increases. As our maximum loss from an option contract is the price but the profit is either infinite or large this benefits the owner hugely.
%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
\subsection{Markov and Wiener Processes}
The next concept we will need to introduce is that of a \emph{Markov process}. This is defined as a stochastic process with the \emph{Markov property}; where the future value of the variable is dependent only on its current value and not its history. We often say that a Markov process is memoryless due too this property. We will examine a very specific type of Markov process, known as the \emph{Wiener process} \citep[Chap. 3]{GSS}.
\begin{defn}[Wiener process]
Let $W=(W_t)_{t \in [0, \infty)}$ be a continuous process. We say that this process is a \textbf{Wiener process} if the following properties hold:
\begin{enumerate}[\textnormal{Property} $1.$]
\item We require that if $\Delta t$ is a small period of time, \label{W1}
\begin{align*}
W_{t + \Delta t} - W_t = \epsilon \sqrt{\Delta t},
\end{align*}
where $\epsilon$ has a standardised normal distribution (i.e. a normal distribution with a mean of zero and standard deviation of one).
\item We require that each increment is independent, so that $W_t$ and $W_s$ are independent for $0 \leq s <t$, \label{W2}
\item The function defined by $t \rightarrow W_t$ is almost surely everywhere continuous (i.e. the probability that $W_t$ is everywhere continuous is one).
\end{enumerate}
\end{defn}
We see this is indeed a Markov process as the attribute described in Property \ref{W2} is precisely the Markov property. Also note that the first property shows that the Wiener process is normally distributed. An interesting implication of the Wiener process is found when we consider these normally distributed Markovian variables.
\begin{prop}
\label{prop1}
Let $X$ and $Y$ be two independent normally distributed variables with means $\mu_X$ and $\mu_Y$ respectively and variances $\sigma^2_X$ and $\sigma^2_Y$ respectively. Then the variable defined as $Z = X+Y$ is normally distributed with mean $\mu_X+\mu_Y$ and variance $\sigma^2_X + \sigma^2_Y$. We will use the notation $A \sim N(\mu, \sigma^2)$ to mean ``A variable $A$ is normally distributed with mean $\mu$ and variance $\sigma^2$''.
\end{prop}
\begin{proof}
%i's here need to be imaginary i's not variable i's
The characteristic function of $X$ and $Y$ are by definition,
\begin{align*}
G_X(s) = \mathbb{E}\left[ \mathrm{e}^{isX}\right] \text{ and } G_Y(s) = \mathbb{E}\left[ \mathrm{e}^{isY}\right].
\end{align*}
We sum these two variables to generate a new variable $Z = X+Y$, which has characteristic function given by \citep{DR},
\begin{align*}
G_Z(s) = \mathbb{E}\left[ \mathrm{e}^{\mathrm{i}sZ}\right] = \mathbb{E}\left[ \mathrm{e}^{\mathrm{i}sX + \mathrm{i}sY)}\right] = \mathbb{E}\left[ \mathrm{e}^{\mathrm{i}sY}\right] \mathbb{E}\left[ \mathrm{e}^{\mathrm{i}sX}\right] = G_X(s)G_Y(s),
\end{align*}
as these variables are independent. Now using the general formula for the characteristic function of a normal distribution from section \ref{pel} we have that,
\begin{align*}
G_z(s) = G_X(s)G_Y(s) &= \exp\left(\mathrm{i}t\mu_x - \frac{\sigma^2_Xs^2}{2}\right)\exp\left( \mathrm{i}t\mu_Y - \frac{\sigma^2_Ys^2}{2}\right) \\
&= \exp\left(\mathrm{i}t(\mu_X + \mu_Y) - \frac{(\sigma_x^2 + \sigma_Y^2)s^2}{2} \right).
\end{align*}
Which is precisely the characteristic function of a normally distributed variable with mean $\mu_X + \mu_Y$ and variance $\sigma^2_X + \sigma^2_Y$. Hence, $W \sim (\mu_X + \mu_Y, \sigma^2_X + \sigma_Y^2)$ mean $W$ has normal distribution with mean $\mu_X + \mu_Y$ and variance $\sigma^2_X + \sigma_Y^2$.
\end{proof}
Consider a variable $X$ that follows a Markov process. If we know that the change in value of a single day is a normal distribution with mean zero and variance one, then we may find the distribution for two days. Due to the Markov property the above proposition applies, and we have that the change in the variable over two days is normally distributed with mean zero and variance two.
Note that we may apply this as many times as we please to find the change in the variable over any period but the variance increases showing the uncertainty of these predictions.
It follows that we may consider the change in the variable $W$ over a relatively large period of time $W$. We denote this $W_T - W_0$ and it can be thought of as the sum of changes in $W$ in $N$ small time intervals of length $\Delta t = T/N$. Hence we have that,
\begin{align}
W_T - W_0 = \sum^N_{i=1} \epsilon_i \sqrt{\Delta t}, \label{wp1}
\end{align}
where each of the $\epsilon_i$ for $i = 1,2,\dots,N$ are normally distributed with mean zero and variance one.
It follows from (\ref{wp1}) that, the mean of $W_T - W_0$ is zero and has variance $N\Delta t = T$. In normal calculus we consider the limit of a discrete process as the changes head toward zero. Similar notations and conventions exist in stochastic calculus, here we will use the notation $\mathrm{d}W$ to refer to the Wiener process $W_{\Delta t + t} - W_t$ in the limit $\Delta t \rightarrow 0$.
%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
\subsection{Generalisation of the Wiener Process and It\^os Process}
Insofar we have discussed the Wiener process, however there is a problem with this. We can see this by defining the \emph{drift rate} and \emph{variance rate}, these are given as the mean change per unit time and variance change per unit time. In our prior discussion of the Wiener process it is obvious that these have been zero and one respectively. This leads to a small issue, the expected value of $W$ at any given time $t$, is equal to its current value. This is clearly an issue as stocks trend upward or downward. We may capture this aspect of stocks by generalising the Wiener process. The \emph{generalised Wiener process} for a variable, here denoted as $x$ is defined as,
\begin{align}
\mathrm{d}x = a \mathrm{d} t + b \mathrm{d} W, \label{gwp}
\end{align}
where $a$ and $b$ are constants and $\mathrm{d}W$ is the Wiener process. Here if the $\mathrm{d}W$ term where removed from (\ref{gwp}), then it would have solution $x = x_0 + at$ for some $x_0$ specified by an initial condition. This is the way we assume the stock will grow, moving at a constant rate $a$. The $\mathrm{d}W$ term adds ``noise'' to this. As stock movement is assumed to be random, this term introduces the randomness in the form of a Wiener process.
Expressing this in discrete terms we have that $\Delta x = a \Delta t + b\epsilon \sqrt{\Delta t}$, with $\epsilon$ as before. Hence we have that the mean is now $a\Delta t$ and variance is $b^2\Delta t$. Again following similar arguments as with the Wiener process, we see that for the continuous process that the mean and variance are now $aT$ and $b^2T$ respectively.
It is easy to see that even this model which allows us to drift the stock either up or down in a given direction is not particularly favourable. It assumes that the drift is constant. This is, in practice, is not the case, a stock may rise at a constant rate but crash the next day. Here we introduce an \emph{It\^o process} to compensate for this.
\begin{defn}[It\^o processes]
An It\^o process is a is a type of generalised Wiener process where the constants $a$ and $b$ are now dependent on the value of the underlying variable $x$ and time $t$. In mathematical terms,
\begin{align*}
\mathrm{d}x = a(x,t)\mathrm{d}t + b(x,t)\mathrm{d}W.
\end{align*}
\end{defn}
It\^o processes address the issues discussed. Following the same structure as previous arguments this has drift rate $a(x,t)$ and variance rate $b(x,t)^2$.
%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
\subsection{It\^os Lemma}
Before we may prove It\^os lemma we must first prove a very important result, often regarded as the corner stone of It\^o calculus, $\mathrm{d}W^2_t = \mathrm{d}t$. This form of expressing this is merely short hand, the exact statement is given and proven below \citep[P. 87-88]{GSS}.
\begin{thrm}[$\mathrm{d}W_t^2 = \mathrm{d}t$]
Given a function $g(X_t)$ we have that,
\begin{align*}
\int^t_{t_0} g(X_t)\mathrm{d}W_t^2 = \int^t_{t_0}g(X_t)\mathrm{d}t.
\end{align*}
\end{thrm}
\begin{proof}
To begin first let us define a new variable $Y$ as,
\begin{align*}
Y &= \int^t_{t_0} g(X_t)\mathrm{d}W_t - \int^t_{t_0}g(X_t)\mathrm{d}t \\
&= \sum g(X_{i-1})\Delta W_{t_i}^2 - \sum g(X_{i-1})\Delta t_i \\
&= \sum g(X_{i-1})(\Delta W_{t_i}^2-\Delta t_i)
\end{align*}
We may then calculate the moments of this new variable.
\begin{align*}
\mathbb{E}\left[ Y\right] &= \lim_{n \rightarrow \infty}\sum^n_{i=1}\mathbb{E}\left[ g(X_{i-1})(\Delta W_{t_i}^2-\Delta t_i)\right] \\
&= \lim_{n \rightarrow \infty}\sum^n_{i=1}\mathbb{E}\left[ g(X_{i-1})\right]\mathbb{E}\left[(\Delta W_{t_i}^2-\Delta t_i)\right],
\end{align*}
as $g(X_{i-1})$ and $(\Delta W_{t_i}^2-\Delta t_i)$ are independent, so the mean of the product is the product of the mean. %%Should i include a proof of this?
Furthermore, as $\Delta t$ is deterministic and merely a constant, from the properties discussed in the preliminaries,
\begin{align*}
\mathbb{E}\left[ Y\right] &= \lim_{n \rightarrow \infty}\sum^n_{i=1}\mathbb{E}\left[ g(X_{i-1})\right](\mathbb{E}\left[\Delta W_{t_i}^2\right] -\Delta t_i).
\end{align*}
Now note that, as $\mathbb{E}\left[W_{t_i}W_{t_{i-1}}\right] = \mathbb{E}\left[(W_{t_i}-W_{t_{i-1}})W_{t_{i-1}} + W_{t_{i-1}}^2\right]$. Splitting this up we have that,
\begin{align*}
\mathbb{E}\left[(W_{t_i}-W_{t_{i-1}})W_{t_{i-1}} + W_{t_{i-1}}^2\right] = \mathbb{E}\left[(W_{t_i}-W_{t_{i-1}})W_{t_{i-1}}\right] + \mathbb{E}\left[W_{t_{i-1}}^2\right]
\end{align*}
The first term consists of two independent variables, thus,
\begin{align*}
\mathbb{E}\left[(W_{t_i}-W_{t_{i-1}})W_{t_{i-1}}\right] = \mathbb{E}\left[(W_{t_i}-W_{t_{i-1}})\right]\mathbb{E}\left[W_{t_{i-1}}\right] = 0,
\end{align*}
as the average of the Wiener process is zero. We then notice the second term may be expressed as,
\begin{align*}
\mathbb{E}\left[W_{t_{i-1}}^2\right] &= \mathbb{E}\left[W_{t_{i-1}}\right]^2 + \mathbb{E}\left[W_{t_{i-1}}^2\right] - 2\mathbb{E}\left[W_{t_{i-1}}\right]^2 + \mathbb{E}\left[W_{t_{i-1}}\right]^2 \\
&= \operatorname{Var}[W_{t_{i-1}}] + \mathbb{E}\left[W_{t_{i-1}}\right]^2.
\end{align*}
The average and variance of the Wiener process are known. Hence, we know that $\mathbb{E}\left[W_{t_{i-1}}^2\right] = t_{i-1}$. Then we have,
\begin{align*}
\mathbb{E}\left[\Delta W_{t_i}^2\right] = \mathbb{E}\left[ (W_{t_i} - W_{t_{i-1}})^2\right] = \mathbb{E}\left[W_{t_{i}}^2\right] + \mathbb{E}\left[W_{t_{i-1}}^2\right] - 2\mathbb{E}\left[W_{t_i}W_{t_{i-1}}\right] = |t_i - t_{i-1}| = \Delta t_i,
\end{align*}
hence the average of our variable $Y$ is zero. Now note that,
\begin{align*}
\operatorname{Var}[Y] = \mathbb{E}\left[Y^2\right]-\mathbb{E}\left[Y\right]^2 = \mathbb{E}\left[ Y^2\right].
\end{align*}
So the second moment is precisely the variance. Considering the second moment we have that,
\begin{align*}
\mathbb{E}\left[ Y^2\right] &= \lim_{n \rightarrow \infty}\mathbb{E}\left[\left( \sum^n_{i=1} g(X_{t_{i-1}})(\Delta W_{t_i}^2-\Delta t)\right)^2 \right] \\
&=\lim_{n\rightarrow\infty}\sum^n_{i=1}\mathbb{E}\left[g(X_{t_{i-1}})^2\right]\mathbb{E}\left[(\Delta W_{t_i} - \Delta t)^2\right] \\
&\qquad + 2 \sum^n_{i=1}\sum_{j<i}\mathbb{E}\left[\Delta W^2_{t_{i}}-\Delta t\right]\mathbb{E}\left[g(X_{t_{i-1}})g(X_{t_{j-1}})(\Delta W_{t_j}^2-\Delta t) \right].
\end{align*}
Note that in the second term in the above we have that $\mathbb{E}\left[\Delta W_{t_i}^2 -\Delta t \right]$ is independent of all the other terms, so we may split up the angular brackets. However $\mathbb{E}\left[\Delta W_{t_i}^2 -\Delta t \right] = 0$ hence there is no contribution from the second term. Now consider $\mathbb{E}\left[(\Delta W_{t_i} - \Delta t)^2\right]$,
\begin{align*}
\mathbb{E}\left[(\Delta W_{t_i} - \Delta t)^2\right] &= \mathbb{E}\left[W^4_{t_i}\right] - 2\Delta t\mathbb{E}\left[W^2_{t_i}\right] + \Delta t^2\\
&= 3\Delta t ^2 - 2 \Delta t^2 + \Delta t^2 = 2\Delta t^2.
\end{align*}
Here we used that as $\Delta W$ is a Gaussian variable with zero mean this means that $\mathbb{E}\left[W^4\right] = 3\mathbb{E}\left[W^2\right]$, this can be seen by direct integration of the Gaussian distribution and is quoted but not derived here. Hence, we now have,
\begin{align*}
\mathbb{E}\left[ Y^2\right] &= \lim_{n \rightarrow\infty}\sum_{i=1}^n \mathbb{E}\left[g(X_{t_{i-1}})^2 2\Delta t^2\right] \\
&= \lim_{n \rightarrow\infty} 2\Delta t\sum_{i=1}^n \mathbb{E}\left[g(X_{t_{i-1}})^2 2\Delta t\right].
\end{align*}
Note that taking the limit $n \rightarrow \infty$ is equivalent to taking $\Delta t \rightarrow 0$, as $\Delta t = (t-t_0)/n$. Hence, $\mathbb{E}\left[Y^2\right] =0$. We have shown $Y$ has mean and variance that are zero. The argument also holds for higher moments as $\Delta t \rightarrow 0$ and these can be shown to be zero as well. We now have a variable with all moments equal to zero and hence $Y \equiv 0$. This proves our result.
\end{proof}
With this corner stone of stochastic calculus we may now prove It\^o's lemma, sometimes referred to as It\^o's equation or It\^o's differentiation rule. Proving this will then allow us to derive the Black-Scholes formula. We follow the method of \citet[P. 95]{GSS}.
\begin{lem}[It\^os Lemma]
\label{itol}
Consider an It\^o process described by,
\begin{align}
\mathrm{d}X_t = A(t,X_t)\mathrm{d}t + B(t,X_t)\mathrm{d}W_t, \label{SDE}
\end{align}
where $W_t$ is the Wiener process. Then, if $g(t,x)$ is a twice-differentiable scalar function of two variables $x,t \in \mathbb{R}$, then,
\begin{align*}
\mathrm{d}g(t,X_t) = \left(\frac{\partial g}{\partial t} + \frac{\partial g}{\partial X}A + \frac{1}{2}\frac{\partial^2 g}{\partial X^2}B^2 \right)\mathrm{d}t + \frac{\partial g}{\partial X}B\mathrm{d}W.
\end{align*}
\end{lem}
\begin{proof}
We do not give a full rigorous proof here as it is beyond the scope of this project, however we may derive this result using results from Riemann calculus. Consider $g(t,X_t)$, then from Taylor's Theorem we know that an approximation for the derivative of $g$of order $\mathcal{O}(dt)$ is given by,
\begin{align*}
\mathrm{d}g(t,X_t) = \frac{\partial g}{\partial t} \mathrm{d}t + \frac{\partial g}{\partial X}\mathrm{d}X_t + \frac{1}{2}\frac{\partial^2 g}{\partial X^2}\mathrm{d}X_t^2.
\end{align*}
We include the final term as when we substitute in (\ref{SDE}) we will obtain a $\mathrm{d}W^2$ term which we know to be $dt$. Substituting (\ref{SDE}) into the above we obtain,
\begin{align*}
\mathrm{d}g(t,X_t) &= \frac{\partial g}{\partial t} \mathrm{d}t + \frac{\partial g}{\partial X}(A(t,X_t)\mathrm{d}t + B(t,X_t)\mathrm{d}W_t) \\
&\qquad + \frac{1}{2}\frac{\partial^2 g}{\partial X^2}(A(t,X_t)\mathrm{d}t + B(t,X_t)\mathrm{d}W_t)^2.
\end{align*}
Recall that $\left<\Delta W_t^2\right> = \Delta t$. Now taking limits we see that $\left<\mathrm{d} W_t^2\right> = \mathrm{d} t$, so $\mathrm{d}W_t$ can be thought of as $\mathcal{O}(\mathrm{d}t)$. Expanding and removing terms using this rule of order greater than $\mathrm{d}t$ we have that,
\begin{align*}
\mathrm{d}g(t,X_t) &= \frac{\partial g}{\partial t} \mathrm{d}t + \frac{\partial g}{\partial X}(A(t,X_t)\mathrm{d}t + B(t,X_t)\mathrm{d}W_t) + \frac{1}{2}\frac{\partial^2 g}{\partial X^2}\left[A(t,X_t)^2\mathrm{d}t^2\right. \\
&\qquad\left.+ B(t,X_t)^2\mathrm{d}W_t^2 + 2A(t,X_t)B(t,X_t)\mathrm{d}t\mathrm{d}W_t^2\right] \\
&= \frac{\partial g}{\partial t} \mathrm{d}t + \frac{\partial g}{\partial X}(A(t,X_t)\mathrm{d}t + B(t,X_t)\mathrm{d}W_t) + \frac{1}{2}\frac{\partial^2 g}{\partial X^2}\left[B(t,X_t)^2\mathrm{d}W_t^2\right].
\end{align*}
Replacing $\mathrm{d}W_t^2$ with $\mathrm{d}t$ we have,
\begin{align*}
\mathrm{d}g(t,X_t) = \left(\frac{\partial g}{\partial t} + \frac{\partial g}{\partial X}A(t,X_t) + \frac{1}{2}\frac{\partial^2 g}{\partial X^2}B(t,X_t)^2 \right) \mathrm{d}t + \frac{\partial g}{\partial X}B(t,X_t)\mathrm{d}W_t,
\end{align*}
as required.
\end{proof}
%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
\subsection{Black-Scholes Model}
As discussed before %%in section generalised weiner process. Refference thsi
it is incorrect to assume that the stock follows a generalized Wiener process. We then introduced an It\^o process to counter this. It remains to be determined what the functions $a(x,t)$ and $b(x,t)$ need to be in our It\^o process. The appropriate assumption is that the expected return (the drift over the stock price) is constant. Thus we must have that the expected return is $\mu S$ for some $\mu \in \mathbb{R}$ \citep[Chap. 12]{JCH}.
Furthermore, we assume that the percentage return's variability over a small time interval, $\Delta t$, is constant and independent of stock price. This means that a buyer is as uncertain of the return (as a percentage) when the stock costs $\$1$ as when the stock costs $\$1000$. This leads to the fact that the the stock price should be proportional to the standard deviation over a small period of time $\Delta t$. This leads to the following final model,
\begin{align}
\mathrm{d}S = \mu S\mathrm{d}t + \sigma S \mathrm{d}W, \label{model}
\end{align}
where the variable $\sigma$ is the volatility of the stock per year and $\mu$ is the expected rate of return on the stock per year. This is the most widely used model for stock behavior. Now using (\ref{model}) and applying It\^o's lemma we obtain, for a function $G(S,t)$ we have the process $G$ follows is given by,
\begin{align}
\mathrm{d}G = \left( \frac{\partial G}{\partial S} + \frac{\partial G}{\partial t} +\frac{1}{2} \sigma^2 S^2 \frac{\partial^2 G}{\partial S^2}\right)\mathrm{d} t + \frac{\partial G}{\partial S} \sigma S \mathrm{d} W.\label{BSM1}
\end{align}
We see from (\ref{model}) that the volatility is a measure of how unsure we are about the path the stock will take. This is because it is multiplying the random component. It can also be viewed as the standard deviation of the lognormal distribution of $S_T$, as we will see in the next section.
%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
\subsection{The Lognormal Property}
\label{lnp}
Consider the equation as described in (\ref{model}). We may use It\^os lemma to derive the process followed by $\log(S)$ \citep[Chap. 13]{JCH}. Let $G= \log(S)$ by applying It\^os lemma, (\ref{itol}), with $X_t =S$ and $g(x,t) = \log(S)$ we obtain,
\begin{align*}
\mathrm{d}G = \left( \mu - \frac{\sigma^2}{2}\right)\mathrm{d}t + \sigma\mathrm{d}W.
\end{align*}
With $\mu$ as, the expected rate of return, and $\sigma$ being the volatility of the stock, these are constant. Thus $G$ follows a generalised Wiener process, with drift rate $\mu - \frac{\sigma^2}{2}$ and variance rate $\sigma$. Therefore, by Proposition \ref{prop1}, between time $0$ and $T$ we see that $G$ has mean $(\mu - \frac{\sigma^2}{2})T$ and variance $\sigma^2 T$. Hence,
\begin{align*}
\log(S_T) - \log(S_0) \sim \phi\left[\left(\mu - \frac{\sigma^2}{2}\right)T, \sigma \sqrt{T} \right] \\
\Rightarrow\log(S_T) \sim \phi\left[ \log(S_0) + \left(\mu - \frac{\sigma^2}{2}\right)T, \sigma \sqrt{T} \right].
\end{align*}
Where $S_t$ is the stock price at time $t$. Then by Definition \ref{lnd} we have that,
\begin{align*}
\mathbb{E}(S_T) &= S_0 \mathrm{e}^{\mu T} \\
\textnormal{Var}(S_T) &= S_0^2\mathrm{e}^{2\mu T}\left(\mathrm{e}^{\sigma^2 T} - 1 \right).
\end{align*}
%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^6
\subsection{The Black-Scholes Differential Equation}
Using this model for stock prices, we may derive the Blakc-Scholes differential equation \citep[Chap. 13]{JCH}. Given that $f$ is an option subject to $S$, then $f$ must be some function of $S$ and $t$. Hence, from (\ref{BSM1}),
\begin{align*}
\mathrm{d}f = \left( \frac{\partial f}{\partial S} + \frac{\partial f}{\partial t} +\frac{1}{2} \sigma^2 S^2 \frac{\partial^2 f}{\partial S^2}\right)\mathrm{d} t + \frac{\partial f}{\partial S} \sigma S \mathrm{d} W.
\end{align*}
The equations (\ref{model}) and the above have discretized versions,
\begin{align}
\Delta S &= \mu S\Delta t + \sigma S \Delta W. \label{Dmodel} \\
\Delta f &= \left( \frac{\partial f}{\partial S} + \frac{\partial f}{\partial t} +\frac{1}{2} \sigma^2 S^2 \frac{\partial^2 f}{\partial S^2}\right)\Delta t + \frac{\partial f}{\partial S} \sigma S \Delta W, \label{DBSM}
\end{align}
over a time interval $\Delta t$. As the Wiener processes contained in $\Delta f$ and $\Delta S$ are the same, it follows that we may construct a portfolio to eliminate this. Such a portfolio should sell an option and buy $\frac{\partial f}{\partial S}$ shares. Then by definition our portfolio, $\Pi$, is,
\begin{align}
\Pi = -f + \frac{\partial f}{\partial S} S. \label{BSE1}
\end{align}
The change in this over $\Delta t$ is,
\begin{align*}
\Delta\Pi = -\Delta f + \frac{\partial f}{\partial S}\Delta S.
\end{align*}
By substituting in (\ref{Dmodel}) and (\ref{DBSM}) we obtain,
\begin{align}
\Delta \Pi = \left(\frac{\partial f}{\partial t} - \frac{1}{2}\frac{\partial^2 f}{\partial S^2}\sigma^2 S^2\right)\Delta t \label{BSE}
\end{align}
Over the time period $\Delta t$ we have eliminated $\Delta W$, so the portfolio must be riskless in this time period and must therefore make the riskfree interest rate. Thus,
\begin{align*}
\Delta \Pi = r\Pi \Delta t
\end{align*}
substituting (\ref{BSE1}) and (\ref{BSE}) into the above, we yield,
\begin{align*}
\left(\frac{\partial f}{\partial t} + \frac{1}{2}\sigma^2 S^2\frac{\partial^2 f}{\partial S^2} \right)\Delta t = r\left( f - \frac{\partial f}{\partial S}S\right)\Delta t
\end{align*}
so that,
\begin{align*}
\frac{\partial f}{\partial t} + rS\frac{\partial f}{\partial S} + \frac{1}{2}\sigma^2 S^2\frac{\partial^2 f}{\partial S^2} = rf. %\label{BS}
\end{align*}
The above is known as the Black-Scholes differential equation. It is solvable for some boundary conditions and unsolvable analytically for others. In particular this is solvable for European options, which we investigate in the next section.
%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
\subsection{Black-Scholes Formula for European Options}
\label{bsedddd}
\begin{thrm}[Black-Scholes Pricing Formula]
The value of a call option, $f_c$, and a put option, $f_p$, are given by the following formulae;
\begin{align*}
f_c &= S_0\Phi(d_1) - K \mathrm{e}^{-rT}\Phi(d_2), \\
f_p &= K\mathrm{e}^{-rT}\Phi(-d_2) - S_0 \Phi(-d_1),
\end{align*}
where
\begin{align}
d_1 &= \frac{\log(S_0/K)+(r+\sigma^2/2)T}{\sigma \sqrt{T}} \label{bsdii1}, \\
d_2 &= \frac{\log(S_0/K)+(r-\sigma^2/2)T}{\sigma \sqrt{T}} = d_1 - \sigma. \sqrt{T}.
\end{align}
\end{thrm}
Before proving this we prove the following claim \citep[P. 310-312]{JCH}.
\begin{claim*}
If $V$ is lognormally distributed and the standard deviation of $\log(V)$ is $\omega$, then,
\begin{align*}
\mathbb{E}(\max(V - K,0)) = \mathbb{E}(V)N(d_1) - KN(d_2),
\end{align*}
where $\mathbb{E}$ denotes the expected value, and we have that,
\begin{align*}
d_1 = \frac{\log\left(\frac{\mathbb{E}(V)}{K}\right) + \frac{\omega^2}{2}}{\omega} \\
d_2 = \frac{\log\left(\frac{\mathbb{E}(V)}{K}\right) - \frac{\omega^2}{2}}{\omega}.
\end{align*}
\end{claim*}
\begin{proof}[Proof of claim]
Define $h(V)$ to be the probability density function of $V$. Then we must have that,
\begin{align}
\mathbb{E}(\max(V-K,0)) &= \int^{\infty}_0 \max(V-K,0)h(V)\mathrm{d}V \\
&=\int^{\infty}_K (V-K)h(V)\mathrm{d}V \label{pc1}.
\end{align}
By assumption the variable $\log(V)$ is normally distributed with standard deviation $\omega$. Then from Definition \ref{lnd} we have that the mean, $m$, is given by,
\begin{align*}
m = \log(\mathbb{E}(V))- \frac{\omega^2}{2}.
\end{align*}
We further define a new variable, $W$, by the following,
\begin{align}
W = \frac{\log(V)-m}{\omega}. \label{trans}
\end{align}
This is the transformation that turns the distribution of $\log(V)$ to the standard normal distribution. Let the probability distribution function of $W$ be $g(W)$, so that,
\begin{align*}
g(W) = \frac{1}{\sqrt{2\pi}}\mathrm{e}^{\frac{-W^2}{2}}.
\end{align*}
Using (\ref{trans}) as a change of variable for (\ref{pc1}) we obtain that,
\begin{align}
\mathbb{E}(\max(V-K,0))&= \int^\infty_{\frac{\log(K)-m}{\omega}} \left(\mathrm{e}^{Q\omega + m}-K\right)h(Q)\mathrm{d}Q \\
&= \int^\infty_{\frac{\log(K)-m}{\omega}} \mathrm{e}^{Q\omega + m}h(Q)\mathrm{d}x - K\int^\infty_{\frac{\log(K)-m}{\omega}}h(Q)\mathrm{d}Q. \label{cp2}
\end{align}
To solve this first consider,
\begin{align*}
\mathrm{e}^{Q\omega +m }h(Q) &= \frac{1}{\sqrt{2\pi}}\exp\left(-\frac{(Q-\omega)^2 + 2m + \omega^2}{2}\right) \\
&=\mathrm{e}^{m + \frac{\omega^2}{2}}h(Q-\omega),
\end{align*}
thus we now have that (\ref{cp2}) is,
\begin{align}
\mathbb{E}(\max(V-K,0))&= \mathrm{e}^{m + \frac{\omega^2}{2}}\int^\infty_{\frac{\log(K)-m}{\omega}}h(Q-\omega)\mathrm{d}Q - K\int^\infty_{\frac{\log(K)-m}{\omega}}h(Q)\mathrm{d}Q. \label{inteqn}
\end{align}
In (\ref{inteqn}) the first integral, the integrand is a normal distribution with a shifted mean. Thus it must be a function of cumulative normal distribution as we are summing the area under the probability distribution function, giving that,
\begin{align*}
\int^\infty_{\frac{\log(K)-m}{\omega}}h(Q-\omega)\mathrm{d}Q &= \left[\Phi(Q-\omega)\right]^\infty_{\frac{\log(K)-m}{\omega}} \\
&= 1 - \Phi\left( \frac{\log(K)-m}{\omega} - \omega\right) \\
&= \Phi\left( \frac{-\log(K) + m}{\omega} + \omega \right).
\end{align*}
Substituting for $m$ we obtain,
\begin{align*}
\Phi\left( \frac{-\log(K) + m}{\omega} + \omega \right) = \Phi\left(\frac{\log\left(\frac{\mathbb{E}(V)}{K}\right) + \frac{\omega^2}{2}}{\omega} \right) = \Phi(d_1).
\end{align*}
Similarly we obtain that the second integral in (\ref{inteqn}) is $\Phi(d_2)$. Thus we have that,
\begin{align*}
\mathbb{E}(\max(V-K,0)) = \mathrm{e}^{m+\frac{\omega^2}{2}}\Phi(d_1) -K\Phi(d_2).
\end{align*}
Substituting for $m$ and we obtain our required result.
\end{proof}
\begin{proof}[Proof of Black-Scholes Formula]
We will prove this for a call option. The proof for a put follows similarly using a similarly proved claim and a similar strategy for this proof.
Consider a call option on a non-dividend paying stock with time of expiry $T$ (initial time $t=0$), strike price $K$, risk-free interest rate $r$, current stock price $S_0$ and volatility $\sigma$. The value of such a call at $T$ in a risk-neutral world would be,
\begin{align}
\mathbb{E}(\max(S_T-K,0)). \label{bsp1}
\end{align}
Note here this is the expected value in a risk neutral world and not necessarily the real world. For the rest of this proof all expected values will be as such.
Then as this is its value at time $T$, its value, $c$, would be (\ref{bsp1}) discounted to the initial time,
\begin{align*}
c = \mathrm{e}^{-rT}\mathbb{E}(\max(S_T-K,0)).
\end{align*}
Then as we have shown the stochastic process underlying the stock is log normally distributed. Then at time $t=T$, we have $S_T$ is log normally distributed, and from Section \ref{lnp} $\mathbb{E}(S_T) = S_0 \mathrm{e}^{rS_T}$, the standard deviation of $\log(S_T)$ is $\sigma\sqrt{T}$. Using our claim,
\begin{align*}
c &= \mathrm{e}^{-rT}\left(S_0 \mathrm{e}^{rT}\Phi(d_1) - \Phi(d_2) \right) \\
&= S_0\Phi(d_1) - K\mathrm{e}^{-rT}\Phi(d_2).
\end{align*}
Further in this case,
\begin{align*}
d_1 = \frac{\log\left(\frac{\mathbb{E}(S_T)}{K}\right) + \frac{\sigma^2T}{2}}{\sigma\sqrt{T}} = \frac{\log(S_0/K)+(r+\sigma^2/2)T}{\sigma \sqrt{T}} \\
d_2 = \frac{\log\left(\frac{\mathbb{E}(S_T)}{K}\right) - \frac{\sigma^2T}{2}}{\sigma\sqrt{T}} = \frac{\log(S_0/K)+(r-\sigma^2/2)T}{\sigma \sqrt{T}}.
\end{align*}
\end{proof}
%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
\subsection{Black-Scholes Model for Other Options}
In the previous section we found an analytic solution to the Black-Scholes differential equation for European options. Due to the inequality constraints for American and averaging for Asian options, it is not possible to solve the Black-Scholes equation analytically for these in general. There does exist an analytic formula for American options if there is only one dividend to pay known as Roll-Geske-Whaley model and for some types of Asian options however these are specialist and shall not be discussed here.
%http://www.m-hikari.com/ijma/ijma-2011/ijma-25-28-2011/elshegmaniIJMA25-28-2011.pdf
%Int. Journal of Math. Analysis, Vol. 5, 2011, no. 26, 1259 - 1265
%Analytical Solution for an Arithmetic Asian Option Using Mellin Transforms
%Zieneb Ali Elshegmani and Rokiah Rozita Ahmed
%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
\subsection{Black-Scholes for Dividend Paying Stocks}
%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
On the date a dividend is paid the stock declines by the amount of the dividend \citep[Chap. 13]{JCH}. Therefore we have that our expected rate of return is now $\mu = r-q$ and the model in (\ref{model}) still holds. The derivation follows similarly and we yield the following equation,
\begin{align}
\frac{\partial f}{\partial t} + (r-q)S\frac{\partial f}{\partial S} + \frac{1}{2}\sigma^2 S^2\frac{\partial^2 f}{\partial S^2} = rf. \label{bsd}
\end{align}
Again this may be solved for European options and the previous statements for Asian and American options are still true, so we may not solve this analytically for these. To solve this notice that the previous probabilistic analysis as seen in Section \ref{bsedddd} still holds if we take $\mu = r-q$ instead of $\mu = r$. So we see that as all of our formulae still hold, we merely replace $r$ with $r-q$. Thus we have the following theorem,
%%% We gonna do this through solving the equation for completeness.
%##################################################################################
%\subsubsection{Reduction to Heat 5Equation With Constant Coefficients}
%We begin with the first transformation, let $u = \log\left(\frac{S}{K}\right)$. Then it is apparent that $\frac{\partial u}{\partial S} = \frac{1}{S}$. Then define $\hat{f}(u,t)= f(S,t)$. We have then that the derivatives are,
%\begin{align*}
%\frac{\partial \hat{f}}{\partial S} &= \frac{1}{S}\frac{\partial \hat{f}}{\partial u} \\
%\frac{\partial^2 \hat{f}}{\partial S^2} &= \frac{1}{S^2}\left(\frac{\partial^2 \hat{f}}{\partial u^2} - \frac{\partial \hat{f}}{\partial u} \right).
%\end{align*}
%When we substitute back into (\ref{bsd}) we obtain,
%\begin{align*}
%\frac{\partial \hat{f}}{\partial t} + (r-q - \frac{1}{2}\sigma^2)\frac{\partial \hat{f}}{\partial S} + \frac{1}{2}\sigma^2\frac{\partial^2 \hat{f}}{\partial S^2} = r\hat{f}.
%\end{align*}
%Noting here that $\hat{f}$ grows exponentially due to the $r\hat{f}$ term we need to remove this. We write our solution in the form,
%\begin{align}
%\hat{f}(u,t) = \mathrm{e}^{-r(T-t)}y(u,t). \label{asump}
%\end{align}
%This has the advantage of normalizing our boundary condition at $t=T$. Substituting into our transformed Black-Scholes equation we obtain,
%\begin{align*}
%\frac{\partial y}{\partial t} + (r-q - \frac{1}{2}\sigma^2)\frac{\partial y}{\partial S} + \frac{1}{2}\sigma^2\frac{\partial^2 y}{\partial S^2} = 0.
%\end{align*}
%Now to get the above into canonical form we perform our penultimate last change of variables. Let,
%\begin{align}
%u' &= u\frac{r-q-\frac{1}{2}\sigma}{\frac{\sigma^2}{2}} \label{utrans}\\
%t' &= (T-t)\frac{(r-q-\frac{1}{2}\sigma^2)^2}{\frac{\sigma^2}{2}} \label{ttrans}
%\end{align}
%Then if we let $y'(u',t') = y(u,t)$ we have
%\begin{align}
%\frac{\partial y'}{\partial t'} = \frac{\partial y'}{\partial u'} + \frac{\partial^2 y'}{\partial u'^2}.
%\end{align}
%Finally to see the equivalence of the Black-Scholes equation and the heat equation let $z = u' + t'$ with $\hat{y}(z,t') = \hat{y}(u' + t',t') = y'(u,t)$ and we finally have the canonical form of the one dimensional heat equation,
%\begin{align}
%\frac{\partial \hat{y}}{\partial t'} = \frac{\partial^2 \hat{y}}{\partial z^2}
%\end{align}
%Further we have that the boundary conditions are, $\max(S-K,0)$ at time $t = T$. It is easily seen that this is now for $t'=0$ and $u'\geq 0$. Hence, we have that $\hat{y}(z,0) = x-K = K(\mathrm{e}^u -1)$ for $z \geq 0$ and $\hat{y}(z,0) = 0$ for $z<0$. We may now solve this to obtain the analytic formula for European options.
%\subsubsection{Solving the Heat Equation}
%The solution of the one dimensional heat equation is known. Using a Greens function it is known that the solution is given by,
%\begin{align}
%\hat{y}(z,t') = \int_{-\infty}^\infty \hat{y}(w,0) \frac{1}{\sqrt{4\pi t'}} \exp\left(-\frac{(z-w)^2}{4t'}\right)\mathrm{d}w. \label{dbs1}
%\end{align}
%where,
%\begin{align*}
%G(z-w;t') = \frac{1}{\sqrt{4\pi t'}} \exp\left(-\frac{(z-w)^2}{4t'}\right)
%\end{align*}
%is our Greens function. Now as $\hat{y}(w,0)$ is zero for all negative $w$ and is $K(\mathrm{e}^u -1)$ for zero or positive $z$ it follows that (\ref{dbs1}) becomes,
%\begin{align}
%\hat{y}(z,t') &= \int_{0}^\infty K\left(\mathrm{e}^u -1\right) \frac{1}{\sqrt{4\pi t'}} \exp\left(-\frac{(z-w)^2}{4t'}\right)\mathrm{d}w \\
%&= \int_{0}^\infty K\mathrm{e}^u \frac{1}{\sqrt{4\pi t'}} \exp\left(-\frac{(z-w)^2}{4t'}\right)\mathrm{d}w - K \int_{0}^\infty \frac{1}{\sqrt{4\pi t'}} \exp\left(-\frac{(z-w)^2}{4t'}\right) \mathrm{d}w \label{bsdd}
%\end{align}
%To solve this we first make the change of variables, $v = \frac{w-z}{\sqrt{2t'}}$. For simplicity let $ a = \frac{\frac{\sigma^2}{2}}{r-q-\frac{1}{2}\sigma}$ then $\mathrm{d}w = \sqrt{2t'}\mathrm{d}q$. Then (\ref{bsdd}) becomes,
%\begin{align*}
%\hat{y}(z,t') = \frac{1}{\sqrt{2\pi}}\int_{-d_2}^\infty K\exp\left(v\sqrt{2t'} + z\right) \exp\left(-\frac{v^2}{2}\right)\mathrm{d}v - K \frac{1}{\sqrt{2\pi}} \int_{-d_2}^\infty \exp\left(-\frac{v^2}{2}\right) \mathrm{d}v
%\end{align*}
%where $d_2$ will be determined later and $v\sqrt{2t'} + z$ is gained from the formula for $u$ and using that at $t=T$ we have $t'=0$ hence, $z = u'$. As this is evaluated at $w$ the exponent is $aw$, rearranging for $v$ we obtain the above.
%We separate these integrals here for simplicity. Note that we have essentially solved the rightmost integral as the integrand is of the form of a normal distribution. Thus the integral is the cumulative distribution and we have that this integral is $K\Phi(d_2)$.
%For the left most integral note that. We may complete the square for the exponents. This yields,
%\begin{align*}
%\frac{1}{\sqrt{2\pi}}\int_{-d_2}^\infty K\exp\left(-\frac{(q-a\sqrt{2t'})}{2} + az + a^2t'\right)\mathrm{d}v.
%\end{align*}
%This inspires a final change of variables, define $v' = q - \sqrt{2t'}a$. Then we have the leftmost term is,
%\begin{align*}
%\exp\left( az + a^2t'\right)\frac{1}{\sqrt{2\pi}}\int_{-d_1}^\infty K\exp\left(-\frac{v'}{2}\right)\mathrm{d}v'
%\end{align*}
%where again $d_1$ will be defined later. Now the integrand is a standadised normal distribution and further our term is now $\exp\left( az + a^2t'\right)K \Phi(d_1)$. Finally we may see that, using the definitions of our changes of variables we may show that $az + a^2 t' = \log(\frac{S}{K}) +(r-q)(T-t)$.
%Thus our solution is,
%\begin{align*}
%\hat{y}(z,t') = K \left( \exp\left(\frac{S}{K} + (r-q)(T-t)\right)\Phi(d_1) - \Phi(d_2) \right),
%\end{align*}
%using (\ref{asump}) we have our final solution is,
%\begin{align*}
%f(S,t) &= K \exp\left(-q(T-t)\right)\Phi(d_1) - K\exp\left( -r(T-t)\right)
%&= K \exp\left(-qT\right)\Phi(d_1) - K\exp\left( -rT\right),
%\end{align*}
%assuming $t=0$. Further through substituting backwards again we obtain that,
%\begin{align*}
%d_1 = \frac{\log\left(\frac{\mathbb{E}(S_T)}{K}\right) - \frac{\sigma^2T}{2}}{\sigma\sqrt{T}} = \frac{\log(S_0/K)+(r-q+\sigma^2/2)T}{\sigma \sqrt{T}}\\
%d_2 = \frac{\log\left(\frac{\mathbb{E}(S_T)}{K}\right) - \frac{\sigma^2T}{2}}{\sigma\sqrt{T}} = \frac{\log(S_0/K)+(r-q-\sigma^2/2)T}{\sigma \sqrt{T}}
%\end{align*}
%The solution for a put is similar. We have proven then,
%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
\begin{thrm}[Black-Scholes Pricing Formula with Dividends]
The value of a call option, $f_c$, and a put option, $f_p$, are given by the following formulae.
\begin{align*}
f_c &= S_0\mathrm{e}^{-qT}\Phi(d_1) - K \mathrm{e}^{-rT}\Phi(d_2), \\
f_p &= K\mathrm{e}^{-rT}\Phi(-d_2) - S_0\mathrm{e}^{-qt} \Phi(-d_1),
\end{align*}
where
\begin{align*}
d_1 &= \frac{\log(S_0/K)+(r-q+\sigma^2/2)T}{\sigma \sqrt{T}}, \\
d_2 &= \frac{\log(S_0/K)+(r-q-\sigma^2/2)T}{\sigma \sqrt{T}} = d_1 - \sigma \sqrt{T}.
\end{align*}
\end{thrm}
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
\section{Monte Carlo Simulation}
%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
\subsection{Monte Carlo Concept}
\label{mcs1}
Monte Carlo is the first numerical method we will look at for finding the price of options. The method is incredibly simple and relies on some of the ideas we built in the previous section. We shall develop these ideas following the same ideas as \citet[Chap. 1]{PG}. The concept is as follows;
\begin{enumerate}
\item Attempt to predict the path of the price of the underlying stock from $t=0$ to $t= T$ where $T$ is the expiry date of the option, \label{a}
\item Evaluate the stock price at $t=T$,
\item Using the previous step, calculate the option price at maturity ($t=T$), \label{c}
\item Repeat (\ref{a}) - (\ref{c}) a statistically significant number of times, \label{d}
\item Calculate the average option price at maturity,
\item Discount the average option price by the interest rate to obtain the option price at $t=0$.
\end{enumerate}
For the first step we use the model for stock behavior we developed in (\ref{model}). We may solve this with an application of It\^os lemma along with the $\mathrm{d}W_t^2 = \mathrm{d}t $ formula. Firstly note that by It\^os lemma we have that,
\begin{align}
\mathrm{d}(\ln(S_t)) = \frac{1}{S_t}\mathrm{d}S_t - \frac{1}{2S_t^2}\mathrm{d}S_t^2 \label{mc1}.
\end{align}
Then by using (\ref{model}) with (\ref{mc1}) and substituting in for $\mathrm{d}S_t$ we have that,
\begin{align*}
\mathrm{d}(\ln(S_t)) &= \frac{1}{S_t}S_t(\mu \mathrm{d}t + \sigma W_t) - \frac{1}{2S_t^2}S_t^2(\sigma^2 \mathrm{d}W_t^2) \\
&= \mu\mathrm{d}t + \sigma \mathrm{d}W_t - \frac{1}{2}\sigma^2\mathrm{d}t.
\end{align*}
Then through exponentiation we have that,
\begin{align}
S_t = S_0 \left[\exp \left( (\mu - \frac{1}{2}\sigma^2) t + \sigma W_t\right) \right] \label{gbms},
\end{align}
where $S_0$ is our initial stock price as before. Now note that as the Wiener process is normally distributed with mean zero and standard deviation $\sqrt{T}$ we may rewrite (\ref{gbms}) as
\begin{align}
S_t = S_0 \left[\exp \left( (\mu - \frac{1}{2}\sigma^2) t + \sigma\sqrt{T} N\right) \right] \label{gbms1},
\end{align}
where $N$ is a standardized normal random variable. Note here that by taking logarithms of the above we may again see that $\ln(S_t)$ is normally distributed.
Using (\ref{gbms1}) we may perform our first step by sampling for $N$ to generate our path. This works well for European options however for Asian options as we are considering the average this may not generate a realistic average.
To compensate for this we split our interval into many different time points $0 = t_0 < t_1 < \dots < t_{n-1}< t_n = T$. We sample at each of these time points to generate a path between $[t_i, t_{i+1}] $ for $i = 0,\dots,n-1$. This allows us to generate more realistic paths.
The next step in the algorithm that was described above is very simple. However the third step changes and varies depending on what kind of option we are considering so we shall see how we adapt these methods for European, American and Asian options.
%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
\subsection{Applications of Monte Carlo}
%########################################################################################################
\subsubsection{Monte Carlo for European Options}
\label{mce}
Monte Carlo is not a particularly useful method for European options. In light of the analytic formulas that we have there is no need to use Monte Carlo for European option. Further as \citet{PG} said ``Monte Carlo is not a competitive method for computing one dimensional integrals'' so a poor convergence rate is expected of Monte Carlo.
For European options the method is fairly simple. For the third step in our algorithm we need only use the payoff function to calculate the value of the option at the termination of each path. We then calculate the average price at the terminal date and discount to the initial time $t=0$. Using the methods we have discussed insofar our algorithm for European options is as follows;
\begin{enumerate}
\item Simulate a path for the interval by taking a random sample from $N$,\label{qqq}
\item Calculate the payoff at maturity $t=T$, using the payofff function for European options $\textnormal{Payoff}(S_{j,T})$, \label{qq}
\item Repeat (\ref{qqq}) - (\ref{qq}) a statically significant number of times to generate a large number, $M$, of different terminal points indexed by $j= 1, \dots, M$,
\item Calculate the average payoff, $\overline{S_T}$, at maturity,
\begin{align*}
\overline{S_T} = \frac{1}{M}\sum^{M}_{j=1}\textnormal{Payoff}(S_{j,T}),
\end{align*}
\item Discount by the risk-free interest rate to the initial time $t=0$ to find the value of the option $V$ such that
\begin{align*}
V = \mathrm{e}^{-rT}\overline{S_T}.
\end{align*}
\end{enumerate}
%########################################################################################################
\subsubsection{Monte Carlo for American Options}
It is possible to adapt Monte Carlo to American options however it is very difficult. The problem arises from the possibility of early exercise, therefore we would need to find an optimal exercise rule\citep[Chap. 8]{PG}. As we will develop other methods to price American options we will not consider this here.
%########################################################################################################
\subsubsection{Monte Carlo for Asian Options}
\label{fdfdfdf}
Monte Carlo is a method we may employ to price Asian options \citep[Chap. 1]{PG}. The difficulty with Asian options is in payoff function for these. For Asian options we have the following payoff function; let $\overline{S}$ be the average price of the option between the initial time and terminal time $[0,T]$. Then,
% Check this.
\begin{align}
\textnormal{Payoff} = \left\{
\begin{array}{ll}
\max{(\overline{S}-K,0)} & \text{Call} \\
\max{(K - \overline{S},0)} & \text{Put}.
\end{array}
\right. \label{mcpa}
\end{align}
The way that the average is calculated yields several different types of Asian options. We will focus on the discrete case where the average is calculated in the following way. Given a set of discrete dates where the price will be monitored at times $t_1, \dots, t_n$, and a strike price $K$, the average of the stock, $\overline{S}$ will be given by,
\begin{align}
\overline{S} = \frac{1}{n} \sum_{i=1}^{n}S_{t_i}. \label{mcava}
\end{align}
Note that $\overline{S}$ is dependent upon the value of $n$ we choose. We shall just use the notation $\overline{S}$ to mean $\overline{S}_n$. There are other ways such as a continuous average, priced through the methods developed by \citet{GY}. We will mainly focus on the above.
Assuming the same model for stock movement as in the European case and partitioning the interval $[0,T]$ into $n$ sub intervals as before we have that the discrete version of (\ref{gbms1}) is given by,
\begin{align}
S_{t_{i+1}} = S_{t_{i}}\exp\left((\mu - \frac{1}{2}\sigma^2 )(t_i - t_{i-1})+ \sigma \sqrt{t_i - t_{i-1}}N_i\right) \label{mc2}
\end{align}
where $i = 1, \dots, n$ and $N_i$ denotes the $i$th sampling of $N$.
Using this our algorithm becomes;
\begin{enumerate}
\item For $i = 1,\dots, n$ generate a $N_i$ as our random component,
\item Use (\ref{mc2}) to calculate the stock price at the next point where the stock price is monitored,
\item Calculate the average using (\ref{mcava}),
\item Calculate the payoff, $C_j$, using (\ref{mcpa}),
\item Repeat the first three steps for $j = 1,\dots,m$ to generate $m$ different paths,
\item Find the average payoff from the option,
\begin{align*}
C = \frac{1}{m}\sum^m_{j=1}\textnormal{Payoff}(S_i).
\end{align*}
\item Discount this average to find the initial price of the option, so that option price is given by $\mathrm{e}^{-rT}C$.
\end{enumerate}
Using this method we may price Asian options \citep[P. 99]{PG}. There is another type of option known as the \emph{Geometric average option}. The only difference between the Asian and geometric average options is the method in which the average is calculated. For the geometric option, under all the previous assumptions, the average is given by,
\begin{align*}
\overline{S} = \left(\prod^n_{i=1} S(t_i) \right)^{\frac{1}{n}}.
\end{align*}
%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
\subsection{Optimization, Bias and Variance}
The natural questions that arise for these methods are that of efficiency. Here we will analyze two different methods of reducing our error. To begin we note that two different factors will cause our answer to be different from the true answer; these are obviously variance and bias. As such we need a ``measure'' to balance these. The appropriate tool for balancing these factors is the \emph{mean square error} \citep[P. 16]{PG}.
\begin{defn}
Let $\hat{x}$ be an estimator of $x$, the \emph{mean square error} of $\hat{x}$ is then,
\begin{align}
\textnormal{MSE}(\hat{x}) &= \mathbb{E}[(\hat{x}-x)^2] \\
&= \mathbb{E}[(\hat{x} - \mathbb{E}[\hat{x}])^2] + (\mathbb{E}[\hat{x}] - x)^2\label{qwewewe}\\
&= \textnormal{Var}(\hat{x}) + \textnormal{Bias}^2(\hat{x}).
\end{align}
In (\ref{qwewewe}) the right hand term on the right hand side is $\textnormal{Bias}^2(\hat{x})$, and the variance is given by the leftmost term on the right hand side.
\end{defn}
As \citet{PG} said ``In applications of Monte Carlo to financial engineering, estimator variance is typically larger than the bias'', we will focus on variance reduction techniques here on in to reduce the mean square error.
\subsubsection{Control Variate Techniques}
\label{mcs31}
Here the idea is to consider a similar option, with similar price, with a closed from formula for the price available \citep[P. 185-186]{PG}. Then we simulate the paths for both options using the same set of random variables. We compute the price from our Monte Carlo method and the analytic solution to calculate the error. Then the correct simulation for our option, without a closed-form pricing formula, would be our simulated value minus some weighting times the error involved.
Mathematically speaking, let $P_A$ be the price of an option and $P_B$ be the price of a different option with closed form solution. We simulate with random variables, $N_1,N_2,\dots,N_n$, to get two prices for the options from simulation say $P_A^{\textnormal{sim}}$ and $P_B^{\textnormal{sim}}$. Then we have our error is given by $E_B = P_B-P_B^{\textnormal{sim}}$, where $P_B$ is found from our closed form formula. Then we have that the correct price for $P_A$ under this simulation, $P^*_A$ would be $P^*_A = P_A^{\textnormal{sim}} - bE_B$.
So we seek the optimal $b$ for this. Consider $Y_1,\dots Y_n$ as $n$ replications of simulation trying to determine $\mathbb{E}[Y_i]$ with each $Y_i$ independent and identically distributed (i.i.d.). Suppose now at each $i$ we compute some $X_i$ with the same random variable, then from the method as above define $Y^*_i$ and $\overline{Y}^*$ to be,
\begin{align*}
Y_i^* &= Y_i - b (X_i - \mathbb{E}[X]),\\
\overline{Y}^* &= \frac{1}{n}\sum_{i=1}^n Y_i - b (X_i - \mathbb{E}[X].
\end{align*}
Then the control variate estimator above is unbiased as,
\begin{align*}
\mathbb{E}[\overline{Y}^*] = \mathbb{E}[Y_i - b (X_i - \mathbb{E}[X])] = \mathbb{E}[\overline{Y}] = \mathbb{E}[Y].
\end{align*}
This is consistent as it holds under the limit $n\rightarrow\infty$ in the definition of $\overline{Y}$. Now we may calculate the variance of each of the $Y^*_i$,
\begin{align}
\textnormal{Var}[Y^*_i] &= \textnormal{Var}[Y_i - b (X_i - \mathbb{E}[X])] \\
&= \sigma^2_Y - 2b\sigma_X\sigma_Y\rho_{X,Y} + b^2\sigma^2_X, \label{mcv1}
\end{align}
where $\sigma^2_X$ is the variance of $X$ and the analogous definition for $\sigma^2_Y$. It is easy to see that this technique reduces the variance if $2b\sigma_X\sigma_Y\rho_{X,Y} > b^2\sigma^2_X$. To optimize $b$ we use (\ref{mcv1}) and upon solving this we see that $b$ is optimized when,
\begin{align*}
b = \frac{\sigma^2_Y}{\sigma^2_X}\rho_{X,Y} = \frac{\textnormal{Covar[X,Y]}}{\textnormal{Var}[X]}.
\end{align*}
In practice it is not guaranteed that all the values required to calculate $b$ are present. As such we use our simulations to yield the estimate, $b^*$,
\begin{align}
b^* = \frac{\sum^n_{i=1}(X_i-\overline{X})(Y_i - \overline{Y})}{\sum^n_{i=1}(X_i-\overline{X})^2}. \label{bstar}
\end{align}
In the case of Asian options we have a closed form solution for the geometric average case but not the arithmetic average case. As such we may apply the above using $Y$ as our Asian option and $X$ as our geometric option. The close correlation of the price of these two types of assets was shown by \citet{BG}.
To see this formula for geometric average options \citep[P. 99-100]{PG} note that the product of lognormal variables is itself lognormal. Furthermore, the geometric average of lognormal variables is itself lognormal. It follows then from (\ref{gbms1}) that,
\begin{align*}
\left( \prod^n_{i=1}S(t_i)\right)^{1/n} = S_0 \exp\left( (r-\frac{1}{\sigma^2})\frac{1}{n}\sum^n_{i=1}t_i + \frac{\sigma}{n}\sum^n_{i=1}W_{t_i}\right).
\end{align*}
Following similar logic to Proposition \ref{prop1} we can see that,
\begin{align*}
\sum^n_{i=1}W_{t_i} \sim N(0,\sum^n_{i=1}(2i-1)t_{n+1-})
\end{align*}
as each $W_{t_i}$ is independent of the next due to the Markov property.
Thus it follows that at time $t=T$ the geometric average follows a process described by geometric Brownian motion with,
\begin{align*}
\mu &= r- \frac{1}{2}\sigma^2 + \frac{1}{2}\hat{\sigma}^2, \\
\sigma^2 &= \hat{\sigma}^2,
\end{align*}
where,
\begin{align*}
\hat{\sigma}^2 = \frac{\sigma^2}{n^2T}\sum^n_{i=1}(2i-1)t_{n+1-i}.
\end{align*}
Thus we may price these options using the Black-Scholes formula with these values of $\mu$ and $\sigma$. Now we may see how we can apply these for Asian options. We may use a geometric option as our option with a closed form solution to help reduce the variance when finding the price of an Asian option and we have the following algorithm to price the option,
\begin{enumerate}
\item Calculate the simulated price of the Asian and geometric option using the usual technique as given in Section \ref{mcs1}, using the same random variables,
\item Calculate the error in the geometric option by using the closed form solution we have for it,
\item Calculate $b^*$ as given in (\ref{bstar}),
\item Calculate the new option price using by subtracting the error multiplied by $b^*$ from the value of the Asian option.
\end{enumerate}
This will become useful and is how we will implement this method when investigating performance later on.
\subsubsection{Antithetic Variates}
\label{mcs32}
We will here explore another variance reduction technique, antithetic variates \citep[P. 205-207]{PG}. In this method for each $Y_i$ generated we create a corresponding $\hat{Y}_i$. Each pair $(Y_i, \hat{Y}_i)$ must be i.i.d.. Define the antithetic estimator as the average of these $2n$ replications, we then have the value of our estimator is of this is given by,
\begin{align*}
\overline{Y}^* &=\frac{1}{2n}\left(\sum^n_{i=1}Y_i + \sum_{i=1}^n\hat{Y}_i\right) \\
&= \frac{1}{n}\sum^n_{i=1}\left(\frac{Y_i + \hat{Y}_i}{2} \right).
\end{align*}
There are several tricks we may employ here to reduce computing time, one such trick is choosing $\hat{N}_i = - N_i$ for Gaussian variables (note that this is non-zero as we apply the payoff to obtain $Y_i$). These tricks are not always readily available thus in the worst case scenario to do this method we will require $2n$ replications. Thus we need to compare the new variance after the variance reduction technique to the variance of if we used a normal method with $2n$ replications. It follows that this method is reduces variance if,
\begin{align*}
\textnormal{Var}[\overline{Y}^*] < \textnormal{Var}\left[\frac{1}{2n}\sum^{2n}_{i=1}Y_i \right]
\end{align*}
thus the above implies that,
\begin{align}
\textnormal{Var}[Y_i + \hat{Y}_i] < 2\textnormal{Var}[Y_i]. \label{anv1}
\end{align}
We note that the left hand side can be written as,
\begin{align}
\textnormal{Var}[Y_i + \hat{Y}_i] &= \textnormal{Var}[Y_i] + \textnormal{Var}[\hat{Y}_i] + 2\textnormal{Covar}[Y_i,\hat{Y}_i] \\
&= 2\textnormal{Var}[Y_i] + 2\textnormal{Covar}[Y_i,\hat{Y}_i], \label{anv2}
\end{align}
as both $Y_i$ and $\hat{Y}_i$ have the same distribution and therefore must have the same variance. So using (\ref{anv1}) with (\ref{anv2}) we see that the condition for antithetic variance reduction to be effective is,
\begin{align*}
\textnormal{Covar}[Y_i,\hat{Y}_i] < 0,
\end{align*}
for each path.
Now we may apply this to the field of option pricing. To do this we will use the trick we discussed earlier. The algorithm is as follows,
\begin{enumerate}
\item Simulate the price to generate the paths in the usual way as given in Section \ref{mcs1} and calculate their payoff, \label{ytre}
\item Use the same random components and choose, for the new paths, that $N_i = -N_i$ where $N_i$ is the $i$th random variable we simulated in step \ref{ytre} and calculate the pay off of these new paths, \label{qwertt}
\item For each path generated in step \ref{ytre} average it with the new path generated in step \ref{qwertt},
\item Take the average of these paths and divide by $n$ to find the approximate value of the option.
\end{enumerate}
This is how the method will be implemented later on and we will explore its comparison to the control variate technique.
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
\section{Tree-based Models}
Tree based models are another way we may price options. Here we assume that the underlying asset follows a \emph{random walk}, i.e. in each time step we assume it may move in a number of directions, each with their own probability of occurring. We begin with the simplest model ,the binomial tree, we develop this to attain some pricing formulae for European and American options \citep[Chap. 11]{PG}.
%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
\subsection{Binomial Trees}
%%%%% NEED TO SORT OUT NOTATION, DON'T THINK ITS INTRODUCED WELL AND SPECIFY MOST OF THE FORMULAE DISCUSSED UP TO THE SECTION MOMENT MATCHIN IS ONLY FOR EUROPEAN OPTIONS.
The development of this model begins by assuming that there are only two directions in which the asset can go; up or down. This model was first introduced by \citet{CRR} and has the advantages of being easy to understand and yielding accurate prices, leading to its wide spread use. We will follow much the same derivation and discussion.
%########################################################################################################
\subsubsection{One-step Model}
\label{onestep}
This model is best introduced through an example over a single time step. Consider a stock that is currently trading at $\$20$, assume that over the period of a month we know that the price will be either $\$18$ or $\$24$. If we sell a European call option that expires at this date with strike price $\$21$ then the value of the option will be $\$3$ if the end price is $\$24$ or $\$0$ if it is $\$18$.
Then we may set up a riskless portfolio of stock. Consider if we sell one option an buy $x$ shares of stock, then the value of the stock will be $24 x$ if the stock ends at $\$24$ and the value of the portfolio will be $24 x - 3$. If the stock drops to $\$18$ then the value of the portfolio is $18x$. We choose $x$ so that the portfolio is riskless, i.e. so that in the event of both outcomes the value of the portfolio is the same. This gives $24x-3=18x \Rightarrow x = 0.5$, and we have that our riskless portfolio is selling one option and $0.5$ shares, with a value of $24*0.5 - 3 = 9$.
A riskless portfolio must earn the riskless interest rate, so if this is value the of the portfolio at the end of the month, we may discount this and obtain the value of the portfolio now. So assuming the risk free interest rate is $12\%$,
\begin{align*}
9*\mathrm{e}^{-\frac{1}{12} \times 0.12} = 8.910
\end{align*}
It follows that, with the price of the stock today being $\$20$, that the value of the portfolio today is $20 \times 0.5 - f$ where $f$ is the price of the option, and hence the price of the option is $\$11.09$. We see that we may use this concept to price options.
\begin{figure}[ht]
\centering
\begin{tikzpicture}[font=\small, scale=0.6]
\draw (0,5) -- (4,7);
\draw (0,5) -- (4,3);
\draw[fill] (0,5) circle [radius=0.07];
\draw[fill] (4,7) circle [radius=0.07];
\draw[fill] (4,3) circle [radius=0.07];
\node[left] at (0,5) {Stock price $S_0$};
\node[right] at (4,7) {Stock price increases to $S_0u$};
\node[right] at (4,3) {Stock price decreases to $S_0d$ };
\end{tikzpicture}
\caption{\label{fig:BTrees1}A one step binomial tree }
\centering
\end{figure}
We may generalise this algebraically. Given that the underlying may increase or decrease by values $u$ and $d$ respectively, where $u-1$ and $1-d$ are the percentage increase and decrease in the value respectively we assign probabilities to these of $q$ and $1-q$. Hence, if the current stock price is $S_0$, at the end of the period the value of the stock will be either $S_0u$ or $S_0d$ with corresponding option payoff value $f_u$ and $f_d$ respectively. This may be seen from Figure \ref{fig:BTrees1}. We also assume that the interest rate is constant. Then we have that the value of our portfolio is,
\begin{align*}
S_0ux - f_u \quad \text{ or } \quad S_0dx-f_d,
\end{align*}
for up and down movements respectively. Equating and solving for $x$ we obtain that $x = \frac{f_u - f_d}{S_0u - S_0d}$. Furthermore we know that the value of the portfolio now is $(S_0ux-f_u)\mathrm{e}^{-rT}$, where $T/1 = \Delta t$ is our time period and the cost of setting up the portfolio now is $S_0x - f$ with $f$ as the true option price. These must be equal at the initial time so equating and $f = S_0x(1-u\mathrm{e}^{-rT}) + f_u\mathrm{e}^{-rT}$. Substituting in the value for our $x$ we obtain,
\begin{align}
f&= S_0 \left( \frac{f_u-f_d}{S_0u - S_0d} \right) (1-u\mathrm{e}^{-rT}) + f_u\mathrm{e}^{-rT} \\
&= \frac{f_u(1-d\mathrm{e}^{-rT} + f_d (u\mathrm{e}^{-rT} -1))}{u-d} \\
&=\mathrm{e}^{-rT}(pf_u + (1-p)f_d), \label{bt1}
\end{align}
where $p = \frac{\mathrm{e}^{rT} - d}{u-d}$.
Here we see that $p$ is our probability $q$. Consider the expected value of $S$ at time $T$, denoted $\mathbb{E}[S_T]$ this is given by $\mathbb{E}[S_T] = qS_0u + (1-q)S_0d$. We know that the investment is constructed to be riskless, hence the expected value of the stock at time $T$ would be given by $\mathbb{E}[S_T] = S_0\mathrm{e}^{rT}$. Thus,
\begin{align}
&\Rightarrow S_0\mathrm{e}^{rT} = qS_0u + (1-q)S_0d \\
&\Rightarrow \mathrm{e}^{rT} = (u-d)q + d \\
&\Rightarrow q = \frac{\mathrm{e}^{rT} - d}{u-d}. \label{bt2}
\end{align}
So we see that $p$ is indeed our probability, referred to and denoted as such herein.
%########################################################################################################
\subsubsection{Two-step Model}
\label{twostep}
Here we consider a two step binomial tree this can bee seen in Figure \ref{fig:BTrees2}. Here our time step is now $\Delta t = \frac{T}{2}$, and thus we now have that equations (\ref{bt1}) and (\ref{bt2}) are now,
\begin{align}
f &=\mathrm{e}^{-r\Delta t}(pf_u + (1-p)f_d) \label{bt3} \\
p &= \frac{\mathrm{e}^{r\Delta t} - d}{u-d}.
\end{align}
Repeated application of (\ref{bt3}) yields,
\begin{align}
f_u = \mathrm{e}^{-r\Delta t}(pf_{uu} +(1-p)f_{ud}), \label{bt4}\\
f_d = \mathrm{e}^{-r\Delta t}(pf_{ud} + (1-p)f_{dd}). \label{bt5}
\end{align}
Then substituting equations (\ref{bt3}), (\ref{bt4}) and (\ref{bt5}) we obtain,
\begin{align}
f = \mathrm{e}^{-r2\Delta t}(p^2f_{uu} + 2p(1-p)f_{ud} + (1-p)^2f_{dd}). \label{bt6}
\end{align}
%########################################################################################################
\begin{figure}[ht]
\centering
\begin{tikzpicture}[font=\small, scale=0.6]
\draw (0,5) -- (4,7);
\draw (0,5) -- (4,3);
\draw (4,3) -- (8,1);
\draw (4,3) -- (8,5);
\draw (4,7) -- (8,9);
\draw (4,7) -- (8,5);
\draw[fill] (0,5) circle [radius=0.07];
\draw[fill] (4,7) circle [radius=0.07];
\draw[fill] (4,3) circle [radius=0.07];
\draw[fill] (8,1) circle [radius=0.07];
\draw[fill] (8,5) circle [radius=0.07];
\draw[fill] (8,9) circle [radius=0.07];
\node[below] at (0,5) {$S_0$};
\node[below] at (4,7) {$S_0u$};
\node[below] at (4,3) {$S_0d$};
\node[below] at (8,5) {$S_0ud$};
\node[below] at (8,9) {$S_0u^2$};
\node[below] at (8,1) {$S_0d^2$};
\end{tikzpicture}
\caption{\label{fig:BTrees2}A two step binomial tree }
\centering
\end{figure}
So we have again calculated the expected value of the tree, discounting to move from our finishing time $T$ to the initial time. Note that this formula only works for European options. This is as the above does not allow for early exercise. We here give an example to give some substance to this idea, this example is biased on one from \citet[P. 250]{PG}.
\begin{ex}
Given that the stock price is at $\$30$ and it may either increase or decrease by $10\%$, then for two, one month periods our tree can be seen in Figure \ref{fig:ex1}. Here we may combine the two middle nodes as $S_0ud = S_0 du$. Consider if we bought a European call option with a strike price of $\$31$. Then we see at all the nodes the option is out of the money, except the ones with values $\$33$ and $\$36.3$, here it has value $\$2$ and $\$5.3$ respectively. We seek to discount our option price to the initial node. We can either do this in two stages using equation (\ref{bt1}) or much more efficiently using equation (\ref{bt6}). Noting that here $u = 1.1, d= 0.9$ and assume that $r=12\%$ we may solve this using either of these methods which we give below,
\begin{figure}[b]
\centering
\begin{tikzpicture}[font=\small, scale=0.6]
\draw (0,5) -- (4,7);
\draw (0,5) -- (4,3);
\draw (4,3) -- (8,1);
\draw (4,3) -- (8,5);
\draw (4,7) -- (8,9);
\draw (4,7) -- (8,5);
\draw[fill] (0,5) circle [radius=0.07];
\draw[fill] (4,7) circle [radius=0.07];
\draw[fill] (4,3) circle [radius=0.07];
\draw[fill] (8,1) circle [radius=0.07];
\draw[fill] (8,5) circle [radius=0.07];
\draw[fill] (8,9) circle [radius=0.07];
\node[below] at (0,5) {$30$};
\node[below] at (4,7) {$33$};
\node[below] at (4,3) {$27$};
\node[below] at (8,5) {$29.7$};
\node[below] at (8,9) {$36.3$};
\node[below] at (8,1) {$26.4$};
\end{tikzpicture}
\caption{A two step binomial tree for a stock starting at $\$30$ }
\label{fig:ex1}
\centering
\end{figure}
\begin{enumerate}
\item Using (\ref{bt1}) to discount to the node with value $\$33$ we have that, $p = \frac{\mathrm{e}^{-0.12\times\frac{1}{12}}-0.9}{1.1-0.9} = 0.45$. Hence,
\begin{align*}
\mathrm{e}^{-\frac{1}{12}\times0.12}(0.45\times 5.3 + 0.55 \times 0) = 2.36.
\end{align*}
Then applying this again to discount back to the initial node we have, that the option has value $\$2.36$ at this node and value $\$0$ at the one below it. Hence,
\begin{align*}
\mathrm{e}^{-\frac{1}{12}\times0.12}(0.45\times 2.36 + 0.55 \times 0) = 1.053.
\end{align*}
Thus we have found the value of our option.
\item Now using the alternate formula given in equation (\ref{bt6}) we have,
\begin{align*}
f = \mathrm{e}^{-2\times 0.12\times \frac{1}{12}}(0.45^2 \times 5.3 + 2\times 0.45 \times 0.55 \times 0 + 0.55^2 \times 0) = 1.053.
\end{align*}
As before.
\end{enumerate}
We see that both of these methods are effective, however the general two step equation is in general easier to use. We now wish to generalise this method for $n \in \mathbb{N}$ time steps.
\end{ex}
\subsubsection{Generalization}
%########################################################################################################
We may further generalize this model to $n\in \mathbb{N}$ time steps, for European options \citep[Chap. 11]{JCH}. In Section \ref{twostep} we found the formula for a two step tree from a one step tree by repeatedly applying (\ref{bt3}). It follows from this that the formula for the value of the option after $j$ upward moves and $i$ downward moves is given by,
\begin{align*}
f_{u^j, d^{i}} = \mathrm{e}^{-r\Delta t}(p f_{u^{j+1}, d^{i}} + (1-p) f_{u^j, d^{i+1}}).
\end{align*}
Noting that our formula for the two step model has coefficients given by pascals triangle, and the powers of $p$ and $(1-p)$ follow this as well we make the following claim.
\begin{claim*}
The general price of a European option for an $n\in\mathbb{N}$ step tree is given by,
\begin{align}
f = \mathrm{e}^{-rn\Delta t}\left[\sum_{j=0}^{n} \binom{n}{j} p^j(1-p)^{n-j}f_{u^j,d^{n-j}} \right] \label{btp1}.
\end{align}
\end{claim*}
\begin{proof}
We prove this by induction. Clearly for $n = 1$ we have,
\begin{align*}
f &= \mathrm{e}^{-r\Delta t}\left[\sum_{j=0}^{1} \binom{1}{j} p^j(1-p)^{1-j}f_{u^j,d^{1-j}} \right] \\
&= \mathrm{e}^{-r\Delta t}\left[\binom{1}{0} p^0(1-p)^{1}f_{u^0,d^{1}} + \binom{1}{1} p^1(1-p)^{0}f_{u^1,d^{0}} \right] \\
&= \mathrm{e}^{-r\Delta t}[pf_u + (1-p)f_d].
\end{align*}
Which is exactly what we found in Section \ref{onestep} and is equation (\ref{bt1}). We can also check this for $n=2$ against equation (\ref{bt6}). We now perform our inductive step. Assume that (\ref{btp1}) hold for some $k \in \mathbb{N}$. Then note that,
\begin{align*}
f &= \mathrm{e}^{-rk\Delta t}\left[\sum_{j=0}^{k} \binom{n}{j} p^j(1-p)^{k-j}f_{u^j,d^{k-j}} \right] \\
&= \mathrm{e}^{-rk\Delta t} \left[ \binom{k}{0}(1-p)^kf_{d^k} + \binom{k}{1}p(1-p)^{k-1}f_{u,d^{k-1}} + \dots \right.\\
&\qquad + \left.\binom{k}{k-1}p^{k-1}(1-p)f_{u,d^{k-1}} + \binom{k}{k}p^kf_{u^k}\right].
\end{align*}
Now using the formula from (\ref{btp1}) we may expand this giving,
\begin{align*}
f &= \mathrm{e}^{-rk\Delta t} \left[ \binom{k}{0}(1-p)^k[pf_{u,d^k}+ (1-p)f_{d^{k+1}}]\mathrm{e}^{-r\Delta t} \right.\\
&\qquad + \binom{k}{1}p(1-p)^{k-1}[pf_{u^2,d^{k-1}}+ (1-p)f_{u,d^{k}}]\mathrm{e}^{-r\Delta t} + \dots \\
&\qquad \left. + \binom{k}{k}p^k[pf_{u^{k+1}}+ (1-p)f_{u^k,d}] \mathrm{e}^{-r\Delta t}\right].
\end{align*}
Here we note that for a term $f_{u^j,d^{k-j+1}}$ for $j = 1,2,\dots, k+1$ the only contributing factors from are $f_{u^j,d^{k-j+1}}$ and $f_{u^{j+1},d^{k-j}}$. These contributing factors are $\binom{k}{j-1}p^j(1-p)^{k-j+1}$ and $\binom{k}{j}p^{j-1}p(1-p)^{k-j}(p-1)$ respectively. Hence we have that the coefficient of $f_{u^j,d^{k-j+1}}$ is given by,
\begin{align*}
(1-p)^{k-j+1}p^j\binom{k}{j-1} + \binom{k}{j}(1-p)^{k-j+1}p^j &= (1-p)^{k-j+1}p^j \left( \binom{k}{j} + \binom{k}{j-1}\right) \\
&= (1-p)^{k-j+1}p^j \binom{k+1}{j}.
\end{align*}
Hence, as this is the coefficient of our next step we have that,
\begin{align*}
f &= \mathrm{e}^{-rk\Delta t}\sum^{k+1}_{j=1} \mathrm{e}^{-r\Delta t}\binom{k+1}{j}(1-p)^{k-j+1}p^j f_{u^j,d^{k-j+1}} \\
&= \mathrm{e}^{-r(k+1)\Delta t}\sum^{k+1}_{j=1} \binom{k}{j}(1-p)^{k-j}p^j f_{u^j,d^{k-j}}.
\end{align*}
We have shown that if it is true for $k$ then it is true for $k+1$. This completes the proof.
\end{proof}
\subsubsection{Finding $u$ and $d$}
We may find these parameters for both European and American options by considering variance \citep[Chap. 11]{JCH}. From Section (2.3) and (2.5) we see that for a single step $\Delta t$ the variance of the discrete version of (\ref{model}) is $\sigma^2 \Delta t$. Furthermore from a one step tree we see that the variance is clearly, $pu^2 + (1-p)d^2 - [pu + (1-d)]^2$. Equating these two we have that,
\begin{align*}
pu^2 + (1-p)d^2 - [pu + (1-d)]^2 = \sigma^2 \Delta t.
\end{align*}
Using the values for $p$ that we found in Section \ref{onestep} we have,
\begin{align*}
\mathrm{e}^{r\Delta t}(u+d) - ud - \mathrm{e}^{2r\Delta t} = \sigma^2 \Delta t.
\end{align*}
Using the series expansion for $\mathrm{e}^x$ and ignoring terms of higher order than $\Delta t$ we have that,
\begin{align*}
u &= \mathrm{e}^{\sigma \sqrt{\Delta t}}, \\
d &=\mathrm{e}^{-\sigma \sqrt{\Delta t}}.
\end{align*}
This is the Cox, Ross, and Rubenstien model.
\subsubsection{Moment Matching and Other Probabilities}
%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Above we have used a risk neutral argument to choose the parameters $u$ and $d$ \citep[Chap. 11]{JCH}. There are many methods to do this. These include the; Cox, Ross and Rubinstein model, the equal probability model, Cox, Ross and Rubinstein with a drift model and the Titan model.
The problem with the Cox, Ross, Rubinstein model is that it does not match the moments of Black-Scholes stochastic differential equation (S.D.E.) \citep{JRB}. So we seek to find out how we would choose the parameters $u$ and $d$ to match these moments. This is due to the approximation we used in the last section ignoring some terms.
As we have three unknowns to find we need three equations to define these. One idea came from Titan this is that we require that the first three moments of the binomial trees and the Black-Scholes S.D.E. This generates a new model not considered here.
Another model is the Cox, Ross, Rubinstein model with a drift. This introduces a new parameter $\eta$ so that the stock may drift. This allows us to model stocks more effectively as stocks tend to drift up or down. In this case we have that,
\begin{align*}
u &= \mathrm{e}^{\eta \Delta t + \sigma \sqrt{\Delta t}} \\
d &= \mathrm{e}^{\eta \Delta t - \sigma \sqrt{\Delta t}}
\end{align*}
with $p$ as before.
Alternatively as \citet{JRP} said we may derive these in an alternative way. Consider (\ref{model}) this implies which,
\begin{align*}
\log\left( \frac{S_T}{S_t}\right) = \mu t + \sigma \sqrt{t}W_t. \label{jr0}
\end{align*}
We may rewrite the above with $W_t = N$ where $N$ is a standardized normal distribution. Further from our discussion on binomial trees for an $n$ step tree we have that, if $u = \mathrm{e}^{u'}$ and $d = \mathrm{e}^{d'},$
\begin{align}
S_T &= S_t \mathrm{e}^{ju' + (n-j)d'} \\
\Rightarrow \log\left( \frac{S_T}{S_t}\right) &= nd' + (u'-d')j. \label{jr-1}
\end{align}
for a tree with $j$ upward movements, with probability,
\begin{align*}
\binom{n}{j}p^j(1-p)^{n-j}.
\end{align*}
It follows that this is binomialy distributed with mean $np$ and variance $np(1-p)$. Then we wish that for large $n$ that these two distributions are equal, causing their moments to be the same. This is a known result in probability \citet{KC}. Thus we have that, if $p$ is a constant independent of $n$, as $n\rightarrow\infty$,
\begin{align}
\frac{j-np}{\sqrt{np(1-p)}} \sim N \label{jr1}
\end{align}
where $\sim$ denotes convergence in the distribution. Now here we may choose $p = \frac{1}{2}$ as the above guarantees the first two moments are equivalent. We require a third equation to solve for our three unknowns. Then we have that (\ref{jr1}) becomes,
\begin{align*}
j \sim \frac{\sqrt{n}}{2}N + \frac{n}{2},
\end{align*}
and substituting this into (\ref{jr-1}) we have,
\begin{align*}
\log\left( \frac{S_T}{S_t}\right) \approx nd' + (u'-d')\left(\frac{\sqrt{n}}{2}N + \frac{n}{2} \right).
\end{align*}
Equating this with (\ref{jr0}) and equating coefficient of $N$ and constants we yield,
\begin{align*}
\left.
\begin{array}{ll}
\hat{\mu} t = nd' + (u'-d')\frac{n}{2} \\
\sigma\sqrt{t} = (u'-d')\frac{\sqrt{n}}{2}
\end{array}
\right\}
\Rightarrow
\left\{
\begin{array}{ll}
u' &= \frac{\hat{\mu} t }{n} + \sigma \sqrt{\frac{t}{n}} \\
d' &= \frac{\hat{\mu} t }{n} - \sigma \sqrt{\frac{t}{n}}
\end{array}
\right.
\end{align*}
where $\mu = r - \sigma^2/2$ thus we have that,
\begin{align*}
u = \mathrm{e}^{u'} = \mathrm{e}^{(r - \sigma^2/2)t/n + \sigma\sqrt{t/n}} \\
d = \mathrm{e}^{d'} = \mathrm{e}^{(r - \sigma^2/2)t/n - \sigma\sqrt{t/n}} \\
p = \frac{1}{2}.
\end{align*}
This is known as the \emph{Equal probability model} or \emph{Jarrow-Rudd} model.
%In so far we have been developing the \citet{CRR} model, through a similar method as in the paper. There is another approach which can be more fruitful to us.
%For the \citet{CRR} model consider the expected return after one step (i.e. in a time interval $\Delta t$). This is obviously, $\mathbb{E}(S_{\Delta t}) = puS_0 + (1-p)dS_0$. Further we may calculate the second moment as well yielding $\mathbb{E}(S^2_{\Delta t}) = pu^2S^2_0 + (1-p)d^2S_0^2$. Then, in our section (reference $E$ equations) we found that for geometric Brownian motion we have that $\mathbb{E}(S_t) = \mathrm{e}^{\mu t}$ and $\mathbb{E}(S^2_t) = S_0^2\mathrm{e}^{2 \mu t + \sigma^2t }$. We may derive the Cox, Ross, Rubinstein model by equating these and solving simultaneously. Observe that by setting $\mathrm{e}^{\mu t} = A$ and $\mathrm{e}^{2\mu t + \sigma^2 t} = B $ then we must have that,
%\begin{align*}
%A &= pu + (1-p)d \\
%B &= pu^2 (1-p)d^2.
%\end{align*}
%For the three unknowns we need an extra equation as discussed before this will be $ud=1$.
%There is a more sophisticated model baised the Cox, Ross, Rubenstein model. As is seen in Cox, Ross, Rubenstein model, the stock follows a random walk however this does not allow the stock to drift as it can be seen to do in reality, as in (ref section on ito processes). We instead now introduce a new drift parameter $\eta$ which will allow the stock to move in a given direction to derive this observe that...
%Another model we may look at follows a similar form of logic. This is the model proposed by (cite titan 93). Here we choose to match the third moment to give three equations instead of setting $ud=1$. This leads us to an alternative model...
\subsection{Trinomial Trees}
Trinomial trees have a near identical concept to that of binomial trees \citep[Chap. 11]{JCH}. The key difference is that instead of having two paths the stock can take at any given node we now have three. This now assumes that the stock can rise or fall by specified amounts or stay the same in any given time step. This can be seen from Figure \ref{fig:TTrees1}. The reason that we write $S_0m$ for the middle branch is as for models such as the Boyles model with a drift we may have that $m\neq 1$.
\begin{figure}[ht]
\centering
\begin{tikzpicture}[font=\small, scale=0.6]
\draw (0,5) -- (4,7);
\draw (0,5) -- (4,3);
\draw (0,5) -- (4,5);
\draw[fill] (0,5) circle [radius=0.07];
\draw[fill] (4,7) circle [radius=0.07];
\draw[fill] (4,3) circle [radius=0.07];
\draw[fill] (4,5) circle [radius=0.07];
\node[below] at (0,5) {$S_0$};
\node[below] at (4,7) {$S_0u$};
\node[below] at (4,3) {$S_0d$};
\node[below] at (4,5) {$S_0m$};
\end{tikzpicture}
\caption{A one step trinomial tree}
\label{fig:TTrees1}
\centering
\end{figure}
%########################################################################################################
\subsubsection{One-step Model}
Following exactly the same logic as before we have that for a one step tree our option price would be given by,
\begin{align*}
f = \mathrm{e}^{-rT}(p_uf_u + p_mf_m + p_df_d),
\end{align*}
where $p_u$ is the probability of going up and the analogous definitions for $p_m$ and $p_d$.
%########################################################################################################
\subsubsection{Generalization}
Here the generalization is very akin to that of the binomial tree case. Here we may work backwards from the end nodes to obtain the price at the initial node as before. The concepts and algorithms for price calculation are similar to the binomial case with the obvious changes for the middle branch.
Here we do not have the availability of a general formula for the European or American case. Thus we must work backwards through the tree to calculate our price.
%########################################################################################################
\subsubsection{Moment Matching and Other Probabilities}
Here we now have six unknowns to solve for. Assuming that we match the first two moments as in \citet{JRB} or the alternative derivation of Cox, Ross, and Rubinstein, we still need four more equations. Using that $p_u + p_d + p_m = 1$, $ud =1 $ and $m=1$ we still require another equation.
The approach proposed by \citet{BO} is to introduce a stretch parameter, $\lambda$, so that $u = \mathrm{e}^{\lambda\sigma\sqrt{t}}$. This can take many values but we shall only consider $\lambda = \sqrt{2}$ due to the equivalence of the trinomial trees method and the explicit finite difference method. Taking this value of $\lambda$ we obtain that,
\begin{align*}
u &= \mathrm{e}^{\lambda\sqrt{2\Delta t}},\\
d &= \mathrm{e}^{-\lambda\sqrt{2\Delta t}}, \\
m &= 1,
\end{align*}
with probabilities,
\begin{align*}
p_u = \left(\frac{\mathrm{e}^{\mu\Delta t/2 - } - \mathrm{e}^{-\sigma\sqrt{\Delta t/2}}}{\mathrm{e}^{\sigma\sqrt{\Delta t/2}}-\mathrm{e}^{-\sigma\sqrt{\Delta t/2}}} \right)^2 \\
p_d = \left(\frac{\mathrm{e}^{\mu\Delta t/2 - } - \mathrm{e}^{-\sigma\sqrt{\Delta t/2}}}{\mathrm{e}^{\sigma\sqrt{\Delta t/2}}-\mathrm{e}^{\sigma\sqrt{\Delta t/2}}} \right)^2 \\
p_m = 1-p_d-p_m.
\end{align*}
There are other methods including equal probabilities or four moment matching trees that we shall not consider here.
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
\section{Finite Difference Methods}
%all this from
\label{fdms}
We may use finite difference methods to solve the Black-Scholes equation and therefore price options. We first recall a few basic results about Taylor series and finite difference methods. Let $f(x)$ be a function that is twice differentiable, we know using Taylors theorem that,
\begin{align}
&f''(x) = \frac{1}{(\Delta x)^2} (f(x + \Delta x) + f(x - \Delta x) - 2f(x)) + \mathcal{O}((\Delta x)^2)\label{C1},\\
&f'(x) = \frac{1}{2\Delta x}(f(x + \Delta x) + f(x - \Delta x)) + \mathcal{O}((\Delta x)^2). \label{FDM1}
\end{align}
These are well known approximation to us, for the single variable case. We know that Equation (\ref{FDM1}) is known as the \emph{central difference approximation}, we are also familiar with the \emph{forward difference approximation} and \emph{backward difference approximation}, which respectively are,
\begin{align*}
f'(x) = \frac{1}{\Delta x}(f(x + \Delta x) -f(x)), \\
f'(x) = \frac{1}{\Delta x}(f(x) -f(x-\Delta x)).
\end{align*}
Now we attempt to extend these methods to functions of two variables, we will follow the same approach as \citet{GDS}.
%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
\subsection{Concept}
We begin by extending our notion of how to approximate a single variable function to two variables \citep[Chap. 17]{JCH}. Let $f$ be a function of two variables $x$ and $y$ as we may not divide the whole infinite plane we must truncate, choosing some $x_{min}$ and $x_{max}$ and similarly for $y$. We then subdivide the truncated $x,y$ plane into rectangles of width $\Delta x = h$ and height $\Delta y = k$, this may be seen in Figure \ref{fig:Subdivision}. Then we may represent a point $(x,y)$ on the corner of any rectangle as $x = ih$ and $y = jk$, for appropriate $i,j \in \mathbb{N}$. Then we denote the value of $f$ at this point as, $f(x,y)=f(ih,ik)=f_{i,j}$. The by (\ref{C1}) we have that,
\begin{align*}
\frac{\partial^2 f}{\partial x^2} &= \frac{f((i+1)h, jk) - 2f(ih,jk)+f((i-1)h,jk)}{h^2} \\
&=\frac{f_{i+1,j}-2f_{i,j}+f_{i-1,j}}{h^2},
\end{align*}
with leading error of order $h^2$. Similarly, we may find the second derivative with respect to $y$ and redefine the forward difference approximation and backward difference approximation in terms of our new notation, for two variables.
\begin{align*}
\frac{\partial^2 f}{\partial y^2} &=\frac{f_{i,j+1}-2f_{i,j}+f_{i,j-1}}{k^2} + \mathcal{O}(k^2), \\
\frac{\partial f}{\partial y} &= \frac{f_{i,j+1}-f_{i,j}}{k} + \mathcal{O}(k), \\
\frac{\partial f}{\partial y} &= \frac{f_{i,j}-f_{i,j-1}}{k} + \mathcal{O}(k).
\end{align*}
\begin{figure}[t]
\centering
\maketitle
\begin{tikzpicture}[font=\small, scale=1.2]
\draw[step=1cm,gray,very thin] (0,0) grid (6,6);
\draw[thick,->] (0,0) -- (6.5,0);
\draw[thick,->] (0,0) -- (0,6.5);
\draw[thick,<->] (-0.5,0) -- (-0.5,1);
\draw[thick,<->] (0,-0.5) -- (1,-0.5);
\node[below] at (3,-1) {$x$};
\node[left] at (-1,3) {$y$};
\node[left] at (-0.5,0.5) {$\Delta y$};
\node[left] at (0.75,-0.75) {$\Delta x$};
\filldraw[black] (3,3) circle (2pt) node[anchor=north west] {$i,j$};
\filldraw[black] (4,3) circle (2pt) node[anchor=north west] {$i+1,j$};
\filldraw[black] (3,4) circle (2pt) node[anchor=north west] {$i,j+1$};
\filldraw[black] (2,3) circle (2pt) node[anchor=north west] {$i-1,j$};
\filldraw[black] (3,2) circle (2pt) node[anchor=north west] {$i,j-1$};
%\draw (7.5,0) -- (7.5, -0.1);
%\node[below] at (7.5,0) {$a$};
%\draw (0,3) -- (-0.1,3);
%\node[left] at (0,3) {$b$};
%\node[right] at (10,4) {$y=\frac{b}{a}x$};
\end{tikzpicture}
\caption{The partitioning of the $x,y$ plane for a two dimensional finite difference method }
\label{fig:Subdivision}
\centering
\end{figure}
For options we are attempting to approximate the derivatives of the price of the option with respect to time and the underlyings' price, as well as the second derivative with respect to the underlyings' price. Following the method we described we first need to partition the plane, here the $(t,S)$ plane. We choose divide $T$ into $N$ equally spaced intervals, so $h = \Delta t = \frac{T}{N}$. Similarly we choose to partition $S$ into $M$ equally spaced intervals so $k = \Delta S = \frac{S}{M}$, indexing $N$ by $i$ and $M$ by $j$. It is important that $M$ be chosen so that $S_0$ is one of the points considered ($i=0$), so we may calculate the value of the option at this point.
\subsection{Terminal and Boundary Conditions}
%We have the framwork required to approximate the derivatives, however to solve this we need to intial and boundry conditions. These will differ for both calls and puts, these can be seen in Tabel \ref{tab:bc's}, for European options. These are fairly obvious for European options and can be seen fairly easily. The value at time $T$ is their value at expariation and as such should be the payoff function.
%\begin{table}[ht]
%\centering
%\begin{tabular}{|l||c|c|c|c|}
%\hline
%Boundry & European Call & European Put \\ \hline
%$t = T$ & $V(S,T) = \max(S-K,0)$ & $V(S,T) = \max(K-S,0)$ \\
%$S=0$ & $V(0,t) = 0$ & $V(0,t) = K\mathrm{e}^{-r(T-t)}$ \\
%$S=\infty$ & $\lim_{s\rightarrow \infty}(V(S,t) - (S-K\mathrm{e}^{-r(T-t)})) =0$ & $\lim_{s \rightarrow \infty}(V(S,t)) =0$ \\\hline
%\end{tabular}
%\caption{\label{tab:bc's} A tabel showing the conditions at each boundry.}
%\end{table}
%For a call option the stike price $K$ i always greater than zero. However the value for the option when the stock price is zero is obviously zero, as clearly the strike price is above zero. Further, for a call, as the stock price tends to infinity, the difference between the value of the option at expariation ((((((((((((((((((((((((((((((((((((((((((THIS LAST ONE COMES FROM THE BLACK-SCHOLES FOMRULAE IN JCHULL CHAP 13 PG 291 EQN 13.20 AND PG 293)))))))))))))))))))))))))))
%For a put option if the stock price is zero then the value of the put is $K$, but we have to discount this by the risk-free intrest rate as this would be its value at expairation and not currently. The second coniditon comes from the fact that as $S$ tends to infinity, $S$ will certainly be greater than $K$. Hence the value of this will be zero, and when discounted to the current tim step is still zero.
%Futher we may entend these ideas to American options to obtain boundry conditions for these. Due to the ability to exercise at anytime we need not discount. Otherwise the arguments are identical and we have the conditions in tabel \ref{tab:bc's2}.
%\begin{table}[ht]
%\centering
%\begin{tabular}{|l||c|c|c|c|}
%\hline
%Boundry & American Call & American Put \\ \hline
%$t = T$ & $V(S,T) = \max(S-K,0)$ & $V(S,T) = \max(K-S,0)$ \\
%$S=0$ & $V(0,t) = 0$ & $V(0,t) = K$ \\
%$S=\infty$ & $\lim_{s\rightarrow \infty}(V(S,t) - (S-K\mathrm{e}^{-r(T-t)})) =0$ & $\lim_{s \rightarrow \infty}(V(S,t)) =0$ \\\hline
%\end{tabular}
%\caption{\label{tab:bc's2} A tabel showing the conditions at each boundry.}
%\end{table}
We have the framework required to approximate the derivatives, however to solve this we need terminal and boundary conditions \citep[P. 78]{MG}. As we are required to truncate the plane the actual boundary conditions at $\infty$ and $0$ are of little use to us. We must therefore examine the boundary conditions for the truncated plane. These will differ for both calls and puts, and for European and American options. %These conditions can be seen in Table for European options.
Say that we truncate the plane so that $S \in [S_{\min}, S_{\max}]$, we still consider $t \in [0,T]$. Then we have the boundary and terminal conditions for a European option, are given in Table \ref{tab:bc'st}, where $V(S,t)$ is the value of the option at stock price $S$ and at time $t$.
\begin{table}[ht]
\centering
\begin{tabular}{|l||c|c|c|c|}
\hline
Boundary & European Call & European Put \\ \hline
$t = T$ & $V(S,T) = \max(S_{\min} + j\Delta S - K,0)$ & $V(S,T) = \max(K-S_{\min} - j\Delta S,0)$ \\
$S=S_{\min}$ & $V(S_{\min},t) = \max(S_{\min}-K,0)$ & $V(S_{\min},t) = K\mathrm{e}^{-r(T-t)}$ \\
$S=S_{\max}$ & $V(S_{\max},t) = \max(S_{\max} - K,0)$ & $V(S_{\max},t) = \max(K-S_{\max},0)$ \\\hline
\end{tabular}
\caption{\label{tab:bc'st} The conditions at each boundary of the gird for a European option }
\end{table}
These boundary and side conditions are found by using the payoff function at the necessary points. The points of interest are at the termination date and the sides as this is where we have truncated our plane to, hence we arrive at the above. We may then extend this idea to American options.
The main difference between the American and European case is the lack of need to discount, as we are able to exercise early. This means that our boundary and terminal conditions for an American option are given by Table \ref{tab:bc'st2}.
\begin{table}[ht]
\centering
\begin{tabular}{|l||c|c|c|c|}
\hline
Boundary & American Call & American Put \\ \hline
$t = T$ & $V(S,T) = \max(S_{\min} + j\Delta S - K,0)$ & $V(S,T) = \max(K-S_{\min} - j\Delta S,0)$ \\
$S=S_{\min}$ & $V(S_{\min},t) = \max(S_{\min}-K,0)$ & $V(S_{\min},t) = K$ \\
$S=S_{\max}$ & $V(S_{\max},t) = \max(S_{\max} - K,0)$ & $V(S_{\max},t) = \max(K-S_{\max},0)$ \\\hline
\end{tabular}
\caption{\label{tab:bc'st2} The conditions at each boundary of the gird for an American option}
\end{table}
%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
\subsection{Implicit Method}
\label{ifdms}
Hence partitioning the plane as we have discussed, we now have that $h = \Delta S$, $k = \Delta t$ and $x = S$. When approximating the Black-Scholes differential equation we use a backwards difference approximation for the $\frac{\partial f}{\partial S}$ derivative, a forward difference approximation for the $\frac{\partial f}{\partial t}$ and a central difference approximation for the $\frac{\partial^2 f}{\partial S^2}$ derivative \citep[Chap, 17]{JCH}. Hence, the Black-Scholes equation becomes,
\begin{align}
\frac{f_{i+1,j} - f_{i,j}}{\Delta t} + r j \Delta S\frac{f_{i,j+1} - f_{i,j-1}}{2\Delta S} + \frac{1}{2} \sigma^2 (j\Delta S)^2\frac{f_{i,j+1} + f_{i,j-1} - 2 f_{i,j}}{\Delta S^2} = rf_{i,j}. \label{IM1}
\end{align}
We may then rearrange this to find coefficients for the $f_{i,j-1}$, $f_{i,j}$ and $f_{i, j+1}$. This yields that (\ref{IM1}) is now,
\begin{align}
a_jf_{i,j-1} + b_jf_{i,j} + c_jf_{i,j+1} = f_{i+1,j} \label{FDMI}
\end{align}
where,
\begin{align*}
a_j &= \frac{1}{2}(r-q)j\Delta t - \frac{1}{2}\sigma^2j^2\Delta t, \\
b_j &= 1 + \sigma^2 j^2 \Delta t + r \Delta t, \\
c_j &= -\frac{1}{2}(r-q)j\Delta t - \frac{1}{2}\sigma^2j^2\Delta t.
\end{align*}
Note that these constants are independent of $S$. This is because the $\Delta S$ and $\Delta S^2$ introduced by the derivatives are the same $\Delta S$ that we are considering the option over in any given triangle. Thus once we discretize the plane we have cancellation resulting in the above relationships.
Now we need to solve this so that we can calculate the value of our option at $t=0$ i.e. $i=0$. For a European option we begin by considering $i = N-1$ in (\ref{FDMI}), for $j = 1,\dots,M-1$. Thus we have $M-1$ simultaneous equations to solve. If we write out a few of these,
\begin{align*}
a_1 f_{N-1,0} + b_1 f_{N-1,1} + c_1f_{N-1,2} &= f_{N,1} \\
a_2 f_{N-1,1} + b_2 f_{N-1,2} + c_2f_{N-1,3} &= f_{N,2} \\
&\vdots \\
a_{M-2}f_{N-1,M-3} + b_{M-2}f_{N-1,M-2} + c_{M-2}f{N-1,M-1} &= f_{N,M-1}\\
a_{M-1}f_{N-1,M-2} + b_{M-1}f_{N-1,M-1} + c_{M-1}f{N-1,M} &= f_{N,M}.
\end{align*}
Note that all the $f_{i,N}$ are known from the boundary conditions. Furthermore the terms $a_1f_{N-1,0}$ and $c_{M-1}f_{N-1,M}$ are known from the boundary conditions. We may then express our system with unknowns on the left and knowns on the right.
\begin{align*}
\begin{pmatrix}
b_{1} & c_{1} &&&& \\
a_{2} & b_{2} & c_2 &&&\\
& a_{3} & b_{3} & c_3 &&\\
&& \ddots & \ddots & \ddots&\\
&&& a_{M-2}&b_{M-2}&c_{M-2} \\
&&&&a_{M-1}&b_{M-1}
\end{pmatrix}
\begin{pmatrix}
f_{N-1,1} \\
f_{N-1,2} \\
f_{N-1,3} \\
\vdots \\
f_{N-1,M-2} \\
f_{N-1,M-1}
\end{pmatrix}
=
\begin{pmatrix}
f_{N,1} \\
f_{N,2} \\
f_{N,3} \\
\vdots \\
f_{N,M-2} \\
f_{N,M-1}
\end{pmatrix}
+
\begin{pmatrix}
-a_1f_{N-1,0} \\
0 \\
0 \\
\vdots \\
0 \\
-c_{M-1}f_{N-1,M}
\end{pmatrix}.
\end{align*}
Allowing $A$ to be the leftmost matrix and $F_i$ to be the leftmost vector so that $F_{i+1}$ is the second leftmost vector and $B_i$ to be the rightmost vector we have that,
\begin{align}
AF_i = F_{i+1} + B_i. \label{MF}
\end{align}
In the matrix form above we have taken $i=N-1$ however we note that this is true for arbitrary $i$.
We may solve these using the tridiagonal matrix algorithm for $N-1$, then the $F_{N-1}$ we have found is used in (\ref{MF}) with $i=N-2$ to find $F_{N-2}$. We may continue doing this until we reach $i=0$ and we may find the value of our option at this point.
For an American option we apply the same method however after each iteration of the method we compare the $F_{i}$ we found to a vector, $P$ the payoff at each $j$,
\begin{align*}
P_j = \left\{
\begin{array}{lll}
j\Delta \max(S - K,0) & \mbox{if call} \\
\max(K-j\Delta S,0) & \mbox{if put}
\end{array}
\right.
\end{align*}
we compare the each of the $j=1,\dots,M-1$ vector and take the maximum. After performing our method once, we obtain
\begin{align*}
F_{i,j} = \max(P_j,F_{i,j})
\end{align*}
where $F_{i,j}$ is the $j$th component of $F_i$.
So we have discussed how to use the implicit method to find the value of an option. There is also the explicit method which has a very different solution method.
%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
\subsection{Explicit Method}
\label{FDME}
The main difference between the implicit and explicit methods is that we assume that the values of $\frac{\partial f}{\partial S}$ and $\frac{\partial^2 f}{\partial S^2}$ at a given point $(i,j)$ is the same as at the point $(i+1,j)$ for the explicit method \citep[Chap. 17]{JCH}. We then have that,
\begin{align*}
\frac{f_{i+1,j} - f_{i,j}}{\Delta t} + r j \Delta S\frac{f_{i+1,j+1} - f_{i+1,j-1}}{2\Delta S} + \frac{1}{2} \sigma^2 (j\Delta S)^2\frac{f_{i+1,j+1} + f_{i+1,j-1} - 2 f_{i+1,j}}{(\Delta S)^2} = rf_{i,j}.
\end{align*}
Again we rearrange the above to obtain,
\begin{align*}
f_{i,j} = a^*_{j}f_{i+1,j-1} + b^*_{j}f_{i+1,j} + c^*_jf_{i+1, j+1}
\end{align*}
where,
\begin{align*}
a^*_j &= \frac{1}{1+r\Delta t}\left( -\frac{1}{2}(r-q)j\Delta t + \frac{1}{2}\sigma^2j^2 \Delta t\right), \\
b^*_j &= \frac{1}{1+r\Delta t}\left( 1-\sigma^2j^2\Delta t\right), \\
c^*_j &= \frac{1}{1+r\Delta t}\left( \frac{1}{2}(r-q)j\Delta t + \frac{1}{2}\sigma^2j^2 \Delta t\right).
\end{align*}
Note here that this is a much simpler method to implement. To go back a time step it is merely a linear combination of the nodes that came before it. Thus, as the value of $f_{N,j}$ are known, we may start with these and work back to $i=0$ for the price of our option.
%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
\subsection{Change of Variables}
% refference lognormal property for this one.
\label{coftt}
We discussed before in Section \ref{lnp} that the Black-Scholes equation is lognomally distributed. This gives us motivation to perform the change of variables, $z = \ln(x)$ \citep[Chap. 17]{JCH}.
We begin this by defining $G = \ln(S)$. Then by computing derivatives we see that the Black-Scholes differential equation becomes,
\begin{align*}
\frac{\partial f}{\partial t} + \left( r - q - \frac{\sigma^2}{2}\right)\frac{\partial f}{\partial G} + \frac{1}{2}\sigma^2 \frac{\partial^2 f}{\partial G^2} = rf.
\end{align*}
We now apply the finite difference methods again. Note that we are now considering equally spaced values of $G$ and not $S$. Then we have that the implicit methods' equation becomes,
\begin{align*}
\frac{f_{i+1,j}-f_{i,j}}{\Delta t} + (r-q-\frac{\sigma^2}{2})\frac{f_{i,j+1}-f_{i,j-1}}{2\Delta G} + \frac{1}{2}\sigma^2 \frac{f_{i,j+1} + f_{i,j-1} - 2 f_{i,j}}{\Delta G^2} = rf_{i,j}.
\end{align*}
So that now, we have
\begin{align*}
&\alpha_jf_{i,j-1} + \beta_j f_{i,j} + \gamma_j f_{i,j+1} = f_{i+1,j},
\end{align*}
where,
\begin{align*}
&\alpha_j = \frac{\Delta t}{2\Delta G}(r-q-\frac{\sigma^2}{2}) - \frac{\Delta t}{2\Delta G^2}\sigma^2, \\
&\beta_j = 1 + \frac{\Delta t}{\Delta G^2 }\sigma^2 + r\Delta t, \\
&\gamma_j = -\frac{\Delta t}{2\Delta G}(r-q-\frac{\sigma^2}{2}) - \frac{\Delta t}{2\Delta G^2}\sigma^2.
\end{align*}
Furthermore the explicit methods' equation becomes,
\begin{align*}
\frac{f_{i+1,j}-f_{i,j}}{\Delta t} + (r-q-\frac{\sigma^2}{2})\frac{f_{i+1,j+1}-f_{i+1,j-1}}{2\Delta G} + \frac{1}{2}\sigma^2 \frac{f_{i+1,j+1} + f_{i+1,j-1} - 2 f_{i+1,j}}{\Delta G^2} = rf_{i,j}.
\end{align*}
Again we now have that,
\begin{align*}
&\alpha_j^*f_{i+1,j-1} + b_j^*f{i+1,j} + c_j^*f{i+1,j+1} = f{i,j}, \\
\end{align*}
where,
\begin{align*}
&\alpha_j^* = \frac{1}{1+ r \Delta t} \left( -\frac{\Delta t}{2\Delta G}(r-q-\frac{\sigma^2}{2}) + \frac{\Delta t}{2\Delta G^2}\sigma^2\right), \\
&\beta_j^* = \frac{1}{1+r\Delta t}\left( 1- \frac{\Delta t}{2\Delta G^2}\sigma^2\right), \\
&\gamma_j^* = \frac{1}{1+ r \Delta t} \left( \frac{\Delta t}{2\Delta G}(r-q-\frac{\sigma^2}{2}) + \frac{\Delta t}{2\Delta G^2}\sigma^2\right).
\end{align*}
Note that this change of variables method has eliminated dependency on $j$ for all our constants. Thus these are constant for each $j$ as well as all $i$. This means that they need only be calculated once, a helpful trait. % from JCHull pg 427. May try to find a better source refference this later.
We now need to consider the boundary conditions. If we transform the variable it follows we should transform the boundary conditions. Hence these are now given in Tables \ref{tab:tbc's} and \ref{tab:tbc's2}.
\begin{table}[ht]
\centering
\begin{tabular}{|l||c|c|c|c|}
\hline
Boundary & European Call & European Put \\ \hline
$t = T$ & $V(S,T) = \max(\mathrm{e}^{S_{\min}} + \mathrm{e}^{j\Delta S} - \mathrm{e}^K,1)$ & $V(S,T) = \max(\mathrm{e}^K-\mathrm{e}^{S_{\min}} - \mathrm{e}^{j\Delta S},1)$ \\
$S=S_{\min}$ & $V(S_{\min},t) = \max(\mathrm{e}^{S_{\min}}-\mathrm{e}^K,1)$ & $V(S_{\min},t) = \mathrm{e}^{-r(T-t)}\mathrm{e}^K$ \\
$S=S_{\max}$ & $V(S_{\max},t) = \max(\mathrm{e}^{S_{\max}} - \mathrm{e}^K,1)$ & $V(S_{\max},t) = \max(\mathrm{e}^K-\mathrm{e}^{S_{\max}},1)$ \\\hline
\end{tabular}
\caption{\label{tab:tbc's} The conditions at each boundary of the transformed gird for European option }
\end{table}
\begin{table}[ht]
\centering
\begin{tabular}{|l||c|c|c|c|}
\hline
Boundary & American Call & American Put \\ \hline
$t = T$ & $V(S,T) = \max(\mathrm{e}^{S_{\min}} + \mathrm{e}^{j\Delta S} - \mathrm{e}^K,1)$ & $V(S,T) = \max(\mathrm{e}^K-\mathrm{e}^{S_{\min}} - \mathrm{e}^{j\Delta S},1)$ \\
$S=S_{\min}$ & $V(S_{\min},t) = \max(\mathrm{e}^{S_{\min}}-\mathrm{e}^K,1)$ & $V(S_{\min},t) = \mathrm{e}^K$ \\
$S=S_{\max}$ & $V(S_{\max},t) = \max(\mathrm{e}^{S_{\max}} - \mathrm{e}^K,1)$ & $V(S_{\max},t) = \max(\mathrm{e}^K-\mathrm{e}^{S_{\max}},1)$ \\\hline
\end{tabular}
\caption{\label{tab:tbc's2} The conditions at each boundary of the transformed gird for an American option}
\end{table}
%########################################################################################################
\subsection{Comparison of Explicit Method to Trinomial Trees}
\label{compttress}
The explicit method is equivalent to the trinomial tree approach. This is due to the mechanics of the method. The explicit method as seen in (\ref{FDME}) gives a relationship between the three values at the next time step $(i+1)\Delta t$ and the value current time step, $i,\Delta t$ \citep[P. 425-427]{JCH}. This can be seen in Figure \ref{fig:EM1}.
\begin{figure}[ht]
\centering
\begin{tikzpicture}[font=\small, scale=0.6]
\draw (0,5) -- (4,7);
\draw (0,5) -- (4,3);
\draw (0,5) -- (4,5);
\draw[fill] (0,5) circle [radius=0.07];
\draw[fill] (4,7) circle [radius=0.07];
\draw[fill] (4,3) circle [radius=0.07];
\draw[fill] (4,5) circle [radius=0.07];
\node[below] at (0,5) {$f_{i,j}$};
\node[above] at (4,7) {$f_{i+1,j+1}$};
\node[below] at (4,3) {$f_{i+1,j-1}$};
\node[below] at (4,5) {$f_{i+1,j}$};
\end{tikzpicture}
\caption{A one step trinomial tree}
\label{fig:EM1}
\centering
\end{figure}
Through this graphic it is easily seen how this can be equivalent to trinomial tree approaches. We can see that the terms $a^*_j$, $b^*_j$ and $c^*_j$ may be seen as follows:
\begin{align*}
\begin{array}{ll}
\left( -\frac{1}{2}rj\Delta t + \frac{1}{2}\sigma^2j^2 \Delta t\right) & \text{Probability that the stock price decreases in time } \Delta t \\
& \text{from }j\Delta S \text{ to } (j-1)\Delta S, \\
& \\
\left( 1-\sigma^2j^2\Delta t\right)& \text{Probability that the stock price is constant in time } \Delta t \\
& \text{remaining at }j\Delta S, \\
& \\
\left( \frac{1}{2}rj\Delta t + \frac{1}{2}\sigma^2j^2 \Delta t\right) & \text{Probability that the stock price increases in time } \Delta t \\
& \text{from }j\Delta S \text{ to } (j+1)\Delta S.
\end{array}
\end{align*}
One may be inclined to say that this would not be possible, as the value of $S_0$ increase by exponentiation of either $u$ or $d$. Thus the partitioning grid would have to be non-linear. We deal with this by performing the change of variables as seen in Section \ref{coftt}. Then we may use a linearly increasing grid.
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
\section{Sensitivities of Financial Derivatives and their Numerical Estimation}
\label{greekssection}
%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
There are a number of derivatives that are of interest to us. These pertain to how the price of the option, $f$, changes with respect to the different factors that affect it \citep[Chap. 15]{JCH}. These changes are termed ``Greeks'' due to the Greek letters that refer to them, the three most common are Delta, Theta and Gamma. These are defined as,
\subsection{Delta, Theta and Gamma}
\label{DTG}
\begin{defn}[Delta, Theta and Gamma]
\textbf{Delta}, $\Delta$, is a measure of the rate of change of the options calculated value, $f$, with respect to the change of the underlying assets’ price, $S$. I.e.
\begin{align*}
\Delta = \frac{\partial f}{\partial S}.
\end{align*}
\textbf{Theta}, $\Theta$, is a measure of $f$ with respect to the passage of time, $t$. This is sometimes referred to as the \textbf{time decay} of the value of the option. I.e.
\begin{align*}
\Theta = - \frac{\partial f}{\partial t}.
\end{align*}
\textbf{Gamma}, $\Gamma$, is a measure of the rate of change of $\Delta$ with respect to the price of the underlying asset. It could be thought of as the Delta of Delta. I.e.
\begin{align*}
\Gamma = \frac{\partial \Delta}{\partial S} = \frac{\partial^2 f}{\partial S^2}.
\end{align*}
\end{defn}
%########################################################################################################
\subsubsection{Relationship Between Delta, Theta and Gamma}
From the Black-Scholes differential equation we can see a useful relationship between Delta, Theta and Gamma \citep[Chap. 359]{JCH}. Delta, Theta, and Gamma are defined as derivatives consider the Black-Scholes differential equation,
\begin{align*}
\frac{\partial f}{\partial t} + rS\frac{\partial f}{\partial S} + \frac{1}{2}\sigma^2 S^2\frac{\partial^2 f}{\partial S^2} = rf.
\end{align*}
Note that every derivative is a Greek. Thus by direct substitution we have,
\begin{align}
\Theta + rS\Delta + \frac{1}{2} \sigma^2 S^2 \Gamma = rf. \label{greeks1}.
\end{align}
\begin{figure}[ht]
\centering
\begin{subfigure}[b]{0.3\textwidth}
\includegraphics[width=\textwidth]{ThetaGammaSmall.png}
\caption{The relationship between Theta and Gamma when their values are small}
\label{fig:TGS}
\end{subfigure}%
\quad %add desired spacing between images, e. g. ~, \quad, \qquad, \hfill etc.
%(or a blank line to force the subfigure onto a new line)
\begin{subfigure}[b]{0.3\textwidth}
\includegraphics[width=\textwidth]{ThetaGammaLarge.png}
\caption{The relationship between Theta and Gamma when their values are large}
\label{fig:TGL}
\end{subfigure}
\caption{The relationship between Gamma and Theta}\label{fig:animals}
\end{figure}
Alternatively if we use $f$ as the value of a portfolio, $\Pi$ and have the analogous definitions Theta, Delta and Gamma the above still holds for the whole portfolio.
When we constructed the ideas for a binomial tree for one step we stipulated that the portfolio should be riskless. Stating that the up and down values of the portfolio should be equal, this is also known as creating a delta-neutral portfolio. That is we have constructed a portfolio where $\Delta = 0$. If this is the case then (\ref{greeks1}) yields,
\begin{align*}
\Theta + \frac{1}{2}\sigma^2S^2\Gamma = r\Pi.
\end{align*}
Notice that if Theta or Gamma is large then the contribution from the right hand side becomes negligible. It follows that if Theta is large and positive then Gamma must be large and negative. This relationship for Theta and Gamma can be seen in Figure \ref{fig:TGL}; furthermore we may see the lack of this correlation in Figure \ref{fig:TGS}. This will become useful to us as we investigate the performance of these methods in later sections.
%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
\subsection{Vega}
The other most used Greek is Vega. This is defined as below.
\begin{defn}[Vega]
\textbf{Vega}, $\nu$, is the derivative of the value of the option with respect to its volatility, $\sigma$. I.e.
\begin{align*}
\nu = \frac{\partial f}{\partial \sigma}.
\end{align*}
\end{defn}
It is easy to see why these are not as popular as some of the other Greeks. Knowing the change of option price with respect to the stock price is obviously more likely to be useful than that of the risk-free interest rate.
These Greeks are seldom exactly calculable, as a general formula for $f$ is not often readily available. As such it is often the case that we must estimate these using numerical methods just as we must solve the Black-Scholes equation using numerical methods. The next section details how we may estimate the Greeks through the methods we have discussed insofar.
%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
\subsection{Estimation of Greeks}
%$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
\subsubsection{Estimation of Greeks Through Black-Scholes}
Due to the existence of the solution of the Black-Scholes differential equation for European options it is simple to differentiate our solution to obtain exact formulae for the calculation of the Greeks \citep[P. 180-181]{NAC}. These can be seen in Table \ref{tab:gbs}. Note here that $\Phi'(x)$ is the derivative of $\Phi(x)$ with respect to $x$,
\begin{table}[ht]
\centering
\begin{tabular}{|l||c|c|}
\hline
Greek & European Call & European Put \\ \hline
Delta & $\mathrm{e}^{-qT}\Phi(d_1)$ & $\mathrm{e}^{-qT}(\Phi(d_1) -1)$ \\
Theta & $-\frac{S_0\Phi'(d_1)\sigma\mathrm{e}^{-qT}}{2\sqrt{T}} $ & $-\frac{S_0\Phi'(d_1)\sigma\mathrm{e}^{-qT}}{2\sqrt{T}} $ \\
&$-rK\mathrm{e}^{-rT}\Phi(d_2) + qS_0\Phi(d_1)\mathrm{e}^{-qT}$&$+rK\mathrm{e}^{-rT}\Phi(-d_2) - qS_0\Phi(-d_1)\mathrm{e}^{-qT}$\\
Gamma & $\frac{\Phi'(d_1)\mathrm{e}^{-qT}}{S_0\sigma\sqrt{T}}$ &$\frac{\Phi'(d_1)\mathrm{e}^{-qT}}{S_0\sigma\sqrt{T}}$ \\
Vega & $S_0\sqrt{T}\Phi'(d_1)\mathrm{e}^{-qT} $ & $S_0\sqrt{T}\Phi'(d_1)\mathrm{e}^{-qT} $ \\
\hline
\end{tabular}
\caption{\label{tab:gbs} The exact formula for the Greeks for European options }
\end{table}
%%Black-scholes book pg 180 Chirss BS and Beyond.
For other types of options, such as American and Asian, as the Black-Scholes equation does not admit an analytic solution calculation of the Greeks through analytic formulae is not possible.
%aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa
\subsubsection{Estimation of Greeks Through Monte Carlo}
\label{hghghghghf}
Estimation of Greeks through Monte Carlo must follow a very different tact as we are numerically finding the value of the option and as such we must needs calculate the Greeks in a different way \citep[Chap. 7]{PG}.
As we shall, see finite difference methods will serve us well and it is through these, that we will calculate the Greeks for Monte Carlo. The parameters that are required to use Monte Carlo are $\sigma, r,S_0,T$ from (\ref{mce}), along with the number of replications. Then we may calculate the Greek Delta in the following method,
\begin{enumerate}
\item Simulate $N$ option prices $C(\sigma, r,S_0,T)$,
\item Simulate a further $N$ option prices slightly incremented $C(\sigma, r,S_0 + \Delta S,T) = C^i$,
\item Calculate the average of these in the usual way $\overline{C}^i$ and $\overline{C}$,
\item Gain a finite difference approximation using these values.
\end{enumerate}
Obviously we may choose to increment $\sigma$, $S_0$ or $T$ to gain our approximations for $\nu,\Delta$ or $\Gamma,$ and $\Theta$ respectively. For Delta, Theta and Vega we may use a central difference approximation, forward difference or backward difference. Each of these has its own advantages and disadvantages but notice that for Gamma, as in Section \ref{fdms}, we need to use a central difference approximation. So we have that Gamma is given by,
\begin{align*}
\Gamma \approx \frac{\hat{C}(\sigma, r,S_0 + \Delta S,T) - 2\hat{C}(\sigma, r,S_0,T) + \hat{C}(\sigma, r,S_0 - \Delta S,T) }{\Delta S^2}.
\end{align*}
For the other Greeks depending on the difference approximation used we obtain different formulae. Consider a forward difference approximation of Delta denoted $\Delta_F$, then we have that,
\begin{align*}
\Delta_F = \frac{\overline{C}(S_0 + \Delta S) - \overline{C}(S_0)}{\Delta S},
\end{align*}
removing the other parameters as we need not consider them. Then it follows that, if $\alpha(S) = \mathbb{E}[\Delta_F]$,
\begin{align*}
\mathbb{E}[\Delta_F] = \Delta S_0^{-1} (\alpha(S_0+\Delta S) - \alpha).
\end{align*}
Then if $\alpha(S_0)$ is twice differentiable at $S_0$, then using a Taylor expansion of order $o(\Delta S^2)$ we have that,
\begin{align*}
\textnormal{Bias}(\Delta_F) = \mathbb{E}[\Delta_F - \alpha'(S_0)] = \frac{1}{2}\alpha''(S_0)h + o(h).
\end{align*}
Following an similar argument we may show that,
\begin{align*}
\textnormal{Bias}[\Delta_C] = \frac{1}{6}\alpha'''(S_0)\Delta S^2 + o(\Delta S^2).
\end{align*}
Then it follows that the variance of a forward difference estimator is,
\begin{align*}
\textnormal{Var}[\Delta_F] = \Delta S^2 \textnormal{Var}[\overline{C}(S_0 + \Delta S_0) - \overline{C}(S_0)].
\end{align*}
Then $\textnormal{Var}[\overline{C}(S_0 + \Delta S_0) - \overline{C}(S_0)] = \frac{1}{n}\textnormal{Var}[C(S_0 + \Delta S_0) - C(S_0)]$ is either $\mathcal{O}(1)$ if we sample using different random numbers. Then,
\begin{align*}
\textnormal{Var}(\Delta_\beta) &= \frac{\sigma^2}{n\Delta S^2} + o(\Delta S^{-2}) \\
\textnormal{Bias}(\Delta_\beta) &= b\Delta S^\beta + o(\Delta S^\beta),
\end{align*}
for some non-zero $b$ and taking $\beta = 1$ for the forward difference case and $\beta = 2$ for the central difference case. Then consider a $\Delta S$ of the form $\Delta_n S = \Delta S_* n ^{-\gamma}$ for some positive $\Delta S_*$ and where $\Delta_n$ is the estimated value of Delta from a Monte Carlo simulation using $n$ paths. It can be shown that the optimal $\gamma$ is $\frac{1}{2\beta + 2}$. Then we have our mean square error for a general Delta is given by,
\begin{align*}
\textnormal{MSE}(\Delta) = b^2 \Delta S_n^{2\beta} + \frac{\sigma^2}{n\Delta S_n^2}.
\end{align*}
Taking the square root of this yields the \emph{root mean square error}, $\textnormal{RMSE}$. This is a good measure of convergence, so taking using $\Delta_n S = \Delta S_* n ^{-\gamma}$ with our value of$\gamma$ we see that,
\begin{align*}
\textnormal{RMSE}(\Delta) = \mathcal{O}\left( n^{-\frac{\beta}{2\beta + 2}}\right).
\end{align*}
So the rates of convergence for the forward and central difference methods are $\mathcal{O}(n^{-1/4})$ and $\mathcal{O}(n^{-1/3})$. Furthermore it can be shown with a little more analysis that the optimal $\Delta S$ in for each of these methods is $\Delta S^{-4} \sim n$ for the forward difference case and $\Delta S^{-6} \sim n$ for the central difference case, from the optimal $\gamma$ and the appropriate values of $\beta$.
%aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa
\subsubsection{Estimation of Greeks Through Trees}
\label{greekstreesss}
We may estimate the Greeks from binomial trees using finite difference methods \citep[P. 397-398]{JCH}. This allows us to evaluate the Greeks at each node. Consider a tree constructed as in Figure \ref{fig:eg1}. Then at the node denoted $S_0u^id^j$ let $f_{i,j}$ be the payoff of the option at that point. Thus at time step $i+j$, and the stock price be $S_{i,j}$. Then we further denote $\Delta_{i,j}$ to be the value of Delta at the node with the stock value $S_0u^id^j$. Then we can estimate Delta at that point using the difference between the two nodes that emanate from our point. Thus we would have,
\begin{align}
\Delta_{i,j} = \frac{f_{i+1,j} - f_{i,j+1}}{S_{i+1,j} - S_{i,j+1}}. \label{deltatrees}
\end{align}
\begin{figure}[ht]
\centering
\begin{tikzpicture}[font=\small, scale=0.4]
\draw (0,5) -- (4,7);
\draw (0,5) -- (4,3);
\draw (4,3) -- (8,1);
\draw (4,3) -- (8,5);
\draw (4,7) -- (8,9);
\draw (4,7) -- (8,5);
\draw[fill] (0,5) circle [radius=0.07];
\draw[fill] (4,7) circle [radius=0.07];
\draw[fill] (4,3) circle [radius=0.07];
\draw[fill] (8,1) circle [radius=0.07];
\draw[fill] (8,5) circle [radius=0.07];
\draw[fill] (8,9) circle [radius=0.07];
\node[below] at (0,5) {$S_0$};
\node[below] at (4,7) {$S_0u$};
\node[below] at (4,3) {$S_0d$};
\node[below] at (8,5) {$S_0ud$};
\node[below] at (8,9) {$S_0u^2$};
\node[below] at (8,1) {$S_0d^2$};
\end{tikzpicture}
\caption{A two step binomial tree with stock prices given at each node }
\label{fig:eg1}
\centering
\end{figure}
Using this method we may also calculate some of the other Greeks. As Gamma is the second derivative of with respect to stock price, we need to go two steps forward to calculate Delta at a given step before the node. Then we may calculate Gamma at the node we are considering. Thus it follows that,
\begin{align*}
\Gamma_{i,j} = \frac{[(f_{i+2,j} - f_{i+1,j+1})/(S_{i+2,j} - S_{i+1,j+1})] + [(f_{i+1,j+1} - f_{i,j+2})/(S_{i+1,j+1} - S_{i,j+2})]}{\frac{1}{2}(S_{i+2,j} - S_{i,j+2})}.
\end{align*}
For Theta if $ud = 1$ is fulfilled then we have that the option price is the same at nodes $f_{i,j}$ and $f_{i+1,j+1}$ the only factor that would have affected the price of the option is time. Thus we may use the finite difference approximation and the definition of Theta to see that,
\begin{align*}
\Theta_{i,j} &= -\frac{f_{i,j} - f_{i+1,j+1}}{2\Delta t} \\
&= \frac{f_{i+1,j+1} - f_{i,j}}{2\Delta t}
\end{align*}
where $\Delta t$ is the time between steps.
If $ud=1$ is not fulfilled then \citet{RUB} found using the relationship between $\Delta, \Theta$ and $\Gamma$ as is (\ref{greeks1}) we must have,
\begin{align}
\Theta_{i,j} = rf_{i,j} - (r-q)S_{i,j}\Delta_{i,j} - \frac{1}{2}\sigma^2 S^2_{i,j}\Gamma_{i,j}. \label{theta}
\end{align}
For Vega we must assume a small change in $\sigma$. Let $f_{i,j}(\sigma) =S_{i,j}$ with $u$ and $d$ as defined by our model. Then $f_{i,j}(\sigma + 0.01) =S_{i,j}$ would be the same with $u$ and $d$ calculated with a very small change. Hence using another finite difference approximation we may calculate Vega in the following way,
\begin{align*}
\nu_{i,j} = \frac{f_{i,j}(\sigma + 0.01) - f_{i,j}(\sigma)}{0.01}.
\end{align*}
Here we have used $0.01$ as a small change however any small value may be taken.
%aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa
\subsubsection{Estimation of Greeks Through Finite Difference Methods}
The idea for approximating the Greeks can be extended from the way they are estimated from trees to finite difference methods in a simple way \citep[P. 430]{JCH}. Once we have calculated all of the values then we may employ exactly the same method as from the tree estimation. For a point $(i,j)$ in the partition, the same formulae hold.
Obviously here we have never stipulated that $ud=1$ need be a condition. As such (\ref{theta}) does not hold. Thus the equations for our Greeks are the same as trees. Though now for us to estimate $\nu$, as $u$ and $d$ are not readily available we must instead calculate the values for $a_j,b_j$ and $c_j$ for a slightly different $\sigma$, say $\sigma + 0.01$. Then our formula for $\nu$ as before holds.
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
\section{Implementation of Numerical Methods and Investigation of their Performance}
In this section we will discuss the implementation of the methods previously described and their performance under certain conditions. To implement these methods MatLab was chosen. This is for a number of reasons; MatLab has many functions built in that we may use, such as the cumulative normal distribution function. Furthermore unlike some programing languages, such as Visual Basic, it allows us to view arrays as vectors and perform vector operations with them. As we will see later on in this section this will be invaluable for calculating option prices using binomial trees or finite difference methods.
%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
\subsection{European Options}
%########################################################################################################
\subsubsection{Black-Scholes}
We have developed an analytic formula for the price of a European option through the solution of the Black-Scholes differential equation. Implementing this is easy as we merely have to calculate $d_1$ and $d_2$, test whether the option is a call or put so that we may use the appropriate formula and finally implement the formula using the command ``normcdf($d_1$)'' to calculate $\Phi(d_1)$ or the required variable we need to take the cumulative normal distribution function we are considering. The code for this may be found in Appendix \ref{App:BS}. As this formula is analytic it is helpful for us to compare the performance of our other methods. This is how this will be used herein.
%########################################################################################################
\subsubsection{Binomial Trees}
\label{imbte}
Before we discuss the implementation of this method it is useful for us to make the following observations.
To begin note that there are some tricks we may employ to make things easier. Notice that for a binomial tree once we have calculated $u$ and $d$ appropriately regardless of the model or parameters entered the method remains the same. As such we, for European options, we use Appendix \ref{App:BTSE} as a hub program that calculates the appropriate $u$ and $d$ and then merely calls Appendix \ref{App:BTECRR} which is our method using the general formula. Alternatively we could calculate this by using our back stepping method as in Appendix \ref{App:BTECRRV}.
Also note that to calculate Greeks we must use vectors. We generate the vector of stock prices at the terminal time, $T$, then shift the vectors so that the values that are needed to calculate the previous nodes' value are aligned. Then we may use our formula and loop back through to obtain our initial option price. Then we can use the values to calculate the value of the Greeks at each node.
We have talked about how Binomial Trees converge to the Black-Scholes formula in the case of European options. We may see this in Figure \ref{fig:Er}. We see that as the number of steps increase the error decreases quickly.
\begin{figure}[ht]
\centering
\begin{subfigure}[b]{0.3\textwidth}
\includegraphics[width=\textwidth]{Error.png}
\caption{A normal plot of the error when using a binomial tree approximation}
\label{fig:Er}
\end{subfigure}%
\quad %add desired spacing between images, e. g. ~, \quad, \qquad, \hfill etc.
%(or a blank line to force the subfigure onto a new line)
\begin{subfigure}[b]{0.3\textwidth}
\includegraphics[width=\textwidth]{LogLogError.png}
\caption{A Log-Log plot of the error when using a binomial tree approximation}
\label{fig:LogEr}
\end{subfigure}
\caption{The error of a binomial tree approximation }\label{fig:errorbt}
\end{figure}
We now investigate the performance of these methods. In Section \ref{greekssection} we discussed the Greek, which are measures of how the price changes with respect to different factors. It therefore makes sense for us to investigate how the price of the option changes with respect to these. This will allow us to evaluate the robustness of our models.
To begin with we consider Delta. We have already discussed how these are estimated in Section \ref{greekstreesss}. The best way to evaluate the performance of binomial trees is to choose a $\varepsilon > 0$ then evaluate the number of steps required to find the price of the option to this level of accuracy. We may do this by comparing the value obtained by a Binomial Tree approach to the value found from the Black-Scholes formula. This allows us to see how the value of Delta affects the efficiency of the method.
For us to be able to do this we need to be able to change Delta. Notice that as $u$ increases, Delta increases. As the numerator increases more than the denominator does. So we affect Delta by increasing $\sigma$ then use the method we just developed to see how Delta affects the number of steps needed. Using this method we generate the graphs as seen in Figures \ref{fig:d1}-\ref{fig:d3}.
\begin{figure}[ht]
\centering
\begin{subfigure}[b]{0.3\textwidth}
\includegraphics[width=\textwidth]{Tressnumberofsteps1.png}
\caption{$r=0.12$, $K=48$, $S_0 = 40$ }
\label{fig:d1}
\end{subfigure}%
\quad %add desired spacing between images, e. g. ~, \quad, \qquad, \hfill etc.
%(or a blank line to force the subfigure onto a new line)
\begin{subfigure}[b]{0.3\textwidth}
\includegraphics[width=\textwidth]{Tressnumberofsteps2.png}
\caption{$r=0.04$, $K=22$, $S_0 = 20$}
\label{fig:d2}
\end{subfigure}
\quad %add desired spacing between images, e. g. ~, \quad, \qquad, \hfill etc.
%(or a blank line to force the subfigure onto a new line)
\begin{subfigure}[b]{0.3\textwidth}
\includegraphics[width=\textwidth]{Tressnumberofsteps3.png}
\caption{$r=0.02$, $K=102$, $S_0 = 100$}
\label{fig:d3}
\end{subfigure}
\caption{The number of steps required to achieve an error tolerance as Delta varies }\label{fig:btdvns}
\end{figure}
Each of these graphs is generated with different $r,S_0$ and $K$. We see that there is no correlation between these and therefore no correlation in general. It is important to understand that we do not affect directly. To cause our change in Delta we use the formula for Delta from Table \ref{tab:gbs}. Notice that from (\ref{bsdii1}) as $\sigma$ decreases $d_1$ increases meaning that $\Phi(d_1)$ increases. This causes Delta to increase. So by changing $\sigma$ is how we affect Delta in Figure \ref{fig:d1} - \ref{fig:d3}.
We see from these graphs that there is no correlation between the value of Delta and the number of steps it requires to calculate the price accurate to a certain value. This shows that the convergence of the method is unaffected by Delta. Computationally speaking this is beneficial. We are able to price an option with the same amount of steps as an equivalent option with a higher Delta. Note that this does not mean that there is no correlation between the price and Delta. Furthermore we may see that as there is no correlation here it follows that there is no correlation between the convergence and Gamma, as Gamma is the derivative of Delta with respect to Theta.
Finally notice that in Figure \ref{fig:d3} all the Deltas are negative. This is because this was done with a put. Notice that from our formula for the estimation of Delta from (\ref{deltatrees}) for a call the numerator is always positive and for a put it is always negative, whilst the denominator is always positive.
%########################################################################################################
\subsubsection{Trinomial Trees}
We shall not consider the performance of these, as they are redundant. This is due to their equivalence to the explicit finite difference method, so all conclusions we draw for the explicit finite difference method apply to trinomial trees. However note that they may be implemented in a very similar way to binomial trees, using vectors and looping to backdate to our initial time, to find the value of the option. The code for this method for European options can be seen in Appendix \ref{App:TTEB}. Notice that again once $u$ and $d$ are calculated the algorithm is identical, so we may use a source program. This can be seen in Appendix \ref{App:TTSE}.
Finally notice that the algorithm for a European and American option is similar, the only difference being the need to see if early exercise is optimal after we backdate a step. Thus these algorithms are near identical as can be seen in Appendices \ref{App:TTEB} and \ref{App:TTAB}. Again, we may use a source program given in Appendix \ref{App:TTSA}.
%########################################################################################################
\subsubsection{Monte Carlo}
For the Monte Carlo methods in Section \ref{mcs1} we have essentially already formed our algorithm. It is merely a case of encoding these using the appropriate functions. This can be seen in Appendix \ref{App:MCE}.
The first question is of convergence. We have remarked that Monte Carlo converges to the analytic solution given by the Black-Scholes formula. In Figure \ref{fig:mce1}-\ref{fig:mce3} we see this behavior. Notice that in Figure \ref{fig:mce1} and \ref{fig:mce2} the convergence is different. This is due to the random component, as the random numbers were different for these sets of simulations, it generates different graphs.
\begin{figure}[ht]
\centering
\begin{subfigure}[b]{0.3\textwidth}
\includegraphics[width=\textwidth]{McError.png}
\caption{Convergence to the value given by the analytic solution for European options}
\label{fig:mce1}
\end{subfigure}%
\quad %add desired spacing between images, e. g. ~, \quad, \qquad, \hfill etc.
%(or a blank line to force the subfigure onto a new line)
\begin{subfigure}[b]{0.3\textwidth}
\includegraphics[width=\textwidth]{McError2.png}
\caption{Convergence to the value given by Black-Scholes with different random numbers}
\label{fig:mce2}
\end{subfigure}
\quad %add desired spacing between images, e. g. ~, \quad, \qquad, \hfill etc.
%(or a blank line to force the subfigure onto a new line)
\begin{subfigure}[b]{0.3\textwidth}
\includegraphics[width=\textwidth]{McError3.png}
\caption{The error between the Black-Scholes value and Monte Carlo estimate}
\label{fig:mce3}
\end{subfigure}
\caption{Convergence of Monte Carlo to the Analytic solution for European options }\label{fig:mcerrrror}
\end{figure}
Now consider antithetic variates. When implementing this we use a small trick. As described in Section \ref{mcs32} we may generate secondary paths by choosing the random numbers to be the negative of the random numbers for the usual set of paths. This saves us a lot computational time and further relaxes the condition that $\textnormal{Cov}[Y_i,\hat{Y}_i]<0$.
In Figure \ref{fig:mcee1} we see the convergence of Monte Carlo to the option price. Note that it is erratic, due to the random numbers; this is why we sample a large amount of times. Figure \ref{fig:mcee2} shows the error of each of the two methods as the number of replications increases. We see here how much the antithetic variance reduction technique reduces our error, this is even more greatly seen in Figure \ref{fig:mcee3}. Here we see how many replications are required to gain a certain error tolerance. Notice that the normal method requires a thousand or more replications to find the correct price at the lower error estimates. Helpfully the antithetic variate technique does not require as many paths, in one case requiring one ninth of the paths of the normal method.
\begin{figure}[ht]
\centering
\begin{subfigure}[b]{0.3\textwidth}
\includegraphics[width=\textwidth]{McAvCovergence.png}
\caption{$r=0.04$, $K=22$ and $S_0 = 20$}
\label{fig:mcee1}
\end{subfigure}
\quad %add desired spacing between images, e. g. ~, \quad, \qquad, \hfill etc.
%(or a blank line to force the subfigure onto a new line)
\begin{subfigure}[b]{0.3\textwidth}
\includegraphics[width=\textwidth]{McAVErrorCompare.png}
\caption{$r=0.12$, $K=48$ and $S_0 = 40$ }
\label{fig:mcee2}
\end{subfigure}%
\quad %add desired spacing between images, e. g. ~, \quad, \qquad, \hfill etc.
%(or a blank line to force the subfigure onto a new line)
\begin{subfigure}[b]{0.3\textwidth}
\includegraphics[width=\textwidth]{ErrorNumberOfsimulationsMc.png}
\caption{$r=0.04$, $K=22$ and $S_0 = 20$}
\label{fig:mcee3}
\end{subfigure}
\caption{Convergence of Monte Carlo methods to the Analytic solution for European options }\label{fig:mcerooror2}
\end{figure}
As stated earlier the requirement for the covariance being less than zero is more relaxed as we do not simulate $2n$ paths for the variance reduced method due to our trick. This is seen as the covariance in the graph above is only slightly smaller than zero with the covariance being negative and of order $10^{-3}$ in general. This shows how beneficial our method can be.
As we mentioned earlier Monte Carlo is not particularly competitive for European options. This can be seen as the error in Figure \ref{fig:mcee2} once below $0.1$ is erratic and in general much higher than in the binomial tree approximation given in Figure \ref{fig:Er}. The key difference here is that we have need to simulate $500-1000$ paths to achieve this in Monte Carlo, and then only the variance reduced method achieved this error tolerance and not for the normal method, however we achieve better than this with a $50$ step binomial tree. This shows us how bad this method is as both of these programs are $\mathcal{O}(n)$ as can be seen in their respective codes, Appendices \ref{App:MCE}, \ref{App:MCEA} and \ref{App:BTGE}. The reason is we needed the Greeks when we generated our binomial tree graph hence the use of this code. We now see how this method compares to finite difference methods.
%########################################################################################################
\subsubsection{Finite Difference Methods}
\label{efdmi}
For the finite difference method we use different approaches for the implicit and explicit methods. For the explicit method, due to the nature of it being equivalent to the trinomial tree approach we may implement using a similar method for the trinomial trees; vector operations. We need only add the appropriate boundary conditions in.
For the implicit method we need another approach. We may use the Thomas Algorithm to help us implement this problem. This has the effect of solving the equations in $\mathcal{O}({MN})$ instead of $\mathcal{O}(M^2N^2)$ required by Gaussian elimination. To see how we implement this consider the Thomas Algorithm. This is a method for solving tridiagonal matrices.
Given the tridiagonal matrix as $A$ in Section \ref{ifdms} multiplied by $X = (x_1,x_2,\dots,x_n)^T$ and a vector of knowns, $D = (d_1, \dots ,d_n)^T$, we perform the following transformation,
\begin{align*}
c'_i &=
\begin{cases}
\begin{array}{lcl}
\cfrac{c_i}{b_i} & ; & i = 1 \\
\cfrac{c_i}{b_i - a_i c'_{i - 1}} & ; & i = 2, 3, \dots, n-1 \\
\end{array}
\end{cases} \\
d'_i &=
\begin{cases}
\begin{array}{lcl}
\cfrac{d_i}{b_i} & ; & i = 1 \\
\cfrac{d_i - a_i d'_{i - 1}}{b_i - a_i c'_{i - 1}} & ; & i = 2, 3, \dots, n. \\
\end{array}
\end{cases}
\end{align*}
This has the bonus of eliminating the $a_i$ so that $a_i = 0 \ \forall i$ and transforming the $b_i$ so that $b_i = 1 \ \forall i$. Thus we may then solve this system of equations through backwards substitution.
\begin{align*}
x_n &= d'_n\\
x_i &= d'_i - c'_i x_{i + 1} \qquad ; \ i = n - 1, n - 2, \ldots, 1
\end{align*}
We apply this to our system as seen in Section \ref{ifdms} by saying $F_i = X$ and $D = F_{i+1}+B_i$. We may now implement this method again using vector operations to reduce computing time. We do this for each iteration step, going backwards until the initial value found at the start is found.
We need to compare the performance of these two methods separately. This is because the explicit method, as seen in Section \ref{compttress}, is equivalent to trinomial trees. The implicit method is very different to trinomial tress using the Thomas Algorithm to solve this. This means that due to the large amount of work required for the implicit method, that is not used in the explicit method, we expect this to require more time to achieve the same error tolerance.
\begin{figure}[ht]
\centering
\begin{subfigure}[b]{0.3\textwidth}
\includegraphics[width=\textwidth]{FDMError.png}
\caption{Convergence to the value given by the analytic solution for European options}
\label{fig:fdm1}
\end{subfigure}%
\quad %add desired spacing between images, e. g. ~, \quad, \qquad, \hfill etc.
%(or a blank line to force the subfigure onto a new line)
\begin{subfigure}[b]{0.3\textwidth}
\includegraphics[width=\textwidth]{FMDsNumberOfParitionslarge.png}
\caption{Number of partitions for a given error tolerance with large $M$ and $n$ and with $M=N$}
\label{fig:fdm2}
\end{subfigure}
\quad %add desired spacing between images, e. g. ~, \quad, \qquad, \hfill etc.
%(or a blank line to force the subfigure onto a new line)
\begin{subfigure}[b]{0.3\textwidth}
\includegraphics[width=\textwidth]{FMDsNumberOfParitionsSmall.png}
\caption{Number of partitions for a given error tolerance with small $M$ and $n$ and with $M=N$}
\label{fig:fdm3}
\end{subfigure}
\caption{Convergence of finite difference methods to Convergence of Monte Carlo to the Analytic solution for European options }\label{fig:animals}
\end{figure}
%his behavouir is characterized in Figures \ref{fig:fdm2} and \ref{fig:fdm3}. We see that for larger value of error the convergence is much quicker. However as the error tolerance decreases the implicit method requires much high $M$ and $N$ to achieve the same tolerance.
We can see this as if we time each method in Appendices \ref{App:FDMIE} and \ref{App:FDMEE} for say $N=M=50$ the implicit method requires twice as much time as the explicit method to achieve the same error tolerance level. We also see that in Figures \ref{fig:fdm2} and \ref{fig:fdm3} the number of rectangles we partition into (choosing $M=N$ so that the $x$-axis is the square root of the number of rectangles used) to gain a given error tolerance is vastly more for lower error tolerances than the explicit method. This is due to the implicit method being $\mathcal{O}(\Delta t, \Delta S^2)$ while the explicit method is $\mathcal{O}(\Delta t^2, \Delta S^2)$.
This shows the behavior we see in Figures \ref{fig:fdm2} and \ref{fig:fdm3}. In the second of these we see for $\Delta t$ and $\Delta S$ larger these have similar convergence, however as $\Delta t \rightarrow 0$ the number of rectangles needed increases greatly. As in the graphs we have taken $M=N$ this causes the vast increase as our error tolerance decreases.
Lastly see that these orders of convergence show that error of the explicit method is less than the error of the implicit method for the same given $M$ and $N$. This can be seen in Figure \ref{fig:fdm1}, showing that in general as we increases the number of rectangles the error of the explicit method has a lower error than the implicit method. It is important to see that all of these result hold for American options which can be seen in the next section.
%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
\subsection{American Options}
%########################################################################################################
We shall only investigate the performance of the method we have discuses for American options. While some that we have discussed may be adapted to American options, we have not looked at how these may be adapted and as such shall not consider them here. We only considering finite difference methods and binomial trees for which we have developed models for American options.
\subsubsection{Binomial Trees}
For an American option our implementation is near identical to the vector implementation of European options seen in Appendix \ref{App:BTECRRV}. As can be seen in Appendix\ref{App:BTSAV} the only minor alteration is having to, after calculating the value at the previous nodes, take the maximum of these and the possible payoff for early exercise. Further again we may use a source program to allow for more user input as seen in Appendix \ref{App:BTSA}. The Greek calculation is again near identical with the only change being taking the maximum of payoff at the current node and the backward calculated value as seen in Appendix \ref{App:BTGA}.
The model is developed in the same way, meaning that all of the conclusions we discussed in Section \ref{imbte} still hold. We shall therefore investigate how this method performs when different Greeks change.
As when we explored Delta in Section \ref{imbte} we may not change Theta directly. Note that for European options from Table \ref{tab:gbs} we see that as $\sigma$ increases Theta decreases. Thus we shall take the same approach changing $\sigma$ slowly to observe the behavior of the price of the option as Theta changes. We first examine the price of American options against the price of European options.
\begin{figure}[ht]
\centering
\begin{subfigure}[b]{0.3\textwidth}
\includegraphics[width=\textwidth]{AmeiricanOptionsVsEuropean.png}
\caption{Option price of equivalent American and European options}
\label{fig:bta1}
\end{subfigure}%
\quad %add desired spacing between images, e. g. ~, \quad, \qquad, \hfill etc.
%(or a blank line to force the subfigure onto a new line)
\begin{subfigure}[b]{0.3\textwidth}
\includegraphics[width=\textwidth]{DeltaAmericanOptions.png}
\caption{Option price compared to the value of Delta}
\label{fig:bta2}
\end{subfigure}
\quad %add desired spacing between images, e. g. ~, \quad, \qquad, \hfill etc.
%(or a blank line to force the subfigure onto a new line)
\begin{subfigure}[b]{0.3\textwidth}
\includegraphics[width=\textwidth]{ThetaAmericanOptions.png}
\caption{Option Price compared to the value of Theta}
\label{fig:bta3}
\end{subfigure}
\caption{American options compared to the Greeks and European options}\label{fig:EvAA}
\end{figure}
From Figure \ref{fig:bta1} we see that the price of an American option is always greater than the equivalent European option. The above is obviously only for a given $r, S_0, T$ and $K$, however this is a trend we expect to see. The reason for this is that an American option is always more valuable than a European option. This is because the owner of an American option has every chance the owner of a European option has to exercise and more. This fairly obviously makes the option more valuable.
We see that there is a definite correlation between Theta and Delta and the price of the option. In Figure we observe that \ref{fig:bta2} the higher Delta is the higher the price of the option. This from the same line of reasoning as to why we expect the volatility to increase the value of an option as discussed in Section \ref{fasp1}. This is because as the rate the stock price changes increases if the stock price goes up we are likely to make a profit. However if it decreases severely our loses are bounded by the price of the option. Hence we want this to be as big as possible. Further note that the gradient of the line is increasing. This means that as Gamma increases the stock price increases as well. The reason for this follows similar logic to that of Gamma.
Furthermore we see in Figure \ref{fig:bta3} that as Theta increases the price of the option decreases. This is because Theta is a measure of the time decay of the option price. As such we wish for this to be as small as possible; the faster the price decays the less the option is worth.
Again note that from our formula for the approximation of these Greeks we see that the above is for calls. As for puts both Theta and Delta would be negative.
\begin{figure}[ht]
\centering
\begin{subfigure}[b]{0.3\textwidth}
\includegraphics[width=\textwidth]{AmericanOptionsVegaB.png}
\caption{The price of an American option compared to the value of Vega}
\label{fig:av1}
\end{subfigure}%
\quad %add desired spacing between images, e. g. ~, \quad, \qquad, \hfill etc.
%(or a blank line to force the subfigure onto a new line)
\begin{subfigure}[b]{0.3\textwidth}
\includegraphics[width=\textwidth]{EuropeanOptionVegaBT.png}
\caption{The price of a European option compared to the value of Vega}
\label{fig:ev1}
\end{subfigure}
\caption{American and European option as Vega varies, if early exercise is optimal }\label{fig:Agdgd}
\end{figure}
In Figures \ref{fig:av1} and \ref{fig:ev1} we see that there is a general correlation between the European option and Vega. In general we have that as Vega increases the price of the option decreases. This is because if $\sigma$ is changing rapidly it does not guarantee that it will be larger, which will increase our price.
The lack of correlation seen in Figure \ref{fig:av1} is not necessarily a trend as seen in Figures \ref{fig:av2} and \ref{fig:ev2}. Here the relationship between Vega and the option price is identical for European and American options. The lack of correlation seen in Figure \ref{fig:av1} is due to early exercise, as it is optimal here, however in Figure \ref{fig:av2} it is not leading to the presence of correlation. As when we exercise early $\sigma$ plays a less significant role, due to the option acting over a shorter time. This is as when we exercise early for a binomial tree we multiply by $u$ and $d$ less, meaning $\sigma$ has less impact when we exercise early. Thus we expect that Vega would have less impact when we exercise early, resulting in our lack of correlation. Alternatively when we do not exercise early, the relationship is near identical to European options.
\begin{figure}[ht]
\centering
\begin{subfigure}[b]{0.3\textwidth}
\includegraphics[width=\textwidth]{AmericanVegaNoEarlyB.png}
\caption{American option compared to Vega when early exercise is never optimal}
\label{fig:av2}
\end{subfigure}%
\quad %add desired spacing between images, e. g. ~, \quad, \qquad, \hfill etc.
%(or a blank line to force the subfigure onto a new line)
\begin{subfigure}[b]{0.3\textwidth}
\includegraphics[width=\textwidth]{EuropeanOptionVegaEarly.png}
\caption{European option compared to Vega where early exercise is no possible}
\label{fig:ev2}
\end{subfigure}
\caption{American and European option as Vega varies, if early exercise is never optimal }\label{fig:Agdgd}
\end{figure}
%########################################################################################################
\subsubsection{Finite Difference Methods}
This is identical in both cases to the European option implementation. Again the only difference between the American and European case is we need to consider early exercise, this can be seen in Appendices \ref{App:FDMIA} and \ref{App:FDMEA}. In the case of the implicit method we need compare the values of $F_i$ we have found to the values that would be given by early exercise at this point.
In the case of the explicit method this is essentially the same. Though in the formation of the problem we never defined a $F_i$ it is helpful to do so in implementing this. As explained in the European case it is advantageous to do this using vectors. The method is identical for the American case however we now only do the comparison as above.
\begin{figure}[ht]
\centering
\begin{subfigure}[b]{0.3\textwidth}
\includegraphics[width=\textwidth]{AvEA.png}
\caption{The value of an American option compared to a European option approximated using the implicit method}
\label{fig:ave1}
\end{subfigure}%
\quad %add desired spacing between images, e. g. ~, \quad, \qquad, \hfill etc.
%(or a blank line to force the subfigure onto a new line)
\begin{subfigure}[b]{0.3\textwidth}
\includegraphics[width=\textwidth]{AvEAExplicit.png}
\caption{The value of an American option compared to a European option approximated using the explicit method}
\label{fig:ave2}
\end{subfigure}
\caption{Price of an American option compared to a European option using different finite difference approximations }\label{fig:Agdgd}
\end{figure}
All of the conclusions we drew from Section \ref{efdmi} still hold. We may see in Figure \ref{fig:ave1} for the implicit method we still have that American options are more valuable than European options. This is again echoed in Figure \ref{fig:ave2} for the explicit method.
Furthermore all of our analysis for European options hold. Note here that we would observe the same behavior w.r.t. Vega here as before due to early exercise.
As these options are so similar to European options, we see that there are not many new conclusions we can draw from these. We now consider a very different type of option; Asian options.
%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
\subsection{Asian Options}
As with American options we only consider here the numerical method that we have applied to Asian options in our development in the theory. Thus we only consider Monte Carlo, so we shall see how the different types of Monte Carlo we developed perform.
%########################################################################################################
\subsubsection{Monte Carlo}
As with the European case the algorithms we are using are known from Section \ref{fdfdfdf}. The implementation the normal, antithetic and control variate techniques can be seen in Appendices \ref{App:MCA}, \ref{App:MCAA} and \ref{App:MCACV} respectively. Thus we need only consider how these perform.
\begin{figure}[ht]
\centering
\begin{subfigure}[b]{0.3\textwidth}
\includegraphics[width=\textwidth]{MCANvA.png}
\caption{Option price of equivalent American and European options}
\label{fig:mcanva}
\end{subfigure}%
\quad %add desired spacing between images, e. g. ~, \quad, \qquad, \hfill etc.
%(or a blank line to force the subfigure onto a new line)
\begin{subfigure}[b]{0.3\textwidth}
\includegraphics[width=\textwidth]{MCANvC.png}
\caption{Option price compared to the value of Delta}
\label{fig:mcanvc}
\end{subfigure}
\quad %add desired spacing between images, e. g. ~, \quad, \qquad, \hfill etc.
%(or a blank line to force the subfigure onto a new line)
\begin{subfigure}[b]{0.3\textwidth}
\includegraphics[width=\textwidth]{MCAAvC.png}
\caption{Option Price compared to the value of Theta}
\label{fig:mcaavc}
\end{subfigure}
\caption{Different Variance reduction techniques and their affect on a Monte Carlo approximation }\label{fig:Agdgd}
\end{figure}
We see in Figures \ref{fig:mcanva} and \ref{fig:mcanvc} that both of these variance reduction techniques improve on the normal Monte Carlo estimates. So the natural question is that of which is best. If one is not always optimal, under what circumstances is each better in.
In Figure \ref{fig:mcaavc} we see that the antithetic variance reduction techniques is much better in the given case. The way in which the antithetic values technique reduces variance is always the same and always improves the precision. However, consider if the price of an Asian option is very close to that of a geometric average option. In this case the control variate will bring the price much closer much quicker as the error correction is nearly exactly accurate. Thus if this was the case we would expect to see the control variate technique out perform the antithetic values technique.
Now we seek to examine the performance of the estimation of Greeks. In Section \ref{hghghghghf} we considered that convergence of two different methods of estimating Delta.
We see in Figure \ref{fig:mcdf} the value of Delta compared to $n$, the number of paths. We showed in Section \ref{hghghghghf} that as the number of paths increases Delta should converge approximately like $n^{-1/4}$. We see this Figure \ref{fig:mcdf} as the estimations of Delta seem to follow roughly the same behavior.
Furthermore we stipulated that for the central difference method the convergence should behave like $n^{-1/6}$. This is demonstrated in Figure \ref{fig:mcdc} where we see the estimation seems to follow the same path roughly. Note that in Figures \ref{fig:mcdf} and \ref{fig:mcdc} are generated with antithetic values. This is as for the option this graph was generated from, the method had reduced variance more showing the correlation more clearly.
\begin{figure}[ht]
\centering
\begin{subfigure}[b]{0.3\textwidth}
\includegraphics[width=\textwidth]{DetlaForward.png}
\caption{Convergence of Delta using a forward difference approximation}
\label{fig:mcdf}
\end{subfigure}%
\quad %add desired spacing between images, e. g. ~, \quad, \qquad, \hfill etc.
%(or a blank line to force the subfigure onto a new line)
\begin{subfigure}[b]{0.3\textwidth}
\includegraphics[width=\textwidth]{DeltaCentral.png}
\caption{Convergence of Delta using a central difference approximation}
\label{fig:mcdc}
\end{subfigure}
\caption{Convergence of Delta using as estimated from a Monte Carlo approximation }\label{fig:Agdgd}
\end{figure}
Here ends our discussion of numerical and analytic methods in option pricing. Having developed our key ideas such as; the Black-Scholes differential equation, Black-Scholes formula for European options, Mont Carlo methods, binomial trees, finite difference methods and the Greeks and their numerical estimation we have implemented these and discussed how they perform under certain circumstances. This has led us to some interesting insights as to which methods are best in given scenarios and how we may improve them. We have further seen the convergence of many methods and convergence of approximations of Greeks with these methods. While our study or investigation has by no means be exhaustive it has given us an interesting insight into the field of option pricing
\newpage
\begin{thebibliography}{999}
%not used
\bibitem[Black and Scholes(1974)]{BS}{Black, F. and Scholes, M. (1973). The Pricing of Options and Corporate Liabilities. \textit{The Journal of Political Economy}. 81 (May-June., 1973), p.637-654.}%notused USEIT
\bibitem[Boyle(1988)]{BO}{Boyle, P.P. (1988). A Lattice Framework for Option Pricing with Two State Variables. \textit{The Journal of Financial and Quantitative Analysis}. 23 (March 1988), p. 1-12.}
\bibitem[Chriss(1997)]{NAC}{Chriss A. N. (1997) \textit{Black-Scholes and Beyond}. 2nd Ed. McGraw-Hill Professional.}
\bibitem[Chung(1974)]{KC}{Chung, K., (1974) \textit{Elementary probability theory with stochastic processes}. Springer-Verlag.}
\bibitem[Cox \textit{et al.}(1979)]{CRR}{Cox, J., Ross, S. and Rubinstein, M. (1979) Option Pricing: A Simplified Approach. \textit{Journal of Financial Economics}. 7 (September 1979). p.229-263}
\bibitem[Durrett(2010)]{DR}{Durrett, R. 2010., \textit{Probability.} [online]. Cambridge University Press. Available from:$<$http://www.myilibrary.com?ID$=281866>$ $6$ November 2014} % think this is for characteristic functions stuff but unsure.
\bibitem[Gardiner(1985)]{GSS}{Gardiner, C.W., \textit{Handbook of stochastic methods} 2nd Ed. Springer}
%not used
\bibitem[Geman and Yor(1993)]{GY}{German H. and Yor, M., Bessel processes, Asian options and perpetuities, \textit{Mathematical Fiance}. 3. (December 1993) p.349-375.}
\bibitem[Gilli(2011)]{MG}{Gilli M.,\textit{Numerical Methods and Optimization in Finance.} 1 Edition. Academic Press.}
\bibitem[Glasserman(2003)]{PG}{Glasserman, P., (2003) \textit{Monte Carlo Methods in Financial Engineering.} 3rd Ed. Springer.}
%not used
%not used
\bibitem[Hull(2005)]{JCH}{Hull J. C. (2005) \textit{Options, Futures, and Other Derivatives}. 6th Ed. Pearson.}
\bibitem[Broadie and Glasserman(2005)]{BG}{Broadie, M., and Glasserman, P. (1997) \textit{Handbook of Risk Management and Analysis}. Wiley, Chichester, England}
%not used
\bibitem[Jarrow and Rudd(1983)]{JRB}{Jarrow R., and Rudd A., (1983) \textit{Option pricing}. 1st Ed. Richard D. Irwin, Inc.}
\bibitem[Jarrow and Rudd(1986)]{JRP}{Jarrow, R., and Rudd, A., (1986). Option Pricing \textit{The Journal of Banking and Finance.} 10 (March 1986), p. 157-161.}
\bibitem[Merton(1973)]{RM}{Merton, R., Theory of Rational Option Pricing, \textit{Bell Journal of Economics \& Management}. 4. (June 1973) p.141-183.}
\bibitem[Rubinstein(1994)]{RUB}{Rubinstein, M.,(1994). Implied Binomial Trees. \textit{Journal of Finance }. 49, pp. 771-818.}
\bibitem[Sinclair(2010)]{ES}{Sinclair E. (2010) \textit{Option Trading: Pricing and Volatility Strategies and Techniques}. 1st Ed. John Wiley \& Sons.}
\bibitem[Smith(1965)]{GDS}{Smith G. D. (1965) \textit{Numerical solutions of partial differential equations}. 3rd Ed. Oxford university press.}
\end{thebibliography}
\newpage
\appendixtitleon
\appendixtitletocon
\begin{appendices}
%\section{Binomial Trees MatLab Code (American CRR)} \label{App:BTACRR}
%\lstinputlisting{BinomialTreesACrr.m}
\section{Black-Scholes}
\subsection{Black-Scholes Model MatLab Code}\label{App:BS}
\lstinputlisting{BlackScholesE.m}
\section{Monte Carlo}
\subsection{European}
\subsubsection{Monte Carlo MatLab Code} \label{App:MCE}
\lstinputlisting{MonteCarloE.m}
\subsubsection{Monte Carlo Antithetic MatLab Code} \label{App:MCEA}
\lstinputlisting{MonteCarloEAnithetic.m}
\subsection{Asian}
\subsubsection{Monte Carlo MatLab Code} \label{App:MCA}
\lstinputlisting{MonteCarloA.m}
\subsubsection{Monte Carlo Antithetic MatLab Code} \label{App:MCAA}
\lstinputlisting{MonteCarloAAntithetic.m}
\subsubsection{Monte Carlo Control Variate MatLab Code} \label{App:MCACV}
\lstinputlisting{MonteCarloAControlVariate.m}
\section{Binomial Trees}
\subsection{European}
\subsubsection{Binomial Trees MatLab Code Using General Formula MatLab Code} \label{App:BTECRR}
\lstinputlisting{BinomialTreesCRRE.m}
\subsubsection{Binomial Trees MatLab Code Using Vectors MatLab Code} \label{App:BTECRRV}
\lstinputlisting{BinomialTreesEVectors.m}
\subsubsection{Source Program Binomial Trees MatLab Code} \label{App:BTSE}
\lstinputlisting{BinomialTrees.m}
\subsubsection{Numerical Estimation of Greeks from a Binomial Tree MatLab Code} \label{App:BTGE}
\lstinputlisting{BinomialTreesCRREGreeks.m}
\subsection{American}
\subsubsection{Binomial Trees MatLab Code MatLab Code} \label{App:BTSAV}
\lstinputlisting{BinomialTreesACrr.m}
\subsubsection{Source Program Binomial Trees MatLab Code} \label{App:BTSA}
\lstinputlisting{BinomialTreesA.m}
\subsubsection{Numerical Estimation of Greeks from a Binomial Tree MatLab Code} \label{App:BTGA}
\lstinputlisting{BinomialTreesAGreeks.m}
\section{Trinomial Trees}
\subsection{European}
\subsubsection{Trinomial Trees MatLab Code (Boyle)} \label{App:TTEB}
\lstinputlisting{TrinomialTreesEBoyle.m}
\subsubsection{Trinomial Trees Source Program (Boyle)} \label{App:TTSE}
\lstinputlisting{TrinomialTreesE.m}
\subsection{American}
\subsubsection{Trinomial Trees MatLab Code MatLab Code (Boyle)} \label{App:TTAB}
\lstinputlisting{TrinomialTreesABoyle.m}
\subsubsection{Trinomial Trees Source Code (Boyle)} \label{App:TTSA}
\lstinputlisting{TrinomialTreesA.m}
\section{Finite Difference Methods}
\subsection{European}
\subsubsection{Implicit Finite Difference Method MatLab Code} \label{App:FDMIE}
\lstinputlisting{FiniteDiffereceMethodsEI.m}
\subsubsection{Explicit Finite Difference Method MatLab Code} \label{App:FDMEE}
\lstinputlisting{FiniteDifferenceMethodsEE.m}
\subsection{American}
\subsubsection{Implicit Finite Difference Method MatLab Code} \label{App:FDMIA}
\lstinputlisting{FiniteDiffereceMethodsAI.m}
\subsubsection{Explicit Finite Difference Method MatLab Code} \label{App:FDMEA}
\lstinputlisting{FiniteDifferenceMethodsAE.m}
%\section{Binomial Trees MatLab Code (American CRR with a drift)} \label{App:BTACRRD}
%\lstinputlisting{BinomialTreesACrrDrift.m}
%\section{Binomial Trees MatLab Code (American CRR different parameters)} \label{App:BTCRR}
%\lstinputlisting{BinomialTreesACrrParameters.m}
%\section{Binomial Trees MatLab Code (European CRR with drift)} \label{App:BTECRRD}
%\lstinputlisting{BinomialTreesCRREDrift.m}
%\section{Binomial Trees MatLab Code (European CRR with drift (Different parameters))} \label{App:q}
%\lstinputlisting{BinomialTreesCRREDiffParameters.m}
%\section{Trinomial Trees MatLab Code (American Boyle different parameters)} \label{App:TTABD}
%\lstinputlisting{TrinomialTreesABoyleParameters.m}
%\section{Trinomial Trees MatLab Code (European Boyle different parameters)} \label{App:TTEBD}
%\lstinputlisting{TrinomialTreesEBoyleParameters.m}
%%
\end{appendices}
\end{document}
