b1<-0.75 a1<-0.0 c<-0.6 eps=rnorm(1:2000) y=rep(0,2000) for (i in 2:2000){ y[i]=b1*y[i-1]+eps[i]+a1*eps[i-1]+c*y[i-1]*eps[i-1] } plot(1:1000,y[1001:2000],'l',lwd=2,ylab='',xlab='')