Deriving the Heston-Hull-White PDE

I'm trying to derive the Heston-Hull-White PDE. The correct backwards PDE is equation (1.3) of this paper on page (2). I will begin deriving the forward PDE, but switching between the two is trivial.

The model I am working with is the Heston-Hull-White model, given below:

$$\mathrm{d}S = rS\mathrm{d}t + \sqrt{v}S\mathrm{d}W_1$$ $$\mathrm{d}v = \kappa (\bar{v}-v)\mathrm{d}t+\omega\sqrt{v}\mathrm{d}W_2$$ $$\mathrm{d}r = \lambda(\theta(t)-r)\mathrm{d}t + \eta\mathrm{d}W_3\text{.}$$

I assume that $$S$$ and $$v$$ have correlation $$\rho_{S,v}$$, $$S$$ and $$r$$ have correlation $$\rho_{S,r}$$, and $$v$$ and $$r$$ have correlation $$\rho_{v,r}$$.

My approach is to apply the Feynman-Kac theorem. This is standard, so I will skip most steps.

1. Let $$h(S(T))$$ be the payoff function of the option. For a vanilla call, $$h(S(T))=S(T)-K\text{.}\tag{1}$$

2. Let $$g(t, S(t), v(t), r(t)) =\tilde{\mathbb{E}}\Big(e^{-\int_{u=t}^{u=T}r(u)\mathrm{d}u}h(X(T))\Big)\tag{2}$$ be the price of the option. My goal is to find the PDE (implied by Feynman-Kac) for $$g$$.

3. $$g$$ is not a martingale, so we want to first make a transformation to get a martingale. If we follow a Black-Scholes example, we may try to do something like this: $$f(t, S(t), v(t), r(t)) = e^{-\int_{u=0}^{u=t}r(u)\mathrm{d}u}g(t, S(t), v(t), r(t))\text{.}\tag{3}$$ so that $$f$$ becomes a martingale. The problem here is that $$r$$ is a random variable, and we can't pull it outside of the expected value. I continue now as if the definition of $$f$$ makes sense.

4. I now apply Itô's lemma and set the coefficient of $$\mathrm{d}t$$ equal to $$0$$. That gives the following PDE for $$f$$: $$f_t + rSf_S +\kappa(\bar{v}-v)f_v+\lambda(\theta(t)-r)f_r + \rho_{S,v}Sv\omega f_{s,v} + \rho_{S, r}\eta S\sqrt{v}f_{S, r} + \rho{v, r}\omega \sqrt{v}\eta f_{v, r}+\frac{1}{2}vS^2f_{S, S} + \frac{1}{2}v\omega^2f_{v, v} + \frac{1}{2}f_{r, r}\eta^{2}=0\text{.}\tag{4}$$

The PDE in (4) is the PDE for $$f$$, but I need the PDE for $$g$$. Following a Black-Scholes example, I get the PDE of $$g$$ by computing the partial derivatives of $$f$$ in terms of $$g$$ using (3).

For example, I compute

$$f_t = e^{-\int_{u=0}^{u=t}r(u)\mathrm{d}u}(-r(t) g + g_t)\text{.}\tag{5}$$

I then substitute this back into equation (4).

Now I need to compute $$f_r$$, but again, $$r$$ is a random variable. Step (3) was indeed a mistake. How do I continue?

(Just what I think is the right start)

The pricing PDE comes out of the dynamics of a self-financing portfolio, $$\Pi$$, hedged against the movements of stock, $$S$$, its volatility, $$v$$, and interest rate $$r$$.

With $$V$$ the target path-independent derivative, $$U$$ vanilla European option (different from $$V$$), and $$P$$ zero-coupon bond, the portfolio would be:

$$\Pi = V+\alpha S + \beta U + \gamma P,$$ $$d\Pi = dV+\alpha dS + \beta dU + \gamma dP,$$

and we would exploit the fact that $$\Pi$$ is riskless (set the integrands of $$dS$$, $$dv$$ and $$dr$$ to $$0$$) and also the fact that it returns $$r$$ :

$$d\Pi = r \Pi dt.$$

We note that, under HW dynamics, $$P(t,T)$$ is just a deterministic function of $$r_t$$ only (so we don't need to worry about the integrated short rate) and it has its own PDE (in $$r$$ variable).

This would be an extension of the detailed derivation of Heston PDE by Rouah here.