程序2

sym.poly<-function(z,k)
{

# Sum of the products of k
# distinct elements of the vector z

if(k==0) {r <- 1}
else if(k==1)
{r <- sum(z)}
else {
r <- 0
for(i in 1:length(z)) {r <- r + z[i]*sym.poly(z[-i],k-1)}
r <- r/k
# Each term appeared k times

}
r
}

sym.poly(c(1,2,3),1)
#6
sym.poly(c(1,2,3),2)
#11
sym.poly(c(1,2,3),3)
#6

roots.to.poly<-function(z)
{
n <- length(z)
p <- rep(1,n)
for(k in 1:n) {p[n-k+1]<-(-1)^k*sym.poly(z,k)}
p <- c(p,1)
p
}

roots.to.poly(c(1,2))
#2-31

round(Re(polyroot(roots.to.poly(c(1,2,3)))),digits=1)
# After this interlude, we can finally construct an MA process and one of its inverse

sn <- 200
k <- 3

ma<-runif(k,-1,1)
# The roots

z <- polyroot(c(1,-ma))
# The inverse of the roots

zi<-1/z
# The polynomial

p <- roots.to.poly(zi)
# The result should be real, but because of rounding errors, it is not.

p <- Re(p)
# We want the constant term to be 1.

p <- p/p[1]
mai<--p[-1]
op<-par(mfrow=c(4,1),mar=c(2,4,1,2)+.1)
x <- arima.sim(list(ma=ma),n)
plot(x,xlab="")
acf(x,main="",xlab="")
lines(0:n,ARMAacf(ma=ma,lag.max=n),lty=2,lwd=3,col='red')
x <- arima.sim(list(ma=mai),n)
plot(x,xlab="")
acf(x,main="",xlab="")
lines(0:n,ARMAacf(ma=mai,lag.max=n),lty=2,lwd=3,col='red')
par(op)