有时 R “不太友好”

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

问题 1: 截距和均值到底用那一个?

拟合ARIMA模型,R处理的是均值还是截距？如果是中心化的AR模型，就无所谓了，如果不是，那就麻烦了。例如，假设x(t)=α+φ*x(t-1)+w(t)是平稳序列。取期望，则有μ=α+φ*μ 或 α=μ*(1-φ)。所以，截距α不是均值μ，除非φ=0。一般来讲，中心化的AR模型，截距和均值是相同的。如下例：

# 生成均值是50的 AR(1) 序列

set.seed(66)      # 设置种子，可以重复验证

x = arima.sim(list(order=c(1,0,0), ar=.9), n=100) + 50  

mean(x) 

  [1] 50.60668   # 样本均值

arima(x, order = c(1, 0, 0))  

  Coefficients:

           ar1  intercept  <--  出问题啦

        0.8971    50.6304  <--  这里要修改

  s.e.  0.0409     0.8365

拟合AR(1)模型是
x(t) = 50.6304 + .8971*x(t-1) + w(t)
而实际的模型是
x(t)−50.6304 = .8971*[x(t-1)−50.6304] + w(t)
或
x(t) = 5.21 + .8971*x(t-1) + w(t)

【5.21 = 50.6304*(1-.8971)】

处理的方法是把截距改成均值即可：

  Coefficients:

           ar1       mean

        0.8971    50.6304

  s.e.  0.0409     0.8365

R在这一点上处理的不好，但也有人说这是仁者见仁。

如果使用 ar.ols，却会得到另一组数据：

ar.ols(x, order=1, demean=F, intercept=T)

  Coefficients:

           1 

      0.9052               <-- 估计 φ

  Intercept: 4.806 (2.167) <-- 这里是截距

说明：arima() 使用 MLE, 而 ar.ols() 使用 OLS 拟合模型，因此模型参数有差异。有一点可以确定，R中对截距的解释是开放的，不是准确的最优解释。

问题 2: 为什么 arima 拟合与差分拟合有差异?  xreg是什么?

用R拟合ARIMA模型，有差分，模型不包含常数项；无差分，默认的方式是拟合均值[细节见 ?arima ]。这样做有什么问题？ 还是用一个例子说明：

arima(x, order = c(1, 1, 0))          #(1)

和下面的模型产生不同的结果

arima(diff(x), order = c(1, 0, 0))    #(2)

因为在模型(1)中，R拟合的是 [ ∇x(s) = x(s)-x(s-1)]
∇x(t)= φ*∇x(t-1) + w(t) 
而模型(2), R拟合的是
∇x(t) = α + φ*∇x(t-1) + w(t).

如果有漂移 (例如, α 不等于零), 两个模型差别明显，使用模型 (1) 将导致拟合不正确，那么后续的预测自然掺不忍睹 (见后面的问题)。

如果 α 不等于零，则必须在模型 (1) 中使用 xreg :

arima(x, order = c(1, 1, 0), xreg=1:length(x))    #(1+)

为什么要这样做？ 选择 xreg = t ， R 用x(t)-β*t替代 x(t) ；也就是它拟合模型 
∇[x(t) - β*t] = φ*∇[x(t-1) - β*(t-1)] + w(t)
简化成 
∇x(t) = α + φ*∇x(t-1) + w(t) 这里 α = β*(1-φ)

下面用随机行走的例子建立ARIMA(1,1,0)模型来说这一点：

set.seed(1)           #设置种子，可以重复验证

v = rnorm(100,1,1)    # v 包含 100 iid N(1,1)

x = cumsum(v)         # x 是一步随机行走

plot.ts(x)            # 生成图形...

 

arima(x, order = c(1, 1, 0))         #(1)

 

   Coefficients:

           ar1

        0.6031

  s.e.  0.0793

 

arima(diff(x), order = c(1, 0, 0))   #(2)

 

  Coefficients:

           ar1  intercept  <-- 记住,这是diff(x)的均值不是截距

       -0.0031     1.1163    

  s.e.  0.1002     0.0897

 

arima(x, order = c(1, 1, 0), xreg=1:length(x))    #(1+)

 

  Coefficients:

            ar1  1:length(x)  <-- 这里是截距，对diff(x)...头痛吗？

        -0.0031      1.1169                     

  s.e.   0.1002      0.0897

解释一下到底发生了什么。产生模型的数据是
x(t) = 1 + x(t-1) + w(t) 
这里 w(t) 是服从正态分布 N(0,1)的白噪声。换一种写法 
[x(t)-x(t-1)] = 1 + 0*[x(t-1)-x(t-2)] + w(t) 
或
∇x(t) = 1 + 0*∇x(t-1) + w(t) 
如果用∇x(t)拟合 AR(1)，那么参数应该是ar1 = 0和intercept = 1。

