Hay bastantes paquetes en R que permiten hacer filtrado, predicción, suavizado y estimación sobre modelos en espacio de estado. Todos tienen sus puntos fuertes y débiles, y es frecuente que un análisis requiera hacer uso de varios de ellos. Nosotros ilustraremos el empleo de dlm
(en la presente sesión) y de KFAS
(en otra posterior).
dlm
es un paquete muy sólido y “todo terreno”, que incluye buena documentación a la que se añade el libro Petris et al.(2009) Dynamic Linear Models in R, Springer Verlag. Ha sido utilizado (al sólo efecto de leer los datos en el objeto Nile
) en una de las sesiones prácticas precedentes. Comenzaremos por cargarlo y leer de nuevo dichos datos.
library(dlm)
data(Nile)
dlm
Un modelo se construye utilizando la función dlm
con argumentos que especifican las matrices intervinientes. En el caso más simple de un modelo invariante en el tiempo (=con todas sus matrices y parámetros constantes) debemos especificar:
m0
, y su matriz de covarianzas, C0
. Si no tenemos información a priori, típicamente señalaremos un vector de ceros y una matriz diagonal con varianzas muy grandes.GG
del estado.FF
.W
.V
.Para especificar el modelo de nivel local,
\[\begin{eqnarray} \alpha_t &=& \alpha_{t-1} + \eta_t \\ y_t &=& \alpha_t + \epsilon_t \end{eqnarray}\]nivel local
haremos lo siguiente:
<- dlm(m0=0,C0=10^7,FF=1,GG=1,W=500,V=15000)
mod1 mod1
## $FF
## [,1]
## [1,] 1
##
## $V
## [,1]
## [1,] 15000
##
## $GG
## [,1]
## [1,] 1
##
## $W
## [,1]
## [1,] 500
##
## $m0
## [1] 0
##
## $C0
## [,1]
## [1,] 1e+07
class(mod1)
## [1] "dlm"
Hemos fijado los dos únicos parámetros —las varianzas de los ruidos— en valores similares a los empleados previamente, a efectos de comparación. Cualesquiera elementos de un modelo se puede ver o alterar directamente, o haciendo uso de funciones de extracción:
V(mod1)
## [,1]
## [1,] 15000
$V mod1
## [,1]
## [1,] 15000
V(mod1) <- 10000
Para hacer tanto predicción una etapa hacia adelante del estado (es decir, para obtener \(a_t = \alpha_{t|t-1}\)) como para obtener valores filtrados (es decir, \(\alpha_{t|t}\)) empleamos la misma función dlmFilter
:
<- dlmFilter(y=Nile,mod=mod1) result1
Examinemos el objeto result1
:
str(result1,max.level=1)
## List of 9
## $ y : Time-Series [1:100] from 1871 to 1970: 1120 1160 963 1210 1160 1160 813 1230 1370 1140 ...
## $ mod:List of 6
## ..- attr(*, "class")= chr "dlm"
## $ m : Time-Series [1:101] from 1870 to 1970: 0 1119 1140 1076 1115 ...
## $ U.C:List of 101
## $ D.C: num [1:101, 1] 3162.3 100 71.6 60 53.9 ...
## $ a : Time-Series [1:100] from 1871 to 1970: 0 1119 1140 1076 1115 ...
## $ U.R:List of 100
## $ D.R: num [1:100, 1] 3162.4 102.4 75 64 58.4 ...
## $ f : Time-Series [1:100] from 1871 to 1970: 0 1119 1140 1076 1115 ...
## - attr(*, "class")= chr "dlmFiltered"
En el componente y
tenemos una copia de la serie original, en m
los valores filtrados \(\alpha_{t|t}\), en a
las predicciones un periodo hacia adelante, \(\alpha_{t|t-1}\) y en f
las predicciones una etapa hacia adelante de las observaciones (es decir, \(y_{t|t-1} = z_t'\alpha_{t|t-1}\)). Podemos representar gráficamente las predicciones (a
ó f
, son iguales en este caso) frente a las observaciones y comparar con lo obtenido anteriormente.
plot(dropFirst(result1$y))
lines(dropFirst(result1$a), col="red")
lines(dropFirst(result1$m), col="blue")
Por lo que hace a las matrices de covarianzas del estado (en este caso, simples varianzas), no aparecen directamente en el objeto result1
devuelto por dlmFilter
. El motivo es que no se emplea el algoritmo de filtro de Kalman tal como lo hemos presentado, sino una variedad de filtro “raíz cuadrada” (que, en el caso particular de dlm
, propaga los componentes de la descomposición en valores singulares de \(P_t\), a partir de los cuáles recomponer ésta). Para recuperar las matrices de covarianzas de los valores filtrados de result1
podemos hacer lo siguiente:
<- rep(0,length(result1$U.C))
Var.m for (i in seq(along=result1$U.C)) {
<- result1$U.C[[i]] %*% result1$D.C[i,]^2 %*% t(result1$U.C[[i]])
Var.m[i] }
Observamos la rápida convergencia al estado estacionario (eliminamos el primer valor, correspondiente a la varianza a priori, para evitar la compresión de la escala).
plot(ts(Var.m[-1],start=1870),main="Varianza de los valores filtrados del estado")
Con modelos de alguna complejidad, la construcción directa de las matrices de transición, etc. puede hacerse larga y proclive a error. dlm
ofrece funciones auxiliares para simplificar el trabajo, junto con la posibilidad de especificar separadamente “trozos” de modelo que luego se unen en un modelo final.
La función dlmModPoly
permite generar con una invocación modelos de nivel local (order=1), de tendencia local lineal (order=2), y de órdenes superiores. Como argumentos opcionales están las (matrices de) covarianzas de los ruidos y la media y matriz de covarianzas iniciales del estado. Por ejemplo, un modelo de tendencia local lineal con observaciones unidimensionales puede especificarse así:
<- dlmModPoly(order=2)
mod mod
## $FF
## [,1] [,2]
## [1,] 1 0
##
## $V
## [,1]
## [1,] 1
##
## $GG
## [,1] [,2]
## [1,] 1 1
## [2,] 0 1
##
## $W
## [,1] [,2]
## [1,] 0 0
## [2,] 0 1
##
## $m0
## [1] 0 0
##
## $C0
## [,1] [,2]
## [1,] 1e+07 0e+00
## [2,] 0e+00 1e+07
Observemos que la especificación que hace difiere un poco de la que puede encontrarse en algunos textos (e.g., Harvey, A.C.(1989) Forecasting structural time series models and the Kalman filter, Cambridge Univ. Press), en que a la primera componente del vector de estado se le atribuye perturbación de varianza cero. Nada impide, si lo deseamos, “parchear” la matriz W
introduciendo una varianza en el lugar correspondiente:
$W[1,1] <- 0.5
mod mod
## $FF
## [,1] [,2]
## [1,] 1 0
##
## $V
## [,1]
## [1,] 1
##
## $GG
## [,1] [,2]
## [1,] 1 1
## [2,] 0 1
##
## $W
## [,1] [,2]
## [1,] 0.5 0
## [2,] 0.0 1
##
## $m0
## [1] 0 0
##
## $C0
## [,1] [,2]
## [1,] 1e+07 0e+00
## [2,] 0e+00 1e+07
Del mismo modo podemos proceder siempre que el resultado de una función constructora no proporcione exactamente lo que queremos: por ejemplo, si queremos correlación contemporánea entre las perturbaciones del vector de estado o de las observaciones.
Podemos también especificar fácilmente un modelo estacional. Por ejemplo, con datos diarios y estacionalidad semanal (= 7 “estaciones”) podríamos hace:
<- dlmModSeas(frequency=7)
mod mod
## $FF
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 1 0 0 0 0 0
##
## $V
## [,1]
## [1,] 1
##
## $GG
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] -1 -1 -1 -1 -1 -1
## [2,] 1 0 0 0 0 0
## [3,] 0 1 0 0 0 0
## [4,] 0 0 1 0 0 0
## [5,] 0 0 0 1 0 0
## [6,] 0 0 0 0 1 0
##
## $W
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 1 0 0 0 0 0
## [2,] 0 0 0 0 0 0
## [3,] 0 0 0 0 0 0
## [4,] 0 0 0 0 0 0
## [5,] 0 0 0 0 0 0
## [6,] 0 0 0 0 0 0
##
## $m0
## [1] 0 0 0 0 0 0
##
## $C0
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 1e+07 0e+00 0e+00 0e+00 0e+00 0e+00
## [2,] 0e+00 1e+07 0e+00 0e+00 0e+00 0e+00
## [3,] 0e+00 0e+00 1e+07 0e+00 0e+00 0e+00
## [4,] 0e+00 0e+00 0e+00 1e+07 0e+00 0e+00
## [5,] 0e+00 0e+00 0e+00 0e+00 1e+07 0e+00
## [6,] 0e+00 0e+00 0e+00 0e+00 0e+00 1e+07
Cuando el intervalo de observación es muy pequeño en relación a la duración de los ciclos estacionales, el empleo de dummies se hace engorroso y es mejor emplear funciones trigonométricas. Este sería el modo de especificar con datos diarios una estacionalidad semanal:
<- dlmModTrig(s=7,q=1)
mod $GG mod
## [,1] [,2]
## [1,] 0.6234898 0.7818315
## [2,] -0.7818315 0.6234898
cos(2*pi/7)
## [1] 0.6234898
sin(2*pi/7)
## [1] 0.7818315
A veces, la estacionalidad no queda bien aproximada por una sinusoide, sino que necesitamos introducir varios armónicos. Haríamos lo siguiente para introducir la frecuencia fundamental y dos armónicos:
<- dlmModTrig(s=7,q=3)
mod mod
## $FF
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 1 0 1 0 1 0
##
## $V
## [,1]
## [1,] 1
##
## $GG
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 0.6234898 0.7818315 0.0000000 0.0000000 0.0000000 0.0000000
## [2,] -0.7818315 0.6234898 0.0000000 0.0000000 0.0000000 0.0000000
## [3,] 0.0000000 0.0000000 -0.2225209 0.9749279 0.0000000 0.0000000
## [4,] 0.0000000 0.0000000 -0.9749279 -0.2225209 0.0000000 0.0000000
## [5,] 0.0000000 0.0000000 0.0000000 0.0000000 -0.9009689 0.4338837
## [6,] 0.0000000 0.0000000 0.0000000 0.0000000 -0.4338837 -0.9009689
##
## $W
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 0 0 0 0 0 0
## [2,] 0 0 0 0 0 0
## [3,] 0 0 0 0 0 0
## [4,] 0 0 0 0 0 0
## [5,] 0 0 0 0 0 0
## [6,] 0 0 0 0 0 0
##
## $m0
## [1] 0 0 0 0 0 0
##
## $C0
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 1e+07 0e+00 0e+00 0e+00 0e+00 0e+00
## [2,] 0e+00 1e+07 0e+00 0e+00 0e+00 0e+00
## [3,] 0e+00 0e+00 1e+07 0e+00 0e+00 0e+00
## [4,] 0e+00 0e+00 0e+00 1e+07 0e+00 0e+00
## [5,] 0e+00 0e+00 0e+00 0e+00 1e+07 0e+00
## [6,] 0e+00 0e+00 0e+00 0e+00 0e+00 1e+07
Una característica que hace el paquete dlm
cómodo para el trabajo es la capacidad de sumar “trozos de modelo”. El operador “+” ha sido sobrecargado para construir modelos a partir de sus partes integrantes. Por ejemplo, si quisiéramos un modelo que incluyera a la vez una tendencia local lineal y estacionalidad, especificada mediante dummies, de periodo 4, podríamos especificarlo así:
<- dlmModPoly(order=2) + dlmModSeas(frequency=4)
mod mod
## $FF
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 0 1 0 0
##
## $V
## [,1]
## [1,] 2
##
## $GG
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 1 0 0 0
## [2,] 0 1 0 0 0
## [3,] 0 0 -1 -1 -1
## [4,] 0 0 1 0 0
## [5,] 0 0 0 1 0
##
## $W
## [,1] [,2] [,3] [,4] [,5]
## [1,] 0 0 0 0 0
## [2,] 0 1 0 0 0
## [3,] 0 0 1 0 0
## [4,] 0 0 0 0 0
## [5,] 0 0 0 0 0
##
## $m0
## [1] 0 0 0 0 0
##
## $C0
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1e+07 0e+00 0e+00 0e+00 0e+00
## [2,] 0e+00 1e+07 0e+00 0e+00 0e+00
## [3,] 0e+00 0e+00 1e+07 0e+00 0e+00
## [4,] 0e+00 0e+00 0e+00 1e+07 0e+00
## [5,] 0e+00 0e+00 0e+00 0e+00 1e+07
Podemos identificar en la matriz de transición los bloques recogiendo las dinámicas de la tendencia local lineal y la estacionalidad. Si quisiéramos un nivel local y estacionalidad anual con datos diarios, podriamos especificar:
<- dlmModPoly(order=1) + dlmModTrig(s=365,q=1)
mod mod
## $FF
## [,1] [,2] [,3]
## [1,] 1 1 0
##
## $V
## [,1]
## [1,] 2
##
## $GG
## [,1] [,2] [,3]
## [1,] 1 0.00000000 0.00000000
## [2,] 0 0.99985184 0.01721336
## [3,] 0 -0.01721336 0.99985184
##
## $W
## [,1] [,2] [,3]
## [1,] 1 0 0
## [2,] 0 0 0
## [3,] 0 0 0
##
## $m0
## [1] 0 0 0
##
## $C0
## [,1] [,2] [,3]
## [1,] 1e+07 0e+00 0e+00
## [2,] 0e+00 1e+07 0e+00
## [3,] 0e+00 0e+00 1e+07
Un modelo de tendencia local lineal con estacionalidad anual, datos diarios y en el que la estacionalidad estuviera representada por la frecuencia fundamental y el primer armónico, podría especificarse así:
<- dlmModPoly(order=2) + dlmModTrig(tau=365.25,q=2)
mod mod
## $FF
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 1 0 1 0 1 0
##
## $V
## [,1]
## [1,] 2
##
## $GG
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 1 1 0.00000000 0.00000000 0.00000000 0.00000000
## [2,] 0 1 0.00000000 0.00000000 0.00000000 0.00000000
## [3,] 0 0 0.99985204 0.01720158 0.00000000 0.00000000
## [4,] 0 0 -0.01720158 0.99985204 0.00000000 0.00000000
## [5,] 0 0 0.00000000 0.00000000 0.99940821 0.03439806
## [6,] 0 0 0.00000000 0.00000000 -0.03439806 0.99940821
##
## $W
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 0 0 0 0 0 0
## [2,] 0 1 0 0 0 0
## [3,] 0 0 0 0 0 0
## [4,] 0 0 0 0 0 0
## [5,] 0 0 0 0 0 0
## [6,] 0 0 0 0 0 0
##
## $m0
## [1] 0 0 0 0 0 0
##
## $C0
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 1e+07 0e+00 0e+00 0e+00 0e+00 0e+00
## [2,] 0e+00 1e+07 0e+00 0e+00 0e+00 0e+00
## [3,] 0e+00 0e+00 1e+07 0e+00 0e+00 0e+00
## [4,] 0e+00 0e+00 0e+00 1e+07 0e+00 0e+00
## [5,] 0e+00 0e+00 0e+00 0e+00 1e+07 0e+00
## [6,] 0e+00 0e+00 0e+00 0e+00 0e+00 1e+07
Obsérvese que aunque los bloques aparezcan desacoplados, nada impide que tomemos el objeto mod
resultante y lo modifiquemos en el sentido que deseemos.
Podemos también al coste de una cierta complicación considerar modelos cuyas matrices de sistema (transición, observación, matrices de covarianzas…) evolucionan a lo largo del tiempo. El procedimiento es el siguiente. Además de los argumentos FF
, V
, GG
y W
, podemos proporcionar a dlm
otros cinco. Los cuatro primeros, son matrices nombradas como las de sistema, con el prefijo J
: JFF
, JV
, JGG
y JW
. El quinto tiene por nombre X
. El convenio es el siguiente: si un elemento de Jxx
es cero, el correspondiente término de xx
es invariante. Si por el contrario un elemento de Jxx
es un entero j
distinto de cero, el elemento correspondiente de xx
en el momento i
se toma de la fila i
y columna j
de X
.
Imaginemos por ejemplo una serie y
que suponemos depende linealmente de otra, z
, pero con coeficientes variables a lo largo del tiempo, siguiente un paseo aleatorio. Podríamos especificar algo como lo siguiente:
<- rnorm(100) # Valores cualesquiera
z <- rnorm(100) # Valores cualesquiera
y <- dlmModPoly(order=1) + dlmModPoly(order=1)
mod mod
## $FF
## [,1] [,2]
## [1,] 1 1
##
## $V
## [,1]
## [1,] 2
##
## $GG
## [,1] [,2]
## [1,] 1 0
## [2,] 0 1
##
## $W
## [,1] [,2]
## [1,] 1 0
## [2,] 0 1
##
## $m0
## [1] 0 0
##
## $C0
## [,1] [,2]
## [1,] 1e+07 0e+00
## [2,] 0e+00 1e+07
Tenemos ya la “base”, que podemos parchear. Añadimos lo que nos falta:
$JFF <- matrix(c(0,1),1,2)
mod$X <- matrix(z,100,1)
mod mod
## $FF
## [,1] [,2]
## [1,] 1 1
##
## $V
## [,1]
## [1,] 2
##
## $GG
## [,1] [,2]
## [1,] 1 0
## [2,] 0 1
##
## $W
## [,1] [,2]
## [1,] 1 0
## [2,] 0 1
##
## $JFF
## [,1] [,2]
## [1,] 0 1
##
## $X
## [,1]
## [1,] -1.919
## [2,] 0.643
## [3,] ...
##
## $m0
## [1] 0 0
##
## $C0
## [,1] [,2]
## [1,] 1e+07 0e+00
## [2,] 0e+00 1e+07
El caso de regresión con parámetros cambiantes es con mucho el más frecuente, y por ello dlm
proporciona una función constructora que simplifica el trabajo. El mismo modelo especificado más arriba puede especificarse así:
<- dlmModReg(X=z,addInt=TRUE,dV=1,dW=c(1,1))
mod.bis mod.bis
## $FF
## [,1] [,2]
## [1,] 1 1
##
## $V
## [,1]
## [1,] 1
##
## $GG
## [,1] [,2]
## [1,] 1 0
## [2,] 0 1
##
## $W
## [,1] [,2]
## [1,] 1 0
## [2,] 0 1
##
## $JFF
## [,1] [,2]
## [1,] 0 1
##
## $X
## [,1]
## [1,] -1.919
## [2,] 0.643
## [3,] ...
##
## $m0
## [1] 0 0
##
## $C0
## [,1] [,2]
## [1,] 1e+07 0e+00
## [2,] 0e+00 1e+07
Si no indicamos dW=c(1,1)
se fijan las varianzas por defecto a cero (parámetros fijos).
El procedimiento de especificación es completamente general y permite hacer cualesquiera matrices variables en el tiempo.
Como ilustración, podemos estudiar la posibilidad de un cambio de régimen en el caudal del Nilo tras el cierre de la primera presa de Assuan en 1898. Comenzamos por definir una variable que recoja el efecto de la intervención,
<- ts(c(rep(0,28),rep(1,72)), start=1870)
x plot(x)
A continuación, especificamos un modelo que podría se el precedente, en que como variable “regresor” tomamos
x
y forzamos a que su coeficiente sea constante:
<- dlmModReg(X=x,addInt=TRUE,dV=1,dW=c(1,0)) mod.interv
A continuación filtramos,
<- dlmFilter(y=Nile,mod=mod.interv)
result2 plot(dropFirst(result2$y), ylim=c(0,1450),
ylab=expression(Hm^3/año), xlab="Año",
main="Caudal anual Nilo y ajustes")
# Ajuste con salto
# lines(dropFirst(result2$a[,1]) + dropFirst(result2$a[,2]), col="red")
lines(dropFirst(result2$f), col="red")
# Ajuste anterior, para comparar
lines(dropFirst(result1$f), col="blue")
# Valor estimado del salvo
lines(dropFirst(mean(Nile[1:28]) + result2$a[,2]), col="green", lwd=2)
legend("bottomleft",
legend=c("Caudal","Nivel local","Id. + Salto","Salto"),
lwd=c(1,1,1,1),
col = c("black","blue","red","green"))
Alternativamente,
library(ggplot2)
##
## Attaching package: 'ggplot2'
## The following object is masked from 'package:dlm':
##
## %+%
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
## ✓ tibble 3.1.6 ✓ dplyr 1.0.7
## ✓ tidyr 1.1.4 ✓ stringr 1.4.0
## ✓ readr 2.1.1 ✓ forcats 0.5.1
## ✓ purrr 0.3.4
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x ggplot2::%+%() masks dlm::%+%()
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
<- data.frame(Nile=Nile,
result Fecha=time(Nile),
predic.ser=result2$f,
salto=mean(Nile[1:28]) + result2$a[,2])[-1,]
head(result)
## Nile Fecha predic.ser salto
## 2 1160 1872 1120.000 1097.75
## 3 963 1873 1146.667 1097.75
## 4 1210 1874 1031.875 1097.75
## 5 1160 1875 1142.143 1097.75
## 6 1160 1876 1153.182 1097.75
## 7 813 1877 1157.396 1097.75
<- result %>% gather('Nile','predic.ser', 'salto', key="Serie", value="Caudal")
result head(result)
## Fecha Serie Caudal
## 1 1872 Nile 1160
## 2 1873 Nile 963
## 3 1874 Nile 1210
## 4 1875 Nile 1160
## 5 1876 Nile 1160
## 6 1877 Nile 813
ggplot(result, aes(x = Fecha, y = Caudal, colour = Serie)) +
geom_line() +
labs(title = "Efecto construcción presa de Asswan", x = "Año", y = "Caudal")
¡Observad que estamos filtrando con las varianzas de los elementos del vector de estado y de la perturbación de observación iguales a 1, salvo la del efecto intervención, algo del todo arbitrario! Este es un ejemplo muy artificial para mostrar cómo invocar las funciones de construcción de modelos y de filtrado.
En sesiones sucesivas examinaremos cómo estimar estos parámetros, que pueden ser conocidos del ingeniero en aplicaciones de diseño, pero no lo son en general en otros casos.
dlm
Hay otros paquetes en R que facilitan el uso y estimación de modelos en espacio de estado mediante el filtro de Kalman. Alguno aparecerá más adelante. Entre otros existen los que se relacionan a continuación:
Seguramente el más elaborado y activamente mantenido. Emplea filtrado secuencial y admite falta parcial de datos, ventajas que examinaremos en otra sesión. También permite especificar distribuciones a priori difusas sobre los parámetros (a diferencia de dlm
, en que recurrimos a la aproximación de emplear varianzas iniciales muy grandes).
Implementa también estimación de modelos generalizados (con distribuciones no gaussianas)
Muy limitado, pero muy rápido. Hay una secuela en el paquete FKF.SP
que emplea también procesado secuencial. Hay una viñeta informativa sobre las diferencias.
Facilita la estimación de modelos con cambio de régimen y modelos en tiempo continuo.
Son paquetes que proporcionan la misma funcionalidad que dlm
o KFAS
pero intentan proporcionar un uso más cómodo.