Ok got another error/s. From the output I can not tell what is cauing it?
Please advise.
Code: Select all
\documentclass{article} % Your input file must contain these two lines
\usepackage{fancyvrb}
\DefineVerbatimEnvironment{code}{Verbatim}{fontsize=\small}
\DefineVerbatimEnvironment{example}{Verbatim}{fontsize=\small}
\usepackage{natbib}
\usepackage{amsmath}
\usepackage{pdflscape}
\usepackage{amsfonts}
\usepackage{bm}
\usepackage{grffile}
\usepackage{graphicx}
\usepackage{multirow}
\usepackage[%
font=small,
labelfont=bf,
figurewithin=section,
tablewithin=section,
tableposition=top
]{caption}
\makeatletter
\def\env@matrix{\hskip -\arraycolsep
\let\@ifnextchar\new@ifnextchar
\array{*\c@MaxMatrixCols l}}
\makeatother
\begin{document} % plus the \end{document} command at the end
\parskip = 1pc %change spacing between paragraphs
\parindent = 0pc %change paragraph indentation
\pagestyle{empty}
\clearpage
\setcounter{page}{1}
\pagestyle{plain}
A close relation between the SABR and heston model exisits [14]
SABR:
CEV = space transformed squared Bessel process
Vol=GBM
Heston:
Asset: GBM
Vol=squared Bessel process
Therefore natural to extrapolate to generalise unbiased schemes for the Heston Stochastic vol model
Heston BK scehem
Beased on Willard [33] the conditional distribution given the terminal volatility and integrated variance in a time interval is log normal.
An acceptance-rejection technique is employed to tsample the variance process
A fourier inversion technique is applied to recover the variance priocess
Application hampered by slow computation speed. Anderson [1] developed 2 efficient variants
TG Truncated Gaussian
QE Quadratic exponential
Both schemes the non-central chi-squared distribution is approximated by a distribution whose moments are matched with the exact distribution
QE does not work well for SABR because the QE scheme is based in a squares Bessel process with reflecting boundary at zero vol, which gives rise to a sub-martingale process. Clearly this is not suited to SABR's asset price dynamics.
An absorbing boundary in required but the absorbing boundary behaviour in non-trivial as the transition density is not known in closed form and the absorbed process does not integrate to unity, and the moments are not known in closed form
Simulation scheme for squared Bessel Processes
Andersen and andreasen [2] Euler as well as Log-Euler MC schemes. Noted a bias in the price of caps floors and swaptions for even an infinite amount of simulations.
Jaekel [19] higher order implicit milstein scheme. Does not perform well when zero is attainable due to the discontinuity of the first derivative of the diffusion coefficient.
Campolieti and Makarov [7] exact scheme based on acceptance-rejection sampling of the Bessel bridge process. The Bessel bridge scehme is however quite complex and the computational time relatively high.
Lord et at [24] consider a Euler scheme in combination with certain rules to deal with the negative paths produced by the Euler scheme. . the4 authors conclude that the computational efficiency of their scheme is far superior to the more complicated schemes.
For certain "relevant" parameter configurations the scheme produces a significant bias for practical sizes of the time steps.
Paper proposes:
Unbiased path simulation scheme for continuous time CEV and SABR models based on Willard's [3] idea of mixing conditional distributions of a stochastic volatility model given the terminal volatility and integrated varaince.
Contribution:
1.) Derivation of conditional distribution of SABR model over a discrete time step, showing that, conditional on the terminal volatility and the integrated variance it is a space treanformed squared Bessel process with a shifted initial condition.
2.) easy to implement algorithm for squared Bessel process with absorption at zero.
3.) Provide a simple approximation formula for the conditional moments of the integrated variance by means of a small disturbance expansion, which facilitates effective sampling from the joint distribution of the terminal volatility and integrated variance
Distribution of CEV process
CEV process
\begin{equation}\label{eq:CEV}
dS_t = \sigma S_t^{\beta}dW_t
\end{equation}
consider invertible transformation
\begin{equation}
X_t = \frac{S_t^{1-\beta}}{(1-\beta)}
\end{equation}
for $\beta \not= 0$
Applying Ito's lemma
\begin{equation}
\begin{split}
dX_t & = (1-\beta)\frac{S_t^{- \beta}}{(1-\beta)} \sigma S_t^\beta dW_{1,t}-\frac{1}{2} \beta (1- \beta) \frac{s_t^{-1-\beta}}{1-\beta} \sigma^2 S_t^{2 \beta} dt \\
& = \sigma dW_t - \frac{\beta \sigma^2}{(2-2\beta)X_t}dt \\
\end{split}
\end{equation}
which is a time-changed Bessel process
A second transformation is then defined as $T_t = X_t^2 $ which results in a time-changes squared Bessel process of dimension $ \delta := \frac{(1-2 \beta)}{(1-\beta)}$ that satisfies the following SDE:
\begin{equation}
dY_t = 2 \sqrt{|Y_t|} \sigma dW_t + \delta \sigma^2 dt
\end{equation}
let $\nu(t)$ be a time-change function so that $\nu(t) = \sigma^2 t $. Then, $Y_t=Z_{\nu(t)}$, where ${Z_t}$ is a $\delta-$di,emsional squared Bessel process, i.e., the strong solution of the SDE:
\begin{equation}\label{eq:solve}
dZ_t = 2 \sqrt{|Z_t|} \sigma dW_t + \delta dt
\end{equation}
with degree of freedom $\delta$. The squared Bessel process is a Markov process and its transition densities are known explicitly.
Then sample random numbers from the analystic transition density in Z-space and to apply an inverse variabe transformation to obtain the random numbers in the original coordinates.
For a standard squared Bessel process, as defined by SDE \ref{eq:solve} the following statements hold:
1.) All solutions to \ref{eq:solve} are non explosive
2.) for $\delta < 2, Z=0$ is an attainable boundary for process \ref{eq:solve}
3.) for $\delta \geq 2$ SDE \ref{eq:solve} has a unique solution and zero is not attainable
4.) for $0<\delta<2$ the SDE \ref{eq:solve} does not have a unique solution, unless the boundary condition is specified for the solution to Eq \ref{eq:solve} at Z=0
5.) for $\delta \leq 0$, thers is a unique strong solution to the SDE \ref{eq:solve} and the boundary condition zero is absorbing
For the latter two cases the transition densities are known:
Transition density for Squared Bessel process $q^{\delta}(t,x,y)$
1.) For $\delta \leq 0$ and for $0 < \delta < 2$ in Eq \ref{eq:solve} but only when the boundary is absorbing at y=0:
\begin{equation}
q^{\delta}(t,x,y) = \frac{1}{2t}(\frac{y}{x})^{\frac{\delta-2}{4}} exp(-\frac{x+y}{2t})I_{|\frac{\delta-2}{2}|}(\frac{\sqrt{xy}}{t}), y>=0, t>0
\end{equation}
2.) For $0 < \delta < 2$ when y=0 is a reflecting boundary:
\begin{equation}
q^{\delta}(t,x,y) = \frac{1}{2t}(\frac{y}{x})^{\frac{\delta-2}{4}} exp(-\frac{x+y}{2t})I_{\frac{\delta-2}{2}}(\frac{\sqrt{xy}}{t}), y>=0, t>0
\end{equation}
Where $I_a(x)$ denotes the Bessel function, defined by
\[
I_a(x):= \sum_{j=0}^{\infty}{\frac{(x/2)^{2j+a}}{j!\Gamma(a+j+1)}}
\]
where the Gamma function $\Gamma(x):= \int_0^{\infty} u^{x-1}e^{-u}du$
\begin{table}[hp]
\centering
\begin{tabular}{|c| |c|}\hline
CEV exponent & Squared Bessel $\delta$\\\hline
$-\infty < \beta < \frac{1}{2}$ & $0< \delta <2$ \\
$\frac{1}{2} \leq \beta < 1$ & $-\infty < \delta \leq 0$\\
$\beta > 1 $& $2< \delta <\infty$ \\
\hline
\end{tabular}
\end{table}
Given the transition density for the squares Bessel processes in X space given in result 2.2 one can easily obtain the density for the CEV process 2.2 in S-space
\[
S_{\Delta}=((1-\beta)\sqrt{|Z_{\nu(\Delta)}|})^{\frac{1}{1-\beta}}
\]
define mapping
\[
h:s\rightarrow((1-\beta)\sqrt{s})^{\frac{1}{1-\beta}}, s \geq 0,
\]
with inverse
\[
h^{-1}:y \rightarrow \frac{y^{2(1-\beta)}}{(1-\beta)^2}, y \geq 0,
\]
So $S_{\Delta} = h(Z_{\nu{\Delta}})$ and $Z_0=h^{-1}(S_0)=S_0^{2(1-\beta)}/(1-\beta)^2$. Then $Z_{\nu{\Delta}}$ has density $q^{\delta}(t,x,y)$ and it follows that the density for $S_{\Delta}$ is given by
\[
p(S|S_0)=q^{\delta}(\nu{\Delta}),Z_0,h^{-1}(s))\frac{dh^{-1}(s)}{dS}
\]
Hence
1.) For $0<\beta< \frac{1}{2}$ with absorbtion at 0 and for $\frac{1}{2}<\beta<1 $:
\[
p(S|S_0) = \frac{1}{\nu{\Delta}}(\frac{S_{\Delta}}{S_0})^{-\frac{1}{2}} exp(-\frac{S_{\Delta}^{2(1-\beta)}+S_0^{2(1-\beta)}}{2(1-\beta)^2 \nu(\Delta)})I_{|\frac{\delta-2}{2}|}(\frac{(S_0 S_{\Delta})^{1-\beta}}{\nu(\Delta)(1-\beta)^2})\frac{s_{\Delta}^{1-2\beta}}{1-\beta},
\]
where $\nu(\Delta) =\sigma^2\Delta$ and $\delta=\frac{1-2 \beta}{1-\beta}$
2.) For $0 < \beta < \frac{1}{2}$ with a reflecting boundary at S=0:
\[
p(S|S_0) = \frac{1}{\nu{\Delta}}(\frac{S_{\Delta}}{S_0})^{-\frac{1}{2}} exp(-\frac{S_{\Delta}^{2(1-\beta)}+S_0^{2(1-\beta)}}{2(1-\beta)^2 \nu(\Delta)})I_{\frac{\delta-2}{2}}(\frac{(S_0 S_{\Delta})^{1-\beta}}{\nu(\Delta)(1-\beta)^2})\frac{s_{\Delta}^{1-2\beta}}{1-\beta},
\]
By integrating these indentities, we find the cumulative distribution functions:
The cumulative distribution function for the CEV price process is given by
1.) For $0<\beta< \frac{1}{2}$ with absorbtion at 0 and for $\frac{1}{2}<\beta<1 $:
\begin{equation}
Pr(S_{\Delta} \leq x| S_0) = 1-\chi^2(a;b,c).
\end{equation}
2.) For $0 < \beta < \frac{1}{2}$ with a reflecting boundary at S=0:
\begin{equation}
Pr(S_{\Delta} \leq x| S_0) = 1-\chi^2(c;2-b,a).
\end{equation}
with the following parameters:
$a=\frac{s+0^{2(1-\beta)}}{(1-\beta)^2\nu(\Delta)}$, $b=\frac{1}{1-\beta}$,$c=\frac{x^{2(1-\beta)}}{(1-\beta)^2\nu(\Delta)}$, $\nu(\Delta)=\sigma^2\Delta$,
and $\chi^2(x;\delta,\lambda)$ is the non-central chi-square distribution for the random variable x with non-centrality parameter $\lambda$ and degree of freedom $\delta.$
The density will not integrate to unity when the boundary is absorbing.
Collary For $0<\beta<1$, the probability of $S_{\Delta}$ given by SDE \ref{eq:CEV} and intital condition $S_0$ reads
\begin{equation}
Pr(S_{\Delta}=0| S_0) =1-\gamma(\frac{1}{2(1-\beta)},\frac{S_0^{2(1-\beta)}}{2(1-\beta)^2 \nu(\Delta)})/\Gamma(\frac{1}{2(1-\beta)}),
\end{equation}
where $\gamma(\alpha,\beta)$ is the lower complete Gamma function and $\Gamma(alpha)$ is the Gamma function
Fig.1 shows two plots comparing the exact cumulative distribution of the CEV process versus the lognormal and normal approximations at T=0.25 for $S_{\Delta}$, given two differnt starting values $S_0=6\%$ and $S_0=2\%$. The model parameters used were $\theta=0.3$ and $\beta=0.4$. The labels 'absorbing' and 'reflecting' distinguish between the CEV process with absorbing and reflecting boundaries repectively. The two matched distributions so not represent accurate approximations of the true distributions of $S_{\Delta}$.
To keep the exposition arbitrage free only the absorbing boundary condition is considered here,
SABR Conditional Distribution
The SABR model
\begin{equation}\label{eq:SABR}
\begin{split}
dS_t & = \sigma_t S_t^{\beta}dW_{1,t}, \\
d \sigma_t & = \alpha \sigma_t dW_{2,t}, \\
\end{split}
\end{equation}
RESULT Conditional in the level of terminal volatility and $\sigma_{\Delta}$ and the integrated variance $\int_0^{\Delta}{\sigma_s ds}$ path
1.) For an invertibale transformant $X(S)=S^{1-\beta}/(1-\beta)$ application of Ito's lemma gives
\begin{equation}
X_{\Delta} = X_0 +\frac{\rho}{\alpha}(\sigma_{\Delta}-\sigma_0)+\sqrt{1-\rho^2}\int_0^{\Delta}{\sigma_s dU_s}-\int_0^{\Delta}{\frac{\beta \sigma_s^2}{(2-2 \beta)X_s}ds},
\end{equation}
Where $U_s$ is a standard brownian motion independent of $W_{2,s}$ in system \ref{eq:}SABR
2.) $Y(S)=S^{2-2\beta}/(1-\beta)^2$ and application of the time change $\nu(t)=(1=\rho^2)\int_0^t{\sigma_s^2ds,}$, $Y_t$ is a squared Bessel process of dimension $\frac{1-2\beta-\rho^2(1-\beta)}{(1-\beta)(1-\rho^2)}$ solvig the SDE:
\begin{equation}
dY_{\nu(t)} = 2 \sqrt{Y_{\nu(t)}}dU_{\nu(t)} + \frac{1-2\beta-\rho^2(1-\beta)}{(1-\beta)(1-\rho^2)}d\nu(t),
\end{equation}
with the inital condition $Y_0 = (\frac{S_0^{1-\beta}}{1-\beta}+\frac{\rho}{\alpha}(\sigma_{\Delta}-\sigma_0))^2$
3.) Let $\tau$ be the stopping time for which process Y hits zero, i.e. $\tau=inf{\nu(s)|Y_{\nu(s)}=0}$, the stopped priocess is then
\begin{equation}
dY_{\nu(t) \wedge \tau} = 2 \int_0^{\nu(t)\wedge \tau}\sqrt{Y_{\nu(t)}}dU_{\nu(t)} + \frac{1-2\beta-\rho^2(1-\beta)}{(1-\beta)(1-\rho^2)}(nu(t)\wedge \tau),
\end{equation}
Presposition (Cumulative Distibution for the SABR process) For some $S_0$, strictly larger than 0, the conditional cumulative distribution of $S_{\Delta}$ with an absorbing boundary at $S_t=0$ given $\sigma_{\Delta}$ and $\int_0^{\Delta}{\sigma_s^2 ds}$ reads
\begin{equation}
P(S_{\Delta} \leq K| S_0 >0, \sigma_{\Delta}, \int_0^{\Delta}{\sigma_s^2 ds})1-\chi^2(a;b;c),
\end{equation}
where
\begin{equation}
\begin{split}
&a=\frac{1}{\nu(\Delta)}(\frac{S_0^{1-\beta}}{(1-\beta)}+\frac{\rho}{\alpha}(\sigma_{\Delta}-\sigma_0))^2 , b=2-\frac{1-2 \beta - rho^2(1-\beta)}{(1-\beta)(1-\rho^2)}, \\ &c =\frac{K^{2(1-\beta)}}{(1-\beta)^2 \nu(\Delta)}, \nu(\Delta)=(1-\rho^2)\int_0^{\Delta}{\sigma_s^2 ds} \\
\end{split}
\end{equation}
$chi^2(x;\delta,\lambda)$ is again the noncentral chi-square distribution function.
REMARK It was argues by andersen[4] that the contimuous time process, $S_t$, will be a martingale with
\[
E[S_{t+\Delta}|S_t]<\infty
\]
the equivalent discrete time process, $\hat{S_t}$, generated by the unbiased simulation scheme may not satisfy the martingale condition
\[
E[\hat{S_{t+\Delta}}|\hat{S_t}] \neq \hat{S_t}
\]
The net drift away from the martingale is visible for parameter sets with small Beta and close to zero rates. However (according to Osterlee) this combination of parameters does not appear in practical applications as it gives rise to impractical implied volatiltiy levels. For practical SABR parameters, the martingale bias is very small and can be controlled by reducing the size of the time step.
The Euler, Milstein and Log-Euler discreteisation schemes are considered. It is noted that although log-euler preserves non-negativity of the asset price process, numerical experiments show that the scheme may become unstable for specific time0step sizes. The instabilities occur beause the diffusion terms below approach infinity quickly when $\hat{S_t}$ reaches zero.
\begin{table}[hp]
\centering
\begin{tabular}{|l| |l| |l|}
\hline
Issue & Stepsize & Euler & Milstein & Log-Euler \\\hline
\multirow{3}{*}{Negativity} & \Delta =0.5 & 85\% & 49\% & 0\% \\
& \Delta =0.25 & 88\% & 56\% & 0\% \\
& \Delta =0.125 & 94\% & 59\% & 0\% \\\hline
\multirow{3}{*}{Infinity} & \Delta =0.5 & 0\% & 34\% & 96\% \\
& \Delta =0.25 & 0\% & 31\% & 96\% \\
& \Delta =0.125 & 0\% & 37\% & 95\% \\\hline
\hline
\end{tabular}
\end{table}
Local Linearisation Shoji and Ozaki [32] are also considered but are not followed up becuase the resultant density is normal giving rise to non-zero probability of becoming negative-valued. Also the drift term diverges for $x_0$ close to zero.
