\documentclass[a4paper]{article}
\usepackage{amsfonts,amssymb,amsmath,amsthm,latexsym,amsbsy,graphicx,float,hyperref,ifthen,color}
\usepackage{makecourse}
\usepackage{tikz}
\allowdisplaybreaks

\newcommand{\mvs}{\vspace{2.5cm}\\}
\newcommand{\mvv}{\vspace{4.0cm}\\}
\newcommand{\expt}{\textrm{E}}
\newcommand{\pr}{\textrm{Pr}}
\newcommand{\logn}{\textrm{ln}}
\newcommand{\var}{\textrm{Var}}
\newcommand{\RR}{\mathrm{I\!R\!}}
\newcommand{\R}{$\mathsf{R\;}$}

\setcounter{section}{4}
\begin{document}
\section{An Introduction to It\^o Calculus}
In this final Section, we introduce a further type of stochastic process known as an It\^o process. We have already met two such processes -- Brownian motion and geometric Brownian motion. Here, we will formulate such models as stochastic differential equations (S.D.E.s) and develop a framework for working with these processes.

\subsection{The It\^o Integral}
Assume that stock price follows a G.B.M. with parameters $\mu$ and $\sigma^2$. Suppose that we are interested in the behaviour of stock price in the interval $[t,t+\Delta t]$. Using the relationship between G.B.M. and B.M., we can express the change in the logarithm of an assets price as
\[
\log(S_{t+\Delta t})-\log(S_{t})=(\mu-\sigma^{2}/2)\Delta t+\sigma (W_{t+\Delta t}-W_{t}) 
\]
where $W_{t}$ is standard Brownian motion. Assuming zero volatility (i.e. $\sigma^{2}=0$) we have 
\begin{eqnarray*}
\frac{\log(S_{t+\Delta t})-\log(S_{t})}{\Delta t}&=& \mu \\
\Rightarrow \frac{d\log(S_{t})}{dt}&=& \mu \qquad \text{upon taking $\Delta t\to 0$.}
\end{eqnarray*} 
In other words, the change in the logarithm of an asset's price is constant. Note that solving the above ordinary differential equation (ODE) system 
gives 
\[
S_t = S_0 e^{\mu t}
\]
which we may have anticipated, since this is also the expectation of $S_t$ when $\sigma$ is non zero. What if we don't assume zero volatility? Clearly,
\[
\frac{\log(S_{t+\Delta t})-\log(S_{t})}{\Delta t}= \mu-\sigma^{2}/2 +\sigma\frac{W_{t+\Delta t}-W_{t}}{\Delta t}\, .
\]
We would like to take the limit (as $\Delta t\to 0$) of both sides and write $dW_{t}/dt$ on the RHS. However, although sample paths of $W_{t}$ are continuous functions of $t$, they are almost everywhere non-differentiable. To see this, note the difference quotient
\[
\frac{W_{t+\Delta t}-W_{t}}{\Delta t}
\]
is not bounded in the neighbourhood of $t$. Intuitively, the variance of the difference quotient explodes as $\Delta t\to 0$. Consequently, $dW_{t}/dt$ has no meaning in the usual sense. In particular, the Riemann integral 
\[
\int_{0}^{t}g(\tau)\frac{dW_{\tau}}{d\tau}d\tau
\]
has no meaning since $dW_{\tau}/d\tau$ is not defined. Therefore, a new definition of a stochastic integral is needed.  

%Let $g(t,X_{t})$ be a function of the random variable $X_{t}$, where $t\in[a,b]$ and $\{X_{t},t\geq 0\}$ is a stochastic process. We define the It\^o stochastic integral for random functions $g$ satisfying
%\begin{equation}\label{ass1}
%\int_{a}^{b}\expt\left(g^{2}(t,X_{t})\right)\, dt < \infty \, .
%\end{equation}
Let $g(t)$ be a function of time, where $t\in[a,b]$. We define the It\^o stochastic integral for functions $g$ satisfying
\[
\int_{a}^{b}g^{2}(t)\, dt < \infty \, .
\]

\subsubsection*{Definition 5.1 (It\^o Stochastic Integral)}
%Partition $[a,b]$ as $a=t_{0}<t_{1}<\ldots <t_{n-1}<t_{n}=b$ and set $\Delta t=t_{i+1}-t_{i}=(b-a)/n$, $\Delta W_{t_{i}}=W_{t_{i+1}}-W_{t_{i}}$. The It\^o stochastic integral of $g$ is defined as 
%\begin{equation}\label{ito}
%\int_{a}^{b}g(t,X_{t})\, dW_{t}=l.i.m._{n\to\infty}\sum_{i=0}^{n-1}g(t_{i},X_{t_{i}})\Delta W_{t_{i}}\, ,
%\end{equation}
%where $l.i.m.$ denotes mean square convergence. In other words, if $G_{n}=\sum_{i=0}^{n-1}g(t_{i},X_{t_{i}})\Delta W_{t_{i}}$ and $\mathcal{I}=\int_{a}^{b}g(t,X_{t})\, dW_{t}$ then $l.i.m._{n\to\infty}G_{n}=\mathcal{I}$ means that
%\[
%\lim_{n\to \infty} \textrm{E}[(G_{n}-\mathcal{I})^{2}] = 0 \, .
%\]
Partition $[a,b]$ as $a=t_{0}<t_{1}<\ldots <t_{n-1}<t_{n}=b$ and set $\Delta t=t_{i+1}-t_{i}=(b-a)/n$, $\Delta W_{t_{i}}=W_{t_{i+1}}-W_{t_{i}}$. The It\^o stochastic integral of $g$ is defined as 
\[
\int_{a}^{b}g(t)\, dW_{t}=l.i.m._{n\to\infty}\sum_{i=0}^{n-1}g(t_{i})\Delta W_{t_{i}}\, ,
\]
where $l.i.m.$ denotes mean square convergence. In other words, if $G_{n}=\sum_{i=0}^{n-1}g(t_{i})\Delta W_{t_{i}}$ and $\mathcal{I}=\int_{a}^{b}g(t)\, dW_{t}$ then $l.i.m._{n\to\infty}G_{n}=\mathcal{I}$ means that
\[
\lim_{n\to \infty} \textrm{E}[(G_{n}-\mathcal{I})^{2}] = 0 \, .
\]

