The following P program, which can be easily ported to R, implements a standard algorithm to simulate a fractional Brownian motion on a grid of time points.

source('stdlib.P')

# returns a matrix with d rows

# each row is a simulated trajectory of fractional Brownian

# motion with parameter H

sim_fBM <- function( grid, H , d ) {

var H2, N, traj, C, i, j, L, k, x;

H2 = 2 * H

N = length( grid ) - 2

traj = matrix( d, N+1 )

C = matrix( N+1, N+1 )

for (i in (1:(N+1))) {

j = i:(N+1)

C[i, j] = 0.5 * (abs(grid[i+1]-grid[j])^H2 +

abs(grid[i]-grid[j+1])^H2 - abs(grid[i] - grid[j])^H2 -

abs(grid[i+1] - grid[j+1])^H2 )

C[j, i] = C[i, j]

}

L = t( chol( C ) )

for ( k in (1:d) ) {

x = rnorm( N+1, 0, 1 )

traj[ k, ] = L %*% x

}

return( traj )

}

fBM_test = function() {

var i, traj;

for ( i in (1:10) ) traj = sim_fBM( seq(0,1,0.001), 0.4, 4 )

return(0)

}

tm = time

call fBM_test()

print( time-tm )

R needs 7.237 seconds and the byte compiler 6.489 corresponding to a speedup factor of 1.12. P's public edition needs 4.028 seconds (speedup factor 1.8) and the Professional Edition 3.8 seconds, that is 1.9 times faster than R.

Even in this example code, where most of computing time is spent in vectorized computations and the internal function which calculates the Cholesky decomposition, P outperforms R substantially.