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')
_______________________________________________
[email protected] mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-finance
-- Subscriber-posting only. If you want to post, subscribe first.
-- Also note that this is not the r-help list where general R questions should
go.
