I tried to create my first PyMC model yesterday, but it seems to be very slow and maybe there are also some convergence issues. Basically, I want to implement a hierarchical Bayesian model that originally looks like this:
where \epsilon_1, \epsilon_2, \kappa_0, P_0 and \nu_0 are known.
In order to implement it in PyMC, I changed from Inv-Wishart to Cholesky decomposition and also used a Normal distribution instead of Uniform(-\infty,\infty). This is my model:
with pm.Model() as model: # Prior for unknown scale parameter tau_sq = pm.InverseGamma('tau_sq', alpha=epsilon_1, beta=epsilon_2) # Cholesky decomposed covariance matrix using LKJCholeskyCov chol, corr, stds = pm.LKJCholeskyCov('chol_cov', n=n_assets, eta=2.0, sd_dist=pm.Exponential.dist(1.0)) # Scale the covariance matrix by tau_sq L = chol * pt.sqrt(tau_sq) # Scaling factor applied here Sigma = pm.Deterministic('Sigma', pt.dot(L, pt.transpose(L))) # Prior for the hierarchical mean mu_0 = pm.Normal('mu_0', mu=0, sigma=10) # Hierarchical prior for mu mu = pm.MvNormal('mu', mu=pt.ones(n_assets) * mu_0, cov=Sigma / kappa_0, shape=n_assets) # Likelihood obs = pm.MvNormal('obs', mu=mu, cov=Sigma, observed=returns_np) # Sampling from the posterior trace = pm.sample()
returns_np contains 200 x 50 stock returns, i.e. 50 stocks, and 200 observations per stock. After maybe 30 minutes it has not finished, and initially I get a lot of messages like this:
DEBUG:pymc.stats.convergence:Energy change in leapfrog step is too large: 658625.0689315796.
What could be the issue with my model/implementation? Or is this expected?