#Sunspot Numbers 1770-1889, Yearly.

sund<- scan("sunspot.dat")         #evtl: sunspot.txt
dim(sund)<- c(120)

postscript(file="SUNsp.ps",height=16,width=8,horizontal=F)
par(mfrow=c(2,1))

mean(sund); sqrt(var(sund))
sunts<- ts(sund,start=1770)              #start im ts plot "xlab" relevant

ts.plot(sunts,xlab="Year",ylab="Sunspot number")
title(main="Yearly Sunspot Numbers, 1770 - 1889",cex=0.9) 
abline(h=c(50,100),lty=2)
#abline(v=1778)

pund<- pmax(sund,1)
lsunts<- ts(log10(pund),start=1770)
 
ts.plot(lsunts,xlab="Year",ylab="log10(Sunspot number)")
abline(h=c(0.5,1,1.5),lty=2)

# Erzeugen der Residuenzeitreihe e[t] aus sunts[1]...sunts[120]
sunar<- ar(sunts,order=2)
# Parameter alpha_1, alpha_2, theta_0 des AR(2)-Modells
al1<- sunar$ar[1]; al2<- sunar$ar[2] ; the<- (1 - al1 -al2)*mean(sunts)
c(al1,al2,the)
et<- 1:120
for (t in 3:120)
{et[t]<- sunts[t]-(al1 * sunts[t-1] + al2 * sunts[t-2] + the)}
acf(et[3:120],lag.max=12,type="corr")

dev.off()

rm(sund,pund,sunts,lsunts,al1,al2,the,et)