# Estimate Ornstein-Uhlenbeck Parameters with Euler Maruyama

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?

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)

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:

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

return advi, trace ,ppc_trace

Any suggestion to increase this outcome? thanks in advance

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
1 Like

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

1 Like

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.

Dear all,

I am working on the same topic too.
But how can I change the prior from univariate ot multivariate?

I can change to MvNormal for my volatility_mu and the obs.

How about the theta and sigma? Both of them will be n x n matrix in multivariate
Any suggestion for the prior?

def sde(x, theta, mu, sigma):
return theta * (mu - x), sigma

import pymc3 as pm
import pymc3.distributions.timeseries as ts
with pm.Model() as wiggins_model:

volatility_theta = pm.Uniform('volatility_theta', lower=0., upper=1., testval=0.5)
volatility_mu = pm.Normal('volatility_mu', mu=-5., sd=.1, testval=-5)
volatility_sigma = pm.Uniform('volatility_sigma', lower=0.001, upper=0.2, testval=0.05)

volatility = ts.EulerMaruyama('volatility',
1.0,
sde,
(volatility_theta, volatility_mu, volatility_sigma),
shape=len(returns),
testval=np.ones_like(returns["change"]))

pm.Normal('obs', mu=0., sd=pm.math.exp(volatility), observed=returns["change"].values)

trace = pm.sample(4000, cores=8, chains=2, tune=3000, random_seed=42)

Thank you very much.

Dear @lucianopaz ,

I read about your post regarding PyMC3 shape handling | Luciano Paz.

I have a question for my projects regarding the shape for prior.
According to Ornstein-Uhlenbeck equation,
dX(t) = theta*(mu-X(t))dt + sigma*dW

I wish to transform from univariate to multivariate.
So,
theta = n x n matrix
mu = n dimensional vector
sigma = n x n matrix

May I know how to set the shape for theta and sigma in prior so that they are in n x n matrix?

Thank you very much