# load data
data(sunspot,package=ts)

# calculate periodogram
n <- length(sunspot.year)
ind <- 1:(ceiling(n/2)-1)
y <- (abs(fft(sunspot.year))^2/n)[-1][ind]
x <- ind/n

# fit the spectrum
sunspot.fit <- gssanova(y~x,family="Gamma",method="u",varht=1)
xx <- seq(0,.5,length=101)
est <- predict(sunspot.fit,data.frame(x=xx),se=TRUE)

# plot the fit
plot(x,y,log="y",xlab="frequency",ylab="spectrum")
lines(xx,exp(est$fit),col=4)
lines(xx,exp(est$fit+1.96*est$se),col=5)
lines(xx,exp(est$fit-1.96*est$se),col=5)