Estou tentando entender como os parâmetros são estimados na modelagem ARIMA / Box Jenkins (BJ). Infelizmente, nenhum dos livros que encontrei descreve o procedimento de estimativa, como o procedimento de estimativa do Log-Likelihood em detalhes. Achei o site / material didático que foi muito útil. A seguir, é apresentada a equação da fonte mencionada acima.
Quero aprender a estimativa do ARIMA / BJ fazendo isso sozinho. Então, usei para escrever um código para estimar o ARMA manualmente. Abaixo está o que eu fiz em R ,
- Simulei o ARMA (1,1)
- Escreveu a equação acima como uma função
- Utilizou os dados simulados e a função otim para estimar os parâmetros AR e MA.
- Também executei o ARIMA no pacote de estatísticas e comparei os parâmetros do ARMA com o que fiz manualmente. Abaixo está a comparação:
** Abaixo estão minhas perguntas:
- Por que há uma pequena diferença entre as variáveis estimadas e calculadas?
- O ARIMA funciona em backcasts R ou o procedimento de estimativa é diferente do descrito abaixo no meu código?
- Eu atribui e1 ou erro na observação 1 como 0, está correto?
- Também existe uma maneira de estimar os limites de confiança das previsões usando o hessian da otimização?
Muito obrigado pela sua ajuda, como sempre.
Abaixo está o código:
## Load Packages
library(stats)
library(forecast)
set.seed(456)
## Simulate Arima
y <- arima.sim(n = 250, list(ar = 0.3, ma = 0.7), mean = 5)
plot(y)
## Optimize Log-Likelihood for ARIMA
n = length(y) ## Count the number of observations
e = rep(1, n) ## Initialize e
logl <- function(mx){
g <- numeric
mx <- matrix(mx, ncol = 4)
mu <- mx[,1] ## Constant Term
sigma <- mx[,2]
rho <- mx[,3] ## AR coeff
theta <- mx[,4] ## MA coeff
e[1] = 0 ## Since e1 = 0
for (t in (2 : n)){
e[t] = y[t] - mu - rho*y[t-1] - theta*e[t-1]
}
## Maximize Log-Likelihood Function
g1 <- (-((n)/2)*log(2*pi) - ((n)/2)*log(sigma^2+0.000000001) - (1/2)*(1/(sigma^2+0.000000001))*e%*%e)
##note: multiplying Log-Likelihood by "-1" in order to maximize in the optimization
## This is done becuase Optim function in R can only minimize, "X"ing by -1 we can maximize
## also "+"ing by 0.000000001 sigma^2 to avoid divisible by 0
g <- -1 * g1
return(g)
}
## Optimize Log-Likelihood
arimopt <- optim(par=c(10,0.6,0.3,0.5), fn=logl, gr = NULL,
method = c("L-BFGS-B"),control = list(), hessian = T)
arimopt
############# Output Results###############
ar1_calculated = arimopt$par[3]
ma1_calculated = arimopt$par[4]
sigmasq_calculated = (arimopt$par[2])^2
logl_calculated = arimopt$val
ar1_calculated
ma1_calculated
sigmasq_calculated
logl_calculated
############# Estimate Using Arima###############
est <- arima(y,order=c(1,0,1))
est
time-series
forecasting
arima
optimization
box-jenkins
previsor
fonte
fonte
g1
+0.000000001
Respostas:
Existe o conceito de probabilidade exata. Requer o conhecimento de parâmetros iniciais, como o primeiro valor do erro MA (uma de suas perguntas). As implementações geralmente diferem quanto à forma como tratam os valores iniciais. O que eu costumo fazer é (o que não é mencionado em muitos livros) também é maximizar também a ML com os valores iniciais.
Por favor, dê uma olhada no seguinte de Tsay, ele não cobre todos os casos, mas foi bastante útil para mim:
http://faculty.chicagobooth.edu/ruey.tsay/teaching/uts/lec8-08.pdf
fonte
Você leu a página de ajuda da
arima
função? Aqui está o trecho relevante:Também relevante é um parâmetro
method=c("CSS-ML", "ML", "CSS")
:Seus resultados não diferem muito daqueles produzidos por
arima
função, então você definitivamente acertou tudo.Lembre-se de que, se você quiser comparar os resultados de dois procedimentos de otimização, verifique se os valores iniciais são os mesmos e o mesmo método de otimização é usado, caso contrário, os resultados podem ser diferentes.
fonte