Wrong means given by `RandomFunction` for `StratonovichProcess`

Question

Hello I was trying to compute mean of multiple Stratonovich integral (for $W$ standard Wiener process).

$$ J_{(1,1)} = \int_0^1 \left(\int_0^s 1\,\circ \mathrm{d}W_t\right) \circ \mathrm{d}W_s $$

Using the following code

proc = StratonovichProcess[{
    \[DifferentialD]x[t] == \[DifferentialD]w[t],
    \[DifferentialD]y[t] == x[t] \[DifferentialD]w[t]
    }, {x[t], y[t]}, {{x, y}, {0, 0}}, 
   t, {w \[Distributed] WienerProcess[]}];

symbolic Mean gives correct result:

Mean[proc[t]]
{0, t/2}

However RandomFunction gives samples with incorrect mean:

samples = RandomFunction[proc, {0, 1, 1/2^10}, 1000]["LastValues"];
Mean[samples]
{0.0270488, 0.0491107} (*Both values random*)

Answer 1

This is a comment to suggest that the issue is not a sampling issue.

Consider taking 10,000 samples rather than just 1,000:

proc = StratonovichProcess[{\[DifferentialD]x[t] == \[DifferentialD]w[t], 
  \[DifferentialD]y[t] == x[t] \[DifferentialD]w[t]}, {x[t], y[t]}, 
  {{x, y}, {0, 0}}, t, {w \[Distributed] WienerProcess[]}];
SeedRandom[12345];
samples = RandomFunction[proc, {0, 1, 1/2^10}, 10000]["LastValues"];

Here, too, the mean of the samples is around 0.5 less than the theoretical mean for the second element:

Mean[proc[1]] // N
(* {0., 0.5} *)
mean = Mean[samples]
(* {-0.0224524, -0.00163318} *)

But the theoretical and sample standard deviations are almost identical:

σ = StandardDeviation[proc[1]] // N
(* {1., 0.707107} *)
s = StandardDeviation[samples]
(* {0.998043, 0.702239} *)

Why would the mean of the second element be consistently off by 0.5 (even when the mean of the WienerProcess is changed) but the standard deviations are nearly identical?

I don't know the answer but it ain't about a finite sample.