library(TSA)
data(hare)
?hare
plot(hare)
acf(sqrt(hare))
pacf(sqrt(hare))
eacf(sqrt(hare),ma.max=7)
m0.hare=arima(sqrt(hare),order=c(2,0,0))
m0.hare
tsdiag(m0.hare,gof.lag=18)
par(mfrow=c(1,1))
qqnorm(m0.hare$residuals);qqline(m0.hare$residuals)
shapiro.test(m0.hare$residuals)

#overfit
m1.hare=arima(sqrt(hare),order=c(3,0,0))
m1.hare
m2.hare=arima(sqrt(hare),order=c(2,0,1))
m2.hare
tsdiag(m1.hare)
# forecasts on the sqrt scale
plot(m1.hare, n.ahead=25,type='o',xlab='Year',ylab='Sqrt(hare)')
abline(h=coef(m1.hare)[names(coef(m1.hare))=='intercept'])

#forecast on the original scale
plot(m1.hare, n.ahead=25,type='o',xlab='Year',ylab='hare',transform=function(x){x^2})
abline(h=(coef(m1.hare)[names(coef(m1.hare))=='intercept'])^2)

#The above approach is not correct because some prediction intervals include
#negative numbers and the square function is not monotone over the whole real line
#An ad hoc method is to first add three times the standard deviation to the
# data before doing a transformation. However, after undoing the transformation
# the prediction intervals may contain negative number, still a problem
# for interpretation. 

hist(sqrt(hare+3*var(hare)^.5))
m3.hare=arima(sqrt(hare+3*var(hare)^.5),order=c(3,0,0))
m3.hare
plot(m3.hare, n.ahead=25,type='o',xlab='Year',ylab='hare',
transform=function(x){x^2-3*var(hare)^.5})