Let’s assume we have two instantaneously correlated processes: S_1, S_2 both driven by Geometric Brownian Motion such that

\left( {\begin{array}{*{20}{c}}
{d{S_1}(t)}\\
{d{S_2}(t)}
\end{array}} \right) = \left( {\begin{array}{*{20}{c}}
{{\mu _1}{S_1}(t)dt}\\
{{\mu _2}{S_2}(t)dt}
\end{array}} \right) + \left( {\begin{array}{*{20}{c}}
{{\sigma _1}{S_1}(t)d{W_1}(t)}\\
{{\sigma _2}{S_2}(t)\left[ {\rho d{W_1}(t) + \sqrt {1 - {\rho ^2}} d{W_2}(t)} \right]}
\end{array}} \right)

This can alternatively be written as a vector-valued process

\begin{equation}
d{\bf{S}}(t) = \mu (t,{\bf{S}}(t))dt + \sigma (t,{\bf{S}}(t))dW(t)
\end{equation}

and we can easily see that

\begin{equation}
d{\bf{S}}(t)|{\bf{S}}(t)\sim N\left( {\left( {\begin{array}{*{20}{c}}
{{\mu _1}{S_1}(t)dt}\\
{{\mu _2}{S_2}(t)dt}
\end{array}} \right),\left( {\begin{array}{*{20}{c}}
{S_1^2(t)\sigma _1^2dt}&{\rho {S_1}(t){S_2}(t){\sigma _1}{\sigma _2}dt}\\
{\rho {S_1}(t){S_2}(t){\sigma _1}{\sigma _2}dt}&{S_2^2(t)\sigma _2^2dt}
\end{array}} \right)} \right)
\end{equation}

If we observe daily values for S_1,S_2 (and their increments dS_1, dS_2), then I can easily use pymc to estimate their parameters (\mu_1, \sigma_1) and (\mu_2, \sigma_2), e.g

```
dt = 1/252
with pm.Model() as model:
mu_param = pm.Uniform('mu_param', -0.5, 0.5)
sigma_param = pm.Uniform('sigma_param', 0.0001, 0.3)
dS = pm.Normal('dS', mu = (data['seriesA'].values*mu_param)*dt,
sigma = data['seriesA'].values*sigma_param*np.sqrt(dt),
observed = data['dS_seriesA'].values)
trace = pm.sample(50000)
```

However, I’d like to know how to use pymc to estimate the correlation parameter \rho from the joint dynamics. I’ve explored MvNormal but I can’t find the right way how to set up the estimation task.

My working ipynb+data is on Github: https://github.com/marekkolman/pymc_demo/blob/master/estimate_params.ipynb

Can someone pls give me a hint how to estimate the parameters in the multidimensional model?