Computing Ornstein-Uhlenbeck Process Trajectories in Ruby
This code is adapted from some Matlab code I found that simulates the OU process exactly. I had to code up a quick Gaussian random number generator in Ruby because I didn’t find a method to handle that. As noted in the comments I used the Box-Muller Transformation but if lots of random numbers are required the Ziggurat Algorithm could be used instead.
#!/usr/bin/env ruby include Math require 'rubygems' require 'gnuplot' #Gaussian random number generator using the Box-Muller Transformation #If this is too slow it can be replaced with the Ziggurat Algorithm def randn z1 = rand z2 = rand rand_normal = sqrt(-2.0*log(z1))*sin(2.0*PI*z2) return rand_normal.to_f end #Enter the parameters for the simulation steps = 400 #start_time = 0 #simulation start time #end_time = 400 #simuation end time dt = 0.01 #time step tau = 0.1 #relaxation time c = 1.0 #diffusion constant x0 = 0.0 x =  #initial value for stochastic variable x mu = 0.0 #mean of stochatic process x y0 = 0.0 y =  #initial value for integral x start_dist = -2.0 #start of OU pdf end_dist = 2.0 #end of OU pdf time = (0..steps).step(dt) i = 0 x = x0 y = y0 time.each do |t| i = i + 1 r1 = randn r2 = randn puts x[i-1] x[i] = x[i-1] * exp(-dt/tau) + sqrt((c*tau*0.5)*(1-(exp(-dt/tau))**2))*r1 y[i] = y[i-1] + x[i-1]*tau*(1-exp(-dt/tau))+sqrt((c*tau**3*(dt/tau-2* \ (1-exp(-dt/tau))+0.5*(1-exp(-2*dt/tau))))-((0.5*c*tau**2)* \ (1-exp(-dt/tau))**2)**2/((c*tau/2)*(1-exp(-2*dt/tau))))*r2+ \ ((0.5*c*tau**2)*(1-exp(-dt/tau))**2)/(sqrt((c*tau/2)* \ (1-(exp(-dt/tau))**2)))*r1 end k = 0 j = (start_dist..end_dist).step(dt) p =  j.each do |l| k = k + 1 p[k] = sqrt((1/tau)/(PI*c))*exp(-(1/tau)*(l-mu)**2/(c)) end
That’s it. Not too difficult to compute with Ruby, but I also did another version in R. I’ll post about that shortly.