Simulation of Brownian Motion If I want to simulate Brownian motion in the Euclidean space I can simulate it by a point that is moving a distance $\epsilon$ in an arbitrary direction then it randomly choose a new direction and moves a distance $\epsilon$ again and so on. The smaller the $\epsilon$ the closer the simulation will be to the real Brownian motion.
How can I simulate Brownian motion in the hyperbolic space (Poincare Disk model for instance)? Does the same work here where I replace the Euclidean distance by the hyperbolic distance? My intuition is yes but when I did the simulation the random walk do not seem to be transient but it should be!
 A: Taking the Poincaré disk model for hyperbolic space as being formed by all complex numbers with modulus smaller than 1, the distance between two points can be calculated as
$$d(z_1,z_2)=\tanh^{-1}\left|\frac{z_1-z_2}{1-z_1\bar{z_2}}\right| \; .$$
The set of points at equal distance $\varepsilon=\tanh(d(z,z_0))$ from a fixed point $z_0$ can be shown to be represented by a circle in the Poincaré disk model with center
$$\frac{(1-\varepsilon^2)z_0}{1-\varepsilon^2|z_0|^2}$$
and radius 
$$\frac{\varepsilon(1-|z_0|^2)}{1-\varepsilon^2|z_0|^2} \; .$$
But this is not the end of the story, you can't just pick points uniformly on this circle. What you want is to pick points such that you can arrive from $z_0$ to them such that the direction of departure from $z_0$ was chosen uniformly. This will induce a particular distribution on that circle.
To figure out this distribution, it is necessary to compute the intersection of an h-line* through $z_0$ for an arbitrary direction and the circle we just computed. From that, we can get the sought distribution. I'll update this later if I get to computing that part.
EDIT: Following the suggestion by @Douglas Zare, this can also be achieved by picking a point at a distance $\varepsilon$ from the origin uniformly, where the deformation of the Poincaré disk model won't be a problem since h-lines passing through the origin are all just diameters. All we have to do then is translate the point with the translation that takes the origin to the point $z_0$ which can be done with a Möbius transformation
$$z \mapsto w=\frac{z+z_0}{\bar{z}_0z+1} \; .$$
This is also easily implemented computationally. Here's how I coded it in R:
    # Brownian motion on Poincaré disk (no comments about crappy code plz kthx)

epsilon=0.01;
z<-0;
path<-c(z);
for (t in 1:10000) { 
    jitter=runif(1,0,1);
    dz=epsilon*complex(1,cos(2*pi*jitter),sin(2*pi*jitter));
    z<-(z+dz)/(Conj(dz)*z+1);
    path<-c(path,z);
}

plot(Re(path),Im(path),type="l",col="blue",asp=1,xlim=c(-1,1),ylim=c(-1,1));
curve(sqrt(1-x^2),-1,1,col="red",add=TRUE);
curve(-sqrt(1-x^2),-1,1,col="red",add=TRUE);

And here is a picture of a run of the code

*By h-line, I mean a line in the hyperbolic plane, which will be either a diameter, either a circle segment orthogonal on the boundary of the Poincaré disk.
A: IIRC, there a discussion of Brownian motion in Feller (An Introduction To Probability Theory and it's Applications). 
A: I think that the start will be to define what is a Brownian motion on certain subset of $\mathbb{R}^3$ in this case. 
The key property of Brownian motion $B_t$ is that it has independent increments, i.e. $B_t-B_s$ is independent of $B_s$. Hence one way of simulating Brownian motion is by summing up independent random variables. For two-dimensional Brownian motion this is achieved by generating $n$ random vectors $X_i$ and then taking
$$\xi_t=\frac{1}{\sqrt{n}}\left(\sum_{i=1}^{[nt]}X_i+(nt-[nt])X_{[nt]+1}\right),$$
where $[\cdot]$ is whole part of a real number. The functional central limit theorem ensures that $\xi_t$ converges to $B_t$ as $T$ goes to infinity, so take large $T$ and you will get a simulation of $B_t$. It is possible also to take only the partial sums $\sum_{i=1}^{[nt]}X_i$, the end result will be the same.
Now to apply this approach to a certain space, you need,  as I said before the definition of the Brownian motion is needed in said space and the corresponding functional central limit theorem. My own impression when writing my thesis (which was about multi-parameter Wiener processes) was that functional central limit theorems go hand in hand with the brownian motion, i.e. if you can define Brownian motion, you can prove the FCLT. 
Before though delving into literature I would advise to find a way how to generate random sample from the space you are intend to work in. This might be the most difficult part in simulation.
