The proof is relatively long, so I focus on displaying the reasoning and major steps.
We work on a Black-Scholes model. Without loss of generality, we focus on an option with strike $P$ to buy at $t_e$ a European call option expiring at $T$, written on a stock $S$. Expectations are always taken with respect to the risk-neutral measure $Q$ unless otherwise stated, and we write $E_t(\cdot):=E(\cdot|\mathscr{F}_t)$.
The value $C_t$ of a European call option is equal to: $$C_t=E_t\left(e^{-r(T-t)}\max\{S_T-K\}\right)$$
The value $O_t$ of an option over a European call option is given by: $$\begin{align} O_t &=E_t\left(e^{-r(t_e-t)}\max\{C_{t_e}-P,0\}\right) \\&=e^{-r(t_e-t)}E_t\left( \max\left\{E_{t_e}\left(e^{-r(T-t_e)}\max\{S_T-K,0\}\right)-P,0\right\}\right) \\&=e^{-r(T-t)}E_t \left(1_{S_{t_e}\geq S^\star}\left(E_{t_e}\left(1_{S_T\geq K}\left(S_T-K\right)\right)-\tilde{P}\right)\right) \\&=e^{-r(T-t)}E_t\left(E_{t_e} \left(1_{S_{t_e}\geq S^\star}1_{S_T\geq K}\left(S_T-K\right)-1_{S_{t_e}\geq S^\star}\tilde{P}\right)\right) \\\tag{1}&=e^{-r(T-t)}\left( E_t\left(1_{S_{t_e}\geq S^\star, S_T\geq K}S_T\right) -E_t\left(1_{S_{t_e}\geq S^\star, S_T\geq K}\right)K -E_t\left(1_{S_{t_e}\geq S^\star}\right)\tilde{P} \right) \end{align}$$ where the last inequality stems from the law of iterated expectations, $\tilde{P}:=e^{r(T-t_e)}P$ is the compounded strike, and $S^\star$ is the value of $s$ that solves the following equation: $$\tag{2}c(s,T-t_e)-\tilde{P}=0$$ where $c$ is the undiscounted Black-Scholes price for a European call option: $$c(s,\tau):=se^{r\tau}\Phi\left(\frac{\ln\frac{s}{K}+\left(r+\frac{\sigma^2}{2}\right)\tau}{\sigma\sqrt{\tau}}\right)-K\Phi\left(\frac{\ln\frac{s}{K}+\left(r-\frac{\sigma^2}{2}\right)\tau}{\sigma\sqrt{\tau}}\right)$$
In practice, the quantity $S^\star$ can be computed by solving numerically Equation $(2)$, for example by Newton-Raphson. Now, the third term in Equation $(1)$ is merely the probability of a log-normal variable being above $S^\star$. By analogy with classical Black-Scholes, this is equal to: $$\tag{3.a}E_t\left(1_{S_{t_e}\geq S^\star}\right)=\Phi\left(\frac{\ln\frac{S_t}{S^\star}+\left(r-\frac{\sigma^2}{2}\right)(t_e-t)}{\sigma\sqrt{t_e-t}}\right)$$
To compute the second term in $(1)$, we define $Z$ and $Y$ as two independent normal random variables with zero mean and unit variance in order to represent the Brownian increments from $t$ to $t_e$ and from there to $T$. Note that: $$\begin{align} 1_{\{S_{t_e}\geq S^\star,\ S_T\geq K\}} &=1_{\left\{S_te^{\left(r-\frac{\sigma^2}{2}\right)(t_e-t)+\sigma \sqrt{t_e-t}Z}\geq S^\star,\ S_te^{\left(r-\frac{\sigma^2}{2}\right)(T-t)+\sigma\left(\sqrt{t_e-t}Z+\sqrt{T-t_e}Y\right)}\geq K\right\}} \\\tag{4} &=1_{\left\{Z\leq \frac{\ln\frac{S_t}{S^\star}+\left(r-\frac{\sigma^2}{2}\right)(t_e-t)}{\sigma\sqrt{t_e-t}},\ X\leq \frac{\ln\frac{S_t}{K}+\left(r-\frac{\sigma^2}{2}\right)(T-t)}{\sigma\sqrt{T-t}}\right\}} \end{align}$$ where $X$ is a third standard normal variable with the following correlation with $Z$: $$\rho:=\frac{\text{Cov}(\sqrt{t_e-t}Z+\sqrt{T-t_e}Y,Z)}{\sqrt{V(\sqrt{T-t_e}Z+\sqrt{t_e-t}Y)V(Z)}}=\sqrt{\frac{t_e-t}{T-t}}$$ Hence the second term in Equation $(1)$ is the cumulative bi-variate normal probability over $Z$ and $X$ parameterised by their correlation $\rho$: $$\begin{align} &E_t\left(1_{S_{t_e}\geq S^\star,\ S_T\geq K}\right) \\[6pt] &\qquad=\Phi_\rho\left(\frac{\ln\frac{S_t}{S^\star}+\left(r-\frac{\sigma^2}{2}\right)(t_e-t)}{\sigma\sqrt{t_e-t}},\frac{\ln\frac{S_t}{K}+\left(r-\frac{\sigma^2}{2}\right)(T-t)}{\sigma\sqrt{T-t}}\right) \end{align}$$
For the first term in Equation $(1)$, we change the measure of the expectation to the stock measure where the numéraire is the stock $S$, which we will write $\mathcal{S}$: $$\begin{align} E_t^\mathcal{Q}\left(1_{S_{t_e}\geq S^\star,\ S_T\geq K}S_T\right) &=E_t^\mathcal{S}\left(e^{r(T-t)}\frac{S_t}{S_T}1_{S_{t_e}\geq S^\star,\ S_T\geq K}S_T\right) \\ &=e^{r(T-t)}S_tE_t^\mathcal{S}\left(1_{S_{t_e}\geq S^\star,\ S_T\geq K}\right) \end{align}$$
The Radon-Nikodym process implied by this change of measure is: $$\begin{align} \left.\frac{d\mathcal{Q}}{d\mathcal{S}}\right|_{\mathscr{F}_t} &=e^{r(T-t)}\frac{S_t}{S_T} \\ &=e^{\frac{\sigma^2}{2}(T-t)-\sigma W^\mathcal{Q}_t} \end{align}$$ The "Novikov process" for this change of measure is therefore $\theta_t:=-\sigma t$. Thus the following process is the Brownian Motion under the stock measure: $$W^\mathcal{S}_t=W^\mathcal{Q}_t-\sigma t$$
This implies that the drift of the stock under this new measure increases by $\sigma^2$. We can leverage Equation $(4)$ but this time on standard normal variables $Z'$, $Y'$ and $X'$ under the stock measure $\mathcal{S}$: $$\begin{align} 1_{\{S_{t_e}\geq S^\star,\ S_T\geq K\}} &=1_{\left\{Z'\leq \frac{\ln\frac{S_t}{S^\star}+\left((r+\sigma^2)-\frac{\sigma^2}{2}\right)(t_e-t)}{\sigma\sqrt{t_e-t}},\ X'\leq \frac{\ln\frac{S_t}{K}+\left((r+\sigma^2)-\frac{\sigma^2}{2}\right)(T-t)}{\sigma\sqrt{T-t}}\right\}} \\ &=1_{\left\{Z'\leq \frac{\ln\frac{S_t}{S^\star}+\left(r+\frac{\sigma^2}{2}\right)(t_e-t)}{\sigma\sqrt{t_e-t}},\ X'\leq \frac{\ln\frac{S_t}{K}+\left(r+\frac{\sigma^2}{2}\right)(T-t)}{\sigma\sqrt{T-t}}\right\}} \end{align}$$ That is: $$\begin{align} &E_t\left(1_{S_{t_e}\geq S^\star,\ S_T\geq K}S_T\right) = \\[6pt] \tag{3.c} & \qquad e^{r(T-t)}S_t \Phi_\rho\left(\frac{\ln\frac{S_t}{S^\star}+\left(r+\frac{\sigma^2}{2}\right)(t_e-t)}{\sigma\sqrt{t_e-t}},\frac{\ln\frac{S_t}{K}+\left(r+\frac{\sigma^2}{2}\right)(T-t)}{\sigma\sqrt{T-t}}\right) \end{align}$$ Defining: $$\begin{align} d_1 & := \frac{\ln\frac{S_t}{K}+\left(r+\frac{\sigma^2}{2}\right)(T-t)}{\sigma\sqrt{T-t}} \\ d_1^\star & := \frac{\ln\frac{S_t}{S^\star}+\left(r+\frac{\sigma^2}{2}\right)(t_e-t)}{\sigma\sqrt{t_e-t}} \\[12pt] d_2 & := d_1-\sigma\sqrt{T-t} \\[18pt] d_2^\star & := d_1^\star-\sigma\sqrt{T-t} \end{align}$$ Then combining Equations $(1)$, $(3.a)$, $(3.b)$ and $(3.c)$, we obtain the desired result: $$O_t=S_t\Phi_\rho\left(d_1^\star,d_1\right)-e^{-r(T-t)} K\Phi_\rho\left(d_2^\star,d_2\right)-e^{-r(t_e-t)}P\Phi(d_2^\star) \quad \square$$
As you can see, it is very similar to the Black-Scholes Equation for a call option:
- The third term of the valuation formula for $O_t$ is equivalent to the second term in the Black-Scholes formula, namely the discounted strike $P$ times the probability that the option on the option will be exercised;
- The first two terms combined look very close to the Black-Scholes value of a call option: this is to be expected, because the option is written on a call option. However, the probabilities $\Phi_\rho(d_1^\star,d_1)$ and $\Phi_\rho(d_2^\star,d_2)$ account for the value of the stock price at $t_e$ and $T$. This is because when you enter into a vanilla call option, you know the stock price at trade date $t$ but not at expiry $T$. In this case, you will enter the option at a future time $t_e>t$, hence in this case you do not know the value of the underlying neither at inception $t_e$ nor at expiry $T$: this heightened uncertainty is captured by the bi-variate normal distribution.