Rob J Hyndman gives great explanation here
(https://robjhyndman.com/hyndsight/estimation/) for reasons why results
from R's arima may differ from other softwares.
@iacobus, to cite one, 'Major discrepancies between R and Stata for
ARIMA'
(https://stackoverflow.com/questions/22443395/major-discrepancies-between-r-and-stata-for-arima),
assign the, sometimes, big diferences from R and other softwares to
different optimization algorithms. However, I think this is overstate
the reason.
I explain better. I fit arima models regularly using |forecast| or
|fable| packages, besides using Stata, Eviews and Gretl. All these
packages, except for R, give very consistent results with each other.
I'll give one example using R, Eviews and Gretl. I'll use "BFGS"
algorithm and observed hessian based standard errors (in Gretl and in
Eviews).
|library(GetBCBData) library(lubridate) library(tsibble)
library(tsbox) library(forecast) library(tidyr) library(dplyr)
#============================================================# #Data
---- #============================================================#
#Brazilian CPI and analytical components ipca <-
gbcbd_get_series(c(433, 4449, 10844, 11428, 27863, 27864),
first.date = dmy("01/01/2004")) ipca <- ipca %>% mutate(series.name
= case_when(id.num == 433 ~ "ipca", id.num == 4449 ~
"administrados", id.num == 10844 ~ "servicos", id.num == 11428 ~
"livres", id.num == 27863 ~ "industriais", id.num == 27864 ~
"alimentos", TRUE ~ series.name)) ipca <- ipca %>% select(data =
ref.date, valor = value, series.name) %>% pivot_wider(names_from =
"series.name", values_from = "valor") ipca_tsb <- ipca %>%
mutate(data = yearmonth(data)) %>% arrange(data) %>% as_tsibble()
ipca_ts <- ipca_tsb %>% ts_ts() ##Eviews and Gretl can easily import
'dta' files ---- ipca %>% foreign::write.dta("ipca.dta")
#============================================================#
#Model ----
#============================================================#
#ARIMA(2,0,1)(2,0,2)[12] modelo <- ipca_ts %>% .[, "servicos"] %>%
Arima(order = c(2, 0, 1), seasonal = c(2, 0, 2), include.mean = T,
method = "ML", optim.method = "BFGS", optim.control = list(maxit =
1000)) #'fable' gives identical results: ipca_tsb%>%
model(ARIMA(servicos ~ 1 + pdq(2, 0, 1) + PDQ(2, 0, 2), method =
"ML", optim.method = "BFGS", optim.control = list(maxit = 1000)))
%>% report() summary(modelo) |*|Series: . ARIMA(2,0,1)(2,0,2)[12]
with non-zero mean Coefficients: ar1 ar2 ma1 sar1 sar2 sma1 sma2
mean 0.7534 0.0706 -0.5705 0.1759 0.7511 0.3533 -0.6283 0.5001 s.e.
NaN NaN 0.0011 NaN NaN NaN NaN 0.1996 sigma^2 = 0.05312: log
likelihood = 1.75 AIC=14.5 AICc=15.33 BIC=45.33 Training set error
measures: ME RMSE MAE MPE MAPE MASE ACF1 Training set -0.006082139
0.2263897 0.1678378 -33.39711 79.74708 0.7674419 0.01342733 Warning
message: In sqrt(diag(x$var.coef)) : NaNs produce|*
|Gretl||output:
https://drive.google.com/file/d/1T_thtM0mRXvlJbPrgkwqlc_tLqCbsXYe/view?usp=sharing|
|Eviews output:
https://drive.google.com/file/d/1Ta8b5vPwftRFhZpb3Pg95aVRfvhl01vO/view?usp=sharing|
Coefficients comparisson:
https://docs.google.com/spreadsheets/d/1TfXmQaCEOtOX6e0foSHI9gTAHJ1vW6Ui/edit?usp=sharing&ouid=104001235447994085528&rtpof=true&sd=true
Both Eviews and Gretl give results that differ after a few decimal
places, which is expected. R, by its turn, gives completlely different
results. Again, all of them used "BFGS" algorithm. Even the standard
erros, which I'm not questioning here, are very similar between Gretl
and Eviews. In this example, the major difference is in AR(1) coefficient.
I know this is just only one example, but I see these kind of divergent
results everyday. Not to mention situations where including a single
observation messes up the entire result, or when R cannot calculate
standard errors ("In sqrt(diag(best$var.coef))").
All that make me wonder: as results may differ greatly from other
software, is |arima| from R truly reliable? Is |optim| the problem or
other intricacies inside |arima| function and its methodology? Or are
all the other softwares making mistakes (converging when they shouldn't,
for example)?
In a summary, the question is: is R's base |arima| (which |Arima|, from
'forecast', and |ARIMA| from 'fable' are based on) really reliable?
Edit: Stata gives the same results as did Gretl and Eviews. Stata
output:
https://drive.google.com/file/d/1TeKfj59aJNjxaWx0Uslke4y8RJWqMxRl/view?usp=sharing
Rodrigo Remédio
rreme...@hotmail.com
||*||*||
______________________________________________
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
