# Spectral Expansions for Asian Options

## Distributional properties

We assume that  under the risk-neutral probability measure, the underlying asset price follows a geometric Brownian motion process:
$$S_t = S_0 e^{\sigma B_t + (r-q-\sigma^2 /2)t}$$

Let $\mathscr{A}_t$ be the continuously sampled arithmetic average price:
$$\mathscr{A}_t=\frac{1}{t}\int_0^t S_u du$$


From the general form of option pricing formula $e^{-rT}\mathbb{E}[(K-\mathscr{A}_T)^+]$, make the following derivation:
$$\left(K-\frac{1}{T}\int^T_0S_udu\right)^+=\frac{4S_0}{\sigma^2T}\left(k-\frac{\sigma^2}{4S_0}\int^{\tau}_0S_0e^{\sigma B_u+(r-q-\sigma^2/2)u}d\left(\frac{\sigma^2 u}{4}\right)\right)^+\\
=\frac{4S_0}{\sigma^2T}\left(k-\int^{\tau}_0e^{2(B_t+\nu t)}dt\right)$$


Define *Brownian exponential functional* $A_{\tau}^{(\nu)}$:
$$A_{\tau}^{(\nu)}:=\int^{\tau}_0e^{2(B_t+\nu t)}dt$$


From the invariance by time reversal of the distribution of a Levy process, the following process $X_t$ has an identity in law, $X_t \overset{(law)}{=}A_{t}^{(\nu)} $:
$$X_t = \int_0^T \exp{\left(2(B_t-B_u)+2\nu (t-u)\right)}du, X_0=0$$


From Ito's fomula, we express the SDE of $X_t$ as:
$$dX_t = \left[2(\nu+1)X_t+1\right]+2X_tdB_t$$
This $X_t$ is a one-dimensional diffusion process on $[0,\infty)$ started at $X_0 = 0$ and with the infinitesimal diffusion $a(x)=2x$, infinitesimal drift $b(x)=2(\nu+1)x+1$, and infinitesimal generator
$$
\begin{align}
(\mathscr{G}f)(x)&=\frac{1}{2}a^2(x)f''(x)+b(x)f'(x)\\&=2x^2f''(x)+[2(\nu+1)+1]f'(x)
\end{align}
$$
This diffusion has scale and speed densities:
$$\mathscr{s}(x)=\exp{\left(-\int\frac{2b(x)}{a^2(x)}dx\right)}=x^{-\nu-1}e^{1/2x}$$
$$\mathscr{m}(x)=\frac{2}{a^2(x)\mathscr{s}(x)}=\frac{1}{2}x^{\nu-1}e^{-1/2x}$$

In addition, this paper also consider *up-and-out puts* on the diffusion X. For $b>k$, consider functions
$$P^{(\nu)}(k,\tau):=\mathbb{E}\left[(k-X_t)^+\right],P_b^{(\nu)}(k,\tau):=\mathbb{E}\left[1_{\{\mathscr{T}_b>\tau\}}(k-X_\tau)^+\right]$$
where $\mathscr{T}_b:=inf{t\ge0:X_t=b}$ is the first hitting time of $b$ and $1_{\{\mathscr{T}_b>\tau\}}=1$ if $\mathscr{T}_b>\tau$. Their strategy is to first calculate the function $P_b^{(\nu)}(k,\tau)$ and then take the limit $b\uparrow\infty$ to recover the function $P^{(\nu)}(k,\tau)$

## Whittaker functions

The Whittaker functions are solutions to the ODE
$$\frac{d^2w}{dz^2}+\left(-\frac{1}{4}+\frac{\kappa}{z}+\frac{1/4-\mu^2}{z^2}\right)w=0$$

There are two linearly independent solutions to the ODE:
$$M_{\kappa,\mu}=\exp{(-z/2)}z^{\mu+1/2}M\left(\mu-\kappa+\frac{1}{2},1+2\mu;z\right),$$
$$W_{\kappa,\mu}=\exp{(-z/2)}z^{\mu+1/2}U\left(\mu-\kappa+\frac{1}{2},1+2\mu;z\right)$$
where
$$M(a,b,z)=\sum_{k=0}^\infty\frac{(a)_k}{(b)_k}\frac{z^k}{k!}$$
$$U(a,b,z)=\frac{\pi}{\sin{b\pi}}\left(\frac{M(a,b,z)}{\Gamma(1+a-b)\Gamma(b)}-z^{1-b}\frac{M(1+a-b,2-b,z)}{\Gamma(a)\Gamma(2-b)}\right)$$

## Up-and-out puts price derivation

From the last section, $P_b^{(\nu)}(k,\tau)$ is what we care the most, we need the following integral to calculate:
$$P_b^{(\nu)}(k,\tau)=\int_0^k (k-y)p_b(\tau;0,y)\mathscr{m}(y)dy$$
where $p_b(t;x,y)$ is the transition probability density w.r.t. the speed measure $\mathscr{m}(y)dy$ of the infinitesimal generator $(\mathscr{G}f)(x)$ started at $x\in[0,b)$ and killed at $b$.

