library(TSA)
data(lynx)
plot(log(lynx))
pacf(log(lynx))
arma.lynx=arima(log(lynx),order=c(3,0,0))
arma.lynx
# R-squared
1-.2669/1.653
data(sunspots)
?sunspots
plot(sunspots)
?aggregate
frequency(sunspots)
sunspots.y=aggregate(sunspots,nfrequency=1)
frequency(sunspots.y)
plot(sunspots.y)
sunspots.y1=aggregate(sunspots,nfrequency=1,FUN=median)*12
lines(sunspots.y1,col='brown')
hist(sunspots.y)
boxcox=BoxCox.ar(sunspots.y+min(sunspots.y[sunspots.y!=0]))
boxcox
# MLE suggest square root transformation
plot(sqrt(sunspots.y))
acf(sqrt(sunspots.y),lag.max=60)
pacf(sqrt(sunspots.y),lag.max=60)
eacf(sqrt(sunspots.y))
arima.s=arima(sqrt(sunspots.y),order=c(3,0,0))
arima.s
polyroot(c(1,-arima.s$coef[-4]))
Mod(polyroot(c(1,-arima.s$coef[-4])))
# stationary
Arg(polyroot(c(1,-arima.s$coef[-4])))
theta from the Arg = theta=0.5698923
period=2*pi/theta
period
coefm.cond.norm=arima.boot(arima.s,cond.boot=F,is.normal=F,B=1000,init=sqrt(sunspots.y))
any(is.na(c(1,2,NA)))
any(is.na(coefm.cond.norm))
hist(coefm.cond.norm[,1],main='AR(1) coefficient')
quantile(coefm.cond.norm[,1],c(.025,0.975)) # bootstrap 95% CI of AR(1)
arima.s
1.3217-1.96*0.0648
1.3217+1.96*0.0648
apply(coefm.cond.norm,2,function(x){quantile(x,c(0.025,0.975))})
hist(coefm.cond.norm[,5],main='noise variane')

period.replace=apply(coefm.cond.norm,1,function(x){
roots=polyroot(c(1,-x[1:3]))
# find the complex root with smalles magnitude
min1=1.e+9
rootc=NA
for (root in roots) {
if( abs(Im(root))<1e-10) next
if (Mod(root)< min1) {min1=Mod(root); rootc=root}
}
if(is.na(rootc)) period=NA else period=2*pi/abs(Arg(rootc))
period
})

hist(period.replace)

#95% CI of the period
quantile(period.replace,c(.025,.975))

set.seed(12345)
phi=.5
acf1.resi=rep(0,1000)
acf1.innov=rep(0,1000)
for (j in 1:1000){

Y=rep(0,100)
e=rnorm(n=100)
Y[1]=rnorm(n=1,mean=0,sd=1/sqrt(1-phi^2))
for (i in 2:100) Y[i]=phi*Y[i-1]+e[i]
#plot(ts(Y))
#plot(ts(e))
ar1=arima(Y,order=c(1,0,0))
lines(ar1$residuals,col='blue')
acf1.resi[j]=acf(ar1$residuals,plot=F)$acf[1]
acf1.innov[j]=acf(e,plot=F)$acf[1]

}

hist(acf1.resi)
hist(acf1.innov)
acf(ar1$residuals)