I'm teaching myself Malliavin calculus and Skorohod integrals and with this kind of math I find myself following the logic through but lacking solid intuition about what is going on.
In particular take the Skorohod integral which generalises the Ito integral to possibly anticipating integrands (it reduces to an Ito integral when the integrand is non-anticipating)
So for example take the Skorohod integral, of the Brownian motion $B(t)$, $B(0)=0$ $$ \int_0^T B(T)\delta B(t) = B^2(T)-T$$ (see for example http://www.nhh.no/Files/Filer/institutter/for/dp/1996/wp0396.pdf)
It is emphasised that $$ \int_0^T B(T)\delta B(t) \neq B(T)\int_0^T\delta B(t)=B(T)\int_0^T dB(t)=B^2(T)$$
and I am trying to understand exactly why this is the case. I appreciate that one can't just view it as a Riemann sum which is what gives the intuition that the above is trying to ward off. So the only other kind of intuition I have is that it is because the integrand and integrator are correlated in the same way that the integrand and integrator are correlated in the Stratonovich integral which means it is not a Martingale. I.e. analogously to $$E[f(B)\circ dB]\neq E[f(B)]E[dB]=E[f(B)]\cdot 0=0\quad\left(=E[f(B)dB]\right)$$
However, in which case I don't understand why the difference between the naive and correct interpretations above is $T$ so that the difference grows independently of the correlation.
In other words I would expect that (given $T_2>T_1$) $$\lim_{(T_2-T_1)\to\infty}\int_0^{T_1}B(T_2)\delta B(t)=B(T_2)\int_0^{T_1}\delta B(t)=B(T_2)\int_0^{T_1}dB(t).$$
Or that I would expect the difference $$\int_0^T B(T)\delta B(t) - B(T)\int_0^T\delta B(t)=\int_0^T\xi(t-T)dt$$
where $\xi(t-T)\to 0$ as $T\to\infty$.
But all of this is completely at odds with the definition of the Skorohod integral above. Can someone explain this in a relatively simple way with intuitive ideas?
A secondary question it raises is how to numerically simulate such an integral (say given a pre computed Wiener process)
EDIT:
So it turns out you can (I believe) write Skorohod integrals in terms of a Wick calculus built on Wick products (https://en.wikipedia.org/wiki/Wick_product). And thus create a Wick-Riemann sum viz. (recalling that $B(0)=0$) $$\int_0^TB(T)\delta B(t)=\int_0^TB(T)\diamond d B(t)=\sum_{i=0}^{n-1}B_n\diamond (B_{i+1}-B_i)\\=\sum_{i=0}^{n-1}\left[B_n\cdot(B_{i+1}-B_i)-\underbrace{\langle B_n\rangle}_{=0}(B_{i+1}-B_i)- B_n\underbrace{\langle(B_{i+1}-B_i)\rangle}_{=0}-2\underbrace{\langle B_n\rangle\langle(B_{i+1}-B_i)\rangle}_{=0}-\langle B_n(B_{i+1}-B_i)\rangle\right]\\=\sum_{i=0}^{n-1}\left[B_n\cdot(B_{i+1}-B_i)-\langle B_n(B_{i+1}-B_i)\rangle\right]\\=B_n^2-\sum_{i=0}^{n-1}\langle B_n(B_{i+1}-B_i)\rangle$$
which is intuitively much closer to where I want to be. This is the naive 'Ito' integral plus what looks a lot like a correlation function (or rather covariance). So I'm more or less on track for my 'skorohod = naive + correlation correction' train of thought.
I'm now in the strange situation where I can recover the result this way, but can't quite connect the dots understanding wise.
So The last part is $$\sum_{i=0}^{n-1}\langle B_n(B_{i+1}-B_i)\rangle=\sum_{i=0}^{n-1}\langle B_nB_{i+1}\rangle-\langle B_nB_i\rangle=\underbrace{\langle B_n^2\rangle}_{=T_n=T} - \langle B_n \underbrace{B_0}_{=0}\rangle$$
meaning we find $$\int_0^TB(T)\diamond d B(t)=B^2(T)-T$$ as expected, but I'm still struggling with the fact that the correlation based correction term doesn't converge to some constant as $T\to\infty$, but just keeps growing. Generalising such that $B(0)\neq 0$ looking at the second to last line it appears that the correction in fact converges $to$ $T$. I think this is the intuition bit I'm lacking. Any thoughts?
For the equality you wrote there is some "intuition" based on the Itô integral. Namely, write this integral as $$ \int_{0}^T B(T) d B(t) = \int_0^T B(t) dB(t) + \int_0^T \big(B(T)- B(t)\big) dB(t). $$
The first integral is a usual Itô integral equal to $(B(T)^2 - T)/2$. The second is a kind of "purely anticipative" integral: the integrand is independent of $\mathcal{F}_t$. By the change of variable $t\to T-t$ we can transform it into $\int_0^T W(t) dW(t)$, where $W(t) = B(T) - B(T-t)$ is a Brownian motion. So again $\int_0^T W(t) dW(t) = (W(T)^2-T)/2 = (B(T)^2 - T)/2$. Therefore, $$ \int_{0}^T B(T) d B(t) = B(T)^2 - T. $$