---
abstract: |
    We develop stock option price approximations for a model which takes both the risk of default and the stochastic volatility into account. We also let the intensity of defaults be influenced by the volatility. We show that it might be possible to infer the risk neutral default intensity from the stock option prices. Our option price approximation has a rich implied volatility surface structure and fits the data implied volatility well. Our calibration exercise shows that an effective hazard rate from bonds issued by a company can be used to explain the implied volatility skew of the option prices issued by the same company. We also observe that the implied yield spread that is obtained from calibrating all the model parameters to the option prices matches the observed yield spread.

    **Key Words:** Option pricing, multiscale perturbation methods, defaultable stocks, stochastic intensity of default, implied volatility skew.
author:
- 'Erhan Bayraktar [^1]'
title: '**Pricing Options on Defaultable Stocks [^2] [^3]**'
---

Introduction
============

Bonds and stocks issued by the same company carry the risk of the company going bankrupt. Although this risk has been accounted for in pricing the bonds, see e.g. [@schoenbucher] for a review of this literature, there is only a limited number of models in equity option pricing, [@carr-linetsky], [@carr-wu], [@hull-nelken] and [@Linetsky], [@merton-jump]. Instead, models are developed to match the implied volatility curve, which is observed daily: stochastic volatility models, see e.g. [@sircar], and jump diffusion models, see e.g. [@cont]. These models account for the default risk implicitly. These approaches could be improved accounting for the risk of default. When a company goes bankrupt and the stock of the company plunges, the out-of-the-money put options are suddenly very valuable. So, when these contracts are priced, the market must also be accounting for these crash scenarios, besides the stochastic volatility and small jump effects, and this is partly why the out of the money put options are worth much more than what the Black and Scholes model, which does not account for bankruptcy, would predict.

In addition, one might argue that the default risk should be priced consistently across the markets. When the likelihood of the default of a company increases (for example when a credit downgrade occurs), it is not unreasonable to expect that besides the increase in the yield spread, the trading volume in the out of the money put options increases and their prices go up. In fact, the more the put option is out of the money, the higher the return it will bring when there is a default, so the higher the prices should be. This is the understanding of [@Linetsky] and the convertible bond pricing literature(see e.g. [@vetzal]).

In this note, we will use an intensity based approach to model the default risk of the company that issues the stock. We will also allow the volatility of the stock price to be stochastic. We observe that the implied yield spread that is obtained from calibrating all the model parameters to the option prices matches the observed yield spread (Figure 15). On the other hand, our calibration exercise shows that an effective hazard rate from bonds issued by a company can be used to explain the implied volatility skew of the option prices issued by the same company.

In our model, the volatility is driven by two processes that evolve on two different time scales, fast and slow. The intensity of default is a function of volatility and an idiosyncratic stochastic component (since the credit market and stock options market can have independent movements), which is driven by two other processes that evolve on fast and slow time scales. Our model is similar to that of [@carr-wu] and the motivation for our model specification can be found on page 4 of that paper. The difference of our model from the Carr-Wu specification is that we also incorporate the market prices of volatility and intensity risks into our model, and we consider the intensity and volatility to be evolving on two different time scales. The purpose of our specification is to be able to obtain asymptotic option price approximations using multiscale perturbation methods developed by [@ronnie-timescale], who used this technique in the context of stochastic volatility models. Our model can be thought of as a hybrid of the stochastic volatility model of [@ronnie-timescale] (whose purpose is to obtain approximations for stock option prices) and the multi-scale intensity model of [@papa] (whose purpose is to price credit derivatives). But one should note that following the Carr-Wu specification we allow the intensity to depend on the volatility, and we also incorporate the market price of default risk into our model. Let us also point out that our default model differs from [@carr-linetsky] and [@Linetsky] since they consider a model in which the intensity of the defaults depends solely on the stock price.

Our price approximation for the option prices has seven parameters (we will refer to it as the 7 parameter model), one of which is the average intensity, that can be calibrated to the observed implied volatility. We compare the performance of this model against 3 simpler models: \[5 parameter model\] a model that does not account for the market price of intensity risk; \[the model of [@ronnie-timescale]\] a model that does not account for the default risk; \[3 parameter model\] a model that does not account for the stochastic volatility effects. When we take the average intensity to be equal to the yield spread of the corporate bond with smallest maturity and calibrate the rest of the parameters to the observed option prices we observe that the seven parameter model not only outperforms the others (which is to be expected) but also that it almost fits the data implied volatility perfectly for longer term options (9 months and over). (In our calibration we pool the options with different strikes and different maturities in the same pool and we have 104 data points, we consider maturities from 3 months to 2 years). The 5 parameter model always outperforms the model of [@ronnie-timescale], which points to the significance of accounting for the default risk. Although, the three parameter model has only two free parameters (the other free parameter is the average intensity of default and is already fixed to be the yield spread of the shortest maturity bond) it fits the implied volatility surface fairly well, considering that we have 104 data points. See Section 3.2.

When we calibrate all parameters of each model, including the average intensity, to the option prices then all the models we consider (except the model of [@ronnie-timescale] since this model does not account for the default risk) fit the the data implied volatility curve quite well (even the 3 parameter model). But when option implied average intensity is plotted against the yield curve of the shortest maturity bond over the entire period of the data (about a year), one observes that all models except the 7 parameter model have unrealistically severe option implied intensities. Only the 7 parameter model almost matches the yield spread of the shortest maturity bond (the maturity is about two years on the last day of our data). Therefore, it might be possible to derive the risk neutral probability of default from the option prices using the 7 parameter model. See Section 3.3. Also we should note that the 7 and 5 parameter models have much richer implied volatility surface structure (see Figures \[fig:s3p\], \[fig:s5p\], \[fig:sfp\], and \[fig:s7p\]) and Section 2.3, which explains why they fit to the option prices better.

We should note that we calibrate the parameters of our model to the daily observed option prices as in [@ronnie-timescale] and [@papa]. [@carr-wu] on the other hand perform a time series analysis to obtain the parameters of their model.

The rest of the paper is organized as follows: In Section 2, we introduce our model specification. In Section \[sec:eop\], we derive an approximation of the prices of the equity options, and in Section 2.2 we derive an approximation to the corporate bond prices, using multi-scale perturbation techniques. In Section 2.3 we analyze the model implied volatility surface structure of the models we consider. In Section 3, we calibrate the 7 parameter model and the other reference models to the data.

A Multiscale Intensity Model
============================

We will use the doubly-stochastic Poisson process framework of [@bremaud], see also [@schoenbucher]. Let $(\Omega,
\mathcal{H},\mathbb{P})$ be a filtered probability space hosting standard Brownian motions $\vec{W}=(W_t^{(0)}, W_t^{(1)}, W_t^{(2)}, W_t^{(3)}, W_t^{(4)})$, whose correlation structure is given by (2.8) and standard Poisson process $N=(N(t))_{t \geq 0}$ independent of $\vec{W}$. These Brownian motion will be used in modeling the the volatility and the intensity (see (2.6)).We model the time of default, which we will denote by $\tau$, as the first time the process the time changed Poisson process $\tilde{N}_t
\triangleq N\left(\int_0^{t}\lambda_s ds\right)$, $t \geq 0$, jumps.

