May be with the code I will have more chances to get an answer:

restart;interface(rtablesize=20):printlevel=100:
                                  10

> with(LinearAlgebra):with(VectorCalculus):with(linalg):with(Student[LinearAlgebra]):
> a1:=1:a2:=2:
> T:=:
> dx:=sqrt((x1-x3)^2+(x2-x4)^2):
> x1dot:=a2/(2*Pi*dx^2)*:
> x2dot:=a1/(2*Pi*dx^2)*:
> dy1:=sqrt((x1-y1)^2+(x2-y2)^2):
> dy2:=sqrt((x3-y1)^2+(x4-y2)^2):
> ydot:=a1/(2*Pi*dy1^2)*+a2/(2*Pi*dy2^2)*+:
> xdot:=:
> zdot:=:
> l:=convert(,Vector):
> g1:=:g2:=:
Lie Bracket Function
> LB:=proc(m,w)
> convert(evalm(multiply(Jacobian(convert(w,Vector),convert(T,list)),m)-multiply(Jacobian(convert(m,Vector),convert(T,list)),w)),Vector[column]):
> end:

> C:=Matrix([g1,g2,LB(l,g1),LB(l,g2),LB(l,LB(l,g1)),LB(l,LB(l,g2)),LB(l,LB(l,LB(l,g1))),LB(l,LB(l,LB(l,g2))),LB(l,LB(l,LB(l,LB(l,g1)))),LB(l,LB(l,LB(l,LB(l,g2))))]):

> Rank(C);

> Basis({g1,g2,LB(l,g1),LB(l,g2),LB(l,LB(l,g1)),LB(l,LB(l,g2)),LB(l,LB(l,LB(l,g1))),LB(l,LB(l,LB(l,g2))),LB(l,LB(l,LB(l,LB(l,g1)))),LB(l,LB(l,LB(l,LB(l,g2))))});

I want to know how can I display the steps for getting either the Rank(C) or Basis(...)

Thanks.

Ok now I am trying with k=1 and it is working. So the problem seems only for k=0.

I don't know why because zeta contains also eta[1](0)[1] for example for k=0 and it is working without any problem.

Alpha and eta are defined exactly the same way, so does somebody see from where the problem is coming from?

For another example, it works fine:

d:=1.0:

Rxi:=sqrt(d)*IdentityMatrix(4):

xi:=WienerProcess(Rxi):

Reta:=sqrt(d)*IdentityMatrix(2):

eta:=Array([WienerProcess(Reta),WienerProcess(Reta),WienerProcess(Reta),WienerProcess(Reta)]):

zeta:=k->Vector[column]([[xi(k+d)[1]-xi(k)[1]],[xi(k+d)[2]-xi(k)[2]],[xi(k+d)[3]-xi(k)[3]],[xi(k+d)[4]-xi(k)[4]],[eta[1](k+d)[1]-eta[1](k)[1]],[eta[1](k+d)[2]-eta[1](k)[2]],[eta[2](k+d)[1]-eta[2](k)[1]],[eta[2](k+d)[2]-eta[2](k)[2]],[eta[3](k+d)[1]-eta[3](k)[1]],[eta[3](k+d)[2]-eta[3](k)[2]],[eta[4](k+d)[1]-eta[4](k)[1]],[eta[4](k+d)[2]-eta[4](k)[2]]]):

ExpectedValue(zeta(0)[1]*zeta(0)[1],timesteps=100,replications=10^4);

for this example, it works. I can calculate the expected value at any time k and compute a matric which has the same looking as sigma5.

I don't know what is the difference between the two :/ 

Thanks again, your help is making me understand a lot of details.

So now I changed everything in terms of wiener process, because reading the help it turns out that it is what I wanted.

I have another question about Finance[ExpectdValue], why

ExpectedValue(zeta(1)[1]);

gives an answer without adding the timesteps and replications?

(zeta is a vector made of wiener processes)

To be sure that I understood:

If I do:

 xi0:=Vector([[0],[0],[0],[0]]):
> Mxi:=Vector([[0],[0],[0],[0]]):
> Rxi:=sqrt(d)*IdentityMatrix(4):
> xi1:=BrownianMotion(xi0,Mxi,Rxi):

> Statistics[ExpectedValue](xi1);

The answer is _X. What does that mean? Why this is not zero?

And shall I take a Brownian Motion or a Wiener Process?

Thank you very much for your time. I think you got the problem. I was trying to calculate the covariance between two floating numbers and it is obviously zero!

I will re-do everything without generating sample paths. I only need gaussian random variables but I did wrong by taking samples.

At least with taking samples I can verify that my formulas are correct, I have more than a hundred lines of code so it was useful to check :)

I will try to do everything symbolically and keep you guys updated.

Thanks angain for the explanations.