arima.sim函数具体是怎么生成一列模拟数据的？

如题，其实就是想问如果不用这个函数，自己写程序来模拟生成一列符合arma(p,q)模型的数据，是怎么生成的。我主要是对产生时间序列的逻辑有些疑惑，求大神解答。

例如：y(t)=a*y(t-1)+error(t)-b*error(t-1)

下面是我猜想的模拟过程
error=rnorm(1000,mean=0,sd=1)


for (i in 1:1000){
if i=1
y(i)=error（i）
else
y(i)=a*y(i-1)+error(i)-b*(y(i-1)-y(target))
}


我主要有两个疑问，标红的第一个地方，表示这列数据的初始值，是直接白噪声生成的么？
标红的第二个地方表示arma模型中的error(t-1)的部分，是不是像我写的这样来的？还是说应该直接写error(i-1)
R是像我这样产生模拟数据的么？

直接输入arima.sim查看代码

suimong 发表于 2015-12-20 02:15
直接输入arima.sim查看代码

不瞒您说 我已经试过了。。我实在没看懂。。。

1. > arima.sim
   
2. function (model, n, rand.gen = rnorm, innov = rand.gen(n, ...),
   
3.     n.start = NA, start.innov = rand.gen(n.start, ...), ...)
   
4. {
   
5.     if (!is.list(model))
   
6.         stop("'model' must be list")
   
7.     if (n <= 0L)
   
8.         stop("'n' must be strictly positive")
   
9.     p <- length(model$ar)
   
10.     if (p) {
    
11.         minroots <- min(Mod(polyroot(c(1, -model$ar))))
    
12.         if (minroots <= 1)
    
13.             stop("'ar' part of model is not stationary")
    
14.     }
    
15.     q <- length(model$ma)
    
16.     if (is.na(n.start))
    
17.         n.start <- p + q + ifelse(p > 0, ceiling(6/log(minroots)),
    
18.             0)
    
19.     if (n.start < p + q)
    
20.         stop("burn-in 'n.start' must be as long as 'ar + ma'")
    
21.     d <- 0
    
22.     if (!is.null(ord <- model$order)) {
    
23.         if (length(ord) != 3L)
    
24.             stop("'model$order' must be of length 3")
    
25.         if (p != ord[1L])
    
26.             stop("inconsistent specification of 'ar' order")
    
27.         if (q != ord[3L])
    
28.             stop("inconsistent specification of 'ma' order")
    
29.         d <- ord[2L]
    
30.         if (d != round(d) || d < 0)
    
31.             stop("number of differences must be a positive integer")
    
32.     }
    
33.     if (!missing(start.innov) && length(start.innov) < n.start)
    
34.         stop(sprintf(ngettext(n.start, "'start.innov' is too short: need %d point",
    
35.             "'start.innov' is too short: need %d points"), n.start),
    
36.             domain = NA)
    
37.     x <- ts(c(start.innov[seq_len(n.start)], innov[1L:n]), start = 1 -
    
38.         n.start)
    
39.     if (length(model$ma)) {
    
40.         x <- filter(x, c(1, model$ma), sides = 1L)
    
41.         x[seq_along(model$ma)] <- 0
    
42.     }
    
43.     if (length(model$ar))
    
44.         x <- filter(x, model$ar, method = "recursive")
    
45.     if (n.start > 0)
    
46.         x <- x[-(seq_len(n.start))]
    
47.     if (d > 0)
    
48.         x <- diffinv(x, differences = d)
    
49.     as.ts(x)
    
50. }
    
51. <bytecode: 0x000000000a197c08>
    
52. <environment: namespace:stats>
    
53. >

复制代码

题主你的问题解决了吗？我遇到了同样的疑惑，求指点