# [R-SIG-Finance] Parameters for equally weighted sums of bivariate normals via MCMC

Alex Bird sunduck at gmail.com
Tue Jan 18 16:15:43 CET 2011

```Hello there,

sorry in advance if this mailing list is not appropriate for my question.
I want to estimate parameters for the following model
Y[t] = 1/2 * X1[t] + 1/2 * X2[t] + e[t]
e[t] ~ NORMAL(0,0.0001)
(X1[t],X2[t]) ~ MULTIVARIATE_NORMAL(mu,varcov)
mu = c(mu1,mu2)
varcov = matrix(c(vol1,rho,rho,vol2),ncol=2)

where the vector of unknown parameters is c(mu1,mu2,rho,vol1,vol2)
and the data is given by Y[] and latent X1[],X2[].

What I tried to do:
1. simulated 1000 (X1[t],X2[t]) with known parameters from the
multivarite normal (as in the R code below).
2. got equally weighted sums of (X1[t],X2[t]) to calculate Y[t]
3. tried to estimate two models:
- first model is where I tried to estimate paramaters just for the
initial (latent) (X1[t],X2[t]) vector
- second model is where I tried to estimate the same vector of
parameters via Y[t] having actual (X1[t],X2[t]) vector as an inital
guess for the latents

I estimated the first model via MLE and MCMC (2000 iterations and
1000 burnin) and the parameters estimated are very similar.
But I was not able to fit the second model even by using 100000 iterations.
The following table shows actual and estimated parameters
mu[1]         mu[2]       rho       vol1       vol2
Actual       0.0005000000  0.0002000000 0.5600000 0.01400000 0.01300000
MLE          0.0001276392  0.0001626320 0.5973657 0.01395924 0.01315293
Model#1 MCMC 0.0001282419  0.0001742142 0.5973248 0.01396082 0.01316128
Model#2 MCMC 0.1106152802 -0.1103177229 0.4215964 0.01464065 0.01429678

What can be wrong with my model/assumptions/bugs_specification/etc?
What can be the reason why MCMC doesn't converge to the actual
parameters?

The code used to simulate the data (X1[t],X2[t] and Y[t]) and to
estimate parameters is below.

Kind regards,
Alex

'
require(dclone)
require(QRMlib)

#initial parameters
vol1<-0.014
vol2<-0.013
rho<-0.56

varcov<-diag(c(vol1,vol2)) %*% matrix(c(1,rho,rho,1),2) %*% diag(c(vol1,vol2))
mu<-c(5e-4,2e-4)

#sampling from bivariate normal
set.seed(1234)
rets<-rmnorm(1000,varcov,mu)
wrets<-1/2*rets[,1]+1/2*rets[,2]

#MLE fit
mlefit<-fit.norm(rets)

#BUGS|model#1 model/data/inits/fit
jfun1<-function(){
for(i in 1:T) {
X[i,1:N] ~ dmnorm(mu[],prec[,])
}
prec[1:N,1:N] ~ dwish(varcov[,],2)
mu[1:N] ~ dmnorm(mu.mu[],mu.prec[,])
vc[1:N,1:N]<-inverse(prec[,])
vol1<-pow(vc[1,1],1/2)
vol2<-pow(vc[2,2],1/2)
rho<-vc[1,2]/(vol1*vol2)
}
jdata1<-list(T=1000,N=2,X=rets,varcov=varcov,mu.prec=diag(1e-6,2),mu.mu=mu)
jinits1<-list(prec=inv(varcov),mu=mu)
jpara1<-c('mu[1]','mu[2]','vol1','vol2','rho')
jfit1<-jags.fit(jdata1,jpara1,jfun1,jinits1,n.chains=1,n.iter=2000)
jsummary1<-summary(window(jfit1,1000))

#BUGS|model#2 model/data/inits/fit
jfun2<-function(){
for(i in 1:T) {
Y[i] ~ dnorm(wmu[i],1E6)
wmu[i]<-1/2*X[i,1]+1/2*X[i,2]
X[i,1:N] ~ dmnorm(mu[],prec[,])
}
prec[1:N,1:N] ~ dwish(varcov[,],2)
mu[1:N] ~ dmnorm(mu.mu[],mu.prec[,])
vc[1:N,1:N]<-inverse(prec[,])
vol1<-pow(vc[1,1],1/2)
vol2<-pow(vc[2,2],1/2)
rho<-vc[1,2]/(vol1*vol2)
}
jdata2<-list(T=1000,N=2,Y=wrets,varcov=varcov,mu.prec=diag(1e-6,2),mu.mu=mu)
jinits2<-list(prec=inv(varcov),mu=mu,X=rets)
jpara2<-c('mu[1]','mu[2]','vol1','vol2','rho')
jfit2<-jags.fit(jdata2,jpara2,jfun2,jinits2,n.chains=1,n.iter=100000)
jsummary2<-summary(window(jfit2,90000))
#here I get the following warnings
#Warning messages:
#1: glm.fit: algorithm did not converge
#2: glm.fit: algorithm did not converge

#comparative statistics
cmprstats<-rbind(
c(mu,rho,vol1,vol2),
c(mlefit\$mu,mlefit\$cor[1,2],sqrt(mlefit\$Sigma[1,1]),sqrt(mlefit\$Sigma[2,2])),
as.numeric(jsummary1[[1]][,1]),
as.numeric(jsummary2[[1]][,1]))
colnames(cmprstats)<-rownames(jsummary1[[1]])
rownames(cmprstats)<-c('Actual','MLE','Model#1 MCMC','Model#2 MCMC')

```