r - 使用 R 的 SARIMAX 模型的数值方法

Question

我的 friend 目前正在使用最大似然估计 (MLE) 方法进行关于时间序列模型参数估计的作业，SARIMAX(季节性 ARIMA 外生)。他使用的数据是2000-2012年的月降雨量，以印度洋偶极子(IOD)指数为外生变量。
这是数据:

    MONTH YEAR   RAINFALL     IOD
1       1 2000 15.3720526  0.0624
2       2 2000 10.3440804  0.1784
3       3 2000 14.6116392  0.3135
4       4 2000 18.6842179  0.3495
5       5 2000 15.2937896  0.3374
6       6 2000 15.0233152  0.1946
7       7 2000 11.1803399  0.3948
8       8 2000 11.0589330  0.4391
9       9 2000 10.1488916  0.3020
10     10 2000 21.1187121  0.2373
11     11 2000 15.3980518 -0.0324
12     12 2000 18.9393770 -0.0148
13      1 2001 19.1075901 -0.2448
14      2 2001 14.9097284  0.1673
15      3 2001 19.2379833  0.1538
16      4 2001 19.6900990  0.3387
17      5 2001  8.0684571  0.3578
18      6 2001 14.0463518  0.3394
19      7 2001  5.9916609  0.1754
20      8 2001  8.4439327  0.0048
21      9 2001 11.8321596  0.1648
22     10 2001 24.3700636 -0.0653
23     11 2001 22.3584436  0.0291
24     12 2001 23.6114379  0.1731
25      1 2002 17.8409641  0.0404
26      2 2002 14.7377067  0.0914
27      3 2002 21.2226294  0.1766
28      4 2002 16.6403125 -0.1512
29      5 2002 10.8074049 -0.1072
30      6 2002  6.3796552  0.0244
31      7 2002 17.0704423  0.0542
32      8 2002  1.7606817  0.0898
33      9 2002  5.3665631  0.6736
34     10 2002  8.3246622  0.7780
35     11 2002 17.8044938  0.3616
36     12 2002 16.7062862  0.0673
37      1 2003 13.5572859 -0.0628
38      2 2003 17.1113997  0.2038
39      3 2003 14.9899967  0.1239
40      4 2003 14.0996454  0.0997
41      5 2003 11.4017542  0.0581
42      6 2003  6.7749539  0.3490
43      7 2003  7.1484264  0.4410
44      8 2003 10.3004854  0.4063
45      9 2003 10.6630202  0.3289
46     10 2003 20.6518764  0.1394
47     11 2003 20.8638443  0.1077
48     12 2003 20.5548048  0.4093
49      1 2004 16.0436903  0.2257
50      2 2004 17.2568827  0.2978
51      3 2004 20.2361063  0.2523
52      4 2004 11.6619038  0.1212
53      5 2004 12.8296532 -0.3395
54      6 2004  8.4202138 -0.1764
55      7 2004 15.5916644  0.0118
56      8 2004  0.9486833  0.1651
57      9 2004  7.2732386  0.2825
58     10 2004 18.0083314  0.3747
59     11 2004 14.4672043  0.1074
60     12 2004 17.3637554  0.0926
61      1 2005 18.9420168  0.0551
62      2 2005 17.0146995 -0.3716
63      3 2005 23.3002146 -0.2641
64      4 2005 17.8689675  0.2829
65      5 2005 17.2365890  0.1883
66      6 2005 14.0178458  0.0347
67      7 2005 12.6925175 -0.0680
68      8 2005  9.3861600 -0.0420
69      9 2005 11.7132404 -0.1425
70     10 2005 18.5768673 -0.0514
71     11 2005 19.6723156 -0.0008
72     12 2005 18.3248465 -0.0659
73      1 2006 18.6252517  0.0560
74      2 2006 18.7002674 -0.1151
75      3 2006 23.4882950 -0.0562
76      4 2006 19.5652754  0.1862
77      5 2006 13.6857590  0.0105
78      6 2006 11.1265448  0.1504
79      7 2006 11.0227038  0.3490
80      8 2006  7.6550637  0.5267
81      9 2006  1.8708287  0.8089
82     10 2006  5.4129474  0.9479
83     11 2006 15.2249795  0.7625
84     12 2006 14.1703917  0.3941
85      1 2007 22.8691932  0.4027
86      2 2007 14.3317829  0.3353
87      3 2007 13.0766968  0.2792
88      4 2007 23.2335964  0.2960
89      5 2007 12.2474487  0.4899
90      6 2007 11.3357840  0.2445
91      7 2007  9.3112835  0.3629
92      8 2007  1.6431677  0.5396
93      9 2007  6.8483575  0.6252
94     10 2007 13.1529464  0.4540
95     11 2007 14.5120639  0.2489
96     12 2007 18.7909553  0.0054
97      1 2008 17.6493626  0.3037
98      2 2008 13.3828248  0.1166
99      3 2008 19.0525589  0.2730
100     4 2008 17.3262806  0.0467
101     5 2008  5.2345009  0.4020
102     6 2008  3.3166248  0.4263
103     7 2008 10.1094016  0.5558
104     8 2008 11.7260394  0.4236
105     9 2008 10.7470926  0.4762
106    10 2008 15.1591557  0.4127
107    11 2008 25.5558213  0.1474
108    12 2008 18.2455474  0.1755
109     1 2009 14.5430396  0.2185
110     2 2009 12.8569048  0.3521
111     3 2009 24.0707291  0.2680
112     4 2009 16.0374562  0.3234
113     5 2009  7.2387844  0.4757
114     6 2009 13.8021737  0.3078
115     7 2009  7.5232972  0.1179
116     8 2009  6.3403470  0.1999
117     9 2009  4.6583259  0.2814
118    10 2009 13.0958008  0.3646
119    11 2009 15.3329710  0.1914
120    12 2009 19.0394328  0.3836
121     1 2010 15.5080624  0.4732
122     2 2010 17.1551742  0.2134
123     3 2010 23.9729014  0.6320
124     4 2010 18.2537667  0.5644
125     5 2010 18.2236111  0.1881
126     6 2010 14.6082169  0.0680
127     7 2010 13.6161669  0.3111
128     8 2010 11.1220502  0.2472
129     9 2010 20.7870152  0.1259
130    10 2010 19.5371441 -0.0529
131    11 2010 24.8837296 -0.2133
132    12 2010 15.5016128  0.0233
133     1 2011 17.3435867  0.3739
134     2 2011 17.6096564  0.4228
135     3 2011 19.0682983  0.5413
136     4 2011 20.4890214  0.3569
137     5 2011 12.0540450  0.1313
138     6 2011 12.5896783  0.2642
139     7 2011  5.0990195  0.5356
140     8 2011  6.5726707  0.6490
141     9 2011  2.5099801  0.5884
142    10 2011 17.6380271  0.7376
143    11 2011 17.5128524  0.6004
144    12 2011 17.2655727  0.0990
145     1 2012 16.6883193  0.2272
146     2 2012 20.8374663  0.1049
147     3 2012 16.7002994  0.1991
148     4 2012 18.7962762 -0.0596
149     5 2012 16.9292646 -0.1165
150     6 2012 11.6490343  0.2207
151     7 2012  6.2529993  0.8586
152     8 2012  5.8991525  0.9473
153     9 2012  7.8485667  0.8419
154    10 2012 12.5817328  0.4928
155    11 2012 24.7770055  0.1684
156    12 2012 23.2486559  0.4899