\subsubsection*{Simple It\^o integrals}

Simple It\^o stochastic integrals can be verified directly from the Definition~5.1. For example taking $g(t)=1$ gives 
\begin{eqnarray}
\int_{a}^{b}dW_t &=& l.i.m._{n\to\infty} \sum_{i=0}^{n-1}\Delta W_{t_{i}} \nonumber \\
		&=& l.i.m._{n\to\infty}(W_{t_{1}}-W_{t_{0}}+W_{t_{2}}-W_{t_{1}}+\ldots +W_{t_{n}}-W_{t_{n-1}}) \nonumber \\
		&=& l.i.m._{n\to\infty}(W_{t_{n}}-W_{t_{0}}) \nonumber \\
		&=& W_{b}-W_{a} \, .	\nonumber
\end{eqnarray}
where the last line follows since $t_0=a$ and $t_n=b$ irrespective of the partition. Similarly, taking $a=0$ and $b=t$ gives 
\[
\int_{0}^{t}dW_{\tau} = W_t - W_0 = W_t
\]
since $W_0=0$ and we used $\tau$ as a dummy variable for time since $t$ is the upper limit of the integral. 
%For any well defined random function $G(t,W_{t})$ satisfying (\ref{ass1}),
%\[
%\int_{a}^{b}dG(t,W_{t})=G(b,W_{b})-G(a,W_{a})\, .
%\]
Now take $g(t)=k$ for some constant $k$ to see that
\[
\int_{0}^{t}kdW_{\tau} = k(W_t - W_0) = k W_t.
\]

\subsubsection*{Comments}
\begin{itemize}
\item Solutions of It\^o integrals are stochastic processes.
\item Note that a general It\^o stochastic integral of the form 
\[
X_{t}=X_0+\int_{0}^{t}g(\tau)dW_{\tau}
\]
can be written in differential form as
\[
dX_{t}= g(t)dW_t \, .
\]
This latter form is called the stochastic differential representation of $X_t$. In the following section we 
will generalise this construct.
\end{itemize}

\subsection{Stochastic Differential Equations (S.D.E.s)}
\subsubsection*{Definition 5.2}
A stochastic process $X_{t}$ is said to satisfy an It\^o stochastic differential equation (S.D.E.), if we may write
\[
dX_{t}=\alpha(t,X_{t})\, dt +\beta(t,X_{t})\, dW_{t}\, .
\]

\subsubsection*{Comments}
\begin{itemize}
\item The corresponding integral representation is
\[
X_{t}=X_{0}+\int_{0}^{t}\alpha(\tau,X_{\tau})\, d\tau + \int_{0}^{t}\beta(\tau,X_{\tau})\, dW_{\tau}\, .
\]
\item Note that the first integral above is the usual Riemann integral and the second is an It\^o stochastic integral. 
\item The $\alpha$ term is known as the \emph{drift} coefficient and $\beta^2$ as the \emph{diffusion} coefficient.
\item If $\alpha(t,X_{t})$ and $\beta(t,X_{t})$ do not depend explicitly on $t$ then the S.D.E. is said to be time-homogeneous.
\end{itemize}

\subsubsection*{Simple S.D.E.s}
\begin{itemize}
\item Suppose that $X_0=0$, $\alpha(t,X_{t})=0$ and $\beta(t,X_{t})=1$. Then 
\[
dX_{t}=dW_{t}.
\]
Integrating both sides of the above between $0$ and $t$ gives
\begin{eqnarray*}
\int_{0}^{t}dX_{\tau} &=& \int_{0}^{t}dW_{\tau}\\
\Rightarrow  X_{t}-X_{0}&=& W_t-W_0
\end{eqnarray*}
for which we see that using $X_0=0$ and $W_0=0$ gives $X_t=W_t$, a standard Brownian motion process. 
\item Take $X_0=x_0$, $\alpha(t,X_{t})=a$ and $\beta(t,X_{t})=b$ where $x_0$, $a$ and $b$ are real numbers. 
The resulting S.D.E. is 
\[
dX_{t}=a\,dt+b\, dW_{t}\, .
\]
Integrating both sides of the above between $0$ and $t$ gives
\begin{eqnarray*}
\int_{0}^{t}dX_{\tau} &=& \int_{0}^t a d\tau+b\int_{0}^{t}dW_{\tau}\\
\Rightarrow  X_{t}-X_{0}&=& at+b(W_t-W_0)
\end{eqnarray*}
for which we obtain
\[
X_{t}=x_0+a\,t+b\, W_{t}
\]
and we recognise $X_t$ as a generalised Brownian motion process.
\item We would like to be able to write down a S.D.E. for the stock price process $S_{t}$, assuming it follows a geometric B.M. We know that $S_{t}=\exp(X_{t})$ where $X_{t}$ is a generalised Brownian motion with initial condition $\log S_0$, drift $\mu-\sigma^{2}/2$ and diffusion coefficient $\sigma^2$ and we know $dX_{t}$. Hence, we just need a way of differentiating $\exp(X_{t})$ in a sensible way. Essentially, we need a \emph{chain rule} for S.D.E.s. This chain rule is the \emph{It\^o Formula}.
\end{itemize}

