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]]
______________________________________________
[email protected] 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.
