Proof behind solution for theta in Hull-White with time-dependent volatility and mean reversion?

Question

I'm studying the following paper on Hull-White model calibration: Hull-White paper

In this paper they study the general form of the HW model with time-dependent mean reversion and volatility: $$dr(t) = (\theta(t) - a(t)r(t))dt + \sigma(t)dW(t)$$ I am trying to prove their formula for the solution $\theta(t)$ which matches a prescribed discount rate at each maturity (on p. 36 of the paper). I have not seen this function derived under such generality.

The solution in the paper is as follows: $$\theta(t) = \frac{\partial}{\partial t} f(0,t) + a(t) f(0,t) + \frac{1}{2} \left(\frac{\partial^2}{\partial t^2} V(0,t) + a(t) \frac{\partial}{\partial t} V(0,t) \right)$$ Where: $$f(0,t) = \text{ instantaneous forward rate at time } t$$ $$V(t,T) = \int_t^T \sigma(u,T)^2 du$$ $$\sigma(u,T) = \sigma(u)B(u,T)$$ $$B(t,T) = E(t) \int_t^T \frac{du}{E(u)}$$ $$E(t) = e^{\int_0^t a(u) du}$$ I imagine one can follow a similar method done in simpler versions of the model and:

1. Compute $r(t)$ via integration

2. Compute the bond prices

3. Take derivative of log of the bond price to compute the forward price

4. Somehow isolate $\theta(t)$ after possibly taking another derivative of the previous equation?

---

Update: Here is what I got so far. Apply Itô's Lemma to the function: $$F(r(t), t) = E(t) r(t)$$ We get: $$F_t = a(t) E(t) r(t) \hspace{10pt} F_r = E(t) \hspace{10pt} F_{rr} = 0$$ Resulting in: $$\begin{align*} dF(t) &= F_t dt + F_r dr \\ &= E(t)\theta(t)dt + E(t)\sigma(t)dW(t) \end{align*}$$ Integrating from $t$ to $s$ one gets: $$\begin{align*} F(s) - F(t) &= \int_t^s E(u)\theta(u)du + \int_t^s E(u)\sigma(u)dW(u) \\ E(s)r(s) &= E(t)r(t) + \int_t^s E(u)\theta(u)du + \int_t^s E(u)\sigma(u)dW(u) \\ r(s) &= \frac{E(t)}{E(s)}r(t) + \frac{1}{E(s)}\int_t^s E(u)\theta(u)du + \frac{1}{E(s)}\int_t^s E(u)\sigma(u)dW(u) \end{align*}$$ We now compute the bond prices by integrating this short rate from $t$ to $T$ and applying the definition of the bond price. First we integrate and make use of Fubini's theorem: $$\begin{align*} \int_t^T r(s) ds &= \int_t^T \frac{E(t)}{E(s)}r(t) ds + \int_t^T \int_t^s \frac{E(u)}{E(s)}\theta(u) du ds + \int_t^T \int_t^s \frac{E(u)}{E(s)}\sigma(u) dW(u) ds \\ &= \int_t^T \frac{E(t)}{E(s)}r(t) ds + \int_t^T \int_u^T \frac{E(u)}{E(s)}\theta(u) ds du + \int_t^T \int_u^T \frac{E(u)}{E(s)}\sigma(u) ds dW(u) \\ &= r(t)\underbrace{E(t) \int_t^T \frac{1}{E(s)} ds}_{B(t,T)} + \int_t^T \theta(u) \underbrace{E(u) \int_u^T \frac{1}{E(s)} ds}_{B(u,T)} du + \int_t^T \sigma(u) \underbrace{E(u) \int_u^T \frac{1}{E(s)} ds}_{B(u,T)} dW(u) \\ &= r(t) B(t,T) + \int_t^T \theta(u) B(u,T) du + \int_t^T \sigma(u) B(u,T) dW(u) \end{align*}$$ Now we can compute the bond prices via: $$\begin{align*} P(t,T) &= \mathbb{E}_t \left[e^{-\int_t^T r(s) ds} \right] \\ &= e^{-r(t)B(t,T) - \int_t^T\theta(u)B(u,T)du} \mathbb{E}_t \left[ e^{-\int_t^T \sigma(u) B(u,T) dW(u)} \right] \\ &= e^{-r(t)B(t,T) - \int_t^T\theta(u)B(u,T)du} e^{\frac{1}{2} \int_t^T \sigma(u)^2 B(u,T)^2 du} \end{align*}$$ Here we used that $r(t)$ is known at $t$ and the formula for the expectation for a lognormal variable. We can now take a log and derivative to produce a formula for the instantaneous forward rate: $$\begin{align*} f(0,T) &= - \frac{\partial}{\partial T} \ln P(0,T) \\ &= \frac{\partial}{\partial T} \left[ r(0) B(0,T) + \int_0^T \theta(u) B(u,T) du - \frac{1}{2} \int_0^T \sigma(u)^2 B(u,T)^2 du \right] \end{align*}$$ Here we will use a couple properties of $B(t,T)$ which follow directly from its definition: $$\frac{\partial}{\partial T} B(t,T) = \frac{E(t)}{E(T)} \hspace{5pt} \text{ and } \hspace{5pt} B(T,T) = 0$$ $$\begin{align*} f(0,T) &= r(0) \frac{E(0)}{E(T)} + \int_0^T \theta(u) \frac{E(u)}{E(T)} du - \frac{1}{2} \int_0^T \sigma(u)^2 2B(u,T)\frac{E(u)}{E(T)} du \end{align*}$$ Taking one more derivative with respect to $T$ will result in two instances of $\theta(t)$ in the formula. I don't see a way to isolate this function by itself to recover the formula cited in the paper above.