\subsubsection{It\^o Formula}
Consider an S.D.E. of the form 
\[
dX_{t}=\alpha(t,X_{t})\, dt +\beta(t,X_{t})\, dW_{t}\, .
\]
Let $G(t,x)$ be a real valued function with continuous partial derivatives
\[
G_{t}=G_{t}(t,x)=\frac{\partial G}{\partial t}\, , \qquad G_{x}=G_{x}(t,x)=\frac{\partial G}{\partial x}\, , \qquad G_{xx}=G_{xx}(t,x)=\frac{\partial^{2}G}{\partial x^{2}}\, ,
\]
then
\begin{eqnarray*}
dG(t,X_{t})&=&\left(G_{t}(t,X_{t})+\alpha(t,X_{t})G_{x}(t,X_{t})+\frac{1}{2}\beta^{2}(t,X_{t})G_{xx}(t,X_{t})\right)\, dt +\nonumber \\ 
& &\,\, \,\, +\, \beta(t,X_{t})G_{x}(t,X_{t})\, dW_{t}\, .
\end{eqnarray*}
\subsubsection*{Sketch proof (not examinable)}
To give an indication of why this is the case, write
\[
dG(t,x)\approx G(t+\Delta t,x+\Delta x)-G(t,x)\, .
\] 
Now take a Taylor series expansion of $G(t+\Delta t,x+\Delta x)$ about $(t,x)$ to give
\begin{eqnarray*}
dG(t,x)&\approx& \Delta t\,G_{t}+\Delta x\,G_{x}+\frac{1}{2}(\Delta t)^{2}\,G_{tt}+\frac{1}{2}\Delta t\Delta x\,G_{tx}+\frac{1}{2}(\Delta x)^{2}\,G_{xx}+\ldots 
\end{eqnarray*}
where $G_{tt}=\partial^{2}G/\partial t^{2}$ and $G_{tx}=\partial^{2}G/\partial t \partial x$. Now replace $\Delta x$ with $\alpha\Delta t+\beta\Delta W$ and $(\Delta x)^{2}$ by $\alpha^{2}(\Delta t)^{2}+2\alpha\beta\Delta t\Delta W+\beta^{2}(\Delta W)^{2}$ to obtain
\begin{eqnarray*}
dG(t,x)&\approx&\Delta t\,G_{t}+(\alpha\Delta t+\beta\Delta W)\,G_{x}+\frac{1}{2}(\Delta t)^{2}\,G_{tt}+\frac{1}{2}(\alpha(\Delta t)^{2}+\beta\Delta t\Delta W)\,G_{tx}+\\
&& \, +\frac{1}{2}(\alpha^{2}(\Delta t)^{2}+2\alpha\beta\Delta t\Delta W+\beta^{2}(\Delta W)^{2})\,G_{xx}+\ldots 
\end{eqnarray*}
Now approximate $(\Delta W)^{2}$ by $\Delta t$ (and note in fact that $\expt[(\Delta W^{2})]=\Delta t$) to write the preceding expression as 
\begin{eqnarray*}
dG(t,x)&\approx&\left(G_{t}+\alpha G_{x}+\frac{\beta^{2}}{2}G_{xx}\right)\Delta t +\beta G_{x}\Delta W+ o(\Delta t)\, .
\end{eqnarray*}
It\^o's formula then follows by letting $\Delta t\to 0$. The above sketch proof is not examinable.

\subsubsection*{Comments}
\begin{itemize}
\item We can write down It\^o's formula for a function of Brownian motion $W_{t}$ by noting that $W_{t}$ satisfies an S.D.E. with $X_{t}=W_{t}$, $\alpha(t,X_{t})=0$ and $\beta(t,X_{t})=1$. Hence for a function $G(t,W_{t})$ we obtain
\begin{equation}\label{ItoBM}
dG(t,W_{t})=\left(G_{t}(t,W_{t})+\frac{1}{2}G_{xx}(t,W_{t})\right)\,dt+G_{x}(t,W_{t})\,dW_{t}\, .
\end{equation}
This is known as the \emph{special case of It\^o formula} and is appropriate for deriving the S.D.E. satisfied by a function of $t$ and $W_t$.
\end{itemize}

\subsubsection*{Example 5.1 (Two applications of It\^o formula)}

\vspace{0.1cm}

\textbf{Deriving the S.D.E. satisfied by G.B.M.}

\vspace{0.2cm}

Consider again the task of formulating geometric Brownian motion (to model stock price $S_{t}$) as an S.D.E. We can write $S_{t}=\exp(X_{t})$ where 
\[
X_{t}=X_{0}+\left(\mu-\frac{\sigma^{2}}{2}\right)\, t +\sigma\,W_{t}
\]
and $X_0=\log(S_0)$. The S.D.E. satisfied by the $X_t$ process is then given by
\[
dX_{t}=\left(\mu-\frac{\sigma^{2}}{2}\right)\, dt +\sigma\,dW_{t}\, .
\]
Hence, we derive an S.D.E. for $S_{t}$ by applying It\^o's formula with $G(t,x)=e^x$. Identify
\[
G_{t}=0\, , \qquad G_{x}=e^{x}\, ,\qquad G_{xx}=e^{x}\, .
\]
Hence
\[
G_{t}(t,X_{t})=0\, , \qquad G_{x}(t,X_{t})=e^{X_{t}}\, ,\qquad G_{xx}(t,X_{t})=e^{X_{t}}
\]
and applying It\^o's formula gives
\begin{eqnarray*}
d\left(e^{X_{t}}\right)&=& \left(\left[\mu-\frac{\sigma^{2}}{2}\right]e^{X_{t}}+\frac{1}{2}\sigma^{2}e^{X_{t}}\right)\,dt +\sigma e^{X_{t}}\,dW_{t}\\
\Rightarrow d\left(e^{X_{t}}\right)&=& \mu e^{X_{t}}\,dt+\sigma e^{X_{t}}\,dW_{t}\, .
\end{eqnarray*}
Now, writing $S_{t}=e^{X_{t}}$, we obtain the usual S.D.E. for geometric Brownian motion as
\begin{eqnarray*}
dS_{t}&=& \mu S_{t}\,dt+\sigma S_{t}\,dW_{t}\, .
\end{eqnarray*}

\vspace{0.1cm}

\textbf{Solving the S.D.E. satisfied by G.B.M.}

\vspace{0.2cm}

