#generate the realization Xt=Z_{t}+theta*Z_{t-1} theta<- -0.8 n<-200 Zt<-vector('numeric',n+1) Zt<-rnorm(n+1) Xt<-vector('numeric',n) for(t in 1:n) { Xt[t]<-Zt[t+1]+theta*Zt[t] } #compute sample mean xbar<-mean(Xt) #compute sample autocorrelation function at lags 0, 1, ..., n/5 ACF<-vector('numeric',n/5+1) rhat0<-sum((Xt-xbar)^2)/n ACF[1]<-1 for(i in 1:(n/5)) { ACF[i+1]<-(sum((Xt[(i+1):n]-xbar)*(Xt[1:(n-i)]-xbar))/n)/rhat0 } # plot realization and ACF par(mfrow=c(2,1)) plot(1:n,Xt,type='l',xlab="index",ylab="x_t",main="realization of 200 MA(1)") plot(0:(n/5),ACF,type='h',xlab="Lag",ylab="ACF",main="Autocorrelation",ylim=c(-0.8,1)) abline(h=0) abline(h=-1.96*sqrt(1+2*ACF[2]^2)/sqrt(n),lty=2,col="blue") abline(h=1.96*sqrt(1+2*ACF[2]^2)/sqrt(n),lty=2,col="blue") points(c(0,1),c(1,theta/(1+theta^2)),pch=21) # Generate an AR(2) using N(0,1) as forcing n<-250 Zt<-vector('numeric',n+1) Zt<-rnorm(n+1) theta1 <- 0.5 theta2 <- -1./3. Xt[1] <- Zt[1] Xt[2] <- Zt[2] - theta1*Xt[1] Xt[3] <- Zt[3] - theta1*Xt[2] - theta2*Xt[1] for(t in 4:250) { Xt[t]<- Zt[t] - theta1*Xt[t-1] - theta2*Xt[t-2] } # now only use the last 200 points; 'burn-in' plot(seq(1:200),Xt[51:250],type="l",lwd=2, xlab="Time (t)", ylab="X_t") library("stats") acf(Xt[51:250],lag.max=50,type="covariance",plot=TRUE)