0

I read the following paper(http://www3.stat.sinica.edu.tw/statistica/oldpdf/A10n416.pdf) where they model the variance-covariance matrix Σ as:

Σ = diag(S)*R*diag(S) (Equation 1 in the paper)

S is the k×1 vector of standard deviations, diag(S) is the diagonal matrix with diagonal elements S, and R is the k×k correlation matrix.

How can I implement this using PyMC ?

Here is some initial code I wrote:

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 )

Thanks

[edit]

I think that can be done using the following:

@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)
akhil
  • 839
  • 3
  • 8
  • 15
  • Please expand on what you mean by "model." That word has many meanings in statistics and science, none of which would seem to apply here. Are you perhaps asking how to *decompose* a covariance matrix into this form? If your question is only about coding an algorithm in PyMC, then please let us know so we can migrate it to the SO community. – whuber Feb 10 '14 at 23:37
  • My question is only about the implementation in PyMC. – akhil Feb 10 '14 at 23:40
  • I think that can be done using the following: @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) – akhil Feb 11 '14 at 00:50

1 Answers1

1

Here is the solution for the benefit of someone who stumbles onto this post:

p=3
prior_mu=np.ones(p)
prior_sdev=np.ones(p)
prior_corr_inv=np.eye(p)


muVector=[10,5,1]
sdevVector=[3,5,10]
corrMatrix=np.matrix([[1,0,-.1],[0,1,.5],[-.1,.5,1]])
cov_matrix=np.diag(sdevVector)*corrMatrix*np.diag(sdevVector)

n_obs=2000
x=np.random.multivariate_normal(muVector,cov_matrix,n_obs)

prior_cov=np.diag(prior_sdev)*np.linalg.inv(prior_corr_inv)*np.diag(prior_sdev)

post_mu = pm.Normal("returns",prior_mu,1,size=p)
post_sdev=pm.Lognormal("sdev",prior_sdev,1,size=p)
post_corr_inv=pm.Wishart("inv_corr",n_obs,prior_corr_inv)

#post_cov_matrix_inv = pm.Wishart("inv_cov_matrix",n_obs,np.linalg.inv(prior_cov))
@pm.deterministic
def post_cov_matrix_inv(post_sdev=post_sdev,post_corr_inv=post_corr_inv,nobs=n_obs):
    post_sdev_inv=(post_sdev)**-1
    return np.diag(post_sdev_inv)*cov2corr(post_corr_inv/nobs)*np.diag(post_sdev_inv)

obs = pm.MvNormal( "observed returns", post_mu, post_cov_matrix_inv, observed = True, value = x )

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

mcmc.sample( 25000, 15000, 1,progress_bar=False )
akhil
  • 839
  • 3
  • 8
  • 15