Compute $\mathbb{E}\left(W_s^2\left(\int_s^t W_r d W_r\right)^2\right)$

Question

Let $W=\left(\Omega, \mathcal{F},\left(\mathcal{F}_t\right)_{t \geq 0},\left(W_t\right)_{t \geq 0}, \mathbb{P}\right)$ be a Brownian motion. I have tried the following, see photo.

The simulation code in R, is the following for anyone interested:

rm(list=ls())
library('e1071')
freq    = 10000
s <- 0.3
t <- 1
true <- numeric(2000)
means2 <- numeric(2000)
means3 <- numeric(2000)
for (i in 1:length(true)) {
  W  = rwiener(end = t, frequency = freq)
 
  true[i] = W[s*freq]^2 * sum( ( 1/freq) * (W[(s*freq):(t*freq)] )^2  )
                               
  means2[i] =  sum( ( 1/freq) * ( W[(s*freq):(t*freq)] - W[s*freq] )^2  )
  means3[i] = sum( ( 1/freq) *( (W[s*freq]-W[(s*freq):(t*freq)] )^2  ) )
}
mean(true)
(t-s)*3*s^2 + s*mean(means2)
(t-s)*3*s^2 + s*mean(means3)
((t-s)*3*s^2 + s*(t^2-s^2)/2 + (t-s)*s)/2
```

Answer 1

You can actually somehow still "apply the conditional Ito isometry", if you use the Tower law, where you condition on $\mathcal{F}_s$, so you do get $E[\int_s^t W_s^2W_r^2dr]$, calculating this is now easy and you obtain the final answer $-(5/2)s^3 + st^2/2+2s^2t$, which is also what the simulation shows.