说明：模型 (1) 给出一个错误答案因为它强迫回归通过原点。但模型 (2) 和模型 (1+)用两种不同的方式给出正确的答案。

问题 3: 为什么做预测，Arima 会给出奇怪的结果?

如果用ARIMA(p,d,q)得到预测结果，那么前面所讲的问题还是要注意。这里使用第三章的全球温度数据来说明。如下所见，第一个拟合模型得出错误结果，第二个模型给出正确的答案。如果使用 sarima 和sarima.for，最好是远离第一个拟合模型，这样就可以避免尴尬的问题。

 u=read.table("/mydata/globtemp2.dat")  # 读入数据

 gtemp=ts(u[,2],start=1880,freq=1)      # 用第二列产生临时变量

 fit1=arima(gtemp, order=c(1,1,1))  

 fore1=predict(fit1, 15)  

 nobs=length(gtemp)               

 fit2=arima(gtemp, order=c(1,1,1), xreg=1:nobs) 

 fore2=predict(fit2, 15, newxreg=(nobs+1):(nobs+15))

 par(mfrow=c(2,1))

 ts.plot(gtemp,fore1$pred,col=1:2,main="WRONG")

 ts.plot(gtemp,fore2$pred,col=1:2,main="RIGHT")

图形如下:

问题 4: tsdiag。Arima 给出错误的p-values

如果ARIMA拟合后使用 tsdiag()诊断，会发现如下是这样的：

p-values for Ljung-Box statistic 给出的p-values 是错误的，因为计算p-values 的自由度是 lag而不是lag - (p+q)。即，程序没有考虑残差来自拟合模型这一事实。至少有一个R的核心开发团队知道这一点。

 

问题5：前移还是后移

 

对时间序列使用延迟函数要小心。lag(x,1) 是前移，lag(x,-1) 是后移。见下例：

x=ts(1:5)

cbind(x,lag(x,1),lag(x,-1))

Time Series:

Start = 0

  End = 6

  Frequency = 1

    x lag(x, 1) lag(x, -1)

0  NA       1          NA

1   1       2          NA

2   2       3           1

3   3       4           2 <-- 这里, x 是 3, lag(x,1) 是 4, lag(x,-1) 是 2

4   4       5           3

5   5      NA           4

6  NA      NA           5

软件默认 lag(x) 是 lag(x,1)。所以对时间序列  x(t)， 
y(t) = lag{x(t)} = x(t+1)，不是 x(t-1)。

如果想用 lag(x,-1)进行回归分析，就要处理两组数据的连接问题。

x = arima.sim(list(order=c(1,0,0), ar=.9), n=20) # 产生20个 AR(1)数据

xb1 = lag(x,-1)

  ##-- 错误的结果:

cor(x,xb1) # 计算相关系数

   [1] 1 ... 尽然是1? 

  lm(x~xb1) # 求回归方程

   Coefficients:

(Intercept)     xb1 

  6.112e-17   1.000e+00

你认为 x(t)=x(t-1)，但软件不是这样想的

    ##-- 正确的结果:

y=cbind(x,xb1) # 首先连接数据

lm(y[,1]~y[,2]) # 求回归

   Coefficients:

(Intercept) y[, 2]

   0.5022 0.7315

  ##-- 或者:

y=ts.intersect(x,xb1) # 用这种方式进行连接

lm(y[,1]~y[,2]) # 求回归方程

Coefficients:

    (Intercept)  y[, 2]

        0.5022   0.7315

 cor(y[,1],y[,2]) # 计算相关系数

[1] 0.842086

顺便提一下，(Intercept) 这里的截距是正确的。

问题 6: ACF和PACF引起的混乱

拟合ARMA 模型，第一步就是观察数据的 ACF 和 PACF 。见下例：模拟 MA(1)：

MA1=arima.sim(list(order=c(0,0,1), ma=.5), n=100)

par(mfcol=c(2,1))

acf(MA1,20) pacf(MA1,20)

输出如下:

发现什么问题吗？首先，两个图形的刻度不同。ACF 范围是 -.2 到 1，而 PACF 是-.2 到 .4。同样ACF的lag起点是0 ；而 PACF 的起点是1。

做如下修改：

par(mfrow=c(2,1)) acf(MA1,20,xlim=c(1,20)) # 把 x-axis 起点设置为1

pacf(MA1,20,ylim=c(-.2,1)) # 调整y-axis范围

输出如下:

 

（注：来自戴申老师的教学博客）

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏6 回帖

关键词：coefficients coefficient EFFICIENT Intercept ARIMA模型 友好

好详细，学习了

对R是大白，学习学习

现象、原因、解决办法，一应俱全

非常好的资料。

现在刚好也在用R做时间序列预估，来学习下~~~

qitelata 发表于 2014-10-27 09:16
现象、原因、解决办法，一应俱全

很牛啊

看看！

顶顶顶

挺难的