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[0] = x0
y[0] = 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.