In this section we look for a more theoretical approach to pricing the option. The underlying model in this section is assumed to be geometric Brownian motion. That is the price of the asset \(S_t\) is given by

\(\frac{\text{d} S_t}{S_t} = (r-d)\text{d}t +\sigma \text{d}W_t\)

where \(r,d\) and \(\sigma\) are constants. We can, for this process, derive an approximating continuous time Markov chain \(X_t\) with Q-matrix \(Q\). See Jump Diffusion Model on how to calculate the Q-matrix for the jump-diffusion model. The geometric Brownian motion model is a special case of the jump-diffusion model when the jump part of the generating matrix has been set to zero. This Markov chain is used to calculate the option price. We start by pricing the asset and then look at pricing the European up-and-in call option. Pricing the option is more difficult and involves finding the density of the stopping time defined by

\[ \tau_H=\inf\{s:X_s=H\}. \]

Price of the Asset

The price of the asset is reasonably easy to calculate using

\[ \mathbb{P}_{X_0}\left[ X_t=x_i \right] = P_0 e^{tQ}(x_i). \]

\(P_0\) is the initial probability distribution and \(Q\) is the Q-matrix for the asset price. The only difficulty is in calculating \(e^{tQ}\). We use that any matrix can be put into Jordan normal form. I.e. there exists a change of basis matrix \(P\) and a Jordan normal matrix \(J\) such that

\[ Q = PJP^{-1}. \]

The Q-matrix has top and bottom row set to zero and is tridiagonal. We assume that the there are only two zero eigenvalues and all the remaining eigenvalues are distinct. With this assumption the Jordan normal Form \(J\) is diagonal. So to exponentiate \(tQ\) we can use

\[ e^{tQ} = Pe^{tJ}P^{-1} \]

and diagonal matrices are easy to exponentiate. This gives a computationally easy method to calculate the price of the asset at any time. The option price pdf is harder to calculate.

Price of the Option

Let \(Pa\) be the payoff of the option, then for \(x>0\)

\[ \mathbb{P}_{X_0}\left[Pa=x\right] = \mathbb{P}_{X_0} \left[ (X_s \geq H \text{ for some } s\in [0,1]) \cap (X_1-K=x) \right]. \]

We can rewrite this as an integral of the asset price and stopping time density using the strong Markov property.

\( \begin{align*} \mathbb{P}_{X_0}\left[Pa=x\right] & = \mathbb{E}_{X_0}\left[ \mathbb{I}_{\tau_H\leq 1} \mathbb{I}_{X_1=x+K} \right] \\ & = \mathbb{E}_{X_0}\left[\mathbb{E}_{X_0}\left[\mathbb{I}_{\tau_H\leq 1}\mathbb{I}_{X_1=x+K}|\mathscr{F}_{\tau_H}\right]\right] \\ & = \mathbb{E}_{X_0}\left[\mathbb{I}_{\tau_H\leq 1}\mathbb{E}_{X_0}\left[\mathbb{I}_{X_{1-\tau_H}=x+K}\circ \Theta_{\tau_H}|\mathscr{F}_{\tau_H}\right]\right] \\ & = \mathbb{E}_{X_0}\left[\mathbb{I}_{\tau_H\leq 1} \mathbb{E}_H\left[ \mathbb{I}_{X_{1-\tau_H}=x+K} \right] \right] \\ & = \int_0^1 \mathbb{P}_H\left[ X_{1-r}=x+K\right] f_{X_0}(r) \; \text{d} r \end{align*} \)

where \(f_{X_0}(r)\) is the density associated to \(\mathbb{P}_{X_0}(\tau_H)\), i.e. \(\mathbb{P}_{X_0}[\tau_H\in A]=\int_A f_{X_0}(r) \; \text{d} r\). We look to find the density of the stopping time \(\tau_H\). We define the Laplace transform of a function \(f:[0,\infty)\rightarrow \mathbb{R}\) to be \(\hat{f}:(0,\infty) \rightarrow \mathbb{R}\) where

\[ \hat{f}(\lambda)=\int_0^\infty e^{-\lambda t} f(t) \; \text{d} t. \]

