I’m still trying to fix the last problem I encountered. I got the idea that I should try to code the likelihood function in an entirely different method as to allow me to check my main code. (BTW, I don’t think what I suspect last time is the problem).

 Since the likelihood function is essentially a multi-dimensional definite integral, I thought I’d try naive Monte Carlo simulation. The actually coding went by quickly, but it didn’t work. It soon became apparent that everything is badly under-flowing. The smallest value R can keep track of is in the region of -10^{-309}, which one can quickly reach when you are multiplying thousands of numbers between 0 and 1 together. The standard tactic to get around it is to log everything. Products become sums and ludicrously small numbers become merely small.

So you end up with code like this:

logdensitySV <- function(y,h,theta){

gamma0 = theta[1];
gamma1 = theta[2];
sigma = theta[3];

h0 = gamma0/(1-gamma1);
T = length(h);

transition = vector("numeric",T);
transition[1] = dnorm((h[1]-gamma0-gamma1*h0)/sigma);
transition[2:T] = dnorm((h[2:T]-gamma0-gamma1*h[1:T-1])/sigma);

condprob = dnorm(y/exp(.5*h));

sum(log(transition))+sum(log(condprob));
}

Which is fine, except when these guys are underflowing too! To consistently get these functions to return a usable value, I have to add about 100 000 to the log value! I need to think carefully about whether doing this will cause an overflow down the line…

Advertisement