#Fitting ARIMA Models to Wolfer sunspot numbers, 1770-1869


wolfer <- scan("https://edoras.sdsu.edu/~babailey/stat496/sun_spot.txt", skip=1)

wolfer <- ts(wolfer, start=1770)
sqrtwolfer <- ts(sqrt(wolfer), start=1770)

# Graph the data.
par(mfrow=c(3,2))
plot(wolfer)
title("Wolfer Sunspot Numbers")
plot(sqrtwolfer)
title("Square Root Transform")
acf(wolfer)
acf(sqrtwolfer)
acf(wolfer, type="partial")
acf(sqrtwolfer, type="partial")


#Fit AR(2) 

fit1 <- arima(sqrtwolfer, order=c(2,0,0))
fit2 <- arima(wolfer, order=c(2,0,0))

# forecast 
pred <- predict(fit1,n.ahead=20)

# graph the forecasts 
par(mfrow=c(1,1))
plot(c(wolfer,(pred$pred)^2), ylim=c(0,180), type="n")
lines(1:100, wolfer, lty=1)
#points(1:100, wolfer, pch=1)
lines(101:120, (pred$pred + 2*pred$se)^2, lty=2)
lines(101:120, (pred$pred - 2*pred$se)^2, lty=2)
points(101:120, (pred$pred)^2, pch=3)
title("95% Forecast intervals")