Estimation ARIMA à la main

J'essaie de comprendre comment les paramètres sont estimés dans la modélisation ARIMA / Box Jenkins (BJ). Malheureusement, aucun des livres que j'ai rencontrés ne décrit en détail la procédure d'estimation telle que la procédure d'estimation de log-vraisemblance. J'ai trouvé le site Web / matériel pédagogique très utile. Voici l'équation de la source référencée ci-dessus.

L L (θ) = - \frac{n}{2} Journal (2 π) - \frac{n}{2} Journal (σ^{2}) - \sum_{t = 1}^{n} \frac{e_{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}$

Je veux apprendre l'estimation ARIMA / BJ en le faisant moi-même. J'ai donc utilisé pour écrire un code pour estimer ARMA à la main. Voici ce que j'ai fait dans , $R$ $R$

J'ai simulé ARMA (1,1)
A écrit l'équation ci-dessus en fonction
Utilisé les données simulées et la fonction optim pour estimer les paramètres AR et MA.
J'ai également exécuté l'ARIMA dans le package de statistiques et comparé les paramètres ARMA à partir de ce que j'ai fait à la main. Voici la comparaison:

Comparaison

** Voici mes questions:

Pourquoi y a-t-il une légère différence entre les variables estimées et calculées?
ARIMA fonctionne-t-il dans les retransmissions R ou la procédure d'estimation est-elle différente de celle décrite ci-dessous dans mon code?
J'ai attribué e1 ou erreur à l'observation 1 à 0, est-ce correct?
Existe-t-il également un moyen d'estimer les limites de confiance des prévisions en utilisant la toile de jute de l'optimisation?

Merci beaucoup pour votre aide, comme toujours.

Voici le code:

## 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 prévisionniste
la source

T

$T$

n

$n$

T

$T$

T = n + 1

$T=n+1$ g1+0.000000001

σ

$\sigma$

J'ai changé l'équation et reflète maintenant ce qui est dans le code. Je ne sais pas comment je pourrais supprimer +0.000000001, car cela entraînera des "NaN" en raison de divisible par 0 et également en raison du problème de journal (0)

prévisionniste

Il y a le concept de vraisemblance exacte. Cela nécessite la connaissance des paramètres initiaux tels que la première valeur de l'erreur MA (une de vos questions). Les implémentations diffèrent généralement en ce qui concerne la façon dont elles traitent les valeurs initiales. Ce que je fais habituellement est (ce qui n'est pas mentionné dans de nombreux livres) est de maximiser également ML par rapport aux valeurs initiales.

Cagdas Ozgenc

Veuillez jeter un coup d'œil à ce qui suit de Tsay, il ne couvre pas tous les cas, mais il m'a été très utile: faculty.chicagobooth.edu/ruey.tsay/teaching/uts/lec8-08.pdf

Cagdas Ozgenc

@CagdasOzgenc, comme indiqué par vos valeurs initiales, est la cause de la différence. Je peux accepter votre commentaire comme réponse si vous postez vos commentaires comme réponse.

prévisionniste

Réponses:

Veuillez jeter un coup d'œil à ce qui suit de Tsay, il ne couvre pas tous les cas, mais m'a été très utile:

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

Cagdas Ozgenc
la source

Avez-vous lu la page d'aide de la arimafonction? Voici l'extrait pertinent:

La vraisemblance exacte est calculée via une représentation en espace d'état du processus ARIMA, et les innovations et leur variance trouvées par un filtre de Kalman. L'initialisation du processus ARMA différencié utilise la stationnarité et est basée sur Gardner et al. (1980). Pour un processus différencié, les composants non stationnaires reçoivent un a priori diffus (contrôlé par kappa). Les observations qui sont toujours contrôlées par l'a priori diffus (déterminé en ayant un gain de Kalman d'au moins 1e4) sont exclues des calculs de vraisemblance. (Cela donne des résultats comparables à arima0 en l'absence de valeurs manquantes, lorsque les observations exclues sont précisément celles abandonnées par la différenciation.)

Un paramètre est également pertinent method=c("CSS-ML", "ML", "CSS"):

Méthode d'ajustement: maximum de vraisemblance ou minimisation de la somme conditionnelle des carrés. La valeur par défaut (à moins qu'il n'y ait des valeurs manquantes) consiste à utiliser la somme conditionnelle des carrés pour trouver les valeurs de départ, puis la probabilité maximale.

Vos résultats ne diffèrent pas beaucoup de ceux produits par la arimafonction, vous avez donc certainement tout bien.

N'oubliez pas que si vous souhaitez comparer les résultats de deux procédures d'optimisation, vous devez vous assurer que les valeurs de départ sont les mêmes et que la même méthode d'optimisation est utilisée, sinon les résultats peuvent différer.

mpiktas
la source