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

[Alex Bird]

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]
without adding any extra noise
  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.

Thanks in advance!

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')