Hi there,
I had a serious problem here . Consider the following Bayesian
model(discretized variance gamma):
#Likelihood
J[i]<-lambda*G[i]+sigma*sqrt(G[i])*rnorm(0,1)
G[i]<-rgamma(1/nu,1/nu)
#Prior:
nu<-rinvgamma(m,M)
# Parameters
lambda=-.04 ; sigam=.38; nu=6.48; m=10,M=10 ; T=5000 (length of data)
An author claimed that he got posterior distribution of nu with standard
deviation .0965
But I could only get posterior sd=.84 with the same setting. I have been working
on my code for quite a while
, but still don't see what is wrong with it. I would really appreciate it
if any body could fit this model and let me
know if you can get it and how ? I attached my R code with data generating and
posterior simulating functions,
I used the slice sampler to do it , you could use your own version of slice
sampler or other generic samplers.
Again , I would really appreciate any help that could save me from the pain.
vg <- function(T=5000,nu=6.48,n.iter=1000,beta=-.04,sig=.38,m=10,M=10){
#------------------- Data Generating -------------------------#
G <- rgamma(T,1/nu,1/nu)
J <- rnorm(T,beta*G,sig*sqrt(G))
#------------------- Gibbs Sampler ---------------------------#
nu.stor <- rep(NA,n.iter)
for(i in 1:n.iter){
G <- MCslice1D(dgamma(x,1/nu,1/nu,log=TRUE)+dnorm(J, beta*x,
sig*sqrt(x),log=TRUE),w=20,m=10,x0=G)
nu <-
MCslice1D(log(dinvgamma(x,m,M))-T*log(x)*1/x-T*log(gamma(1/x))+(1/x-1)*sum(log(G))-sum(G)/x,
w=3,m=12, x0=nu)
nu.stor[i] <- nu
}
return(list(nu=nu.stor))
}
