# clean up workspace
rm(list=ls(all=TRUE))
set.seed(123456789)

#simulation and graphing
y=arima.sim(n=1000,list(order=c(1,0,1),ar=0.8,ma=0.8),n.start=100)
y1=arima.sim(list(ar=0.8),n=100,n.start=100)
y2=arima.sim(list(ar=-0.8),n=100,n.start=100)
e=arima.sim(list(order=c(0,0,0)),n=100,n.start=100)
y2=arima.sim(list(order=c(2,0,2),ar=c(1,-0.16),ma=c(1,0.24)),n=100,n.start=100)

op=par(ask=F)
par(mfrow=c(3,1))
ts.plot(y1,main="y_t=0.8y_[t-1]+e_t")
ts.plot(y2,main="y_t=-0.8y_[t-1]+e_t")
ts.plot(e,main="whote noise")
par(op)


#Asymptotic distribution
nobs=500
nrep=100
b0=0.8
b=numeric(nrep)
for (i in 1:nrep)
{  y=arima.sim(list(ar=0.8),n=nobs,n.start=100)
   b[i]=ar(y,order=1,demean=F)$ar
 }  
 
 bt=sqrt(nobs)*(b-b0)
 hist(bt)  
 
 summary(bt)  
 qqnorm(bt)
 
 mmARMA11 = function(z) {
 # y = b y[-1] + e + th*e[-1]
 n = length(z);
 acvz=acf(z,type="covariance",plot=F)$acf
 g0=acvz[1]
 g1=acvz[2]
 g2=acvz[3]
 b=g2/g1
 k=(g0-b*g1)/(g1-b*g0)
 
 th1=(k+sqrt(k*k-4*(1-b)))/2
 th2=(k-sqrt(k*k-4*(1-b)))/2
 th=th1*(k<0)+th2*(k>=0)
 list(b=b,th=th)
 }
 
 b.mm=numeric(nrep)
 th.mm=numeric(nrep)
 b.mle=numeric(nrep)
 th.mle=numeric(nrep)
 p=numeric(nrep)
 for (i in 1:nrep)
 {
 z=arima.sim(n=nobs,list(order=c(1,0,1),ar=0.8,ma=0.6),n.start=100)
 mm.est= mmARMA11(z)
 b.mm[i]= mm.est$b
 th.mm[i]=mm.est$th
 mle.est= arima(z,order=c(1,0,1))
 b.mle[i]= mle.est$coef[1]
 th.mle[i]= mle.est$coef[2]
 p[i]=Box.test(mle.est$resid)$p.value
 
 }
 b0=0.8
 th0=0.6
 effb.mm=var(b.mm-b0)
 effb.mle=var(b.mle-b0)
 effth.mm=var(th.mm-th0)
 effth.mle=var(th.mle-th0)
 sum(p>0.05)