The stock price dynamics in our framework is the solution of the stochastic differential equation $$d\tilde{S}_t=\tilde{S}_t\left(r_t dt+\sigma_t
dW^{(0)}_t-d\left(\tilde{N}_t-\int_0^{t \wedge \tau}\lambda_u
du\right)\right), \quad \tilde{S}_0=x,$$ in which we take the interest rate $(r_t)_{t \geq 0}$ to be deterministic as Carr and Wu (2006) and Linetsky (2006) do. At the time of default the stock price jumps down to zero. Note that the discounted stock price is a martingale under the measure $\mathbb{P}$. Let $(\mathcal{F}_t)_{t \geq 0}$ be the natural filtration of $\vec{W}$. The price of a European call option with maturity $T>0$ is equal to $$\label{eq:price-option}
\begin{split}
C(t)=\mathbb{E}\left[B(t,T)(\tilde{S}_T-K)^+1_{\{\tau>T\}}
\bigg|\mathcal{F}_t \vee \sigma\{\tilde{N}_s; 0 \leq s \leq t\}\right]=\mathbb{E}\left[B(t,T) e^{-\int_t^{T}\lambda_s
ds}(S_T-K)^+\bigg|\mathcal{F}_t\right],
\end{split}$$ in which $B(t,T)$ is the price of the bond that matures at time $T$, whereas the price of a European put option is given by $$\label{eq:price-option-put}
\begin{split}
P(t)&=\mathbb{E}\left[B(t,T)((K-\tilde{S}_T)^+1_{\{\tau>T\}}+K1_{\{\tau
\leq T\}}) \bigg|\mathcal{F}_t \vee \sigma\{\tilde{N}_s; 0 \leq s \leq t\}\right]
\\&=KB(t,T)-\mathbb{E}\left[B(t,T)e^{-\int_t^{T}\lambda_s
ds}\min(K,S_T)\bigg|\mathcal{F}_t\right],
\end{split}$$ in which $S=(S_t)_{t \geq 0}$ is the solution of $$\label{eq:dyn-stock}
dS_t=S_t(r+\lambda_t)+S_t \sigma_t dW^{(0)}_t, \quad S_0=x.$$ From (\[eq:price-option\]) and (\[eq:dyn-stock\]) we see that pricing a call option on a defaultable stock is equivalent to pricing a call option in a fictitious market that does not default but whose interest rate is stochastic (also see [@Linetsky]). Using (\[eq:price-option\]) and (\[eq:price-option-put\]), the put-call parity is satisfied, $$C(t)-P(t)=S_t-B(t,T)K,$$ as expected.

We will use a model specification similar to that of [@carr-wu] and assume the following dynamics for the volatility $\sigma_t$ and the intensity $\lambda_t$: $$\label{eq:model}
\begin{split}
\sigma_t&=\sigma(Y_t,Z_t), \quad \lambda_t=\beta \sigma_t^2+f(Q_t, U_t),
\end{split}$$ in which $$\begin{split}
dY_t&=\left(\frac{1}{\epsilon}(m-Y_t)-\frac{v \sqrt{2}}{\sqrt{{\epsilon}}}\Lambda(Y_t,Z_t)\right)dt+\frac{v \sqrt{2}}{\sqrt{\epsilon}}dW_t^{(1)}, \quad Y_0=y,
\\ dZ_t&=\left(\delta c(Z_t)-\sqrt{\delta}g(Z_t)\Gamma(Y_t,Z_t)\right)dt+\sqrt{\delta}g(Z_t)dW_t^{(2)}, \quad Z_0=z,
\\dQ_t&=\left(\frac{1}{\epsilon}(\tilde{m}-Q_t)-\frac{{\tilde{v}}\sqrt{2}}{\sqrt{{\epsilon}}}{\widetilde{\Lambda}}(Q_t,U_t)\right)dt+\frac{\tilde{v}
\sqrt{2}}{\sqrt{\epsilon}}dW_t^{(3)}, \quad Q_0=q,
\\  dU_t&=\left(\delta {\tilde{c}}(U_t)-\sqrt{\delta}{\tilde{g}}(Z_t){\widetilde{\Gamma}}(Q_t,U_t)\right)dt+\sqrt{\delta}{\tilde{g}}(U_t)dW_t^{(4)}, \quad U_0=u.
\end{split}$$ Here, ${\Lambda}$ and ${\Gamma}$ are the market prices of volatility risk, ${\widetilde{\Lambda}}$ and ${\widetilde{\Gamma}}$ are the market prices of intensity risk. We will assume the following correlation structure for the Brownian motions in our model: $$\begin{split}
E[W^{(0)}_t,W^{(i)}_t]&=\rho_i t ,\quad i=1,  2,\\
E[W^{(1)}_t,W^{(2)}_t]&=\rho_{12} t, \quad E[W^{(3)}_t,W^{(4)}_t]=\rho_{34} t, \\
E[W^{(i)}_t,W^{(3)}_t]&=0, \quad i=0, 1, 2, \quad  E[W^{(i)}_t,W^{(4)}_t]=0, \quad i=0, 1, 2,
\end{split}$$ for all $t \geq 0$. This model captures the covariation of the implied volatility and the spread of risky bonds, and it also accommodates the fact that the credit market and the stock market can show independent movements (see Section 1.1 of [@carr-wu] for further motivation). One should note that when $\beta=0$, our model specification for the intensity coincides with that of [@papa]. On the other hand, when $\lambda=0$ our model specification becomes that of [@ronnie-timescale]. We assume that the market price of risks are bounded functions as in [@ronnie-timescale].

As in [@ronnie-timescale] and [@papa] we assume that both $\sigma$ and $f$ are bounded, smooth, and strictly positive functions. We also assume that the differential equation defining $Z$ and $U$ have unique strong solutions. We should also note that we think of $\varepsilon>0$, $\delta>0$ as small constants. Hence, the process $Y$ and $Q$ are fast mean reverting and $Z$ and $U$ evolve on a slower time scale. See [@ronnie-timescale] for an exposition of multi-scale modeling in the context of stochastic volatility models.

