Filtrado de Kalman a partir de primeros principios
F. Tusell
2021-11-02
Datos del caudal del Nilo
La serie de datos de caudal del Nilo es un ejemplo clásico en los textos de series temporales (incluido Durbin-Koopman, que los emplea extensivamente). Son datos del caudal anual del Nilo aguas abajo de Asswan durante un siglo. Sólo a efectos de poder hacer uso de estos datos sin teclearlos cargaremos el paquete dlm: el resto de este guión hace sólo uso de R “básico”, por motivos didácticos.
library(dlm)
data(Nile)
str(Nile)## Time-Series [1:100] from 1871 to 1970: 1120 1160 963 1210 1160 1160 813 1230 1370 1140 ...class(Nile)## [1] "ts"Los datos están en formato ts, el más rudimentario para almacenar series temporales: periodicidad fija, sin años bisiestos, pero suficiente para datos anuales como éstos. Veámoslos:
plot(Nile)
Observemos la rotulación automática de los ejes. La función genérica plot “sabe” lo que debe de hacer con un objeto de esta clase.
Estimemos las funciones ACF y PACF para estos datos:
acf(Nile)
pacf(Nile)
La ACF sugiere cierta memoria; la PACF en cambio no sustenta la inclusión de parámetros autorregresivos. Los datos, por otra parte, no parecen estacionarios en media: los más modernos parecen fluctuar en torno a un nivel más bajo que los precedentes.
Modelo de nivel local
Entre los modelos estructurales el más simple en que podemos pensar es el modelo de nivel local (o paseo aleatorio con ruido), que se especifica así:
\[\begin{eqnarray*} \alpha_{t+1} &=& \alpha_t + \eta_t \\ y_t &=& \alpha_t + \epsilon_t \end{eqnarray*}\]
La primera es la ecuación de estado, la segunda la ecuación de observación. En este caso, tanto la matriz de transición como la matriz de observación son unidad y de dimensión 1 x 1 . Podemos fácilmente programar una función de predicción utilizando el filtro de Kalman así:
kalman1 <- function(y,x0,P0,Q,H) {
n <- length(y)
x <- P <- rep(0,n)
x[1] <- x0
P[1] <- P0
for (i in 1:n) {
inno <- y[i] - x[i]
F <- P[i] + H
K <- P[i] / F
P[i+1] <- P[i]*(1-K) + Q
x[i+1] <- x[i] + K*inno
}
return(list(x=x,P=P))
}Las expresiones en el bucle transcriben casos particulares de las fórmulas de las ecuaciones del filtro de Kalman:
| Expresión | Caso particular… | …de la fórmula general |
|---|---|---|
inno <- y[i] - x[i] |
\(v_t = y_t - a_t\) | \(v_t = y_t - Z_ta_t\) |
F <- P[i] + H |
\(F_t = P_t + H_t\) | \(F_t = Z_tP_tZ'_t + H_t\) |
K <- P[i] / F |
\(K_t = P_t(P_t + H_t)^{-1}\) | \(K_t = Z_tP_t(Z_tP_tZ'_t + H_t)^{-1}\) |
P[i+1] <- P[i]*(1-K) + Q |
\(P_{t+1} = P_t(I-K_t) + Q_t\) | \(P_{t+1} = T_tP_tT'_t(I-K_t) + R_rQ_tR'_t\) |
x[i+1] <- x[i] + K*inno |
\(a_{t+1} = a_t + K_t v_t\) | \(a_{t+1} = T_ta_t + K_t v_t\) |
Las simplificaciones hechas consisten en omitir \(Z_t = T_t = R_t = 1\) y en hacer uso del hecho de que la dimensión \(1\times 1\) de las matrices involucradas permite sustituir algunas operaciones matriciales por cocientes: por ejemplo, \(K_t = P_t(P_t + H_t)^{-1}\) se transforma en \(K_t = P_t/(P_t + H_t)\) cuando todas las matrices involucradas son de dimensión \(1\times 1\).
Para invocar la función anterior hemos de especificar la media y matriz de covarianzas del estado inicialmente. Podemos o no tener información al respecto. Si no tenemos, una posibilidad es dar un vector de ceros y una matriz de covarianzas diagonal con valores muy grandes de las varianzas, reflejando una gran incertidumbre. Tenemos que dar también valores para los dos únicos parámetros, las varianzas de los ruidos (como veremos más adelante, en general los parámetros son objeto de estimación; de momento les damos dos valores meramente plausibles). Invocamos entonces la función y visualizamos los resultados así:
result <- kalman1(y=Nile,x0=0,P0=10^7,Q=100,H=1000)
result## $x
## [1] 0.0000 1119.8880 1140.8981 1072.5576 1117.3935 1130.1262 1138.6439
## [8] 1048.1009 1098.0019 1172.0872 1163.3808 1117.7860 1068.3441 1079.6051
## [15] 1056.4702 1046.6157 1023.2136 1065.5727 993.5544 983.9490 1026.1075
## [22] 1046.0701 1090.3569 1106.4699 1145.2454 1176.2471 1188.0672 1145.3644
## [29] 1133.1089 1036.0934 983.1176 953.6388 883.4957 898.7607 880.9951
## [36] 832.3683 854.9619 810.9367 867.4165 916.7425 930.8602 903.8824
## [43] 855.8263 747.8108 768.3938 750.4571 850.2914 917.7517 894.5854
## [50] 859.3069 848.9581 827.0867 831.9261 840.5911 846.3748 806.2904
## [57] 816.7481 797.0947 796.7990 862.5012 834.5397 820.0756 832.2122
## [64] 835.6669 864.9338 897.1003 897.0732 876.7917 912.7787 874.4763
## [71] 820.8567 774.4286 793.7640 798.6906 783.3753 788.1367 856.1791
## [78] 857.2114 861.7469 858.0331 866.6692 833.5293 810.6932 818.0703
## [85] 880.7276 890.7969 916.5166 884.2285 894.7028 916.3956 889.0030
## [92] 924.3926 919.4237 914.4465 983.4858 964.1735 905.2326 908.9519
## [99] 857.3651 818.6341 797.3906
##
## $P
## [1] 1.000000e+07 1.099900e+03 6.237868e+02 4.841556e+02 4.262162e+02
## [6] 3.988440e+02 3.851240e+02 3.780430e+02 3.743332e+02 3.723744e+02
## [11] 3.713359e+02 3.707841e+02 3.704905e+02 3.703342e+02 3.702510e+02
## [16] 3.702067e+02 3.701831e+02 3.701705e+02 3.701638e+02 3.701603e+02
## [21] 3.701584e+02 3.701574e+02 3.701568e+02 3.701565e+02 3.701564e+02
## [26] 3.701563e+02 3.701563e+02 3.701562e+02 3.701562e+02 3.701562e+02
## [31] 3.701562e+02 3.701562e+02 3.701562e+02 3.701562e+02 3.701562e+02
## [36] 3.701562e+02 3.701562e+02 3.701562e+02 3.701562e+02 3.701562e+02
## [41] 3.701562e+02 3.701562e+02 3.701562e+02 3.701562e+02 3.701562e+02
## [46] 3.701562e+02 3.701562e+02 3.701562e+02 3.701562e+02 3.701562e+02
## [51] 3.701562e+02 3.701562e+02 3.701562e+02 3.701562e+02 3.701562e+02
## [56] 3.701562e+02 3.701562e+02 3.701562e+02 3.701562e+02 3.701562e+02
## [61] 3.701562e+02 3.701562e+02 3.701562e+02 3.701562e+02 3.701562e+02
## [66] 3.701562e+02 3.701562e+02 3.701562e+02 3.701562e+02 3.701562e+02
## [71] 3.701562e+02 3.701562e+02 3.701562e+02 3.701562e+02 3.701562e+02
## [76] 3.701562e+02 3.701562e+02 3.701562e+02 3.701562e+02 3.701562e+02
## [81] 3.701562e+02 3.701562e+02 3.701562e+02 3.701562e+02 3.701562e+02
## [86] 3.701562e+02 3.701562e+02 3.701562e+02 3.701562e+02 3.701562e+02
## [91] 3.701562e+02 3.701562e+02 3.701562e+02 3.701562e+02 3.701562e+02
## [96] 3.701562e+02 3.701562e+02 3.701562e+02 3.701562e+02 3.701562e+02
## [101] 3.701562e+02Lo que dicha función produce son valores medios estimados \(a_t = E[\alpha_{t|t-1}]\) (una modificación sencilla proporcionaría \(\alpha_{t|t}\)).
Podemos observar algunas cosas interesantes. Una de ellas es que la varianza estimada del estado se torna constante después de algunas iteraciones. Se dice que el filtro ha alcanzado el estado estacionario (steady state). Esto ocurre muy generalmente con modelos invariantes en el tiempo (las matrices y parámetros de las ecuaciones de observación y estado no dependen de \(t\)) y puede aprovecharse para diseñar algoritmos rápidos (alcanzada la convergencia, una parte relativamente costosa de la recursión —la actualización de la matriz de covarianzas del estado— puede obviarse).
Podemos investigar el efecto de modificar los parámetros. Veamos que ocurre si modificamos (incrementándola) la varianza del ruido de observación. Observemos que la función kalman1 definida devuelve x como un vector. Hemos de convertirlo en una serie temporal de la misma clase que Nile para representarlo en la misma escala.
result2 <- kalman1(y=Nile,x0=0,P0=10^7,Q=500,H=1000)
Nile.ll <- ts(result$x,start=1871)
Nile.ll2 <- ts(result2$x,start=1871)
plot(Nile)
lines(Nile.ll,col="red")
lines(Nile.ll2,col="blue")
Al aumentar la varianza del ruido de observación “nos creemos menos” las medidas (los valores de \(y_t\)). El efecto es incorporar menos de ellas al estado, cuya trayectoria es así más suave.
Un efecto contrario se logra reduciendo el valor de la varianza \(\sigma_\epsilon\) (en H) o incrementando el valor de \(\sigma_\eta\) (en Q). Examinemos lo que ocurre en este último caso:
result <- kalman1(y=Nile,0,10^7,Q=20,H=1000)
Nile.ll3 <- ts(result$x,start=1871)
plot(Nile)
lines(Nile.ll,col="red")
lines(Nile.ll2,col="blue")
lines(Nile.ll3,col="green")
Es equivalente aumentar la varianza en la ecuación de estado o reducirla en la ecuación de observación: ambas cosas tienen por efecto incorporar “más” de las observaciones a la estimación del estado. De hecho, podríamos comprobar con facilidad que lo que realmente importa es la razón \(\sigma_\eta/\sigma_\epsilon\), que podemos ver como una relación señal-ruido.
Datos perdidos
El filtro del Kalman es tolerante a pérdida de datos. En el presente caso, univariante, \(y_t\) es totalmente observado o totalmente ausente. Es fácil modificar la función kalman1 de manera que se propaguen los valores del periodo anterior aunque no haya nueva observación:
kalman2 <- function(y,x0,P0,Q,H) {
n <- length(y)
x <- P <- rep(0,n)
x[1] <- x0
P[1] <- P0
for (i in 1:n) {
if (!is.na(y[i])) {
inno <- y[i] - x[i]
F <- P[i] + H
K <- P[i] / F
} else {
inno <- 0
K <- 0
}
x[i+1] <- x[i] + K * inno
P[i+1] <- P[i]*(1-K) + Q
}
return(list(x=x,P=P))
}Veamos cuál es el efecto de tener una racha de observaciones perdidas:
datos <- Nile
datos[70:76] <- NA
result <- kalman2(y=datos,0,10^7,Q=5000,H=100000)
Nile.ll.NA <- ts(result$x,start=1871)
plot(datos)
lines(Nile.ll.NA,col="red")
Como puede verse, el filtro extrapola horizontalmente (mismo nivel) cuando faltan datos. (En otros modelos, la extrapolación será de otro tipo, por ejemplo una tendencia lineal o de otro orden.) Es también interesante ver lo que ocurre con la varianza estimada del estado:
par(mfrow=c(2,1))
par(mar=c(3,3,0,0))
plot(datos,col="red")
plot(ts(result$P[-1],start=1872),type="l")
En el caso multivariante, si parte del vector \(y_t\) no está disponible, podemos modificar la ecuación de observación \(y_t = Z_ty_t + \epsilon_t\) de manera que se tomen en consideración sólo las filas correspondientes a valores observados: recordemos qje \(Z_t\) puede variar y acomodar así dimensiones cambiantes del vector de observación \(y_t\). El filtro permite aprovechar aquélla información que esté disponible aunque no sea completa.
Hay algo adicional que podemos observar. Estimemos con los parámetros anteriores el valor del estado y su matriz de covarianzas, tanto con los datos originales como con otros completamente diferentes.
result <- kalman1(y=Nile,0,10^7,Q=500,H=15000)
result2 <- kalman1(y=rnorm(100),0,10^7,Q=500,H=15000)
cbind(result$P,result2$P)## [,1] [,2]
## [1,] 10000000.000 10000000.000
## [2,] 15477.534 15477.534
## [3,] 8117.513 8117.513
## [4,] 5767.119 5767.119
## [5,] 4665.565 4665.565
## [6,] 4058.681 4058.681
## [7,] 3694.356 3694.356
## [8,] 3464.282 3464.282
## [9,] 3314.311 3314.311
## [10,] 3214.525 3214.525
## [11,] 3147.221 3147.221
## [12,] 3101.408 3101.408
## [13,] 3070.027 3070.027
## [14,] 3048.442 3048.442
## [15,] 3033.550 3033.550
## [16,] 3023.255 3023.255
## [17,] 3016.129 3016.129
## [18,] 3011.190 3011.190
## [19,] 3007.766 3007.766
## [20,] 3005.391 3005.391
## [21,] 3003.743 3003.743
## [22,] 3002.598 3002.598
## [23,] 3001.804 3001.804
## [24,] 3001.253 3001.253
## [25,] 3000.870 3000.870
## [26,] 3000.604 3000.604
## [27,] 3000.420 3000.420
## [28,] 3000.291 3000.291
## [29,] 3000.202 3000.202
## [30,] 3000.140 3000.140
## [31,] 3000.098 3000.098
## [32,] 3000.068 3000.068
## [33,] 3000.047 3000.047
## [34,] 3000.033 3000.033
## [35,] 3000.023 3000.023
## [36,] 3000.016 3000.016
## [37,] 3000.011 3000.011
## [38,] 3000.008 3000.008
## [39,] 3000.005 3000.005
## [40,] 3000.004 3000.004
## [41,] 3000.003 3000.003
## [42,] 3000.002 3000.002
## [43,] 3000.001 3000.001
## [44,] 3000.001 3000.001
## [45,] 3000.001 3000.001
## [46,] 3000.000 3000.000
## [47,] 3000.000 3000.000
## [48,] 3000.000 3000.000
## [49,] 3000.000 3000.000
## [50,] 3000.000 3000.000
## [51,] 3000.000 3000.000
## [52,] 3000.000 3000.000
## [53,] 3000.000 3000.000
## [54,] 3000.000 3000.000
## [55,] 3000.000 3000.000
## [56,] 3000.000 3000.000
## [57,] 3000.000 3000.000
## [58,] 3000.000 3000.000
## [59,] 3000.000 3000.000
## [60,] 3000.000 3000.000
## [61,] 3000.000 3000.000
## [62,] 3000.000 3000.000
## [63,] 3000.000 3000.000
## [64,] 3000.000 3000.000
## [65,] 3000.000 3000.000
## [66,] 3000.000 3000.000
## [67,] 3000.000 3000.000
## [68,] 3000.000 3000.000
## [69,] 3000.000 3000.000
## [70,] 3000.000 3000.000
## [71,] 3000.000 3000.000
## [72,] 3000.000 3000.000
## [73,] 3000.000 3000.000
## [74,] 3000.000 3000.000
## [75,] 3000.000 3000.000
## [76,] 3000.000 3000.000
## [77,] 3000.000 3000.000
## [78,] 3000.000 3000.000
## [79,] 3000.000 3000.000
## [80,] 3000.000 3000.000
## [81,] 3000.000 3000.000
## [82,] 3000.000 3000.000
## [83,] 3000.000 3000.000
## [84,] 3000.000 3000.000
## [85,] 3000.000 3000.000
## [86,] 3000.000 3000.000
## [87,] 3000.000 3000.000
## [88,] 3000.000 3000.000
## [89,] 3000.000 3000.000
## [90,] 3000.000 3000.000
## [91,] 3000.000 3000.000
## [92,] 3000.000 3000.000
## [93,] 3000.000 3000.000
## [94,] 3000.000 3000.000
## [95,] 3000.000 3000.000
## [96,] 3000.000 3000.000
## [97,] 3000.000 3000.000
## [98,] 3000.000 3000.000
## [99,] 3000.000 3000.000
## [100,] 3000.000 3000.000
## [101,] 3000.000 3000.000Los valores de P que obtenemos son idénticos. ¿Es casualidad?
Problemas
Modificar la función
kalman1de modo que haga filtrado en lugar de predicción un periodo hacia adelante.Justificar el último resultado obtenido (la estimación de la matriz de covarianzas del vector de estado no depende de las medidas \(y_t\)).