The problem
I'm trying to perform a Bayesian approach to the Maximum Likelihood Estimation procedure of Wecker and Ansley (1983). To this end, I need to compute the full likelihood of the data given the model, which after some algebra turns out to have a covariance matrix whose elements are given by:
$$\Sigma_{j,k}=\mathbb{E}[Z^{m}(x_j)Z^{m}(x_k)]+\delta_{j,k}\sigma^2,$$
with $\delta_{j,k}$ the Kronecker delta and where
$$Z^{m}(x_i) = \lambda^{1/2}\sigma \int_{x_1}^{x_i}\frac{(x_i-h)^{m-1}}{(m-1)!}dW(h),$$
and $\lambda^{1/2},\sigma$ are constants and $W(h)$ is a Wiener process with unit dispersion parameter. Here is where I'm stuck (I'm fairly new to stochastic calculus): I need to obtain this expected value. The authors say this matrix should be of the form $\Sigma = \sigma^2\Lambda$, where the $\Lambda$ matrix depends only on $\lambda$, but I need the full functional form of the $\Sigma$ matrix in order to write the full likelihood!
What I have tried
The whole problem is to obtain
$$\mathbb{E}[Z^{m}(x_j)Z^{m}(x_k)]=\lambda \sigma^2 \mathbb{E}\left[\int_{x_1}^{x_j}\frac{(x_j-h)^{m-1}}{(m-1)!}dW(h)\int_{x_1}^{x_k}\frac{(x_k-h)^{m-1}}{(m-1)!}dW(h)\right],$$
and I know, by Ito's isometry that, for $s<t$: $$\mathbb{E}\left[\int_0^s f(u,W(u))dW(u) \int_0^{t}g(v,W(v))dW(v)\right]=\mathbb{E}\left[\int_0^s f(t,W(t))g(t,W(t))dt\right],$$
so, can I say (note the lower integration limit) that, for $x_1<x_j<x_k$,
$$\mathbb{E}\left[\int_{x_1}^{x_j}\frac{(x_j-h)^{m-1}}{(m-1)!}dW(h)\int_{x_1}^{x_k}\frac{(x_k-h')^{m-1}}{(m-1)!}dW(h')\right]=\mathbb{E}\left[\int_{x_1}^{x_j}\frac{(x_j-t)^{m-1}}{(m-1)!}\frac{(x_k-t)^{m-1}}{(m-1)!}dt\right]?$$ Any help will be very much appreciated!
In full generality, for deterministic functions $B$ and $C$, and for $a\leqslant b$, $a\leqslant c$, $$E\left(\int_a^bB(h)\mathrm dW(h)\cdot\int_a^cC(h)\mathrm dW(h)\right)=\int_a^{\min\{b,c\}}B(h)C(h)\mathrm dh.$$