Let \(R(\lambda)\) be the Laplace transform of \(P(t)=\left(\mathbb{P}_i[X_t=j[\right)_{ij}\).

\(\begin{align*} R(\lambda) & = \int_0^\infty e^{-\lambda t}P(t) \; \text{d} t \\ & = \int_0^\infty e^{-(\lambda\mathbb{I}-Q)t} \; \text{d} t \\ & = \left( \lambda \mathbb{I} -Q \right)^{-1} \end{align*}\)

So we can find \(R\) by inverting \(\lambda \mathbb{I}-Q\). Let \(R(\lambda)=\left(r_{ij}\right)_{ij}\), by using the Jordan normal form we can write explicitly that

\[ r_{ij}=\sum_{k=1}^n \frac{p_{kj}p^{(-1)}_{ik}}{\lambda-j_k}. \]

with \(P=(p_{ij})\), \(P^{-1}=(p_{ij}^{(-1)})\) and \(j_k\) are the eigenvalues of \(Q\). This follows since \((\lambda\mathbb{I}-Q)^{-1}=P(\lambda\mathbb{I}-J)^{-1}P^{-1}\). If \(k\neq 1,n\) then \(p_{k1}=p^{(-1)}_{kn}=p_{1k}=p^{(-1)}_{nk}=0\). Hence

\[ r_{ij}=\sum_{k=2}^{n-1} \frac{p_{kj}p^{(-1)}_{ik}}{\lambda-j_k}. \]

Now we note that can write \(r_{ij}\) as

\[ r_{ij} = \frac{1}{\prod_{k=2}^{n-1} (\lambda -j_k)} \left( \sum_{k=2}^{n-1} p_{kj}p_{ik}^{(-1)} \left( \lambda^{n-3} -\sum_{s\neq k} j_s \lambda^{n-4} + \sum_{t,s\neq k} j_sj_t\lambda^{n-5} -\dots \right) \right). \]

Hence if \(i\neq j\) then \(r_{ij}\sim o(-2)\) and if \(i=j\) then \(r_{ij}\sim o(-1)\), since \(\sum_{k=1}^{n} p_{kj}p_{ik}^{(-1)}=\sum_{k=2}^{n-1} p_{kj}p_{ik}^{(-1)}=\delta_{ij}\). With the following theorem we can find the Laplace transform of the stopping time.

Theorem: Let \(f_{ij}\) be the density of the stopping time \(\tau_j\) where the Markov chain is started at \(i\), then \(\hat{f}_{ij}(\lambda) = \frac{r_{ij}(\lambda)}{r_{jj}(\lambda)}\).

Proof: By the strong Markov property:

\[ p_{ij}(t) = \int_0^t f_{ij}(s)p_{jj}(t-s) \; \text{d} s. \]

Defining the convolution \((f*g)(t)=\int_0^t f(s)g(t-s) \; \text{d} s\) we have that

\[ p_{ij}(t) = (f_{ij}*p_{jj})(t). \]

Taking Laplace transforms of each side:

\[ \hat{p}_{ij}(\lambda) = (\widehat{f_{ij}*p_{jj})}(\lambda) = \hat{f}_{ij}(\lambda) \hat{p}_{jj} (\lambda). \]

\( \square \)

Hence we can write

\( \begin{align*} \hat{f}_{ij}(\lambda) & = \frac{ \sum_{k=2}^{n-1} \frac{p_{kj}p_{ik}^{(-1)}}{\lambda-j_k} }{ \sum_{k=2}^{n-1} \frac{p_{kj}p_{jk}^{(-1)}}{\lambda-j_k}} \\ & = \frac{ \sum_{k=2}^{n-1}p_{kj}p_{ik}^{(-1)} \prod_{s\neq k}(\lambda-j_s) }{ \sum_{k=2}^{n-1} p_{kj}p_{jk}^{(-1)} \prod_{s\neq k} (\lambda-j_s) }. \end{align*} \)

Let \(z(\lambda)\) be the denominator in the Laplace transform of the stopping time density,

\[ z(\lambda)=\sum_{k=2}^{n-1} p_{kj}p_{jk}^{(-1)} \prod_{s\neq k} (\lambda-j_s). \]

Notice that

\[ z(j_m) = p_{mj}p_{jm}^{(-1)} \prod_{s\neq m}\left(j_m-j_s\right). \]

If we arrange the eigenvalues in order then we notice that \(z\) switches sign as we move to the next eigenvalue. We have \(n-2\) distinct eigenvalues, hence we must have at least \(n-3\) zero's of \(z\). As \(z\) is an \(n-3\) order polynomial we can factorise

\[ z(\lambda)=c\sum_{k=1}^{n-3} (\lambda-b_k), \]

where \(b_k\) are the zero's of \(z\). Furthermore, since the order of \(f\) is -1,

\[ \hat{f}_{ij}(\lambda) = \sum_{k=1}^{n-3} \frac{a_k}{\lambda-b_k} \]

Now we use that the inverse Laplace transform of \(\frac{A}{\lambda + \gamma}\) is \(Ae^{-\gamma t}\) to show that

\[ f_{ij}(t) = \sum_{k=1}^{n-3} a_ke^{b_k t} \]

We can now give a formula for the probability density function for the option price for \(x>0\)

\[ \mathbb{P}_{X_0}\left[ Pa=x \right] = \sum_{k=1}^n \sum_{i=1}^{n-3} \frac{a_i p_{Hk}p^{(-1)}_{k(x+K)}e^{j_k}}{b_i-j_k}(e^{b_i-j_k}-1). \]

This gives us a method for calculating the option price pdf where we only need numerically estimate the zero's of \(z\).