Answer 1

We assume that the process $\{r(t), \, t \ge 0\}$ satisfies an SDE of the form \begin{align*} dr(t) = \big( \theta(t) - a(t) r(t) \big)dt + \sigma(t) dW_t, \quad t > 0, \end{align*} where $\{W_t, \, t \ge 0 \}$ is a standard Brwonian motion. Note that \begin{align*} d\left(e^{\int_0^t a(u) du}r(t) \right) &=a(t) e^{\int_0^t a(u) du}r(t) dt + e^{\int_0^t a(u) du} dr(t)\\ &= \theta(t) e^{\int_0^t a(u) du} dt + \sigma(t)e^{\int_0^t a(u) du}dW_t. \end{align*} Then, for $v \ge t$, \begin{align*} r(v) = r(t) e^{-\int_t^v a(u) du} + \int_t^v \theta(s)e^{-\int_s^v a(u)du} ds+\int_t^v \sigma(s)e^{-\int_s^v a(u)du} dW_s. \end{align*} For $0 \le t \le T$, let \begin{align*} P(t, T) = E\left(e^{-\int_t^T r_s ds}\,|\, \mathcal{F}_t \right), \end{align*} where $E$ is the expectation operator under the risk-neutral measure and $\mathcal{F}_t$ is the information set up to time $t$, be the value of a zero-coupon bond with maturity $T$ and unit notional amount. Moreover, let \begin{align*} B(t, T) &= \int_t^T e^{-\int_t^v a(u) du} dv,\\ \sigma(t, T) &= \sigma(t) B(t, T),\\ V(t, T) &= \int_t^T \sigma^2(u, T) du. \end{align*} From the above \begin{align*} P(t, T) &= E\left(e^{-\int_t^T r_v dv}\,|\, \mathcal{F}_t \right) \\ &=E\left(e^{-\int_t^T \left(r(t) e^{-\int_t^v a(u) du} + \int_t^v \theta(s)e^{-\int_s^v a(u)du} ds\right)dv - \int_t^T \left(\int_t^v \sigma(s)e^{-\int_s^v a(u)du} dW_s\right) dv}\,|\, \mathcal{F}_t \right) \\ &=e^{-r(t) B(t, T)}E\left(e^{-\int_t^T \left(\int_t^v \theta(s)e^{-\int_s^v a(u)du} ds\right)dv - \int_t^T \left(\int_t^v \sigma(s)e^{-\int_s^v a(u)du} dW_s\right) dv}\,|\, \mathcal{F}_t \right) \\ &=e^{-r(t) B(t, T)}E\left(e^{-\int_t^T \left(\int_s^T \theta(s)e^{-\int_s^v a(u)du} dv\right)ds - \int_t^T \left(\int_s^T \sigma(s)e^{-\int_s^v a(u)du} dv\right) dW_s}\,|\, \mathcal{F}_t \right) \\ &=e^{-r(t) B(t, T)}E\left(e^{ -\int_t^T \theta(s)B(s, T)ds - \int_t^T \sigma(s)B(s, T) dW_s}\,|\, \mathcal{F}_t \right) \\ &=e^{-r(t) B(t, T)-\int_t^T \theta(s)B(s, T)ds + \frac{1}{2}V(t, T)}. \end{align*} Then, \begin{align*} f(0, t)&=-\frac{\partial}{\partial t}\ln P(0, t)\\ &=r(0)\frac{\partial}{\partial t}B(0, t) + \theta(t)B(t, t) + \int_0^t\theta(s)\frac{\partial}{\partial t}B(s, t)ds - \frac{1}{2}\frac{\partial}{\partial t}V(0, t)\\ &=r(0)e^{-\int_0^t a(u) du} + \int_0^t\theta(s)e^{-\int_s^t a(u) du}ds- \frac{1}{2}\frac{\partial}{\partial t}V(0, t). \end{align*} and \begin{align*} \frac{\partial}{\partial t} f(0, t) &= -a(t)r(0)e^{-\int_0^t a(u) du} + \theta(t) - a(t)\int_0^t \theta(s)e^{-\int_s^t a(u) du}ds-\frac{1}{2}\frac{\partial^2}{\partial t^2}V(0, t)\\ &=-a(t)\left(f(0, t) + \frac{1}{2}\frac{\partial}{\partial t}V(0, t) \right) + \theta(t) -\frac{1}{2}\frac{\partial^2}{\partial t^2}V(0, t). \end{align*} That is, \begin{align*} \theta(t) = \frac{\partial}{\partial t} f(0,t) + a(t) f(0,t) + \frac{1}{2} \left(\frac{\partial^2}{\partial t^2} V(0,t) + a(t) \frac{\partial}{\partial t} V(0,t) \right). \end{align*}