> library(tseries)
> x <- arima.sim(n = 63, list(ar = c(0.8897, -0.4858), ma = c(-0.2279, 0.2488)),
+           sd = sqrt(0.1796))
> arma(x)

Call:
arma(x = x)

Coefficient(s):
      ar1        ma1  intercept  
  0.32929    0.31320    0.01793

如果产生MA或AR也还有其他的方法，比如
w=rnorm(500,0,1)                                              # 500 N(0,1) variates
v=filter(w,sides=2,rep(1/3,3))                          # moving average v(t)=(w(t+1)+w(t)+w(t-1))/3
par(mfrow=c(2,1))
plot.ts(w,main="white noise")
plot.ts(v,ylim=c(-3,3), main="moving average")
以及：
w = rnorm(550,0,1)                                      # 50 extra to avoid startup problems
x = filter(w, filter=c(1,-.9), method="recursive")[-(1:50)] #x[3]==x[2]-0.9*x[1]+w[53]
plot.ts(x, main="autoregression")

用R来确定arma的p，q值呢

library(tseries)
x <- arima.sim(n = 63, list(ar = c(0.8897, -0.4858), ma = c(-0.2279, 0.2488)),
          sd = sqrt(0.1796))
ARMA.Order <- function(x,max) {
  armax <- mamax <- max
  require(foreach)
  armafit <- function(x,a,b){
    res.arma <- arma(x,order=c(a,b))
    summary(res.arma)$aic
  }
  res.aic <- foreach(a=1:armax,.combine="rbind") %:% foreach(b=1:mamax,.combine="rbind") %do% {
    aic <- armafit(x,a,b)
    cbind(aic,a,b)
  }
  out <- res.aic[which.max(res.aic[,1]),]
  cat("最大 AIC 是",out[1],"\n AR 阶是",out[2],"\n MA 阶是",out[3],"\n")
  a <- out[2];b <- out[3]
  cat("a =",a,"b =",b,"\n")
  arma.maxaic <- arma(x,order=c(a,b))
  return(arma.maxaic)
}
ARMA.Order(x,max=3)



最大 AIC 是 63.0571
AR 阶是 1
MA 阶是 1
a = 1 b = 1

Call:
arma(x = x, order = c(a, b))

Coefficient(s):
      ar1        ma1  intercept  
  0.41560    0.20058   -0.07382  

警告信息：
1: In arma(x, order = c(a, b)) : Hessian negative-semidefinite
2: In sqrt(diag(object$vcov)) : 产生了NaNs
3: In sqrt(diag(object$vcov)) : 产生了NaNs

学习了！不错！