[R] Really need help here

[Tong Wang]

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

}