Computing Ornstein-Uhlenbeck Process Trajectories in Ruby

    
  #!/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