Random walks down Wall Street, Stochastic Processes in Python

stochasticprocesses.lua

--[[
--   Random walks down Wall Street, Stochastic Processes in Python
--   http://www.turingfinance.com/random-walks-down-wall-street-stochastic-processes-in-python/
--]]

local gnuplot = require 'gnuplot'

local model_parameters = {
  all_s0 = 1000, -- Starting asset value
  all_time = 800, -- Amount of time to simulate for
  all_delta = 0.00396825396, -- Rate of time e.g. 1/12 = monthly
  all_sigma = 0.125, -- Volatility of the stochastic process
  gbm_mu = 0.058, -- Annual drift factor for geometric Brownian motion
  lambda = 0.00125, -- Probability of jump happening at each point in time
  jumps_sigma = 0.001, -- Volatility of the jump size
  jumps_mu = -0.2, -- Average jump size
  cir_a = 3, -- Rate of mean reversion for Cox Ingersoll Ross
  cir_mu = 0.5, -- Long run average interest rate for Cox Ingersoll Ross
  all_r0 = 0.5, -- Starting interest rate value
  cir_rho = 0.5, -- Correlation between the Wiener processes of the Heston model
  ou_a = 3, -- Rate of mean reversion for Ornstein Uhlenbeck
  ou_mu = 0.5, -- Long run average interest rate for Ornstein Uhlenbeck
  heston_a = 0.25, -- Rate of mean reversion for volatility in the Heston model
  heston_mu = 0.35, -- Long run average volatility for the Heston model
  heston_vol0 = 0.06125 -- Starting volatility value for the Heston model
}

-- Converts a sequence of log returns into normal returns and then computes a price sequence given a starting price
local convert_to_prices = function(params, log_returns)
  local returns = torch.exp(log_returns)
  local price_sequence = torch.Tensor(returns:size(1) + 1) -- Sequence of prices starting with params.all_s0
  price_sequence[1] = params.all_s0
  for i = 1, returns:size(1) do
    price_sequence[i + 1] = price_sequence[i] * returns[i] -- Add the price at t-1 x return at t
  end
  return price_sequence
end

-- Plots a stochastic process' returns
local plot_stochastic_process = function(process, title)
  local x = torch.linspace(1, process:size(1), process:size(1))
  gnuplot.figure()
  gnuplot.plot({'', x, process, '-'})
  gnuplot.title(title)
  gnuplot.xlabel('Time')
  gnuplot.ylabel('Simulated Asset Price')
  gnuplot.axis({0, model_parameters.all_time, '', ''})
end


-- (Standard) Brownian motion stochastic process/Wiener process
local brownian_motion_log_returns = function(params)
  local sqrt_delta_sigma = math.sqrt(params.all_delta) * params.all_sigma
  return torch.Tensor(params.all_time):normal(0, sqrt_delta_sigma)
end

local brownian_motion_levels = function(params)
  return convert_to_prices(params, brownian_motion_log_returns(params))
end


-- Geometric Brownian motion stochastic process
local geometric_brownian_motion_log_returns = function(params)
  local wiener_process = brownian_motion_log_returns(params)
  local sigma_pow_mu_delta = (params.gbm_mu - 0.5 * math.pow(params.all_sigma, 2)) * params.all_delta
  return wiener_process + sigma_pow_mu_delta
end

local geometric_brownian_motion_levels = function(params)
  return convert_to_prices(params, geometric_brownian_motion_log_returns(params))
end


-- Produces a sequence of jump sizes which represent a jump diffusion process
local jump_diffusion_process = function(params)
  local s_n, time = 0, 0
  local small_lambda = -1 / params.lambda
  local jump_sizes = torch.zeros(params.all_time)
  while s_n < params.all_time do
    s_n = s_n + small_lambda * math.log(torch.uniform(0, 1))
    for j = 1, params.all_time do
      if time * params.all_delta <= s_n * params.all_delta and s_n * params.all_delta <= j * params.all_delta then
        jump_sizes[j] = jump_sizes[j] + torch.normal(params.jumps_mu, params.jumps_sigma)
        break
      end
    end
    time = time + 1
  end
  return jump_sizes
end

