Can we transform $\int_\rho^1 (W_t - W_{t-\rho}) \,dW_t$ to make its law $\rho$-invariant? I just bumped into the stochastic integral
$$
\int_\rho^1 (W_t - W_{t-\rho}) \,dW_t
$$
where $0 < \rho < 1$ is a constant and $W$ is a standard Wiener process. It would be nice if we have a closed-form representation of the cumulative distribution function, but for my work it suffices to show that there is some simple transform to make its law $\rho$-invariant, that is, to make its law not dependent on $\rho$.
Does it appear possible to any one? Thanks for any help in advance!

Additional Information: I guessed so because when $\rho$ is large enough, $\int_\rho^1 (W_t - W_{t-\rho})\,dW_t \approx W_\rho (W_1 - W_\rho)$. I tried to multiply it by $\frac{1}{\sqrt{\rho(1-\rho)}}$, but in some simple Monte Carlo experiments the distribution of the stochastic integral doesn't seem to be $\rho$-invariant.
 A: (This is not a complete answer, rather an attempt to describe the distribution of the integral.)

The integral $\int_\rho^t W(t) dW(t)$ is simply $\tfrac{1}{2} ((W(1))^2 - 1) - \tfrac{1}{2} ((W(\rho))^2 - \rho)$. Let us handle the other term: $$I_\rho = \int_\rho^1 W(t-\rho) dW(t).$$
Let us approximate $W(t-\rho)$ by an elementary process $X_n(t)$, which is equal to $0$ for $t \in [0, \rho)$, and to $W_{k \rho/n - \rho}$ for $t \in [\tfrac{k}{n} \rho, \tfrac{k+1}{n} \rho)$, $k = n, n + 1, n + 2, \ldots$ Clearly,
$$ I = \lim_{n \to \infty} \int_\rho^{K(n) \rho / n} X_n(t) dW(t) ,$$
where $K(n) = \lfloor n/\rho \rfloor$. Furthermore,
$$ I_n := \int_\rho^{K(n) \rho / n} X_n(t) dW(t) = \sum_{k = n}^{K(n)} W(\tfrac{k}{n} \rho - \rho) (W(\tfrac{k + 1}{n} \rho) - W(\tfrac{k}{n} \rho)) $$
By scaling,
$$ I_n \stackrel{D}{=} \frac{\rho}{n} \sum_{k = n}^{K(n)-1} W(k - n) (W(k + 1) - W(k)) . $$
Write $Y_k = W(k+1) - W_k$, so that $Y_k$ is an i.i.d. sequence of standard Gaussian random variables. Then
$$ I_n \stackrel{D}{=} \frac{\rho}{n} \sum_{k = n}^{K(n)-1} \sum_{j = 0}^{k - n - 1} Y_j Y_k. $$
It follows that $I_n$ is a certain quadratic form, evaluated on a standard Gaussian vector. If $\lambda_{n,j}$ denote the eigenvalues of the corresponding (symmetric) matrix
$$ A_n = [\tfrac{\rho}{2 n} \mathbb{1}_{|j-k| \geqslant n}]_{j,k = 1}^{K(n)}, $$
then
$$ I_n \stackrel{D}{=} \sum_{j = 1}^{K(n)} \lambda_{n,j} Z_j^2, $$
with $Z_j$ another i.i.d. sequence of standard Gaussian random variables.

Numerically, the eigenvalues $\lambda_{n,j}$ (arranged in the decreasing order of their absolute values) appear to converge to a certain limit $\lambda_j$, which essentially depends on $\rho$. This would mean that
$$ I \stackrel{D}{=} c + \sum_{j = 1}^{\infty} \lambda_j (Z_j^2 - 1), $$
with $\lambda_j$ (and possibly $c$) essentially dependent on $\rho$. This, in turn, implies that there is no simple transform of $I$ that would make its distribution invariant with respect to $\rho$.

The numbers $\lambda_j$ are the eigenvalues of the integral operator with kernel $K(s, t) = \tfrac{\rho}{2} \mathbb{1}_{|s - t| > \rho}$ (acting on $L^2([0,1])$). However, I do not know if this helps in establishing a closed-form expression for $\lambda_j$.
A: From $\int_0^1W_t\,dW_t={1\over 2}(W_t^2-t)$.
One can deduce with the change of variable $s=t-\rho$, that your integral is equal to:
$$
{1\over 2}(W_1^2-W_\rho^2-W_{1-\rho}^2)
$$