This prob density admit a spectral representation
$$p_b(t;x,y)=\sum_{n=1}^\infty e^{-\lambda_n \tau}\varphi_n(x)\varphi_n(y)$$
where $\{\lambda_n,\varphi_n\}_{n=1}^\infty$ are eigenvalues and (normalized) eigenfunctions of the Strum-Liouville problem:
$$-(\mathscr{G}u)(x)\equiv \lambda u(x)$$
with the boundary conditions at 0 and $b$:
$$\underset{x\downarrow0}{lim}\frac{u'(x)}{\mathscr{s}(x)}=0, u(b)=0$$

This linear ODE turns out to be the solution of *Whittaker's form of the confluent hypergeometric equation*, with the solution:
$$\psi(x,\lambda)=(2x)^{(1-\nu)/2}e^{1/4x}W_{(1-\nu)/2,(1/2)\sqrt{\nu^2-2\lambda}}\left(\frac{1}{2x}\right)$$

From Theorem 4.1 in Kent (1980), for fixed $x > 0$ the zeros of $\psi(x,\lambda)$ are simple and positive, it can be written as a canonical product
$$\psi(x,\lambda)=\prod_{n=1}^\infty\left(1-\frac{\lambda}{\lambda_{n,x}}\right)$$
by finding root $\psi(x,\lambda)=0$, we solve the eigenvalues $\lambda_n=\lambda_{n,b}$ as the zeros of $\psi(b,\lambda)$. Consider the two cases: $0<\lambda\le\nu^2/2$ and $\lambda>\nu^2/2$

### Approximation for $0<\lambda\le\nu^2/2$

By setting $\lambda = (\nu^2-q^2)/2$, we are equivalently solving $q\in[0,|\nu|)$. To get precise numerical values of $q_{n,b}$, we need to find the roots of $W_{(1-\nu)/2,q/2}\left(\frac{1}{2x}\right)$ numerically. However, for large $b$ we can get estimates by using the Whittaker function
asymptotics for $\mu>0$ and $z > 0$:
$$W_{\kappa,\mu}(z)\sim\frac{\Gamma(2\mu)}{\Gamma(1/2+\mu-\kappa)}z^{1/2-\mu}e^{-z/2}$$

This gives
$$\psi(b,(\nu^2-q^2)/2)\sim\frac{\Gamma(q)}{\Gamma((\nu+q)/2)}\left(\frac{1}{2b}\right)^{(\nu-q)/2}$$

For $\nu<0$,  the reciprocal of the gamma function $1/\Gamma((\nu+q)/2)$ has zeros $(\nu+q)/2=-n+1, n = 1,...,[|\nu|/2+1]$. Thus, for $\nu < 0$ and *large enough* b, the total number of zeros in $[0,|\nu|)$ is equal to $N_\nu(b)=[|\nu|/2+1]$ and 
$$q_{n,b}=(|\nu|-2n+2), n=1,...,[|\nu|/2+1]$$

For $\nu>0$, $N_\nu(b)\equiv0, \forall b>0$

Then the non-normalized eigenfunctions are:
$$\psi(x,\lambda_{n,b})=(2x)^{(1-\nu)/2}e^{1/4x}W_{(1-\nu)/2,q_{n,b}/2}\left(\frac{1}{2x}\right)$$

### Approximation for $\lambda>\nu^2/2$

Setting $\lambda = (\nu^2+p^2)/2$. We can obtain estimates that improve with increasing n by using the following estimate of the Whittaker function with purely imaginary second index $\mu=i\rho, \rho >0, \kappa\in\mathbb{R}$ and fixed $z>0$:
$$W_{\kappa,i\rho}=\sqrt{2z}e^{\pi\rho/2}\rho^{\kappa-1/2}\cos{\left(\rho\ln{\left(\frac{z}{4\rho}\right)}+\rho-\left(\kappa-\frac{1}{2}\right)\frac{\pi}{2}\right)}\times[1+O(\rho^{-1})]$$

In this case, the root is related to the cosine function, let
$$\frac{p}{2}\left(-\ln{(4bp)}+1\right)+\frac{1}{4}\pi\nu=\frac{\pi}{2}+n\pi $$
then the solution of following equation gives the estimate of $p_{n,b}$
$$p_{n,b}[\ln{(4bp_{n,b})}-1]=2\pi\left(n+\frac{\nu}{4}-\frac{1}{2}\right)$$

Then the non-normalized eigenfunctions are:
$$\psi(x,\lambda_{N_\nu(b)+n,b})=(2x)^{(1-\nu)/2}e^{1/4x}W_{(1-\nu)/2,ip_{n,b}/2}\left(\frac{1}{2x}\right)$$

### Normalization of the eigenfunctions

The normalizing factor is related to the other linearly independent solution to the ODE. By calculating the Wronskian, define
$$\eta_{n,b}^{(\nu)}:=\frac{\partial}{\partial p}\left.W_{(1-\nu)/2,ip_{n,b}/2}\left(\frac{1}{2b}\right)\right|_{p=p_{n,b}}$$
$$\xi_{n,b}^{(\nu)}:=-\frac{\partial}{\partial q}\left.W_{(1-\nu)/2,q_{n,b}/2}\left(\frac{1}{2b}\right)\right|_{q=q_{n,b}}$$

The norms of the non-normalized eigenfunctions are:

For $\lambda\le\nu^2/2,\nu<0$ and $N_\nu>0$, the norms are given by:
$$\frac{1}{||\psi(\bullet,\lambda_{n,b})||^2}=2^\nu q_{n,b}\frac{\Gamma(\tfrac{1}{2}(\nu+q_{n,b}))}{\eta_{n,b}^{(\nu)}\Gamma(1+q_{n,b})}M_{(1-\nu)/2,q_{n,b}/2}\left(\frac{1}{2b}\right), n=1,...,N_\nu(b)$$

For $\lambda\in(\nu^2/2,\infty)$, the norms are given by:
$$\frac{1}{||\psi(\bullet,\lambda_{N_\nu(b)+n,b})||^2}=2^\nu p_{n,b}\frac{\Gamma(\tfrac{1}{2}(\nu+ip_{n,b}))}{\xi_{n,b}^{(\nu)}\Gamma(1+ip_{n,b})}M_{(1-\nu)/2,ip_{n,b}/2}\left(\frac{1}{2b}\right),n=1,2,...$$

### Spectral representation for $P_b^{(\nu)}(k,\tau)$

The transition probability density takes the form
$$p_b(t;x,y)=\sum_{n=1}^\infty e^{-t\lambda_{n,b}}\frac{\psi(x,\lambda_{n,b})\psi(y,\lambda_{n,b})}{||\psi(\bullet,\lambda_{n,b})||^2}$$

Define $c(\lambda)$ as
$$
\begin{align}
c(\lambda)&=\int_0^k(k-y)\psi(y,\lambda)\mathscr{m}(y)dy\\
&=2^{-(1+\nu)/2}k^{(\nu+3)/2}e^{-1/4k}W_{-(\nu+3)/2,\sqrt(\nu^2-2\lambda)/2}\left(\frac{1}{2k}\right)
\end{align}$$

Then the transformed payoff $P_b^{(\nu)}(k,\tau)$ have the spectral representation:
$$P_b^{(\nu)}(k,\tau)=\sum_{n=1}^\infty e^{-\tau\lambda_{n,b}}\frac{c(\lambda_{n,b})}{||\psi(\bullet,\lambda_{n,b})||^2}$$


## Proposition 1

The function $P_b^{(\nu)}(k,\tau)$ is given by the following series:
$$\begin{align}
P_b^{(\nu)}(k,\tau)\\
=&\sum_{n=1}^\infty e^{-\tau(\nu^2+p^2_{n,b})/2}p_{n,b}\frac{\Gamma(\tfrac{1}{2}(\nu+ip_{n,b}))}{4\xi_{n,b}^{(\nu)}\Gamma(1+ip_{n,b})}\\
&\times(2k)^{(\nu+3)/2}e^{-1/4k}W_{-(\nu+3)/2,(ip_{n,b})/2}\left(\frac{1}{2k}\right)M_{(1-\nu)/2,ip_{n,b}/2}\left(\frac{1}{2b}\right)\\
&+\sum_{n=1}^{N_\nu(b)} e^{-\tau(\nu^2-q^2_{n,b})/2}q_{n,b}\frac{\Gamma(\tfrac{1}{2}(\nu+q_{n,b}))}{4\eta_{n,b}^{(\nu)}\Gamma(1+q_{n,b})}\\
&\times(2k)^{(\nu+3)/2}e^{-1/4k}W_{-(\nu+3)/2,q_{n,b}/2}\left(\frac{1}{2k}\right)M_{(1-\nu)/2,q_{n,b}/2}\left(\frac{1}{2b}\right)
\end{align}$$


```python

```

fff8de9adaef2ce77fa2a8b8ab9450b7fbfbe57c
ipynb
Jupyter Notebook
ExactAsian.ipynb
yogaball/ASP2021
4985159730bccff36ec07d9e95ee384d4d7d5a09
MIT
ExactAsian.ipynb
yogaball/ASP2021
4985159730bccff36ec07d9e95ee384d4d7d5a09
MIT
ExactAsian.ipynb
yogaball/ASP2021
4985159730bccff36ec07d9e95ee384d4d7d5a09
MIT
Qwen/Qwen-72B
1. YES
2. YES