We can solve the above S.D.E. explicitly (in the sense that we can obtain a closed form expression for $S_{t}$). Take the S.D.E. satisfied by G.B.M. and write
\begin{eqnarray*}
\frac{1}{S_t}dS_{t}&=& \mu \,dt+\sigma \,dW_{t}\, ,\\
\Rightarrow \int_{0}^{t}\frac{1}{S_\tau}dS_{\tau}&=& \int_{0}^{t}\mu \,d\tau+\sigma \,\int_{0}^{t}dW_{\tau}\, ,
\end{eqnarray*}
which suggests that the solution involves $\log(S_t)$. We can check this hunch via It\^o formula. 
Apply the It\^o formula to $\log(S_{t})$ by taking $G(t,x)=\log(x)$ and identify
\[
G_{t}(t,S_{t})=0\, ,\qquad G_{x}(t,S_{t})=\frac{1}{S_{t}}\, ,\qquad G_{xx}(t,S_{t})=-\frac{1}{S_{t}^{2}}\, .
\]
Hence,
\begin{eqnarray*}
d\left(\log(S_{t})\right)&=&\left(\mu S_{t}\frac{1}{S_{t}}-\frac{1}{2}S_{t}^{2}\frac{\sigma^{2}}{S_{t}^{2}}\right)\, dt +\sigma S_{t}\frac{1}{S_{t}}\,dW_{t}\\
&=&\left(\mu-\frac{1}{2}\sigma^{2}\right)\,dt+\sigma\,dW_{t}\, . 
\end{eqnarray*}
Our hunch is just about correct, except for the extra term involving $\sigma^2$ appearing on the RHS. Nevertheless, 
we may still proceed by integrating both sides of the preceding equation to give
\begin{eqnarray*}
\int_{0}^{t}d\left(\log(S_{\tau})\right)&=& \int_{0}^{t}\left(\mu-\frac{1}{2}\sigma^{2}\right)\,d\tau +\int_{0}^{t}\sigma\,dW_{\tau}\\
\Rightarrow \log\left(\frac{S_{t}}{S_{0}}\right)&=&  \left(\mu-\frac{1}{2}\sigma^{2}\right)\, t+\sigma\,W_{t}\\
\Rightarrow S_{t}&=&S_{0}\exp\left\{ \left(\mu-\frac{1}{2}\sigma^{2}\right)\, t+\sigma\,W_{t}\right \}
\end{eqnarray*}
which is exactly what we would expect.

