Proceso AR (1) con errores de medición heteroscedastic

1. El problema

Tengo algunas mediciones de una variable $y_t$ , donde $t=1,2,..,n$ , para el cual tengo una distribución $f_{y_t}(y_t)$ obtenida a través de MCMC, que por simplicidad supondré que es una gaussiana de media $\mu_t$ y varianza $\sigma_t^2$ .

Tengo un modelo físico para esas observaciones, digamos $g(t)$ , pero los residuos $r_t = \mu_t-g(t)$ parecen estar correlacionados; en particular, tengo razones físicas para pensar que un proceso $AR(1)$ será suficiente para tener en cuenta la correlación, y planeo obtener los coeficientes del ajuste a través de MCMC, para lo cual necesito la probabilidad . Creo que la solución es bastante simple, pero no estoy muy seguro (parece tan simple, que creo que me falta algo).

2. Derivando la probabilidad

Un proceso media cero se puede escribir como: X t = ϕ X t - 1 + ε t , ( 1 ) donde asumiré ε t ∼ N ( 0 , σ 2 w ) . Los parámetros a estimar son, por lo tanto, θ = { ϕ , σ 2 w } (en mi caso, también tengo que agregar los parámetros del modelo g ( t $AR(1)$

$$X_t = \phi X_{t-1}+\varepsilon_t,\ \ \ (1)$$ $\varepsilon_t\sim N(0,\sigma_w^2)$$\theta = \{\phi,\sigma_w^2\}$$g(t)$, pero ese no es el problema). Sin embargo, lo que observo es la variable donde supongo que η t ∼ N ( 0 , σ 2 t ) , y se conocen los σ 2 t (los errores de medición) . Como X t es un proceso gaussiano, R t también lo es. En particular, sé que X 1 ∼ N ( 0 , σ 2 w $$R_t = X_t+\eta_t,\ \ \ (2)$$ $\eta_t\sim N(0,\sigma_t^2)$$\sigma_t^2$$X_t$$R_t$ por lo tanto, R 1 ∼ N ( 0 , σ 2 w / [ 1 - ϕ 2 ] + σ 2 t ) . El próximo desafío es obtener R t | R t - 1 para t ≠ 1 . Para derivar la distribución de esta variable aleatoria, tenga en cuenta que, utilizando la ecuación. ( 2 ) Puedo escribir X $$X_1 \sim N(0,\sigma_w^2/[1-\phi^2]),$$ $$R_1 \sim N(0,\sigma_w^2/[1-\phi^2]+\sigma_t^2).$$ $R_t|R_{t-1}$$t\neq 1$$(2)$ Usando la ecuación. (2), y usando la definición de la ecuación. (1), puedo escribir, R t = X t + η t =ϕ X t - 1 + ε t + η t . Usando la ec. (3)en esta última expresión, entonces, obtengo, R $$X_{t-1} = R_{t-1}-\eta_{t-1}.\ \ \ (3)$$ $(2)$$(1)$ $$R_{t} = X_t+\eta_t = \phi X_{t-1}+\varepsilon_{t}+\eta_t.$$ $(3)$ tanto, R t | R t - 1 = ϕ ( r t - 1 - η t - 1 ) + ε t + η t , y, por lo tanto, R t | R t - 1 ∼ N $$R_{t} = \phi (R_{t-1}-\eta_{t-1})+\varepsilon_{t}+\eta_t,$$ $$R_t|R_{t-1} = \phi (r_{t-1}-\eta_{t-1})+\varepsilon_{t}+\eta_t,$$ Finalmente, puedo escribir la función de probabilidad como L ( θ ) = f R 1 ( R 1 = r 1 ) n ∏ t = 2 f R t | R t - 1 ( R t = r $$R_t|R_{t-1} \sim N(\phi r_{t-1},\sigma_w^2+\sigma_t^2-\phi^2\sigma^2_{t-1}).$$ donde f ( ⋅ ) son las distribuciones de las variables que acabo de definir, es decir, que definen σ ′ 2 = σ 2 w / [ 1 $$L(\theta) = f_{R_1}(R_1=r_1) \prod_{t=2}^{n} f_{R_{t}|R_{t-1}}(R_t=r_t|R_{t-1}=r_{t-1}),$$ $f(\cdot)$$\sigma'^2 = \sigma_w^2/[1-\phi^2]+\sigma_t^2,$ y definiendoσ2(t)=σ 2 w +σ 2 t -ϕ2σ 2 t - 1 , fRt| Rt-1(Rt=rt|Rt-1=rt-1)= $$f_{R_1}(R_1=r_1) = \frac{1}{\sqrt{2\pi \sigma'^2}}\text{exp}\left(-\frac{r_1^2}{2\sigma'^2}\right),$$ $\sigma^2(t) = \sigma_w^2+\sigma_t^2-\phi^2\sigma^2_{t-1}$ $$f_{R_{t}|R_{t-1}}(R_t=r_t|R_{t-1}=r_{t-1})=\frac{1}{\sqrt{2\pi \sigma^2(t)}}\text{exp}\left(-\frac{(r_t-\phi r_{t-1})^2}{2\sigma^2(t)}\right)$$

3. Preguntas

1. ¿Está bien mi derivación? No tengo ningún recurso para comparar que no sean simulaciones (que parecen estar de acuerdo), ¡y no soy un estadístico!
2. $MA(1)$$ARMA(1,1)$$ARMA(p,q)$

1. Estás en el camino correcto, pero has cometido un error al derivar la distribución de $R_t$ dado $R_{t-1}$: la media condicional no es $\phi r_{t-1}$. Sus$\phi \widehat{x}_{t-1}$, dónde $\widehat{x}_{t-1}$ es tu mejor estimación de $X$del periodo anterior El valor de$\widehat{x}_{t-1}$ incluye información de observaciones anteriores, así como $r_{t-1}$. (Para ver esto, considere una situación donde$\sigma_w$ y $\phi$son insignificantes, por lo que efectivamente está estimando una media fija. Después de muchas observaciones, su incertidumbre sobre$X$ será mucho más pequeño que $\sigma_{\eta}$.) Esto puede ser confuso al principio, porque observa $R$ y no $X$. Eso solo significa que se trata de un modelo de espacio de estado .

2. Sí, hay un marco muy general para usar modelos lineales-gaussianos con observaciones ruidosas, llamado filtro de Kalman . Esto se aplica a cualquier cosa con una estructura ARIMA y muchos más modelos también. Tiempo variable$\sigma_{\eta}$está bien para el filtro de Kalman, siempre que no sea estocástico. Los modelos con, por ejemplo, volatilidad estocástica , necesitan métodos más generales. Para ver cómo se deriva el filtro de Kalman, pruebe Durbin-Koopman o el capítulo 3 de Harvey . En la notación de Harvey, tu modelo tiene$Z=1$, $d=c=0$, $H_t = \sigma_{\eta,t}^2$, $T=\phi$, $R=1$ y $Q=\sigma^2_w$.

Honestamente, debe codificar esto en ERRORES o STAN y no preocuparse por eso desde allí. A menos que sea una pregunta teórica.