Diagnosing Numerical Problems in Ornstein–Uhlenbeck MLE

# How to fix my Ornstein-Uhlenbeck parameter MLE in Python?


# How to fix my Ornstein-Uhlenbeck parameter MLE in Python?












I am trying to fit time-series data into an Ornstein-Uhlenbeck process. Here is my code so far:

```
# source for computation: https://arxiv.org/pdf/1411.5062.pdf
import math
from math import sqrt, exp, log  # exp(n) == e^n, log(n) == ln(n)
import scipy.optimize as so
import numpy as np

def __compute_log_likelihood(params, *args):
    '''
    Compute the average Log Likelihood, this function will by minimized by scipy.
    Find in (2.2) in linked paper

    returns: the average log likelihood from given parameters
    '''
    # functions passed into scipy's minimize() needs accept one parameter, a tuple of
    #   of values that we adjust to minimize the value we return.
    #   optionally, *args can be passed, which are values we don't change, but still want
    #   to use in our function (e.g. the measured heights in our sample or the value Pi)

    theta, mu, sigma = params
    X, dt = args
    n = len(X)

    sigma_tilde_squared = sigma ** 2 * (1 - exp(-2 * mu * dt)) / 2 * mu

    summation_term = 0

    for i in range(1, len(X)):
        summation_term += (X[i] - X[i - 1] * exp(-mu * dt) - theta * (1 - exp(-mu * dt))) ** 2

    summation_term = -summation_term / (2 * n * sigma_tilde_squared)

    log_likelihood = (-log(2 * math.pi) / 2) + (-log(sqrt(sigma_tilde_squared))) + summation_term

    return -log_likelihood
    # since we want to maximize this total log likelihood, we need to minimize the
    #   negation of the this value (scipy doesn't support maximize)

def estimate_coefficients_MLE(X, dt):
    '''
    Estimates Ornstein-Uhlenbeck coefficients (θ, µ, σ) of the given array
    using the Maximum Likelihood Estimation method

    input: X - array-like data to be fit as an OU process
    returns: θ, µ, σ, Total Log Likelihood
    '''

    bounds = ((0, None), (None, None), (0, None))  # theta > 0, mu ∈ ℝ, sigma > 0
    mu_init = np.mean(X)
    result = so.minimize(__compute_log_likelihood, (1e-6, 1e-6, 1e-6), args=(X, dt), bounds=bounds)
    theta, mu, sigma = result.x
    max_log_likelihood = -result.fun  # undo negation from __compute_log_likelihood
    return theta, mu, sigma, max_log_likelihood
```

But when I simulate an OU process with the following:

```
# simulate Ornstein-Uhlenbeck Process
import numpy as np
import matplotlib.pyplot as plt
t_0 = 0 # define model parameters
t_end = 2
length = 1000
theta = 1.1
mu = 0
sigma = 0.3
t = np.linspace(t_0,t_end,length) # define time axis
dt = np.mean(np.diff(t))

y = np.zeros(length)
y0 = np.random.normal(loc=0.0,scale=1.0) # initial condition
drift = lambda y,t: theta*(mu-y) # define drift term, google to learn about lambda
diffusion = lambda y,t: sigma # define diffusion term
noise = np.random.normal(loc=0.0,scale=1.0,size=length)*np.sqrt(dt) #define noise process
# solve SDE
for i in range(1,length):
    y[i] = y[i-1] + drift(y[i-1],i*dt)*dt + diffusion(y[i-1],i*dt)*noise[i]

plt.plot(t,y)
plt.show()
```

Then fit the data (stored in y) using my function with:

```
theta, mu, sigma, max_ll = estimate_coefficients_MLE(y, 1/len(y))
```

I either get a "Value Error: math domain error" or my coefficients are very off. If someone could point me in the right direction, I would be very grateful, there is a lack of resources online on this subject.

## Answer by Shile Wen (score 1, accepted)

https://quant.stackexchange.com/a/59591

Adding on to Japser's answer, to address the division by 0, we can set a very small value for the lower bound for mu and sigma (e.g. 1x10^-5). To see the algorithm in action, see this

## Answer by Japser Wang (score 2)

https://quant.stackexchange.com/a/59571

That is because of `sigma_tilde_squared == 0`, You could add 0.01 at the add to avoid it `== 0`