\subsubsection*{Example 5.2}
Using It\^o's formula, solve
\[
\int_{0}^{t}W_{\tau}\,dW_{\tau}\, .
\]
\subsubsection*{\emph{Solution}}
\emph{We would expect the solution to the above integral to contain a $W_{t}^{2}$ term. Therefore we find $d(W_{t}^{2})$ via It\^o's formula for functions of Brownian motion, given by equation (\ref{ItoBM}).}

\emph{Set $G(t,x)=x^{2}$ and calculate
\[
G_{t}(t,W_{t})=0\,,\qquad G_{x}(t,W_{t})=2W_{t}\,,\qquad G_{xx}(t,W_{t})=2\, .
\]
Hence, we obtain
\begin{eqnarray*}
d(W_{t}^{2})&=&dt+2W_{t}\,dW_{t}\\
\Rightarrow W_{t}\,dW_{t}&=& \frac{1}{2}d(W^{2}_{t})-\frac{1}{2}dt\\
\Rightarrow \int_{0}^{t}W_{\tau}\,dW_{\tau}&=& \frac{1}{2}\int_{0}^{t}d(W^{2}_{\tau})-\frac{1}{2}\int_{0}^{t}d\tau \\
&=&\frac{1}{2}W_{t}^{2}-\frac{1}{2}t\, .
\end{eqnarray*}}

\subsubsection*{Example 5.3}
Show that $Y_{t}=t^{3}+t^{2}W_{t}^{4}$ is an It\^o process by writing it in the form
\[
dY_{t}=\alpha_{t}\,dt+\beta_{t}\,dW_{t}
\]
for suitable choices of the processes $\alpha_{t}$ and $\beta_{t}$.
\subsubsection*{\emph{Solution}}
\emph{We use the special case of It\^o's formula (\ref{ItoBM}) with $G(t,x)=t^{3}+t^{2}x^{4}$. Hence
\[
G_{t}(t,W_{t})=3t^{2}+2tW_{t}^{4}\,,\,\, G_{x}(t,W_{t})=4t^{2}W_{t}^{3}\,,\,\,G_{xx}(t,W_{t})=12t^{2}W_{t}^{2}
\]
and we get
\[
d(t^{3}+t^{2}W_{t}^{4})=dY_{t}=\left(3t^{2}+2tW_{t}^{4}+6t^{2}W_{t}^{2}\right)\,dt+4t^{2}W_{t}^{3}\,dW_{t} 
\]
which is in the form required.}

\subsection{Solving S.D.E.s numerically}
Consider the task of generating a skeleton path on $[0,T]$ of a process satisfying an S.D.E. of the form given by 
\[
dX_{t}=\alpha(t,X_{t})\, dt +\beta(t,X_{t})\, dW_{t}\, .
\]
This would typically involve solving the S.D.E. analytically to obtain the transition density, and then simulating using this density. For example, a Gaussian transition density is obtained for generalised B.M. and a log-normal transition density is obtained for geometric B.M.  However, for all but the most trivial of cases, solving the S.D.E. is impossible. We therefore seek to approximate a skeleton sample path of $X_{t}$ by a single realisation of a numerical solution; although many types of numerical solution exist, only the simplest scheme, the Euler-Maruyama method, is considered here. For a small time increment $\Delta t$, the first order Euler-Maruyama approximation of the above S.D.E. is
\[
\Delta X_{t}  = \alpha(t,X_{t})\, \Delta t \, + \, \beta(t,X_{t} )\, \Delta W_{t} \, ,
\]
where $\Delta W_t = (W_{t+\Delta t} - W_t) \sim \textrm{N}(0,\Delta t)$. A sample path is then constructed by dividing the time interval $[0,T]$ into $n$ equidistant points, $0=t_{0}<t_{1}<\ldots <t_{n}=T$ (so that $\Delta t=t_{i+1}-t_{i}$ and $\Delta W_{t_{i}}=W_{t_{i}+\Delta t}-W_{t_{i}}$). Denoting the numerical solution at times $t_{0},\ldots ,t_{n}$ by $X_{0},\ldots ,X_{n}$ and using the equation directly above yields the recursion,
\[
X_{i}=X_{i-1} \, + \, \alpha(t_{i-1},X_{i-1})\Delta t \, + \, \beta(t_{i-1},X_{i-1})\, \Delta W_{t_{i-1}} \, , \, \, X_{t_{0}}=X_{0}\, , \, i=1,\ldots ,n \, .
\]
\newpage
Algorithmically:-
\begin{enumerate}
\item Initialise with $X_{0}$. Put $i:=1$
\item Simulate $X_{i}\sim \textrm{N}(X_{i-1}+\alpha(t_{i-1},X_{i-1})\Delta t\,,\,\beta^{2}(t_{i-1},X_{i-1})\Delta t)$
\item If $t_{i}=T$, stop otherwise put $i:=i+1$ and go to step 2.
\end{enumerate}
As an example, consider the S.D.E. formulation of a geometric Brownian motion process (denoted by $X_t$ here),
\[
dX_{t}=\mu X_{t}\,dt+\sigma X_{t}\, dW_{t}\, .
\]
To make the distinction between the exact solution and an approximate solution (obtained via Euler-Maruyama), consider a time interval $[t,t+\Delta t]$ and suppose that we have $X_t=x_t$ at time $t$. The analytic solution over this time interal is
\[
X_{t+\Delta t}=x_t \exp\left\{ (\mu-\sigma^2/2)\Delta t + \sigma (W_{t+\Delta t}-W_t)\right\}
\]
for which we see
\[
X_{t+\Delta t}|X_t = x_t \sim LN\left(\log x_t + (\mu-\sigma^2/2)\Delta t \,,\, \sigma^2 \Delta t\right).
\]
The Euler-Maruyama approximation gives
\[
X_{t+\Delta t}=x_t + \mu x_t \Delta t + \sigma x_t (W_{t+\Delta t} - W_t)
\]
for which we see
\[
X_{t+\Delta t}|X_t = x_t \sim N\left(x_t + \mu x_t \Delta t\,,\, \sigma^2 x_t^2 \Delta t\right).
\]
Although beyond the scope of the course, taking the power series characterisation of the exponential function and truncating terms can be used to obtain the Euler-Maruyama approximation from the analytic solution.

The following suite of \R functions generate a numerical solution to the S.D.E. formulation of geometric Brownian motion: 

\begin{verbatim}
itosim=function(T=20,dt=0.01,x0=40,afun=alpha,bfun=beta) 
{
 n=T/dt
 xvec=vector("numeric",n+1)
 xvec[1]=x0
 for(i in 2:(n+1))
 {
  t=i*dt
  xvec[i]=xvec[i-1]+afun(xvec[i-1],t)*dt+bfun(xvec[i-1],t)*rnorm(1,0,sqrt(dt))
 }
 xvec
}

alpha=function(x,t)
{ 
 mu=0.1
 mu*x
}

beta=function(x,t)
{ 
 sigma=0.2
 sigma*x
}
\end{verbatim}
Note that this flexible setup allows us to simulate from any (univariate) S.D.E. simply by re-writing the \texttt{alpha} and \texttt{beta} functions. The main \texttt{itosim} function remains unchanged. We plot the numerical solution with
\begin{verbatim}
plot(ts(itosim(),start=0,deltat=0.01))
\end{verbatim}

\begin{figure}[h]
\begin{center}
\includegraphics[width=4cm,angle=270]{sde1-eps-converted-to.pdf}
\end{center}
\caption{2 simulated realisations of a geometric B.M. with $x_{0}=40$, $\mu=0.1$ and $\sigma=0.2$ found by numerically solving the S.D.E.}\label{fig:fig_sde1}
\end{figure}

\subsubsection*{Comment}
\begin{itemize}
\item Taking $\Delta t$ smaller and smaller gives increased accuracy at greater computational expense.
\end{itemize}


\subsubsection{Assessing the accuracy of the Euler-Maruyama scheme (not examinable)}

Suppose that we're interested in the accuracy of the numerical scheme at some (maturity) time $T$. Let $X_T^{\Delta t}$ denote the 
Euler-Maruyama approximation of $X_T$ at time $T$, based on a discretisation (step-size) of $\Delta t$.


The numerical scheme is \emph{strongly convergent} if 
\[
\lim_{\Delta t\to 0} \textrm{E}\left(|X_T-X_T^{\Delta t}|\right)=0.
\]
The numerical scheme is \emph{weakly convergent} if 
\[
\lim_{\Delta t\to 0} \big|\textrm{E}\{g(X_T)\}-\textrm{E}\{g(X_T^{\Delta t})\}\big|=0
\]
for polynomials $g(\cdot)$.


The Euler-Maruyama scheme is strongly and weakly convergent subject to suitable conditions on the drift and diffusion 
coefficient. We may also wonder how good the approximation is for a particular $\Delta t$. The numerical scheme is 
strongly convergent with order $\delta$ if 
\[
\textrm{E}\left(|X_T-X_T^{\Delta t}|\right)\leq C_T(\Delta t)^{\delta}
\]
where the constant $C_T$ depends on $T$ and the considered SDE. The numerical scheme is weakly convergent with order $\delta$ if 
\[
\big|\textrm{E}\{g(X_T)\}-\textrm{E}\{g(X_T^{\Delta t})\}\big|\leq C_T^g(\Delta t)^{\delta}
\]
where $C_T^g$ depends on $T$, $g$ and the considered SDE. If a numerical scheme is convergent with order $\delta$ and we make the 
step size $k$ times smaller, the approximation error will decrease by a factor $k^\delta$.


The Euler-Maruyama scheme is weakly convergent with order 1 and strongly convergent with order 0.5. Therefore, the order equal
1 means that if we want to decrease the error 100 times, we have to make the step 100 times smaller. The order equal to
0.5 means that if we want to decrease the error 100 times, we have to make the step $100^2 = 10000$ times smaller. 
Consequently, the Euler-Maruyama scheme is useful for pricing options whose payoff is not path-dependent. Where the 
payoff is path-dependent, a higher order scheme (e.g. the Milstein scheme) should be used.

\newpage

\subsection{Models of Interest Rate}
Interest rate derivatives are instruments whose payoffs are dependent on the level of interest rates. Such derivatives can be more difficult to value than equity or exchange rate derivatives because
\begin{itemize}
\item the behaviour of an individual interest rate is more complicated than that of a stock price,
\item interest rates are used for discounting as well as for defining the payoff from the derivative.
\end{itemize}
In this Section, we will consider several models which provide a description of how short-term interest rate changes through time. The short rate $r_{t}$ is the rate that applies to an infinitesimally short period of time at time $t$.

We will describe the short rate as an It\^o process of the form
\[
dr_{t}=\alpha(r_{t})\,dt+\beta(r_{t})\,dW_{t}
\]
where the drift and diffusion coefficients are independent of time. We will consider three models
\begin{enumerate}
\item $\alpha(r_{t})=\mu r_{t}$, $\beta(r_{t})=\sigma r_{t}$ (Rendleman and Bartter model),
\item $\alpha(r_{t})=a(b-r_{t})$, $\beta(r_{t})=\sigma $ (Vasicek model),
\item $\alpha(r_{t})=a(b-r_{t})$, $\beta(r_{t})=\sigma\sqrt{r_{t}}$ (Cox, Ingersoll and Ross model). 
\end{enumerate}

\subsubsection{Rendleman and Bartter Model}
In the Rendleman and Bartter model, the short rate is governed by the S.D.E.
\[
dr_{t}=\mu r_{t}\,dt+\sigma r_{t}\,dW_{t}\, .
\]
This means that $r_{t}$ follows a geometric Brownian motion and is the same type of process that we assumed for stock price in Section 1. This assumption is a good starting point; for example, the model ensures that interest rates are positive. However, some key properties are also lacking, such as mean reversion.

\subsubsection*{Definition 5.3 (mean reversion)}
When interest rates are high, the economy tends to slow down and there is a low demand for funds from borrowers. As a result, rates go down. When rates are low, there is a high demand from borrowers and rates tend to rise. This is known as \emph{mean reversion}. The effect is negative drift when $r_{t}$ is high and positive drift when $r_{t}$ is low. Overall, interest rates appear to be pulled back to some long-run average level over time. 

\subsubsection*{Comment}

The Rendleman and Bartter model does not incorporate mean reversion since plainly,
\[
\expt\left(r_{t}\right)=r_{0}\exp\left(\mu t\right)\to \infty \qquad \text{as $t\to \infty$}
\]
for $\mu$ positive. If $\mu$ is negative,
\[
\expt\left(r_{t}\right)=r_{0}\exp\left(\mu t\right)\to 0 \qquad \text{as $t\to \infty$}.
\]
Neither scenario is compatible with the definition of mean reversion.

\subsubsection{Vasicek Model}
Here, the short rate is modelled by
\[
dr_{t}=a(b-r_{t})\,dt+\sigma\,dW_{t}\, ,
\]
where $a,b$ and $\sigma$ are positive real numbers. Note that the model incorporates mean reversion and $r_{t}$ is pulled to a level $b$ at rate $a$. The process is also known as a Gaussian Ornstein-Uhlenbeck process.
\begin{figure}[h]
\begin{center}
\includegraphics[width=4cm,angle=270]{sde2-eps-converted-to.pdf}
\end{center}
\caption{2 simulated skeleton paths from the Vasicek model with $b=0.1$, $\sigma=0.02$ and  (a) $a=1$, and (b) $a=0.1$}\label{fig:fig_sde2}
\end{figure}

To see that the Vasicek model is mean reverting, we apply It\^o's formula with $G(t,x)=xe^{at}$. We obtain
\[
G_{t}=xae^{at}\, , \qquad G_{x}=e^{at}\, ,\qquad G_{xx}=0\, .
\]
Hence,
\[
G_{t}(t,r_{t})=ar_{t}e^{at}\, , \qquad G_{x}(t,r_{t})=e^{at}\, ,\qquad G_{xx}(t,r_{t})=0
\]
Thus,
\begin{eqnarray*}
d\left(r_{t}e^{at}\right) &=& \left(ar_{t}e^{at}+a(b-r_{t})e^{at}\right)dt+\sigma e^{at}dW_{t}\\
&=& abe^{at}dt+\sigma e^{at}dW_{t}\\
\Rightarrow r_{t}e^{at}-r_{0} &=& \int_{0}^{t}abe^{a\tau}d\tau +\sigma\int_{0}^{t}e^{a\tau}dW_{\tau}\\
\Rightarrow r_{t}e^{at} &=& r_{0}+\left[be^{a\tau}\right]_{0}^{t}+\sigma\int_{0}^{t}e^{a\tau}dW_{\tau}\\
\Rightarrow r_{t} &=& r_{0}e^{-at}+b(1-e^{-at})+\sigma e^{-at}\int_{0}^{t}e^{a\tau}dW_{\tau}
\end{eqnarray*}
To proceed, we need (the first part) of the following \emph{result} for square-integrable functions $g(t)$.

\subsubsection*{Result 5.1}

\vspace{0.2cm}

\begin{itemize}
\item[(i)] $\expt\left(\int_{0}^{t}g(\tau)dW_{\tau}\right)=0$
\item[(ii)]$\expt\left(\left[\int_{0}^{t}g(\tau)dW_{\tau}\right]^{2}\right)=\int_{0}^{t}g^{2}(\tau)d\tau$. 
\end{itemize}
Note that part (ii) is known as the \emph{It\^o isometry}. Justification of this result is left as an exercise. Combining the 
two results and noting that It\^o integrals involve combinations of normal random variables leads to part (iii) that
\begin{itemize}
\item[(iii)] 
\[
\int_{0}^{t}g(\tau)dW_{\tau} \sim N\left(0\,,\, \int_{0}^{t}g^{2}(\tau)d\tau\right).
\]
\end{itemize}

Now, returning to the Vasicek model and using the fact that $\expt\left(\int_{0}^{t}g(\tau)dW_{\tau}\right)=0$, we have 
\[
\expt(r_{t})=r_{0}e^{-at}+b(1-e^{-at})
\]
which is just $b$ in the limit as $t\to\infty$. The Vasicek model is mean reverting.

As an interesting exercise, we may also find the variance of $r_t$. We have using (ii) of Result 5.1 that
\begin{eqnarray*}
Var(r_t) &=& \sigma^{2}e^{-2at}\int_{0}^{t}e^{2a\tau}d\tau\\
&=& \frac{\sigma^2}{2a}\left(1-e^{-2at}\right).
\end{eqnarray*}
Part (iii) of Result 5.2 then gives
\[
r_t|r_0 \sim N\left(r_{0}e^{-at}+b(1-e^{-at})   \,,\, \frac{\sigma^2}{2a}\left(1-e^{-2at}\right)  \right).
\]
This distribution can be sampled recursively to give a skeleton path of the process (that is, a realisation of the process at discrete times). Two such realisations can be seen in Figure~\ref{fig:fig_sde2}. The effect of mean reversion is apparent, as is the effect of $a$, since the larger the value of $a$ results in stronger reverting behaviour. It should also be clear that the Vasicek model allows for negative interest rates, which is an undesirable property.    

\subsubsection{Cox, Ingersoll and Ross Model}
In the model proposed by Vasicek, the short term rate can become negative. Whilst this is not considered impossible by some (in practical terms), others consider it to be unsatisfactory. Cox, Ingersoll and Ross propose an alternative model where rates are always nonnegative,
\[
dr_{t}=a(b-r_{t})\,dt+\sigma\sqrt{r_{t}}\,dW_{t}\, ,
\]
where $a,b$ and $\sigma$ are positive real numbers. The process is mean reverting and the standard deviation of the change in interest rate in a short time period is proportional to $\sqrt{r_{t}}$. Consequently, as short rate increases, so does the standard deviation. Although the S.D.E. can be solved analytically, this task is beyond the scope of the course. Therefore, skeleton paths were generated using the Euler-Maruyama scheme and plotted in Figure~\ref{fig:fig_sde3}. 
\begin{figure}[h]
\begin{center}
\includegraphics[width=4cm,angle=270]{sde3-eps-converted-to.pdf}
\end{center}
\caption{2 simulated skeleton paths from the CIR model with $b=0.1$, $\sigma=0.02$ and  (a) $a=1$, and (b) $a=0.1$}\label{fig:fig_sde3}
\end{figure}
Finally, note that the process is non-negative. To see this intuitively (rather than formally), note that the infinitesimal variance $\sigma^2 r_t$ tends to zero as the process approaches zero (from above). We expect this to be the case for non-negative processes, otherwise fluctuations around zero due to the non-negligible stochastic noise term near zero will send the process negative.


\newpage

\subsection{Stochastic Volatility Models (not examinable)}
We have seen in Section 3 that it may be unreasonable to assume that the volatility of a stock is constant over time. Here, we use the S.D.E. formalism to construct a class of models whose volatility is also a stochastic process. Such models are called \emph{stochastic volatility (S.V.) models}. We will consider models of the form
\begin{eqnarray*}
dS_{t}&=& (\mu-q)S_{t}\,dt+S_{t}\sqrt{V_{t}}\,dW_{1,t}\\
dV_{t}&=& a(b-V_{t})\,dt+\sigma\sqrt{V_{t}}\,dW_{2,t}
\end{eqnarray*}
where $S_{t}$ is an asset's price, $V_{t}$ is a volatility process, $\mu,q,a,b,\sigma$ are constants and $W_{1,t}$, $W_{2,t}$ are uncorrelated Brownian motions. Clearly, $V_{t}$ is of C.I.R. type and $S_{t}$ (conditional on $V_{t}$) is of G.B.M. type. Note that we can write this model in the form
\[
dX_{t}=\alpha(X_{t})\,dt+\beta(X_{t})\,dW_{t}
\]
by setting
\[
X_{t}=\left(\begin{array}{c}
S_{t}\\
V_{t}
\end{array}\right), \, \, \alpha(X_{t})=\left(\begin{array}{c}
(\mu-q)S_{t}\\
 a(b-V_{t})
\end{array}\right), \, \, \beta(X_{t})=\left(\begin{array}{cc}
S_{t}\sqrt{V_{t}} & 0\\
0 & \sigma\sqrt{V_{t}}
\end{array}\right), \, \, dW_{t}=\left(\begin{array}{c}
dW_{1,t}\\
dW_{2,t}
\end{array}\right)\, .
\]

By adopting the notation of Section~5.3 for a partition of $[0,T]$, the Euler approximation for this model is given by
\begin{eqnarray*}
S_{i}&=& S_{i-1}+ (\mu-q)S_{i-1}\,\Delta t+S_{i-1}\sqrt{V_{i-1}}\,\Delta W_{1,t}\\
V_{i}&=& V_{i-1}+a(b-V_{i-1})\,\Delta t+\sigma\sqrt{V_{i-1}}\,\Delta W_{2,t}
\end{eqnarray*}
for $i=1,\ldots ,n$ with $\Delta W_{1,t}\sim\textrm{N}(0,\Delta t)$ and $\Delta W_{2,t}\sim\textrm{N}(0,\Delta t)$ uncorrelated. Since the volatility process is independent of the price process, we simulate a skeleton path from the S.V. model with the following algorithm:
\begin{enumerate}
\item Initialise with $S_{0},V_{0}$.
\item For $i=1,\ldots ,n$, draw $V_{i}\sim\textrm{N}(V_{i-1}+a(b-V_{i-1})\,\Delta t\, , \, \sigma^{2}V_{i-1}\Delta t)$.
\item For $i=1,\ldots ,n$, draw $S_{i}\sim\textrm{N}(S_{i-1}+ (\mu-q)S_{i-1}\,\Delta t\, , \, S_{i-1}^{2}V_{i-1}\Delta t)$.
\end{enumerate}
\begin{figure}[h]
\begin{center}
\includegraphics[width=4cm,angle=270]{sde4-eps-converted-to.pdf}
\end{center}
\caption{1 simulated skeleton path from the S.V. model with $\mu=0.1$, $q=0.0$, $a=0.2$, $b=0.1$ and $\sigma=0.02$. Plot (a) shows the $S_{t}$ process and plot (b) shows the $V_{t}$ process.}\label{fig:fig_sde4}
\end{figure}
Figure~\ref{fig:fig_sde4} shows a simulated skeleton path from the S.V. model. It is easy to see the effect of the volatility realisation on the stock price realisation.

Naturally, the S.V. model can be used to price options. Since typically no closed form pricing formulae exist when it is assumed that stock price follows the S.V. model, the Monte Carlo methods of Section 3 must be used.

\subsubsection*{Example 5.4}
Suppose that stock price follows a S.V. model with parameters $q=0.0$, $a=0.2$, $b=0.1$ and $\sigma=0.02$. Suppose further that the initial price of stock is $\pounds 40$, the risk free interest rate is $r=0.05$, and we wish to find the fair price of the Asian call option with strike price $K=50$ and maturity $T=1$ year.

We estimate the fair price via Monte Carlo by performing the following steps:
\begin{enumerate}
\item Simulate daily prices $S^{d}_{0},S^{d}_{1},\ldots ,S^{d}_{252}$ in the risk-neutral world. For the S.V. model, this means setting $\mu=r$.
\item Calculate the payoff at time $T=1$:
\[
\left(\sum_{i=1}^{252}\frac{S^{d}_{i}}{252}-K\right)^{+} \, . 
\]
\item Repeat steps 1 and 2 to get $N$ sample values of the payoff.
\item Calculate the mean of these sample payoffs to get an estimate of the expected payoff.
\item Discount the expected payoff at the risk-free interest rate ($r=0.05$ here) to get an estimate of the value of the option.
\end{enumerate}

\clearpage

The following \R function performs the algorithm:

\begin{verbatim}
sv=function(T=1,dt=1/252,s0=c(50,0.1),r=0.05,a=0.2,b=0.1,sigma=0.02,k=50,N=1000) 
{
 n=T/dt
 sdt=sqrt(dt)
 s=matrix(0,ncol=2,nrow=n+1)
 payoff=vector("numeric",len=N)
 for(j in 1:N)
 { 
  s[1,]=s0
  for(i in 2:(n+1))
  {
   s[i,2]=s[i-1,2]+a*(b-s[i-1,2])*dt+sigma*sqrt(s[i-1,2])*rnorm(1,0,sdt)
  }
  for(i in 2:(n+1))
  {
   s[i,1]=s[i-1,1]+r*s[i-1,1]*dt+sqrt(s[i-1,2])*s[i-1,1]*rnorm(1,0,sdt)
  }
  payoff[j]=max(mean(s[,1])-k,0)
 }
 exp(-r*T)*mean(payoff)
}
\end{verbatim}
A single execution of this function gave an estimate of $\pounds 3.80$ for the fair price. 

\subsection{Stochastic Volatility Models with Jumps (not examinable)}
Whilst the stochastic volatility remedies the constant volatility assumption of G.B.M., it 
does not allow for any discontinuities in the price process. Such jumps can correspond to a crash 
in the market (e.g. Black Monday October 1987, 9/11 etc). We therefore extend the S.V. model of the previous 
Section to incorporate jumps in both prices and volatilities. The discretised model is 
\begin{eqnarray*}
S_{i}&=& S_{i-1}+ (\mu-q)S_{i-1}\,\Delta t+S_{i-1}\sqrt{V_{i-1}}\,\Delta W_{1,t}+Z^{s}_{i}J_{i}\\
V_{i}&=& V_{i-1}+a(b-V_{i-1})\,\Delta t+\sigma\sqrt{V_{i-1}}\,\Delta W_{2,t}+Z^{v}_{i}J_{i}
\end{eqnarray*}
defined for $i=1,\ldots ,n$. Now $J_{i}\sim Bernouilli(\Delta t)$ represent the jump times 
and $Z^{s}_{i}$ and $Z^{v}_{i}$ are the jump sizes. Typically $Z^{v}_{i}$ follows an \emph{Exponential} 
distribution and $Z^{s}_{i}$ follows a \emph{Normal} distribution.  

Consider the S\&P500 index data of Section~3. The returns are plotted overleaf as well as the estimated volatility process (using a Bayesian estimation method beyond the scope of this course). Clearly, fitting the above model to this dataset gives a jump in volatility on Black Monday.

\begin{figure}[h]
\begin{center}
\includegraphics[width=63mm,angle=270]{retsp-eps-converted-to.pdf}
\includegraphics[width=63mm,angle=270]{volsp-eps-converted-to.pdf}
\end{center}
\caption{S\&P500 index data (left) and filtered volatilities (right) --- 2.5 and 97.5 percentiles (red lines) and 50 percentiles (black line).}
\end{figure}

\end{document}

\newpage

\subsection*{Section quiz}

\begin{itemize}
\item[\textbf{Q1}] Identify the processes that satisfy each of the following S.D.E.s:
\begin{itemize}
\item[\textbf{(a)}] $dX_t = 2X_t dW_t$
\item[\textbf{(b)}] $dY_t = -Y_t dt + dW_t$
\item[\textbf{(c)}] $dZ_t = 3dt -2dW_t$
\item[\textbf{(d)}] $dA_t = e^2 dW_t$
\end{itemize}

\item[\textbf{Q2}] For the processes that satisfy each of the following S.D.E.s, state which are non-negative:
\begin{itemize}
\item[\textbf{(a)}] $dX_t = 2X_t dW_t$
\item[\textbf{(b)}] $dY_t = -Y_t dt + dW_t$
\item[\textbf{(c)}] $dZ_t = 3dt -2dW_t$
\item[\textbf{(d)}] $dA_t = e^2 dW_t$
\end{itemize}

\item[\textbf{Q3}] For the processes that satisfy each of the following S.D.E.s, state which are mean reverting:
\begin{itemize}
\item[\textbf{(a)}] $dX_t = 2X_t dW_t$
\item[\textbf{(b)}] $dY_t = -Y_t dt + dW_t$
\item[\textbf{(c)}] $dZ_t = 3dt -2dW_t$
\item[\textbf{(d)}] $dA_t = e^2 dW_t$
\end{itemize}


\end{itemize}