set.seed(4532)
N<-10000
alpha<-2.5
sigal<-1
X<-rep(0,N)
X[1]<-0
for(istep in 2:N){
	x<-X[istep-1]
	Z<-abs(rnorm(1,x,sigal))
	al<-min(1,dgamma(Z,alpha)/dgamma(x,alpha))
	u<-runif(1)
	if(u < al){
		X[istep]<-Z
	}else{
		X[istep]<-x
	}
}

postscript("MH-Gamma1.ps")
plot(X[1:200],type="s",xlab="t",ylab="X",lwd=2)
dev.off()

postscript("MH-Gamma2.ps")
hist(X,br=seq(0,12,by=0.25),xlab="t",main=" ",col="black",border="white",ylim=range(0,800))
xv<-seq(0,12,by=0.01)
yv<-dgamma(xv,alpha)
lines(xv,yv*N*0.25,col="red",lwd=2)
dev.off()


postscript("MH-Gamma-ACF1.ps")
acf(X,lag.max=100,main=" ")
dev.off()


#########################################

set.seed(4532)
N<-10000
alpha<-2.5
sigal<-0.1
X<-rep(0,N)
X[1]<-0
for(istep in 2:N){
	x<-X[istep-1]
	Z<-abs(rnorm(1,x,sigal))
	al<-min(1,dgamma(Z,alpha)/dgamma(x,alpha))
	u<-runif(1)
	if(u < al){
		X[istep]<-Z
	}else{
		X[istep]<-x
	}
}

postscript("MH-Gamma3.ps")
plot(X[1:200],type="s",xlab="t",ylab="X",lwd=2)
dev.off()

postscript("MH-Gamma4.ps")
hist(X,br=seq(0,12,by=0.25),xlab="t",main=" ",col="black",border="white",ylim=range(0,800))
xv<-seq(0,12,by=0.01)
yv<-dgamma(xv,alpha)
lines(xv,yv*N*0.25,col="red",lwd=2)
dev.off()

postscript("MH-Gamma-ACF2.ps")
acf(X,lag.max=100,main=" ")
dev.off()

#########################################

set.seed(48532)
N<-10000
alpha<-2.5
sigal<-3
X<-rep(0,N)
X[1]<-0
for(istep in 2:N){
	x<-X[istep-1]
	Z<-abs(rnorm(1,x,sigal))
	al<-min(1,dgamma(Z,alpha)/dgamma(x,alpha))
	u<-runif(1)
	if(u < al){
		X[istep]<-Z
	}else{
		X[istep]<-x
	}
}

postscript("MH-Gamma5.ps")
plot(X[1:200],type="s",xlab="t",ylab="X",lwd=2)
dev.off()

postscript("MH-Gamma6.ps")
hist(X,br=seq(0,12,by=0.25),xlab="t",main=" ",col="black",border="white",ylim=range(0,800))
xv<-seq(0,12,by=0.01)
yv<-dgamma(xv,alpha)
lines(xv,yv*N*0.25,col="red",lwd=2)
dev.off()

postscript("MH-Gamma-ACF3.ps")
acf(X,lag.max=100,main=" ")
dev.off()

#########################################

set.seed(48532)
N<-20000
alpha<-2.5
sigal<-0.1
X<-rep(0,N)
X[1]<-20
for(istep in 2:N){
	x<-X[istep-1]
	Z<-abs(rnorm(1,x,sigal))
	al<-min(1,dgamma(Z,alpha)/dgamma(x,alpha))
	u<-runif(1)
	if(u < al){
		X[istep]<-Z
	}else{
		X[istep]<-x
	}
}


postscript("MH-Gamma7.ps")
plot(X[1:1000],type="s",xlab="t",ylab="X",lwd=2,ylim=range(0,25))
abline(h=qgamma(c(0.025,0.975),alpha),lty=2,lwd=2,col="red")
dev.off()

postscript("MH-Gamma8.ps")
plot(X[1:10000],type="s",xlab="t",ylab="X",lwd=2,ylim=range(0,25))
abline(h=qgamma(c(0.025,0.975),alpha),lty=2,lwd=2,col="red")
dev.off()

postscript("MH-Gamma9.ps")
plot(X[1:N],type="s",xlab="t",ylab="X",lwd=2,ylim=range(0,25))
abline(h=qgamma(c(0.025,0.975),alpha),lty=2,lwd=2,col="red")
dev.off()

postscript("MH-Gamma10.ps")
plot(15001:N,X[15001:N],type="s",xlab="t",ylab="X",lwd=2,ylim=range(0,10))
abline(h=qgamma(c(0.025,0.975),alpha),lty=2,lwd=2,col="red")
dev.off()