Я прочитал следующую статью (http://www3.stat.sinica.edu.tw/statistica/oldpdf/A10n416.pdf) где они моделируют дисперсионно-ковариационную матрицу Σ как:

Σ = diag (S) * R * diag (S) (уравнение 1 в статье)

S - вектор стандартных отклонений k × 1, diag (S) - диагональная матрица с диагональными элементами S, а R - матрица корреляции k × k.

Как я могу реализовать это с помощью PyMC?

Вот некоторый исходный код, который я написал:

import numpy as np
import pandas as pd
import pymc as pm

k=3
prior_mu=np.ones(k)
prior_var=np.eye(k)
prior_corr=np.eye(k)
prior_cov=prior_var*prior_corr*prior_var

post_mu = pm.Normal("returns",prior_mu,1,size=k)
post_var=pm.Lognormal("variance",np.diag(prior_var),1,size=k)
post_corr_inv=pm.Wishart("inv_corr",n_obs,np.linalg.inv(prior_corr))


post_cov_matrix_inv = ???

muVector=[10,5,-2]
varMatrix=np.diag([10,20,10])
corrMatrix=np.matrix([[1,.2,0],[.2,1,0],[0,0,1]])
cov_matrix=varMatrix*corrMatrix*varMatrix

n_obs=10000
x=np.random.multivariate_normal(muVector,cov_matrix,n_obs)
obs = pm.MvNormal( "observed returns", post_mu, post_cov_matrix_inv, observed = True, value = x )

model = pm.Model( [obs, post_mu, post_cov_matrix_inv] )
mcmc = pm.MCMC()

mcmc.sample( 5000, 2000, 3 )

Спасибо

[редактировать]

Я думаю, что это можно сделать с помощью следующего:

@pm.deterministic
def post_cov_matrix_inv(post_sdev=post_sdev,post_corr_inv=post_corr_inv):
    return np.diag(post_sdev)*post_corr_inv*np.diag(post_sdev)