European Option Price {#sec:eop}
---------------------

From an application of the Feynman-Kac Theorem it follows that the price of a call option satisfies the following partial differential equation $$\label{eq:option-pricing-pde}
\begin{split}
{\mathcal{L}}^{{\epsilon},\delta}C^{{\epsilon},\delta}(t,x,y,z,q,u)&=0, \quad t <T,
\\ C^{{\epsilon},\delta}(T,x,y,z,q,u)&=(x-KB(t,T))^+,
\end{split}$$ where the operator ${\mathcal{L}}^{{\epsilon},\delta}$ is given by $${\mathcal{L}}^{{\epsilon},\delta}\triangleq \frac{1}{{\epsilon}}
{\mathcal{L}}_0+\frac{1}{\sqrt{{\epsilon}}}{\mathcal{L}}_1+{\mathcal{L}}_2+\sqrt{\delta}{\mathcal{M}}_1+\delta {\mathcal{M}}_2+
\sqrt{\frac{\delta}{{\epsilon}}}{\mathcal{M}}_3,$$ in which $$\label{eq:operators}
\begin{split}
{\mathcal{L}}_0 & \triangleq(m-y)\frac{\partial}{\partial y}+v^2\frac{\partial^2}{\partial y^2} +(\tilde{m}-q)\frac{\partial}{\partial q}+\tilde{v}^2\frac{\partial^2}{\partial q^2},\\
{\mathcal{L}}_1 &\triangleq  (\rho_1\sigma(y,z)v\sqrt{2})x\frac{\partial^2}{\partial x \partial y}-{\Lambda}(y,z) v \sqrt{2} \frac{\partial}{\partial y}-{\widetilde{\Lambda}}(q,u){\tilde{v}}\sqrt{2} \frac{\partial}{\partial q},\\
{\mathcal{L}}_2 &\triangleq \frac{\partial}{\partial t}+\frac{1}{2}\sigma^2(y,z)x^2\frac{\partial}{\partial x^2}+(\beta \sigma^2(y,z)+f(q,u))x\frac{\partial}{\partial x}-(\beta \sigma^2(y,z)+f(q,u))\cdot,\\
{\mathcal{M}}_1 &\triangleq  \rho_2\sigma(y,z)g(z)x\frac{\partial^2}{\partial x \partial z}-g(z){\Gamma}(y,z)\frac{\partial}{\partial z}-{\tilde{g}}(z){\widetilde{\Gamma}}(q,u)\frac{\partial}{\partial u},\\
{\mathcal{M}}_2 &\triangleq  c(z)\frac{\partial}{\partial z}+\frac{1}{2}(g(z))^2\frac{\partial^2}{\partial z^2}+\tilde{c}(z)\frac{\partial}{\partial z}+\frac{1}{2}({\tilde{g}}(z))^2\frac{\partial^2}{\partial z^2},\\
{\mathcal{M}}_3 &\triangleq \rho_{12}v\sqrt{2}g(z)\frac{\partial^2}{\partial y
\partial z}+\rho_{34}\tilde{v}\sqrt{2}{\tilde{g}}(z)\frac{\partial^2}{\partial q
\partial z}.
\end{split}$$ We have denoted the price of a call option by $C^{{\epsilon},\delta}$ to emphasize the dependence on the parameters ${\epsilon}$ and $\delta$.

Here, the operator ${\mathcal{L}}_2$ is the Black-Scholes operator corresponding to interest rate level $\beta \sigma^2(y,z)+f(q,u)$ and volatility $\sigma(y,z)$. $(1/{\epsilon}) {\mathcal{L}}_0$ is the infinitesimal generator of the Ornstein-Uhlenbeck process $(Y,Q)$. ${\mathcal{L}}_1$ contains the mixed partial derivative due to the correlation between the Brownian motions driving the stock price and the fast volatility factor. ${\mathcal{M}}_1$ contains the mixed partial derivative due to the correlation between the Brownian motions driving the stock price and the slow volatility factor. $\delta {\mathcal{M}}_2$ is the infinitesimal generator of the two-dimensional process $(Z,U)$. ${\mathcal{M}}_3$ contains the mixed partial derivatives due to the correlation between the Brownian motions driving the fast and slow factors of the volatility and the intensity.

**Approximating the Price Using the Multiscale Perturbation Method**

We will first consider an expansion of the price $C^{{\epsilon},\delta}$ in powers of $\sqrt{\delta}$, $$\label{eq:delta-expansion}
C^{{\epsilon},\delta}=C_{0}^{{\epsilon}}+\sqrt{\delta}C_1^{{\epsilon}}+\delta
C_2^{{\epsilon}}+\cdots,$$ and then consider the expansion of each term in (\[eq:delta-expansion\]) in powers of $\sqrt{{\epsilon}}$, i.e., $$\label{eq:expansion-eps}
C_k^{{\epsilon}}=C_{0,k}+\sqrt{{\epsilon}}C_{1,k}+{\epsilon}C_{2,k}+ \cdots, \quad
\text{for all $k \in \mathbb{N}$}.$$ It can be seen from (\[eq:option-pricing-pde\]) and (\[eq:delta-expansion\]) that $C^{{\epsilon}}_0$ and $C^{{\epsilon}}_1$ satisfy $$\label{eq:C-0-eps}
\begin{split}
\left(\frac{1}{{\epsilon}}{\mathcal{L}}_0+\frac{1}{\sqrt{{\epsilon}}}{\mathcal{L}}_1+{\mathcal{L}}_2\right)C^{{\epsilon}}_0&=0,
\\ C^{{\epsilon}}_0(T,x,y,z)&=(x-K)^+, \quad \text{and}
\end{split}$$ $$\label{eq:C-1-eps}
\begin{split}
\left(\frac{1}{{\epsilon}}{\mathcal{L}}_0+\frac{1}{\sqrt{{\epsilon}}}{\mathcal{L}}_1+{\mathcal{L}}_2\right)C^{{\epsilon}}_1&=-\left({\mathcal{M}}_1+\frac{1}{\sqrt{{\epsilon}}}{\mathcal{M}}_3\right)C_{0}^{{\epsilon}},
\\ C_1^{{\epsilon}}(T,x,y,z)&=0.
\end{split}$$ Equating the $1/{\epsilon}$ and $1/\sqrt{{\epsilon}}$ terms in (\[eq:C-0-eps\]), we obtain $$\label{eq:eps-sqrt-eps}
{\mathcal{L}}_0 C_{0,0} =0, \quad {\mathcal{L}}_0 C_{1,0}+ {\mathcal{L}}_1 C_{0,0}=0.$$ Since the first equation in (\[eq:eps-sqrt-eps\]) is a homogeneous equation in $y$ and $q$, we take $C_{0,0}=C_{0,0}(t,x,z,u)$ to be independent of $y$ and $q$. (The $y$ and $q$ dependent solutions have exponential growth at infinity ([@sircar-proof]).) Similarly, from the second equation, it can be seen that $C_{1,0}$ can be taken to be independent of $y$ and $q$, since ${\mathcal{L}}_1 C_{0,0}=0$. When we equate the order one terms and use the fact that ${\mathcal{L}}_1
C_{1,0}=0$, we get $$\label{eq:Poisson-eq-1}
{\mathcal{L}}_0C_{2,0}+{\mathcal{L}}_2 C_{0,0}=0,$$ which is known as a Poisson equation (for $C_{2,0}$). (See [@sircar].) The solvability condition of this equation requires that $$\label{eq:Bs-eqn-average-int-rate}
\left<{\mathcal{L}}_2\right>C_{0,0}=0, \quad \text{in which}$$ $$\left<{\mathcal{L}}_2\right> \triangleq \frac{\partial}{\partial t}+\frac{1}{2}\left<\sigma^2(y,z)\right>x^2\frac{\partial}{\partial x^2}+(\beta \left<\sigma^2(y,z)\right>+\left<f(q,u)\right>)x\frac{\partial}{\partial x}-(\beta \left<\sigma^2(y,z)\right>+\left<f(q,u)\right>)\cdot,$$ where $\left<\sigma^2(y,z)\right>$, which is only a function of $z$, is the average (expectation) with respect to the invariant distribution of $Y$, which is Gaussian with mean $m$ and standard deviation $\nu$, and $\left<f(q,u)\right>$, which is only a function of $u$, is the average with respect to the invariant distribution of $Q$, which is Gaussian with mean $\tilde{m}$ and standard deviation $\tilde{\nu}$. All of this is equivalent to saying that the average $\left<{\mathcal{L}}_2\right>$ of the Black-Scholes operator ${\mathcal{L}}_2$ is with respect to the invariant distribution of $(Y,Q)$, whose density is given by $$\begin{gathered}
\Psi(y,q) \triangleq \frac{1}{2\pi v\tilde
v}\exp\left\{-\frac{1}{2}\left[\left(\frac{y-m}{v}\right)^2+\left(\frac{q-\tilde{m}}{\tilde{v}}\right)^2\right]\right\}.\end{gathered}$$

The solution to (\[eq:Bs-eqn-average-int-rate\]) is the Black-Scholes call option price when the volatility and the interest rate are taken to be $\left<\sigma^2(y,z)\right>$ and $\beta \left<\sigma^2(y,z)\right>+\left<f(q,u)\right>$, $$\label{BS-formula}
\begin{split}
C_{0,0}(t,x,z,u)&= C_{\text{BS}}\left(x; \sqrt{\left<\sigma^2(y,z)\right>},
\beta \left<\sigma^2(y,z)\right>+\left<f(q,u)\right>;K B(t,T),T-t\right) \\ & \triangleq x
N(\tilde{d}_1)-KB(t,T)\exp\left(-(\beta \left<\sigma^2(y,z)\right>+\left<f(q,u)\right>)(T-t)\right)N(\tilde{d}_2),
\end{split}$$ where $$\label{eq:d-1-tilde}
\tilde{d}_{1,2}=\frac{\log(x/[K B(t,T)])+(\beta \left<\sigma^2(y,z)\right>+\left<f(q,u)\right>\pm\frac{1}{2}
\left<\sigma^2(y,z)\right>)(T-t)}{\sqrt{\left<\sigma^2(y,z)\right>
(T-t)}}.$$

From (\[eq:Poisson-eq-1\]) and (\[eq:Bs-eqn-average-int-rate\]), we obtain that $$C_{2,0}=-{\mathcal{L}}^{-1}_0 \left({\mathcal{L}}_2-\left<{\mathcal{L}}_2\right>\right)C_{0,0}.$$ Equating the order $\sqrt{{\epsilon}}$ terms in (\[eq:C-0-eps\]) gives $${\mathcal{L}}_0 C_{3,0}+ {\mathcal{L}}_1 C_{2,0}+ {\mathcal{L}}_2 C_{1,0}=0,$$ which is again a Poisson equation for $C_{3,0}$ and its solvability condition requires that $$\label{eq:av-l2c-10}
\left<{\mathcal{L}}_2\right>C_{1,0}=-\left<{\mathcal{L}}_1 C_{2,0}\right>=\left<{\mathcal{L}}_1 {\mathcal{L}}_0^{-1}({\mathcal{L}}_2-\left<{\mathcal{L}}_2\right>)\right>C_{0,0}.$$ Observe that we have $C_{1,0}(T,x,z,u)=0$. We can compute $C_{1,0}$ explicitly. To facilitate this calculation we will need to introduce two new functions. Let $\phi(y,z)$ and ${\tilde{\phi}}(q,u)$ be the solutions of $${\mathcal{L}}_0\phi(y,z,q,u)=\sigma^2(y,z)-\left<\sigma^2(y,z)\right> \quad \text{and} \quad {\mathcal{L}}_0 {\tilde{\phi}}(y,z,q,u)=f(q,u)-\left<f(q,u)\right>,$$ respectively. Now, one can easily deduce that $${\mathcal{L}}_0^{-1}({\mathcal{L}}_2-\left<{\mathcal{L}}_2\right>)C_{0,0}=\phi\cdot\left(\frac{1}{2}x^2 \frac{\partial C_{0,0}}{\partial x^2}+\beta x \frac{\partial C_{0,0}}{\partial x}-\beta C_{0,0}\right)+{\tilde{\phi}}\cdot\left(x \frac{\partial C_{0,0}}{\partial x}-C_{0,0}\right).$$ On the other hand, using (\[eq:av-l2c-10\]) it can be seen that $$\label{eq:l-2-c-1-0}
\begin{split}
&\left<{\mathcal{L}}_2\right>C_{1,0}=
\left(\frac{1}{\sqrt{2}}v\rho_1\langle \sigma \phi_y\rangle\right)x\frac{\partial}{\partial x}(x^2\frac{\partial^2C_{0,0}}{\partial x^2})
+\beta(\sqrt{2}v\rho_1\langle \sigma \phi_y\rangle)x^2\frac{\partial^2 C_{0,0}}{\partial x^2}+(\sqrt{2}v\rho_1\langle \sigma \tilde{\phi}_y\rangle)x^2\frac{\partial^2 C_{0,0}}{\partial x^2}
\\&- v \sqrt{2}\left<\phi_y {\Lambda}\right>\left(\frac{1}{2}x^2 \frac{\partial C_{0,0}}{\partial x^2}+\beta x \frac{\partial C_{0,0}}{\partial x}-\beta C_{0,0}\right)-v \sqrt{2}\left<{\tilde{\phi}}_y {\Lambda}\right>\left(x \frac{\partial C_{0,0}}{\partial x}-C_{0,0}\right)
\\&-{\tilde{v}}\sqrt{2}\left<{\widetilde{\Lambda}}\phi_q\right>\left(\frac{1}{2}x^2 \frac{\partial C_{0,0}}{\partial x^2}+\beta x \frac{\partial C_{0,0}}{\partial x}-\beta C_{0,0}\right)-{\tilde{v}}\sqrt{2}\left<{\widetilde{\Lambda}}{\tilde{\phi}}_q\right>\left(x \frac{\partial C_{0,0}}{\partial x}-C_{0,0}\right)
\end{split}$$ in which $\phi_y$ and ${\tilde{\phi}}_y$ are the derivatives of $\phi$ and ${\tilde{\phi}}$, respectively, with respect to the $y$ variable. The solution of (\[eq:l-2-c-1-0\]) can be explicitly computed as $$C_{1,0}=-(T-t)\text{RHS},$$ in which $\text{RHS}$ is the right-hand-side of (\[eq:l-2-c-1-0\]), using the fact that $\left<{\mathcal{L}}_2\right>$ and $x^n\partial^n/\partial x^n$ commute for $n\in \mathbb{N}$.

Next, we consider (\[eq:C-1-eps\]) together with (\[eq:expansion-eps\]). Equating $1/{\epsilon}$ and $1/\sqrt{{\epsilon}}$ order terms gives $$\label{eq:1-over-epsilon}
{\mathcal{L}}_0 C_{0,1}=0, \quad {\mathcal{L}}_0 C_{1,1}=0.$$ The first equation in (\[eq:1-over-epsilon\]) implies that $C_{0,1}$ is independent of $y$ and $q$, and the second equation can be obtained by observing that ${\mathcal{M}}_3 C_{0,0}={\mathcal{L}}_1 C_{0,1}=0$, which implies that $C_{1,1}$ is independent of $y$ and $q$ as well.

Matching the first order terms gives $$\label{eq:poisson-eqn-2}
{\mathcal{L}}_0 C_{2,1}+ {\mathcal{L}}_2 C_{0,1}=-{\mathcal{M}}_1 C_{0,0},$$ where we have used the fact that ${\mathcal{M}}_3 C_{1,0}={\mathcal{L}}_1 C_{1,1}=0$. Equation (\[eq:poisson-eqn-2\]) is a Poisson equation and its solvability condition gives $$\label{eq:c10}
\left<{\mathcal{L}}_2\right>C_{0,1}=-\left<{\mathcal{M}}_1\right> C_{0,0}, \quad C_{0,1}(T,x,y,z)=0.$$ Using the facts that $$\begin{split}
\frac{\partial}{\partial z}C_{0,0}&=(T-t)\left<\sigma^2(y,z)\right>_{z}\left( \frac{1}{2}x^2 \frac{\partial^2 C_{0,0}}{\partial x^2}+\beta x \frac{\partial}{\partial x}C_{0,0}-\beta C_{0,0}\right),
\\ \frac{\partial}{\partial u}C_{0,0}&=(T-t)\left<f(q,u)\right>_u\left(x \frac{\partial C_{0,0}}{\partial x}-C_{0,0}\right),
\end{split}$$ in which $$\left<\sigma^2(y,z)\right>_{z} \triangleq \frac{\partial}{\partial z}\left<\sigma^2(y,z)\right>, \quad \left<f(q,u)\right>_u= \frac{\partial}{\partial u} \left<f(q,u)\right>_u,$$ and the definition of the operator ${\mathcal{M}}_1$, we can derive the solution to (\[eq:c10\]) as $$\begin{split}
C_{0,1}(t,x,y,z)&=\frac{(T-t)^2}{2}\rho_2 \left<\sigma(y,z)\right> g(z)
\left<\sigma^2(y,z)\right>_{z} \left(x \frac{\partial}{\partial x}\left(x^2  \frac{\partial^2 C_{0,0}}{\partial x^2}\right)+\beta x^2  \frac{\partial^2 C_{0,0}}{\partial x^2}\right)
\\&-\frac{(T-t)^2}{2}g\cdot \left<\Gamma\right> \cdot \left<\sigma^2\right>_z \cdot \left(\frac{1}{2}x^2 \frac{\partial^2 C_{0,0}}{\partial x^2}+\beta x \frac{\partial}{\partial x}C_{0,0}-\beta C_{0,0}\right)
\\&-\frac{(T-t)^2}{2} \left<f\right>_u \cdot {\tilde{g}}\cdot \left<{\widetilde{\Gamma}}\right>\cdot \left(x \frac{\partial C_{0,0}}{\partial x}-C_{0,0}\right)
\end{split}$$ again using the fact that $\left<{\mathcal{L}}_2\right>$ and $x^n\partial^n/\partial x^n$ commute for $n\in \mathbb{N}$.

Our approximation to the price of the call option then becomes $$\label{eq:app-opt-price}
\begin{split}
\widetilde{C^{{\epsilon},\delta}} &\triangleq C_{0,0}-(T-t)\left(V_1^{{\epsilon}} x \frac{\partial}{\partial x}\left(x^2 \frac{\partial C_{0,0}}{\partial x^2}\right)+V_{2}^{{\epsilon}}x^2 \frac{\partial^2 C_{0,0}}{\partial x^2}+V_{3}^{{\epsilon}}\left(x \frac{\partial C_{0,0}}{\partial x}-C_{0,0}\right)\right)
\\&+(T-t)^2\left(V_1^{\delta} x \frac{\partial}{\partial x}\left(x^2 \frac{\partial C_{0,0}}{\partial x^2}\right)+V_{2}^{\delta}x^2 \frac{\partial^2 C_{0,0}}{\partial x^2}+V_{3}^{\delta}\left(x \frac{\partial C_{0,0}}{\partial x}-C_{0,0}\right)\right),
\end{split}$$ in which $$\begin{split}
V_{1}^{{\epsilon}}&=\sqrt{\frac{{\epsilon}}{2}}v \rho_1 \left<\sigma \phi_y\right>,
\\ V_{2}^{{\epsilon}}&=\sqrt{2 {\epsilon}}\left(\beta v \rho_1 \left<\sigma \phi_y\right>+ v \rho_1 \left<\sigma {\tilde{\phi}}_y\right>
-\frac{ v \left<{\Lambda}\phi_y\right>}{2}-\frac{ {\tilde{v}}\left<{\widetilde{\Lambda}}\phi_q\right>}{2}\right),
\\ V_3^{{\epsilon}}&=-\sqrt{2 {\epsilon}}\left(v \left<{\Lambda}\phi_y\right> \beta+ v \left<{\Lambda}{\tilde{\phi}}_y\right>+
{\tilde{v}}\left<{\widetilde{\Lambda}}\phi_q\right> \beta+ {\tilde{v}}\left<{\widetilde{\Lambda}}{\tilde{\phi}}_q\right>\right),
\\ 2V_{1}^{\delta}&= \sqrt{\delta}\rho_2 \left<\sigma\right> g \left<\sigma^2\right>_z,
\\ 2V_{2}^{\delta}&=\sqrt{\delta}\rho_2 \cdot \left<\sigma\right> \cdot g \cdot \left<\sigma^2\right>_z \cdot \beta - \sqrt{\delta}\frac{g \cdot \left<{\Gamma}\right>\cdot \left<\sigma^2\right>_z}{2},
\\2V_3^{\delta}&=-\sqrt{\delta}g \cdot \left<{\Gamma}\right>\cdot \left<\sigma^2\right>_z \cdot \beta  -\sqrt{\delta}\left<f\right>_u \cdot \tilde{g} \cdot  \left<{\widetilde{\Gamma}}\right>.
\end{split}$$ When we calibrate (daily) to the observed option prices (\[eq:app-opt-price\]) will be our starting point, and we will determine the 7 unknowns, $$\label{eq:def-lambda}
\bar{\lambda}(z,u) \triangleq \beta \left<\sigma^2(y,z)\right>+\left<f(q,u)\right>,$$ which appears in $C_{0,0}$, $V_1^{{\epsilon}}$, $V_2^{{\epsilon}}$, $V_3^{{\epsilon}}$, $V_1^{\delta}$, $V_2^{\delta}$ and $V_3^{\delta}$ from the observed option prices. Alternatively, we will obtain $\bar{\lambda}$ from the yield spread of defaultable bonds (either by using the yield spread of the shortest maturity bond directly or by calibrating the model implied bond price to the yield spread curve) and calibrate rest of the parameters to the observed option prices. In either case, the average volatility term in (\[BS-formula\]) will be calculated from the observed stock prices.

Note that when $\beta=0$ and $f=0$ we obtain the approximation formula in [@ronnie-timescale] as expected. Moreover, as in [@ronnie-timescale], the accuracy of our approximation can be shown to be $$|C^{{\epsilon},\delta}-\widetilde{C^{{\epsilon},\delta}}| \leq C
({\epsilon}|\log({\epsilon})|+\delta),$$ for a constant $C>0$. The proof can be first done for smooth pay-off functions by using a higher order approximation for $C^{{\epsilon},\delta}$ and deriving an the expression for the residual using the Feynman-Kac theorem. The proof for the call option and other non-smooth pay-offs can be performed by generalizing the regularization argument in [@sircar-proof]. This is where we use the boundedness assumption on the functions $f$ and $\sigma$.

The approximation for the put option can be similarly derived. Note that the put-call parity relationship holds between the approximate call price and the approximate put price.

**Note:**

In the calibration exercises below, we will compare the performance of \[the seven parameter model\] in (\[eq:app-opt-price\]) with \[the model of [@ronnie-timescale]\], which can be obtained in our framework by setting the intensity $\lambda$ in (\[eq:model\]) to be zero. In this case stochastic volatility is the only source of the implied volatility smile/smirk. We will also consider the performance of two other models: \[The 5 parameter model\] is obtained when we set the market risk functions to be equal to zero, i.e., ${\Lambda}={\widetilde{\Lambda}}={\Gamma}={\widetilde{\Gamma}}=0$. In this case the approximation formula becomes $$\label{eq:comp-mod-1}
\begin{split}
\widetilde{C^{{\epsilon},\delta}} &\triangleq
C_{0,0}-(T-t)\left[{\widetilde{V}}_1^{{\epsilon}}x \frac{\partial}{\partial x} \left(x^2 \frac{\partial^2 C_{0,0}}{\partial x^2}\right)+{\widetilde{V}}_2^{{\epsilon}}x^2 \frac{\partial^2 C_{0,0}}{\partial x^2}\right]
\\&+(T-t)^2\left[{\widetilde{V}}_1^{\delta}x \frac{\partial}{\partial x} \left(x^2 \frac{\partial^2 C_{0,0}}{\partial x^2}\right)+ {\widetilde{V}}_2^{\delta}x^2 \frac{\partial^2 C_{0,0}}{\partial x^2}\right],
\end{split}$$ and the constants $\bar{\lambda}$ (the average default intensity that appears in the formula for $C_{0,0}$), $\beta$, ${\widetilde{V}}_{1}^{{\epsilon}}$, ${\widetilde{V}}_2^{{\epsilon}}$ and ${\widetilde{V}}_1^{\delta}$ are to be calibrated to the observed option prices.

\[The three parameter model\] is obtained when we set ${\Lambda}={\widetilde{\Lambda}}={\Gamma}={\widetilde{\Gamma}}=0$ and also take $\sigma_t=\sigma$, for some positive constant $\sigma$. In this case the approximation formula becomes $$\label{eq:comp-mod-2}
\widetilde{C^{{\epsilon},\delta}} \triangleq
C_{0,0}+\left(-(T-t)V^{{\epsilon}}+(T-t)^2 V^{\delta}\right)x^2
\frac{\partial^2 C_{0,0}}{\partial x^2},$$ and the constants $\bar{\lambda}$ (the average default intensity that appears in the formula for $C_{0,0}$), $V^{{\epsilon}}$, $V^{\delta}$, which are, again, to be calibrated to the observed option prices. In this last model the implied volatility skew is driven only by the default intensity.

Instead of calibrating $\bar{\lambda}$ to the option prices along with the other parameters, one might prefer to obtain it from the yield spread of the corporate bond prices.

Bond Price
----------

Under the zero recovery assumption, the defaultable bond price is given by $$B^{{\epsilon},\delta}(t,Y_t,Z_t,Q_t,U_t;T)=B(t,T)\mathbb{E}\left[\exp\left(-\int_t^{T}\left(\beta \sigma^2(Y_s,Z_s)+f(Q_s,U_s)\right)dt\right)\bigg|Y_t,Z_t,Q_t,U_t\right].$$ Here, we assume the interest rate to be constant because we are interested in bonds with short maturities. The bond price satisfies $$B^{{\epsilon},\delta}(t,x,y,z,q,u)=:B(t,T){\widetilde{B}}^{{\epsilon},\delta}(t,x,y,z,q,u),$$ in which ${\widetilde{B}}^{{\epsilon},\delta}$ satisfies $$\label{eq:bond-pricing-pde}
\begin{split}
{\widetilde{\mathcal{L}}}^{{\epsilon},\delta}{\widetilde{B}}^{{\epsilon},\delta}(t,x,y,z,q,u)&=0, \quad t <T,
\\ {\widetilde{B}}^{{\epsilon},\delta}(T,x,y,z,q,u)&=1,
\end{split}$$ where the operator ${\widetilde{\mathcal{L}}}^{{\epsilon},\delta}$ is given by $${\widetilde{\mathcal{L}}}^{{\epsilon},\delta}\triangleq \frac{1}{{\epsilon}}
{\mathcal{L}}_0+\frac{1}{\sqrt{{\epsilon}}}{\widetilde{\mathcal{L}}}_1+{\widetilde{\mathcal{L}}}_2+\sqrt{\delta}{\widetilde{\mathcal{M}}}_1+\delta {\mathcal{M}}_2+
\sqrt{\frac{\delta}{{\epsilon}}}{\mathcal{M}}_3,$$ in which ${\mathcal{L}}_0$, ${\mathcal{M}}_1$ and ${\mathcal{M}}_3$ are as in (\[eq:operators\]) and $$\begin{split}
{\widetilde{\mathcal{L}}}_1 &\triangleq -{\Lambda}(y,z) v \sqrt{2} \frac{\partial}{\partial y}-{\widetilde{\Lambda}}(q,u){\tilde{v}}\sqrt{2} \frac{\partial}{\partial q},\\
{\widetilde{\mathcal{L}}}_2 &\triangleq \frac{\partial}{\partial t}-(\beta \sigma^2(y,z)+f(q,u))\cdot,\\
{\widetilde{\mathcal{M}}}_1 &\triangleq  -g(z){\Gamma}(y,z)\frac{\partial}{\partial z}-{\tilde{g}}(z){\widetilde{\Gamma}}(q,u)\frac{\partial}{\partial u},
\end{split}$$

Let us obtain an expansion of the price ${\widetilde{B}}^{{\epsilon},\delta}$ in powers of $\sqrt{\delta}$, $$\label{eq:delta-expansion-bond}
{\widetilde{B}}^{{\epsilon},\delta}={\widetilde{B}}_{0}^{{\epsilon}}+\sqrt{\delta}{\widetilde{B}}_1^{{\epsilon}}+\delta
{\widetilde{B}}_2^{{\epsilon}}+\cdots.$$ We will then consider the expansion of each term in (\[eq:delta-expansion-bond\]) in powers of $\sqrt{{\epsilon}}$, i.e., $$\label{eq:expansion-eps-bond}
{\widetilde{B}}_k^{{\epsilon}}={\widetilde{B}}_{0,k}+\sqrt{{\epsilon}}{\widetilde{B}}_{1,k}+{\epsilon}{\widetilde{B}}_{2,k}+ \cdots, \quad
\text{for all $k \in \mathbb{N}$}.$$ The function ${\widetilde{B}}_{0,0}(t,z,u)$ can be determined from $$\left<{\widetilde{\mathcal{L}}}_2\right>{\widetilde{B}}_{0,0}=0, \quad {\widetilde{B}}_{0,0}(T,z,u)=1,$$ as $${\widetilde{B}}_{0,0}(t,z,u)=\exp(-\bar{\lambda}(z,u) (T-t)),$$ in which $\bar{\lambda}$ is as in (\[eq:def-lambda\]). On the other hand, ${\widetilde{B}}_{1,0}$ can be determined using $$\left<{\widetilde{\mathcal{L}}}_2\right>{\widetilde{B}}_{1,0}=\left<{\widetilde{\mathcal{L}}}_1 {\mathcal{L}}_0^{-1}\left({\widetilde{\mathcal{L}}}_2-\left<{\widetilde{\mathcal{L}}}_2\right>\right)\right>{\widetilde{B}}_{0,0}, \quad {\widetilde{B}}_{1,0}(T,z,u)=0,$$ the solution of which is $${\widetilde{B}}_{1,0}(t,z,u)=L(z,u) (T-t) \exp \left(-\bar{\lambda}(T-t)\right),$$ in which $$L(z,u)=v \sqrt{2} \left(\beta \left<\phi_y \Lambda\right>+\left<{\tilde{\phi}}_y {\Lambda}\right>\right)+{\tilde{v}}\sqrt{2} \left(\beta \left<\phi_q {\widetilde{\Lambda}}\right>+\left<{\tilde{\phi}}_q {\widetilde{\Lambda}}\right>\right).$$ The function ${\widetilde{B}}_{0,1}$ can be determined from $$\left<{\mathcal{L}}_2\right>{\widetilde{B}}_{0,1}=-\left<{\widetilde{\mathcal{M}}}_1\right> {\widetilde{B}}_{0,0}, \quad {\widetilde{B}}_{0,1}(T,x,y,z)=0.$$ The solution of this equation can be written as $${\widetilde{B}}_{0,1}(t,z,u)=-\tilde{L}(z,u)\frac{(T-t)^2}{2}\exp(-\bar{\lambda}(z,u)(T-t)),$$ in which $$\tilde{L}(z,u)=g \left<{\Gamma}\right>\left(\bar{\lambda}\right)_z+{\tilde{g}}\left<{\widetilde{\Gamma}}\right>\left(\bar{\lambda}\right)_q.$$

It can be shown that $$\widehat{B}^{{\epsilon},\delta}(t,z,u;T)=B(t,T)\exp\left(\bar{\lambda}\cdot (T-t)\right)\left(1+L\cdot(T-t)-\tilde{L}\cdot\frac{(T-t)^2}{2}\right),$$ satisfies $|\hat{B}^{{\epsilon},\delta}-{\widetilde{B}}^{{\epsilon},\delta}| \leq C
({\epsilon}+\delta)$, for some constant $C \geq 0$, as in [@papa]. Since very short term maturity options are not available, we will not use this formula to calibrate $\bar{\lambda}$ from the yield spread data. Instead, in Section 3.1, we will take $\bar{\lambda}$ to be equal to the yield spread of the shortest maturity bond.

The Approximate Implied Volatility {#sec:calibration}
----------------------------------

The implied volatility of the approximate pricing formula (\[eq:app-opt-price\]), which we will denote by $I$, is implicitly defined as $$\label{eq:implied-vol}
 C_{\text{BS}}(x; I,r;K,T-t)=C^{{\epsilon},\delta}(t,x,z).$$ In order to understand how a typical implied volatility surface of our model looks like, we expand the implied volatility $I$ in terms of ${\epsilon}$ and $\delta$ as $$I=I_0+\sqrt{{\epsilon}}I^{{\epsilon}}_1+\sqrt{\delta}I^{\delta}_1+\cdots$$ We then apply the Taylor expansion formula to the left-hand-side of (\[eq:implied-vol\]) and use the price formula (\[eq:app-opt-price\]) on the right-hand-side and obtain $$\label{eq:imp-vol-match}
\begin{split}
 &C_{\text{BS}}(x;I_0,r;K,T-t)+(\sqrt{{\epsilon}}I^{{\epsilon}}_1+\sqrt{\delta}I^{\delta}_1)
 \frac{\partial C_{\text{BS}}}{\partial \sigma
 }(x;I_0,r;K,T-t)+\cdots\\
 &=C_{0,0}-(T-t)\left(V_1^{{\epsilon}} x \frac{\partial}{\partial x}\left(x^2 \frac{\partial C_{0,0}}{\partial x^2}\right)+V_{2}^{{\epsilon}}x^2 \frac{\partial^2 C_{0,0}}{\partial x^2}+V_{3}^{{\epsilon}}\left(x \frac{\partial C_{0,0}}{\partial x}-C_{0,0}\right)\right)
\\&+(T-t)^2\left(V_1^{\delta} x \frac{\partial}{\partial x}\left(x^2 \frac{\partial C_{0,0}}{\partial x^2}\right)+V_{2}^{\delta}x^2 \frac{\partial^2 C_{0,0}}{\partial x^2}+V_{3}^{\delta}\left(x \frac{\partial C_{0,0}}{\partial x}-C_{0,0}\right)\right)+\cdots.
\end{split}$$ Matching the order one terms, we get $I_0$ as the solution of $$\label{eq:defn-I-0}
C_{\text{BS}}(x;I_0,0;K B(t,T),T-t)=C_{\text{BS}}\left(x;\sqrt{\left<\sigma^2\right>},\bar{\lambda};K B(t,T),T-t\right).$$ Note that $I_0$ is already a non-linear function of the log moneyness, i.e. $\log(K/x)$, and of the time to maturity, i.e. $T-t$. In fact, $$x\sqrt{T-t}N'(d_1)\frac{\partial I_0}{\partial
K}=B(t,T)(N(d_2)-N(\widetilde{d}_2))<0,$$ i.e. $I_0$ is a strictly decreasing function of $K$. Here, $$d_{1,2} \triangleq \frac{\log(x/[K B(t,T)]) \pm \frac{1}{2}I_0^2\cdot(T-t)}{I_0
\sqrt{T-t}}.$$ The first term in the implied volatility expansion, $I_0$, already carries the features of a typical smile curve and captures the fact that in the money call contracts or out of the money put contracts have larger implied volatilities.

To match a given implied volatility curve, one will need to use $\sqrt{\delta}$ and $\sqrt{{\epsilon}}$ terms in the expansion of the implied volatility. Each of these terms will have one free parameter that can be used to calibrate the model to a given implied volatility curve. Matching the $\sqrt{{\epsilon}}$ and $\sqrt{\delta}$ terms in (\[eq:imp-vol-match\]), we obtain $$\label{eq:implied-vol-components}
\begin{split}
\sqrt{{\epsilon}}I_1^{{\epsilon}}&=-(T-t)\left(V_1^{{\epsilon}} x \frac{\partial}{\partial x}\left(x^2 \frac{\partial C_{0,0}}{\partial x^2}\right)+V_{2}^{{\epsilon}}x^2 \frac{\partial^2 C_{0,0}}{\partial x^2}+V_{3}^{{\epsilon}}\left(x \frac{\partial C_{0,0}}{\partial x}-C_{0,0}\right)\right)\left( \frac{\partial
C_{\text{BS}}}{\partial \sigma
 }\right)^{-1}  \quad \text{and} \quad
\\ \sqrt{\delta}I_1^{\delta}&= (T-t)^2\left(V_1^{\delta} x \frac{\partial}{\partial x}\left(x^2 \frac{\partial C_{0,0}}{\partial x^2}\right)+V_{2}^{\delta}x^2 \frac{\partial^2 C_{0,0}}{\partial x^2}+V_{3}^{\delta}\left(x \frac{\partial C_{0,0}}{\partial x}-C_{0,0}\right)\right)\left( \frac{\partial
C_{\text{BS}}}{\partial \sigma
 }\right)^{-1}.
\end{split}$$ To derive these terms more explicitly we will use the derivatives of the Black-Scholes function $C_{BS}$ on the left hand side of equation (\[eq:defn-I-0\]): $$\begin{gathered}
\frac{\partial C_{BS}}{\partial \sigma}=(T-t)I_0 x^2\frac{\partial^2C_{BS}}{\partial x^2}\\
\left(x\frac{\partial}{\partial x}\right)\frac{\partial}{\partial
\sigma}C_{BS}=\left(1-\frac{d_1}{I_0 \sqrt{T-t}}\right)\frac{\partial}{\partial
\sigma}C_{BS}\quad \text{and,}\\
\frac{\partial C_{BS}}{\partial \sigma}=\frac{x\exp
(-d_1^2/2)\sqrt{T-t}}{\sqrt{2\pi} }.\label{eq:vega}\end{gathered}$$ The corresponding derivatives of $C_{0,0}$ can be obtained by replacing $d_1$ by $\tilde{d}_1$ and $I_0$ by $\sqrt{\left<\sigma^2\right>}$. Now, we can obtain the correction terms in the implied volatilities as $$\begin{split}
\sqrt{{\epsilon}}I_1^{{\epsilon}} &= -V_{1}^{{\epsilon}}\exp\left(\frac{-\tilde{d_1}^2+d_1^2}{2}\right)\left(1-\frac{\tilde{d_1}}{\sqrt{\left<\sigma^2\right>\cdot(T-t)}}\right)
-V_{2}^{{\epsilon}}\frac{1}{\left<\sigma^2\right>}\exp\left(\frac{-\tilde{d_1}^2+d_1^2}{2}\right)
\\&-V_{3}^{{\epsilon}}\sqrt{2 \pi \cdot (T-t)}\exp\left(\frac{d_1^2}{2}\right)\left(\frac{K}{x}B(t,T)N(\tilde{d}_2)\right),
\\ \sqrt{{\epsilon}}I_1^{\delta} &= (T-t)V_{1}^{\delta}\exp\left(\frac{-\tilde{d_1}^2+d_1^2}{2}\right)
\left(1-\frac{\tilde{d_1}}{\sqrt{\left<\sigma^2\right>\cdot(T-t)}}\right)
+(T-t)V_{2}^{\delta}\frac{1}{\left<\sigma^2\right>}\exp\left(\frac{-\tilde{d_1}^2+d_1^2}{2}\right)
\\&+(T-t)V_{3}^{\delta}\sqrt{2 \pi \cdot (T-t)}\exp\left(\frac{d_1^2}{2}\right)\left(\frac{K}{x}B(t,T)N(\tilde{d}_2)\right).
\end{split}$$

Typical implied volatility surfaces corresponding to the models/approximate prices in (\[eq:comp-mod-2\]); (\[eq:comp-mod-1\]); and the model of [@ronnie-timescale], which is the same as setting $\bar{\lambda}=0$ in (\[eq:comp-mod-1\]); and (\[eq:app-opt-price\]) are given respectively by Figures \[fig:s3p\], \[fig:s5p\], \[fig:sfp\] and \[fig:sfp\]. Although all of the models are able to produce the well-known properties of the implied volatility surface of the observed option prices (a-the skewedness (out of the money put options have more implied volatility then in the money ones) and b- decrease in the degree of skewedness as the time to maturity grows, c-smile (the minimum of the implied volatility is at the forward stock price $S_t B(t,T)$)) to some degree, the models in (\[eq:comp-mod-1\]) and (\[eq:app-opt-price\]) are able to produce surfaces with richer structures that looks more like the real implied volatility surfaces.

Calibration to Data
===================

Data Description
----------------

We calibrated (daily) our model’s implied volatility Ford’s options over the period of 2/2/2005-1/9/2006. The sources of data are:

-   The stock price data is obtained from finance.yahoo.com.

-   The stock option data is from OptionMetrics under WRDS database, which is the same database [@carr-wu] use. It is contained in their Volatility Surface file, which contains the interpolated volatility surface using a methodology based on kernel smoothing algorithm. The file contains information on standardized options, with expirations 91, 122, 152, 182, 273, 365, 547 and 730 calender days and different strikes. There are 130 different call (put) options on each day. Exchange traded options on individual stocks are American-style and therefore the price reflects an early exercise premium. Option metrics uses binomial tree to back out the option implied volatility that explicitly accounts for this early exercise premium. We adopt the practice employed by [@carr-wu] for estimating the market prices of European options: We take the average of two implied volatilities at each strike and convert them into out-of-the-money European option prices using the Black-Scholes formula.

-   The treasury yield curve data for each day is also obtained from OptionMetrics. It is derived from the LIBOR rates and settlement prices of CME Eurodollar futures. (Recall that we need the bond price $B(t,T)$ for short time to maturities $T-t$.)

-   Ford’s bond data is obtained from a Bloomberg terminal.

We will consider the performance of the three models given in (\[eq:app-opt-price\]), (\[eq:comp-mod-1\]) and (\[eq:comp-mod-2\]), in fitting the implied volatility curve of the observed option data. In the price formula (\[eq:app-opt-price\]), we have seven parameters to calibrate: $\bar{\lambda}$, $V_i^{{\epsilon}}$, $V^{\delta}_i$, $i \in \{1,2,3\}$; on the other hand (\[eq:comp-mod-1\]) we have five parameters to calibrate: $\bar{\lambda}$, $\tilde{V}^{{\epsilon}}_i$, $\tilde{V}^{\delta}_i$, $i \in \{1,2\};$ and finally, in (\[eq:comp-mod-2\]) we have 3 parameters to calibrate: $\bar{\lambda}$, $V^{{\epsilon}}$, $V^{\delta}$. We will also consider the approximation formula when when we set $\bar{\lambda}=0$ (this is the approximation is the same as the one in [@ronnie-timescale]), and test the performance of the other models against it. The average volatility term $\left<\sigma^2(y,z)\right>$ term that appear in our approximation formula will be estimated using the stock price data as in [@ronnie-timescale].

Following the suggestion of [@cont] (see page 439), we perform the calibration on a particular day by minimizing $$\label{eq:approx-impl}
\begin{split}
\sum_i &\left(I^{\text{observed}}(T_i,K_i)-I^{\text{model}}(T_i,K_i; \text{model parameters})\right)^2
\\ &\hspace{0.5in} \approx\sum_i \frac{\left(O^{\text{observed}}(T_i,K_i)-O^{\text{model}}(T_i,K_i; \text{model parameters})\right)^2}{\text{vega}^2(T_i,K_i)},
\end{split}$$ in which $O^{\text{observed}}(T_i,K_i)$ is the market price of the out of the money option (put or call) with strike price $K_i$, and maturity $T_i$, where as $O^{\text{model}}(T_i,K_i; \text{model parameters})$ is the model implied price of the same option contract. The term $\text{vega}(T_i,K_i)$ is computed using (\[eq:vega\]) with $I_0(T_i,K_i)$ as the market implied volatility.

We will group our calibration scheme into two depending on whether the average hazard rate $\bar{\lambda}$, defined in (\[eq:def-lambda\]) (which enters into the approximate price formulas through $C_{0,0}$ in (\[BS-formula\]) is inferred from the yield spread of the bond prices or it is inferred as a result of calibration to the market option prices. These two groups of calibration schemes will be the topic of the following two subsections.

When the Hazard Rate $\bar{\lambda}$ is Inferred From the Bond Yield Spread
---------------------------------------------------------------------------

We fix the average hazard rate $\bar{\lambda}$ that appears in our price approximation formulas (\[eq:app-opt-price\]), (\[eq:comp-mod-1\]) and (\[eq:comp-mod-2\]), to be the yield spread of the shortest maturity corporate bond (of the same company), which is 2 years in our calibration exercise. We pool the options of maturities $91, 122, 152, 182, 273, 365, 547, 730$ together and use (\[eq:approx-impl\]) to calibrate our model. Each maturity has 13 strike prices available, therefore our data set 104 observed option prices available. After we calibrate our models and the model of [@ronnie-timescale], which is obtained by setting $\bar{\lambda}=0$ in (\[eq:comp-mod-1\]), we plot their implied volatilities against the implied volatilities of the data. These are illustrated in Figures \[fig:9m\], \[fig:1y\], \[fig:1.5y\] and \[fig:2y\]. Each figure is for a different maturity.

As expected the 7 parameter model of (\[eq:app-opt-price\]) fits the data implied volatility better. (Observe that since we fixed $\bar{\lambda}$ to be the yield spread of the shortest maturity bond there are only 6 parameters to calibrate). But the 7 parameter model does not only outperform the other models but also fits the data implied volatility almost perfectly. 7 parameter model outperforms the other models by a lot for very small and very large strikes. Another interesting observation is that the 5 parameter model of (\[eq:comp-mod-1\]) always outperforms (again since the $\bar{\lambda}$ is set to be the yield spread of the shortest maturity bond) the model of [@ronnie-timescale] which is the same as (\[eq:comp-mod-1\]) except that the $\bar{\lambda}$ is set to zero instead of being set to the yield spread of the shortest maturity bond. This shows the significance of taking the yield spread of the corporate bonds into account in pricing options. Note that although the model of (\[eq:comp-mod-2\]) had only two free parameters, it performed fairly well (considering that we had 104 data points available). Also note that, the performance of the 7 parameter model deteriorates as the time to maturity becomes smaller. This is because for the shorter maturity options one should use shorter term yield spreads.

When the Hazard Rate $\bar{\lambda}$ is Inferred from Market Option Prices Directly
-----------------------------------------------------------------------------------

In this section, we will analyze whether it would be possible to infer the risk neural expected time of default if we calibrate the models of (\[eq:app-opt-price\]), (\[eq:comp-mod-1\]) and (\[eq:comp-mod-2\]) to the data implied volatility. That is along with the other parameters, if we also let $\bar{\lambda}$ be decided by the prices of the options would it be possible to say something about the yield spread of the corporate bonds.

Since each model has one extra free variable, namely $\bar{\lambda}$, to fit to the data implied volatility, not surprisingly, the performance of each model enhances, see Figures \[fig:9m-iter\], \[fig:1y-iter\], \[fig:1.5y-iter\] and \[fig:2y-iter\]. But the enhancement is not an ordinary one: all of our three models fit the data exceptionally well, even the three parameter model in (\[eq:comp-mod-2\]). The improvement in the most advanced model is the smallest. But then when we compare the data implied yield spread versus the yield spread of the shortest maturity bond, see Figures \[fig:intensity-3p\], \[fig:intensity-5p\] and \[fig:intensity-7p\], the story is a little bit different. The implied yield spread, $\bar{\lambda}$ of the 3 parameter model of (\[eq:comp-mod-2\]) is much larger than that of the yield spread of the shortest maturity corporate bond. Since, the 3 parameter model does not account for volatility, the implied yield spread is drastically severe. For the 5 parameter model the difference seems to be less severe, since the model accounts for the stochastic volatility effects, but the difference between the yield spread curves is still noticeable. On the other hand, the implied yield spread for the 7 parameter model of (\[eq:app-opt-price\]) seems to fit the yield spread of the corporate bond almost perfectly. The difference between the 7 parameter and the 5 parameter values is that the former one accounts for the market price of intensity risk.

[xx]{}

Ayache, E., Forsyth, P.  Vetzal, K. 2003. The valuation of convertible bonds with credit risk, [*Journal of Derivatives*]{} [**11 (Fall)**]{}: 9–29.

Brémaud, P. 1981. , Springer, New York.

Carr, P.  Linetsky, V. 2006. A jump to default extended cev model: An application of bessel processes, [*Finance and Stochastics*]{} [**10**]{}: 303–330.

Carr, P.  Wu, L. 2006. Stock options and credit default swaps: A joint framework for valuation and estimation, [*NYU preprint, available at [http://faculty.baruch.cuny.edu/lwu]{}*]{} .

Cont, R.  Tankov, P. 2004. , Chapman & Hall, Boca Raton, FL.

Fouque, J.-P., Papanicolaou, G.  Sircar, K. R. 2000. , Cambridge University Press, New York.

Fouque, J. P., Papanicolaou, G., Sircar, R.  Solna, K. 2003a. Multiscale stochastic volatility asymptotics, [*SIAM J. Multiscale Modeling and Simulation*]{} [**2 (1)**]{}: 22–42.

Fouque, J.-P., Papanicolaou, G., Sircar, R.  Solna, K. 2003b. Singular perturbations in option pricing, [*SIAM Journal of Applied Mathematics*]{} [**63 (5)**]{}: 1648–1665.

Hull, J., Nelken, I.  White, A. 2004. Merton’s model, credit risk, and volatility skews, [*Journal of Credit Risk*]{} [**1 (1)**]{}: 1–27.

Linetsky, V. 2006. Pricing equity derivatives subject to bankruptcy, [*Mathematical Finance*]{} [**16 (2)**]{}: 255–282.

Merton, R. 1976. Option pricing when underlying stock returs are discontinuous, [ *Journal of Financial Economics*]{} [**3**]{}: 125–144.

Papageorgiou, E.  Sircar, R. 2006. Multiscale intensity based models for single name credit derivatives, [*to appear in Applied Mathematical Finance*]{}, available at http://www.princeton.edu/ sircar/

Schönbucher, P. J. 2003. , Wiley, New York.

[^1]: Department of Mathematics, University of Michigan, Ann Arbor, MI 48109, USA, email: erhan@umich.edu

[^2]: This work is supported in part by the National Science Foundation, under grant DMS-0604491.

[^3]: This work was inspired by the SAMSI workshops on Financial Mathematics, Statistics and Econometrics (Fall 2005, Spring 2006 North Carolina). The author wishes to thank the organizers for the travel grant to participate in this stimulating event. I also would like to thank Bo Yang for his research assistance and the two anonymous referees and an anonymous associate editor for their valuable suggestions.
