Estimate Ornstein-Uhlenbeck Parameters with Euler Maruyama


#1

Hello,
I want estimate the parameters of 2D Ornstein Uhlenbeck Process, my data are the following: Schermata%20del%202018-06-11%2022-35-41
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:


Schermata%20del%202018-06-11%2023-05-15
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():
    advi = pm.ADVI()

    approx = pm.fit(20000, method='advi', obj_n_mc=20 ,obj_optimizer = pm.adagrad(learning_rate=0.9))
    trace = approx.sample(draws=500)
    ppc_trace = pm.sample_ppc(trace, model=model)

    return advi, trace ,ppc_trace


Schermata%20del%202018-06-12%2017-26-34
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.