Estimate Ornstein-Uhlenbeck Parameters with Euler Maruyama

[Hamalle]

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?

[junpenglao]

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)

[Hamalle]

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

Any suggestion to increase this outcome? thanks in advance

[junpenglao]

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

[hstrey]

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

[Hamalle]

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.

[GMCobraz_T]

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.

[GMCobraz_T]

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