library("polynom")
source("http://www.math.pku.edu.cn/teachers/xirb/Courses/TimeSeries/code/atsa02.r")
demo.ar.roots()

a = 0.85
a=-0.85
sigma2 = 1
f = function(lambda){sigma2*(1+a^2-2*a*cos(lambda))^(-1)/(2*pi)}
x = arima.sim(1000,model=list(ar=a,sd=sqrt(sigma2)))
par(mfrow=c(2,2))
plot(x)
acf(x)
acf(x,type="partial")
spec.ar(x,order=1)
z = seq(0,0.5,by=0.001)
lines(z,f(z*2*pi)*(2*pi),type="l",col="blue")

plot(x[1:100],type="l")
#spectrum(x,span=10)

a1 = 0.75
a2 = -0.5
sigma2 = 1
f = function(lambda){
	z = complex(real=0, imaginary=lambda)
	tmp = 1-a1*exp(z)-a2*exp(2*z)
	return(sigma2*(abs(tmp))^(-2)/(2*pi))
	}
x = arima.sim(1000,model=list(ar=c(a1,a2),sd=sqrt(sigma2)))
par(mfrow=c(2,2))
plot(x)
acf(x)
acf(x,type="partial")
spec.ar(x,order=2)
z = seq(0,0.5,by=0.001)
lines(z,f(z*2*pi)*(2*pi),type="l",col="blue")

plot(x[1:100],type="l")