(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$.
