El $P$ dinámica de los activos subyacentes son:
\begin{align*}
dS=S(\mu dt+\sigma dB_t)
\end{align*}
Que tiene la siguiente solución debajo de los $\mathcal{P}$ dinámica:
\begin{align*}
S_t=S_0 e^{(r-\frac{\sigma^2}{2})t+\sigma W_t}
\end{align*}
Donde $W_t$ es el equivalente de martingala con respecto a la original movimiento browniano geométrico. Definir $Y_t=\int_0^t S_u du$, entonces, de acuerdo a Feynman-Kac el valor de replicar la cartera está dada por
\begin{align*}
V_t&=e^{-r(T-t)}\frac{1}{T^2}\mathbb{E}_t^\mathcal{P}[(Y_t+\int_t^T S_udu)^2)\\
&=\frac{e^{-r(T-t)}}{T^2} \mathbb{E}_t^\mathcal{P}[Y_t^2 + 2Y_t \int_t^T S_u du + (\int_t^T S_u du)^2]\\
\end{align*}
Así:
\begin{align*}
\mathbb{E}_t^\mathcal{P}[\int_t^T S_u du]&=S_t \mathbb{E}_t^\mathcal{P}[\int_t^T \frac{S_u}{S_t} du]\\
&=S_t \int_t^T e^{(r-\frac{\sigma^2}{2})(u-t)}\mathbb{E}^\mathcal{P}[e^{\sigma (W_u-W_t)}] du\\
\end{align*}
Desde $W_u-W_t$ siguiente $\mathcal{N}(0, u-t)$ bajo $P$. De acuerdo a la m.g.f.
\begin{align*}
&=S_t \int_t^T e^{i(u-t)}du\\
&=\frac{S_t}{r}. (e^{r(T-t)}-1)\\
&=\frac{x}{r}. (e^{r(T-t)}-1)
\end{align*}
Permite calcular ahora la expectativa de la plaza de la integral:
\begin{align*}
\mathbb{E}_t^\mathcal{P}[(\int_t^T S_u du)^2]&=S_t^2 \mathbb{E}_t^\mathcal{P}[\int_t^T \frac{S_u}{S_t}du \int_t^T \frac{S_v}{S_t}dv ]\\
&= S_t^2 \int_t^T \int_t^T \mathbb{E}^\mathcal{P}[\frac{S_u}{S_t} \frac{S_v}{S_t}]dv du\\
\end{align*}
Le permite enfocarse en: $\mathbb{E}^\mathcal{P}[\frac{S_u}{S_t} \frac{S_v}{S_t}]$, permite assum $t \leq v \leq u$
\begin{align*}
\mathbb{E}^\mathcal{P}[\frac{S_u}{S_t} \frac{S_v}{S_t}&=\mathbb{E}^\mathcal{Q}[e^{(r-\frac{\sigma^2}{2})(u-t)+\sigma(W_u-W_t))}e^{(r-\frac{\sigma^2}{2})(v-t)+\sigma(W_v-W_t))}]\\
&=e^{(r-\frac{\sigma^2}{2})(u-t)+(r-\frac{\sigma^2}{2})(v-t)}\mathbb{E}^\mathcal{Q}[e^{\sigma(W_u-W_v)+2\sigma (W_v-W_t)}]\\
&=e^{(r-\frac{\sigma^2}{2})(u-t)+(r-\frac{\sigma^2}{2})(v-t)}\mathbb{E}^\mathcal{Q}[e^{\sigma(W_u-W_v)}]\mathbb{E}^\mathcal{Q}[e^{2\sigma(W_u-W_t)}]\\
&=e^{(r-\frac{\sigma^2}{2})(u-t)+(r-\frac{\sigma^2}{2})(v-t)}e^{\frac{\sigma^2}{2}(u-v)}e^{2\sigma^2(v-t)}\\
&= e^{ur}e^{u(r+\sigma^2)}e^{-t(2r+\sigma^2)}
\end{align*}
Así:
\begin{align*}
\mathbb{E}_t^\mathcal{P}[(\int_t^T S_u du)^2]&=2S_t^2 \int_t^T \int_t^u e^{ur}e^{u(r+\sigma^2)}e^{-t(2r+\sigma^2)}dv du\\
& = \frac{2x^2}{r+\sigma^2}(\frac{1}{2r+\sigma^2}e^{(2r+\sigma^2)(T-t)}-\frac{1}{r}e^{r(T-t)}+\frac{r+\sigma^2}{(2r+\sigma^2)r})
\end{align*}
Mezcla todo junto, se obtiene:
\begin{align*}
V(t,x,y)&=\frac{y^2}{T^2}e^{rt-rT}+\frac{1}{T^2}\frac{2xy}{r}(1-e^{rt-rT})\\
& + \frac{1}{T^2}\frac{2x^2}{r+\sigma^2}(\frac{1}{2r+\sigma^2}e^{(r+\sigma^2)(T-t)}-\frac{1}{r}+\frac{r+\sigma^2}{(2r+\sigma^2)r}e^{-r(T-t)})
\end{align*}