ARIMA 的最大似然估计 - 与使用 arima() 计算的值不匹配答案

【问题标题】：Maximum Likelihood estimate of ARIMA-not matching the value computed using arima()ARIMA 的最大似然估计 - 与使用 arima() 计算的值不匹配
【发布时间】：2014-04-09 15:02:18
【问题描述】：

我是 r 的新手。我正在使用 optim 函数来获得 arima 函数的最大似然估计，假设残差是正态分布的。我已经对数据进行了一次差异化以使其静止。我写了以下来计算可能性：-

kings <- scan("http://robjhyndman.com/tsdldata/misc/kings.dat",skip=3)

arima1<-function(a=length(kings))
{e<-array(0,dim=a-2);
e[1:2]=0
y=diff(kings)



likelihood<-function(AR,e,y)
{for(i in 3:41)
{e[i]<-sum(y[i],-AR[1],-(AR[2]*y[i-1]),-(AR[3]*y[i-2]),-(AR[4]*e[i-1]),-(AR[5]*e[i-2]))
}
-sum(-(a-1)*(log(AR[6]*(2*22/7)^.5)),-(sum(e^2)/(2*(AR[6])^2)))
}

optim(par<-c(0,0,1,0,1,14),likelihood,y=y,e=e,control=list(maxit=1000))

}
arima1()

但是使用函数arima(y,order

【问题讨论】：

你应该自己做这个功能吗？ R中有一些包可以为ARIMA模型做MLE，比如TSA（时间序列分析）。
我的一个朋友建议我先手动执行该功能，不要采用自动化功能。他说这是了解系统工作原理的最佳方式。

标签： r statistics

【解决方案1】：

无法理解您的日志功能。此外，我认为您对 arima 和 optim 函数的期望过高。这里有几个cmets：

出于您的目的（学习），我将使用 arima.sim 函数生成一个系列。这将使您知道optim 和arima 应该给出的值。
arima 函数有时会出错。在连续两次运行中，这就是我电脑中的结果：

y = arima.sim(n=1000,list(ar=c(0.6, -0.2), ma=c(-0.7, 0.5)),mean=0.3)
有马(y,order=c(2,0,2))$coef
ar1：0.3489 ar2：0.1060 ma1：-0.4318 ma2：0.1887 截距：0.4919

y = arima.sim(n=1000,list(ar=c(0.6, -0.2), ma=c(-0.7, 0.5)),mean=0.3)
有马(y,order=c(2,0,2))$coef
ar1:0.6100 ar2:-0.1672 ma1:-0.6663 ma2:0.4816 截距：0.4199强>

请注意，第一个结果的 ar2 甚至没有相同的符号，而且到处都是蹩脚的。但第二个非常适合。相同的模拟函数参数，但结果不同。

请注意，时间序列上的值越少，就越难获得准确的结果。这是我电脑上的一些运行。

y = arima.sim(n=42,list(ar=c(0.6, -0.2), ma=c(-0.7, 0.5)),mean=0.3)
有马(y,order=c(2,0,2))$coef
ar1:-1.1102 ar2:-0.4528 ma1:1.2052 ma2:0.9999 截距：0.2353260强>
y = arima.sim(n=42,list(ar=c(0.6, -0.2), ma=c(-0.7, 0.5)),mean=0.3)
有马(y,order=c(2,0,2))$coef
ar1:0.8623 ar2:-0.3935 ma1:-1.0745 ma2:0.9999 截距：0.4109
y = arima.sim(n=42,list(ar=c(0.6, -0.2), ma=c(-0.7, 0.5)),mean=0.3)
有马(y,order=c(2,0,2))$coef
ar1:-0.2749 ar2:0.4170 ma1:-0.0078 ma2:-0.0586 截距：0.3737

所有这些结果都很糟糕。

4.我发现您的代码有点难以阅读/理解。这是我在我的机器上运行的（有结果）。

library(stats4)

kings <- scan("http://robjhyndman.com/tsdldata/misc/kings.dat",skip=3)
y = diff(kings)
Y_0 = y[c(-1,-2)] # this is y(t)
Y_1 = y[c(-1,-length(y))] # this is y(t-1)
Y_2 = y[1:(length(y)-2)] # this is y(t-2)

logARMA22 <- function ( ar1=0.1, ar2=-0.2, ma1=-0.7, ma2=0.1, alpha=0.42, sigma=14)
{
    E_0 = array(0, dim=length(Y_0))
    E_1 = array(0, dim=length(Y_0))
    E_2 = array(0, dim=length(Y_0))

    for ( i in 1:length(E_0) )
    {   
        E_0[i] = Y_0[i] - ar1*Y_1[i] - ar2*Y_2[i]
                 - ma1*E_1[i] - ma2*E_2[i] - alpha
        if ( (i+1) <= length(E_1) )
            E_1 [i+1] = E_0[i]
        if ( (i+2) <= length(E_2) )
            E_2 [i+2] = E_0[i]
    }   

    # e^2
    E2 = (Y_0-alpha-ar1*Y_1-ar2*Y_2-ma1*E_1-ma2*E_2)^2
    res = suppressWarnings( - sum(log( (1/(sqrt(2*pi)*sigma)) * exp(-(E2)/(2*sigma^2)))) )
    return(res)
}

print(arima(y, order=c(2,0,2))$coef)
print(mle(logARMA22))

#Read 42 items
#         ar1          ar2          ma1          ma2    intercept 
# 0.087452103 -0.122367077 -0.740473745  0.006492026  0.372799202 
#
#Call:
#mle(minuslogl = logARMA22)
#
#Coefficients:
#        ar1         ar2         ma1         ma2       alpha       sigma 
#-0.07772158 -0.48516681 -0.50804965  0.03678942  0.10385160 15.33830148

请注意，mle 只是 optim 的前端（因此您实际上在这里调用了 optim）。结果不相等可能是由于样本的大小。本例中的对数似然函数取自“Henrik Madsen, 2008, Time Series Analysis. equation 6.72 (page 167)”。 this link 也可能会有所帮助。

希望对你有帮助

【讨论】：

【解决方案2】：

arima() 函数不使用最大似然法来拟合模型。

通过 ARIMA 过程的状态空间表示以及卡尔曼滤波器发现的创新及其方差来计算确切的可能性。

【讨论】：