-- Geometric Brownian motion jump diffusion/Merton jump diffusion stochastic process
local geometric_brownian_motion_jump_diffusion_log_returns = function(params)
  local jump_diffusion = jump_diffusion_process(params)
  local geometric_brownian_motion = geometric_brownian_motion_log_returns(params)
  return jump_diffusion + geometric_brownian_motion
end

local geometric_brownian_motion_jump_diffusion_levels = function(params)
  return convert_to_prices(params, geometric_brownian_motion_jump_diffusion_log_returns(params))
end


-- Heston stochastic volatility process
local cox_ingeroll_ross_heston = function(params)
  local sqrt_delta_sigma = math.sqrt(params.all_delta) * params.all_sigma
  local brownian_motion_volatility = torch.Tensor(params.all_time):normal(0, sqrt_delta_sigma)
  local volatilities = torch.Tensor(params.all_time)
  volatilities[1] = params.heston_vol0
  for i = 2, params.all_time do
    local drift = params.heston_a * (params.heston_mu - volatilities[i - 1]) * params.all_delta
    local randomness = math.sqrt(volatilities[i - 1]) * brownian_motion_volatility[i - 1]
    volatilities[i] = volatilities[i - 1] + drift + randomness
  end
  return brownian_motion_volatility, volatilities
end

local heston_construct_correlated_path = function(params, brownian_motion_one)
  local sqrt_delta = math.sqrt(params.all_delta)
  local brownian_motion_two = torch.Tensor(params.all_time)
  for i = 1, params.all_time do
    local term_one = params.cir_rho * brownian_motion_one[i]
    local term_two = math.sqrt(1 - math.pow(params.cir_rho, 2)) * torch.normal(0, sqrt_delta)
    brownian_motion_two[i] = term_one + term_two
  end
  return brownian_motion_one, brownian_motion_two
end

local heston_model_levels = function(params)
  local brownian, brownian_motion_market, cir_process
  brownian, cir_process = cox_ingeroll_ross_heston(params)
  brownian, brownian_motion_market = heston_construct_correlated_path(params, brownian)

  local heston_market_price_levels = torch.Tensor(params.all_time)
  heston_market_price_levels[1] = params.all_s0
  for i = 2, params.all_time do
    local drift = params.gbm_mu * heston_market_price_levels[i - 1] * params.all_delta
    local vol = cir_process[i - 1] * heston_market_price_levels[i - 1] * brownian_motion_market[i - 1]
    heston_market_price_levels[i] = heston_market_price_levels[i - 1] + drift + vol
  end
  return heston_market_price_levels, cir_process
end


-- Cox-Ingeroll-Ross stochastic process
local cox_ingeroll_ross_levels = function(params)
  local brownian_motion = brownian_motion_log_returns(params)
  local levels = torch.Tensor(params.all_time)
  levels[1] = params.all_r0
  for i = 2, params.all_time do
    local drift = params.cir_a * (params.cir_mu - levels[i - 1]) * params.all_delta
    local randomness = math.sqrt(levels[i - 1]) * brownian_motion[i - 1]
    levels[i] = levels[i - 1] + drift + randomness
  end
  return levels
end


-- Ornstein-Uhlenbeck stochastic process
local ornstein_uhlenbeck_levels = function(params)
  local brownian_motion_returns = brownian_motion_log_returns(params)
  local ou_levels = torch.Tensor(params.all_time)
  ou_levels[1] = params.all_r0
  for i = 2, params.all_time do
    local drift = params.ou_a * (params.ou_mu - ou_levels[i - 1]) * params.all_delta
    local randomness = brownian_motion_returns[i - 1]
    ou_levels[i] = ou_levels[i - 1] + drift + randomness
  end
  return ou_levels
end


-- Plot all processes
plot_stochastic_process(brownian_motion_levels(model_parameters), 'Brownian')
plot_stochastic_process(geometric_brownian_motion_levels(model_parameters), 'Geometric Brownian')
plot_stochastic_process(geometric_brownian_motion_jump_diffusion_levels(model_parameters), 'Merton Jump Diffusion')
plot_stochastic_process(heston_model_levels(model_parameters), 'Heston')
plot_stochastic_process(cox_ingeroll_ross_levels(model_parameters), 'Cox-Ingeroll-Ross')
plot_stochastic_process(ornstein_uhlenbeck_levels(model_parameters), 'Ornstein-Uhlenbeck')