# Estimate Ornstein-Uhlenbeck Parameters with Euler Maruyama

#1

Hello,
I want estimate the parameters of 2D Ornstein Uhlenbeck Process, my data are the following:
I’m wondering if the attempt to estimate the parameters of this model with Euler-Maruyama class is correct.
Given the OU sde

dX_t = \Theta(\mu - X_t)dt + \sigma dW_t

I’ve built the following model:

with pm.Model() as model:

Theta = pm.Flat('Theta', shape=2)
sigma  = pm.Flat('sigma', shape=2, testval=1)
u = pm.Uniform('u', -1, 1, shape=2)

def lin_sde(x, Theta, u):
return Theta * (u-x), sigma

with model:
# hidden states following a OU SDE
X = EulerMaruyama('X', dt,  lin_sde, (Theta, u), shape=(len(eweValence), 2), testval=np.zeros(shape=(len(eweValence),2)))
Z = pm.Normal('Z', mu=X, sd=sigma, shape=(len(eweValence), 2), observed = obs)
inference_result = Inference('MCMC')


the Inference function is:

   def MCMC():
start = pm.find_MAP(vars=[X], method="L-BFGS-B")

step = pm.NUTS(scaling=start)
trace = pm.sample(2000, step=step)

step = pm.NUTS(scaling=trace[-1], gamma=.25)
trace = pm.sample(1000, step, start=trace[-1], progressbar=True, random_seed=SEED)
return trace


The trace result and sample_ppc() are the following:

the first plot is the real data.
Is this a good attempt to estimate the parameters? and \theta should not it be positive? How can I proof the model’s correctness?

#2

The PPC is actually looks not bad, but the trace is clearly not converge yet.
A few suggestions:

1. Use more informative priors. Flat priors are usually really bad ideas. For example, if you are expecting \theta to be positive, you should use a prior distribution with only positive support.
2. Try using the default sampling kwarg, instead of initialising using MAP: trace = pm.sample(1000, tune=2000)

#3

Thank you very much for the precious help @junpenglao, I followed your suggestions and I used a Beta distribution for \theta, a Normal distribution for u and I have increased the number of the observed data.
The result is the following:

It seems much better than before.
What I want to do now is to make the inference with ADVI but this is my poor result using the following code:

  def ADVI():

trace = approx.sample(draws=500)
ppc_trace = pm.sample_ppc(trace, model=model)



Any suggestion to increase this outcome? thanks in advance

#4

ADVI ignore the posterior correlation, which in many cases gives sub-optimal results.
As general tips, I would suggest the following:

1. Not using ADVI. NUTS is more or less the best thing we have for models it can be sampled. If you want a VI inference that scale a bit better FullRankADVI is the first thing to try
2. Make sure ADVI/fullrankADVI did converge. For example, you can find in this post some discussion about making sure the fitting converge: Poor Accuracy of ADVI for Linear Regression

#5

I am not sure why you use an approximate method for Ornstein-Uhlenbeck. The probabilistic solutions is known analytically (G. Uhlenbeck and I. Ornstein, Phys. Rev. 36, 823 (1930)). I have developed an analytical method to solve the OU likelihood (https://arxiv.org/abs/1805.05977) but I also wrote a OU time-series class for pymc3 (https://github.com/hstrey/Ornstein-Uhlenbeck-Bayesian). This method can easily be extended to two dimensions.

Hope that helps

#6

Thank you very much, you are right. I used an approximate method for have a scalable model, but for this particular problem is not the best. In the next days I will check your paper and your pymc class.