The urgency and the vague description of your problem strongly suggest
that this is homework. This list is not for homework---see the posting
guide at the bottom of every message. Nonetheless since I know this
problem reasonably well I will offer some comments.
QRMlib is a package created to accompany a book. If you read that book
you would see that it fits the generalized hyperbolic to data using the
EM algorithm. If you have QRMlib you have an implementation of the EM
algorithm.
Also why write code to simulate from the generalized hyperbolic (y in
your simulation function below) when you have QRMlib and ghyp, both of
which have functions for simulating from the generalized hyperbolic?
Your code is pretty difficult to follow, with random indenting and zero
comments. The structure of the iteration is totally confused as well.
Not too many marks if you handed something like this in to me to grade.
David Scott
On 21/09/2010 5:32 p.m., snes1...@hotmail.com wrote:
I created a EM algorithm for Generalized hyperbolic distribution.
I want to estimate mutheldaplus, sigmatheldaplus, betasigmaplus in my code.
After getting use these value , then my iteration have to be begin of this code.
But I can not to do iteration part.
Can you help me use my code and get iteration ?
Do know any useful code for EM algorithm for Generalized Hyperbolic
library(QRMlib)
library(ghyp)
############ simulation part
simulation<-function(n,lambda,mu,thelda,gamma,sigma,beta){
set.seed(235)
chi<-thelda^2
psi<-gamma^2
W<- rGIG(n, lambda, chi, psi);
Z<- rnorm(n,0,1);
y<-mu + beta * W + sqrt(W) * Z *gamma;
for (i in 1:n){
theldastar<-rep(0,n)
zi<-rep(0,n)
ti<-rep(0,n)
muthelda<-mu
gammathelda<-thelda*gamma
sigmathelda<-(thelda^2)*sigma
betathelda<-(thelda^2)*sigma*beta
lambdastar<-lambda-0.5
theldastar[i]<-sqrt(1+((y[i]-muthelda)/sigmathelda)^2)
gammastar<-sqrt((gammathelda^2)+((betathelda/sigmathelda)^2))
klambda1<-besselM3(lambdastar+1, x=2, logvalue=FALSE)
klambda<-besselM3(lambdastar,x=2,logvalue=FALSE)
klambda2<-besselM3(lambdastar-1,x=2,logvalue=FALSE)
zi[i]<-((theldastar[i]*klambda1*(theldastar[i]*gammastar))/(gammastar*klambda*theldastar[i]*gammastar))
ti[i]<-((gammastar*klambda2*(theldastar[i]*gammastar))/(theldastar[i]*klambda*theldastar[i]*gammastar))
zimean<-sum(zi)/n
timean<-sum(ti)/n
mutheldaplus<-(zimean*(1/n)* sum((ti[i]*y[i])-mean(y)))/((zimean*timean)-1)
betatheldaplus<- sum(y[i]- mutheldaplus)/(n*zimean)
sigmatheldaplus<-((1/n)*sum((ti[i]*((y[i]-mutheldaplus)^2))-(2*betatheldaplus*(y[i]-mutheldaplus))-((betatheldaplus^2)*zi[i])))
print(muthelda)
print(mutheldaplus)
print(betathelda)
print(betatheldaplus)
print(sigmathelda)
print(sigmatheldaplus)
return(ti)
}
}
a<-simulation(20000,-0.5,0,1,1,1,0)
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
--
_________________________________________________________________
David Scott Department of Statistics
The University of Auckland, PB 92019
Auckland 1142, NEW ZEALAND
Phone: +64 9 923 5055, or +64 9 373 7599 ext 85055
Email: d.sc...@auckland.ac.nz, Fax: +64 9 373 7018
Director of Consulting, Department of Statistics
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
