# coefficient de Lyapounov
gamma.garch11<- function (alpha,beta) {
gam<-integrate(function (x) log(alpha*x^2+beta)*dnorm(x), -Inf, Inf)$value
gam
}
# trouve le alpha tel que gamma.garch11(alpha,beta)=0
objf <- function(alpha, beta)gamma.garch11(alpha,beta)^2     
solve.gamma.0<- function(beta,alphainit,alphamin,alphamax)
     {
valinit<-alphainit
res <- nlminb(valinit,objf, lower=alphamin,
    upper=alphamax, beta=beta)
alpha<-res$par
alpha
}
#alphainit<-0.8; alphamin<-0.5; alphamax<-2; beta<-0.5

npt<-200
abs1<-rep(0,npt)
ord1<-seq(from=0,to=1,length.out=npt)
mom<-integrate(function (x) log(x^2)*dnorm(x), -Inf, Inf)$value
alphamax<-exp(-mom)
abs2<-abs1
for(i in 1:npt){
abs1[i]<-solve.gamma.0(ord1[i],(1-ord1[i]+alphamax)/2,1-ord1[i],alphamax)
abs2[i]<-(-ord1[i]+sqrt(3-2*ord1[i]^2))/3
}
oldpar<-par(las=1)
plot(c(0.13,alphamax),c(0.035,1),type="n",
xlab=expression(alpha[1]),ylab=expression(beta[1]),
xaxp=c(0,4,10),yaxp=c(0,1,5))
segments(0,1,1,0)
lines(abs1,ord1)
lines(abs2,ord1)
text(0.2,0.2,"1")
text(0.65,0.2,"2")
text(1.2,0.2,"3")
text(2.,0.5,"4")
par(oldpar)