为此，他使用了 R，因为它具有用于分析 SARIMAX 模型的包。到目前为止，他一直在使用带有季节性 ARIMA 订单 (1,0,1) 的 TSA 包的 arimax() 函数做得很好。

所以在这里我附上他的语法:

#Importing data
data=read.csv("C:/DATA.csv", header=TRUE)
rainfall=data$RAINFALL
exo=data$IOD

#Creating the suitable model of data that is able to be read by R with ts() function
library(forecast)
rainfall_ts=ts(rainfall, start=c(2000, 1), end=c(2012, 12), frequency = 12)
exo_ts=ts(exo, start=c(2000, 1), end=c(2012, 12), frequency = 12)

#Fitting SARIMAX model with seasonal ARIMA order (1,0,1) & estimation method is MLE (or ML)
library(TSA)
model_ts=arimax(log(rainfall_ts), order=c(1,0,1), seasonal=list(order=c(1,0,1), period=12), xreg=exo_ts, method='ML')

结果如下:

> model_ts

Call:
arimax(x = log(rainfall_ts), order = c(1, 0, 1), seasonal = list(order = c(1, 
    0, 1), period = 12), xreg = exo_ts, method = "ML")

Coefficients:
         ar1      ma1    sar1     sma1  intercept     xreg
      0.5730  -0.4342  0.9996  -0.9764     2.6757  -0.4894
s.e.  0.2348   0.2545  0.0018   0.0508     0.1334   0.1489

sigma^2 estimated as 0.1521:  log likelihood = -86.49,  aic = 184.99

虽然他声称语法有效，但他的讲师期望更多。
理论上，因为他使用了 MLE，他已经证明了对数似然函数的一阶导数给出了隐式解。这意味着估计过程不能用 MLE 分析完成，所以我们需要
继续我们使用数值方法来完成它。

所以这是我 friend 的讲师的期望。他希望他至少可以说服他估计确实需要用数字来完成
如果是这样，他也许可以向他展示 R 使用的方法(数值方法，例如 Newton-Raphson、BFGS、BHHH 等)。

但这里的事情是 arimax() 函数没有给我们选择数值方法来选择是否需要以数值方式执行估计，如下所示:

model_ts=arimax(log(rainfall_ts), order=c(1,0,1), seasonal=list(order=c(1,0,1), period=12), xreg=exo_ts, method='ML')

上面的 'method' 是估计方法，可用的方法是 ML, CSS, 和 CSS-ML 。很明显，上面的语法不包含数值方法，这就是问题。

那么有没有可能的方法来知道 R 使用什么数值方法？或者我的 friend 只是在不依赖 arimax() 函数的情况下构建了自己的程序？

如果我的代码中有任何错误，请告诉我。我也为任何语法或词汇错误道歉。英语不是我的母语。

Answer 1

一些建议:

使用每种方法估计模型:ML , CSS , CSS-ML .参数估计是否一致？

您可以查看arimax()的源代码输入 arimax , View(arimax)或 getAnywhere(arimax)在控制台中。

或者，您可以通过在 model_ts=arimax(...) 行之前放置调试项目符号来进行调试。然后采购或debugSource() -ing 你的脚本。然后您可以进入arimax功能并查看/验证自己使用哪种优化方法 arimax使用。