implement quasi-bayesian maximum likelihood estimation for normal mixtures
Helena,
First, this isn't actually a finance question (although your data or goals may
be related to finance), so it might be more appropriate to r-help than here.
I'm not familiar with the quasi-Bayesian mle, but I am familiar with Bayesian
operations in R. I'll first recommend the three books I'm aware of:
Jim Albert "Bayesian Computation with R"
Jean-Michel Marin and Christian Robert "Bayesian Core"
David Ardia "Financial Risk Management with Bayesian Estimation of GARCH Models"
The first book above has examples that are probably the closest to what you're
looking to do.
Once you've specified a mixture model, then in the Bayesian framework sampling
from that posterior distribution to get an estimate for some probability is
quite straightforward. A simple conditioning step is to use the prior
*observed* mean (or other observed moment) to further condition your mixture
model (or MCMC sampled distribution or gamma distribution or GPD distribution,
or whatever). It seems that your approach is missing this conditioning step.
Again, I'm not familiar with the quasi-Bayesian mle, but it seems to me that
with more modern fully Bayesian techniques available, you might get better
results by going all the way to a full Bayesian sampling method (with the added
bonus that there is quite a lot of R code available already).
Regards,
- Brian
Helena Richter wrote:
Hi, as you can see in the topic, I am trying to fit a normal mixture distribution with the approach suggested by Hamilton (1991). Since I couldn't find any existing packages including the quasi-bayesian mle, I have to write my own function. Unfortunately, I have absolutely no experience in doing this. If you're not familiar with the QB-MLE, I attached the formula as pdf. The idea is to extend the usual MLE with prior beliefs about the values sigma_n and sigma_b. My priors are already included in the code below. I intend to try a mixture of two normal distributions with same mean, and variances 1 and 5 as starting values. This is what I've done so far:
> R <-read.table("C:\\...\\rendite.txt", header=F)
> qbmle <- function(p, data){
mu <- mean(data); (-sum(log(p[1]/p[2]*exp(-0.5*(data-mu)^2/p[2]^2)+(1-p[1])/p[3]*exp(-0.5*data^2/p[3]^2)))-2.772*log(p[2]^2)-2.772*log(p[3]^2) - 2.772/p[2]^2 - 13.86/p[3]^2 )}
> start <-c(0.9, 1, 5) > out <- nlm(qbmle, start, data=R)
The result is: error in nlm(...): non-finite value for nlm, plus a lot of warnings, and the following output:
> out
$minimum [1] -27513.60 $estimate [1] 3.478212e+04 -2.146767e+03 -3.806269e-02 $gradient [1] -5.971628e-02 1.939856e-03 -2.946156e+02 $code [1] 5 $iterations [1] 49 So, what did I do wrong? How can I implement any non-negative constraints, and a restriction for p to be between 0 and 1? I'm sorry to bother you with such a beginners question and am very helpful for any remarks. I don't have to use the qb-mle so if you think there's a better way to do the estimation tell me. Thanks a lot, Helena
Brian G. Peterson http://braverock.com/brian/ Ph: 773-459-4973 IM: bgpbraverock