spectrum(sunspots,spans=10,xlim=c(0,1),main="10:Notquite)
abline(v=1:3/10,lty=3)
2007年2月24日星期六
程序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)
{
# 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)
2007年2月20日星期二
程序1
op<-par(mfrow=c(4,2),mar=c(2,4,4,2))
n <- 200
for(i in 1:4) {
x <- NULL
while(is.null(x)) {
model<-list(ar=rnorm(1))
try(x<-arima.sim(model,n))
}
acf(x, main=paste( "ARMA(1,0)", "AR:", round(model$ar,digits=1) ))
points(0:50, ARMAacf(ar=model$ar,lag.max=50), col='red')
pacf(x,main="")
points(1:50,
ARMAacf(ar=model$ar,lag.max=50,pacf=T),
col='red')
}
par(op)
n <- 200
for(i in 1:4) {
x <- NULL
while(is.null(x)) {
model<-list(ar=rnorm(1))
try(x<-arima.sim(model,n))
}
acf(x, main=paste( "ARMA(1,0)", "AR:", round(model$ar,digits=1) ))
points(0:50, ARMAacf(ar=model$ar,lag.max=50), col='red')
pacf(x,main="")
points(1:50,
ARMAacf(ar=model$ar,lag.max=50,pacf=T),
col='red')
}
par(op)
订阅:
博文 (Atom)