Szacowanie ARIMA ręcznie

Próbuję zrozumieć, w jaki sposób parametry są szacowane w modelowaniu ARIMA / Box Jenkins (BJ). Niestety żadna z książek, które napotkałem, nie opisuje szczegółowo procedury szacowania, takiej jak procedura szacowania wiarygodności logarytmicznej. Znalazłem stronę internetową / materiały dydaktyczne, które były bardzo pomocne. Poniżej znajduje się równanie ze źródła wymienionego powyżej.

L. L. (θ) = - \frac{n}{2)} \log (2) π) - \frac{n}{2)} \log (σ^{2)}) - \sum_{t = 1}^{n} \frac{{mi}_{t}^{2)}}{2) σ^{2)}}

$LL(\theta)=-\frac{n}{2}\log(2\pi) - \frac{n}{2}\log(\sigma^2) - \sum\limits_{t=1}^n\frac{e_t^2}{2\sigma^2}$

Chcę poznać oszacowanie ARIMA / BJ, wykonując to samodzielnie. Użyłem więc do napisania kodu do ręcznego oszacowania ARiMR. Poniżej jest to, co zrobiłem w , $R$ $R$

Symulowałem ARiMR (1,1)
Napisałem powyższe równanie jako funkcję
Wykorzystano symulowane dane i funkcję optymalizacyjną do oszacowania parametrów AR i MA.
Uruchomiłem również ARIMA w pakiecie statystyk i porównałem parametry ARMA z tego, co zrobiłem ręcznie. Poniżej znajduje się porównanie:

Porównanie

** Poniżej znajdują się moje pytania:

Dlaczego istnieje niewielka różnica między zmiennymi szacowanymi i obliczanymi?
Czy ARIMA działa w backcastach R lub czy procedura szacowania różni się od opisanej poniżej w moim kodzie?
Przypisałem e1 lub błąd przy obserwacji 1 jako 0, czy to prawda?
Czy istnieje również sposób oszacowania granic ufności prognoz przy użyciu hessian optymalizacji?

Dziękuję bardzo za pomoc jak zawsze.

Poniżej znajduje się kod:

## 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

— Synoptyk
źródło

T

$T$

n

$n$

T

$T$

T = n + 1

$T=n+1$ g1+0.000000001

σ

$\sigma$

Zmieniłem równanie i teraz odzwierciedla to, co jest w kodzie. Nie jestem pewien, jak mogłem usunąć +0.000000001, ponieważ spowoduje to „NaNs” z powodu podzielności przez 0, a także z powodu problemu z logiem (0)

— prezenter

Istnieje koncepcja dokładnego prawdopodobieństwa. Wymaga znajomości początkowych parametrów, takich jak wartość pięści błędu MA (jedno z twoich pytań). Wdrożenia zwykle różnią się pod względem sposobu traktowania wartości początkowych. To, co zwykle robię, to (o czym nie wspomina wiele książek) jest również maksymalizacja ML względem wartości początkowych.

— Cagdas Ozgenc

Proszę spojrzeć na następujących od TSAY, nie wszystkie przypadki, ale był bardzo pomocny dla mnie: faculty.chicagobooth.edu/ruey.tsay/teaching/uts/lec8-08.pdf

— Cagdas Ozgenc

@CagdasOzgenc wskazane przez ciebie wartości początkowe są przyczyną różnic. Mogę zaakceptować twój komentarz jako odpowiedź, jeśli opublikujesz swój komentarz jako odpowiedź.

— prezenter

Odpowiedzi:

Proszę spojrzeć na następujące z Tsay, nie obejmuje wszystkich przypadków, ale było dla mnie bardzo pomocne:

http://faculty.chicagobooth.edu/ruey.tsay/teaching/uts/lec8-08.pdf

— Cagdas Ozgenc
źródło

Czy czytałeś stronę pomocy arimafunkcji? Oto odpowiedni fragment:

Dokładne prawdopodobieństwo jest obliczane na podstawie reprezentacji w przestrzeni stanu procesu ARIMA, a także innowacji i ich wariancji znalezionych przez filtr Kalmana. Inicjalizacja zróżnicowanego procesu ARMA wykorzystuje stacjonarność i jest oparta na Gardner i in. (1980). W przypadku zróżnicowanego procesu niestacjonarne elementy są wcześniej rozpraszane (kontrolowane przez kappa). Obserwacje, które są nadal kontrolowane przez wcześniejsze rozproszenie (określone przez zysk Kalmana co najmniej 1e4) są wyłączone z obliczeń prawdopodobieństwa. (To daje wyniki porównywalne do arima0 przy braku brakujących wartości, gdy wykluczone są obserwacje dokładnie pomijane przez różnicowanie.)

Istotny jest również parametr method=c("CSS-ML", "ML", "CSS"):

Metoda dopasowania: maksymalne prawdopodobieństwo lub minimalizacja warunkowej sumy kwadratów. Domyślnie (chyba że brakuje wartości) jest użycie warunkowej sumy kwadratów do znalezienia wartości początkowych, a następnie maksymalnego prawdopodobieństwa.

Twoje wyniki nie różnią się tak bardzo od tych wyprodukowanych przez arima funkcję, więc na pewno wszystko masz dobrze.

Pamiętaj, że jeśli chcesz porównać wyniki dwóch procedur optymalizacyjnych, musisz upewnić się, że wartości początkowe są takie same i stosowana jest ta sama metoda optymalizacji, w przeciwnym razie wyniki mogą się różnić.

— mpiktas
źródło