Ejemplos varios (I)
F. Tusell
Esta versión generada el: 2021-11-18
Ejemplo de benchmarking (conciliación de series con intervalos de observación diferentes, de las que una es la serie patrón)
Paro registrado vs. paro encuestado
Leemos primero datos del paro registrado (mensuales):
library(zoo) ##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericsystem(command="head -n 10 ParoReg.txt", intern=TRUE)## [1] " BE_24_15_3" "" "1939:09 632510"
## [4] "1939:10 608040" "1939:11 464140" "1939:12 441388"
## [7] "1940:01 479338" "1940:02 486834" "1940:03 474129"
## [10] "1940:04 462540"temp <- read.table(file="ParoReg.txt",skip=2,header=FALSE,dec=",")
head(temp)## V1 V2
## 1 1939:09 632510
## 2 1939:10 608040
## 3 1939:11 464140
## 4 1939:12 441388
## 5 1940:01 479338
## 6 1940:02 486834tail(temp)## V1 V2
## 890 2013:10 4811383
## 891 2013:11 4808908
## 892 2013:12 4701338
## 893 2014:01 4814435
## 894 2014:02 4812486
## 895 2014:03 4795866ParoReg <- zoo(temp[,2],order.by=as.yearmon(1939 + 8/12 + seq(0,894)/12))
head(ParoReg)## sep 1939 oct 1939 nov 1939 dic 1939 ene 1940 feb 1940
## 632510 608040 464140 441388 479338 486834tail(ParoReg)## oct 2013 nov 2013 dic 2013 ene 2014 feb 2014 mar 2014
## 4811383 4808908 4701338 4814435 4812486 4795866Ahora los de paro EPA, procedentes del INE.
system("head -n 10 ParadosEPA.csv", intern=TRUE)## [1] "Encuesta de Población Activa"
## [2] "Parados"
## [3] "Parados por sexo y grupo de edad. Valores absolutos y porcentajes respecto del total de cada sexo"
## [4] "Ambos sexos;Valor absoluto;2014T3;2014T2;2014T1;2013T4;2013T3;2013T2;2013T1;2012T4;2012T3;2012T2;2012T1;2011T4;2011T3;2011T2;2011T1;2010T4;2010T3;2010T2;2010T1;2009T4;2009T3;2009T2;2009T1;2008T4;2008T3;2008T2;2008T1;2007T4;2007T3;2007T2;2007T1;2006T4;2006T3;2006T2;2006T1;2005T4;2005T3;2005T2;2005T1;2004T4;2004T3;2004T2;2004T1;2003T4;2003T3;2003T2;2003T1;2002T4;2002T3;2002T2;2002T1;Total;"
## [5] "5.427,7;5.622,9;5.933,3;5.935,6;5.943,4;6.047,3;6.278,2;6.021,0;5.824,2;5.731,0;5.667,9;5.287,3;4.998,0;4.844,2;4.921,2;4.702,2;4.585,4;4.655,3;4.617,7;4.335,0;4.121,4;4.139,6;4.018,2;3.206,8;2.600,7;2.385,7;2.190,5;1.942,0;1.806,2;1.773,2;1.863,2;1.819,4;1.766,9;1.834,4;1.942,8;1.860,3;1.783,5;1.969,1;2.121,3;2.176,9;2.199,8;2.247,6;2.309,8;2.276,7;2.247,5;2.216,0;2.328,5;2.232,4;2.196,0;2.103,3;2.152,8;"Un poco difícil de leer:
temp <- read.csv(file="ParadosEPA.csv",skip=4,header=FALSE,dec=",",sep=";",stringsAsFactors=FALSE)
head(temp)## V1 V2 V3 V4 V5 V6 V7 V8 V9
## 1 5.427,7 5.622,9 5.933,3 5.935,6 5.943,4 6.047,3 6.278,2 6.021,0 5.824,2
## V10 V11 V12 V13 V14 V15 V16 V17 V18
## 1 5.731,0 5.667,9 5.287,3 4.998,0 4.844,2 4.921,2 4.702,2 4.585,4 4.655,3
## V19 V20 V21 V22 V23 V24 V25 V26 V27
## 1 4.617,7 4.335,0 4.121,4 4.139,6 4.018,2 3.206,8 2.600,7 2.385,7 2.190,5
## V28 V29 V30 V31 V32 V33 V34 V35 V36
## 1 1.942,0 1.806,2 1.773,2 1.863,2 1.819,4 1.766,9 1.834,4 1.942,8 1.860,3
## V37 V38 V39 V40 V41 V42 V43 V44 V45
## 1 1.783,5 1.969,1 2.121,3 2.176,9 2.199,8 2.247,6 2.309,8 2.276,7 2.247,5
## V46 V47 V48 V49 V50 V51 V52
## 1 2.216,0 2.328,5 2.232,4 2.196,0 2.103,3 2.152,8 NANotad las contorsiones que es preciso realizar con los puntos y comas:
temp <- gsub("\\.","",temp)
temp <- gsub("\\,",".",temp)
temp <- as.numeric(temp[-length(temp)])
head(temp)## [1] 5427.7 5622.9 5933.3 5935.6 5943.4 6047.3tail(temp)## [1] 2216.0 2328.5 2232.4 2196.0 2103.3 2152.8length(temp)## [1] 51ParoEPA <- zoo(rev(temp),order.by=as.yearqtr(2002 + seq(0,50)/4))
head(ParoEPA)## 2002 Q1 2002 Q2 2002 Q3 2002 Q4 2003 Q1 2003 Q2
## 2152.8 2103.3 2196.0 2232.4 2328.5 2216.0tail(ParoEPA)## 2013 Q2 2013 Q3 2013 Q4 2014 Q1 2014 Q2 2014 Q3
## 6047.3 5943.4 5935.6 5933.3 5622.9 5427.7Observemos que el fechado de las dos series de paro es inconsistente:
head(ParoReg)## sep 1939 oct 1939 nov 1939 dic 1939 ene 1940 feb 1940
## 632510 608040 464140 441388 479338 486834head(ParoEPA)## 2002 Q1 2002 Q2 2002 Q3 2002 Q4 2003 Q1 2003 Q2
## 2152.8 2103.3 2196.0 2232.4 2328.5 2216.0Ahora dataremos cada periodo en el último día para poder alinear las series. (¡Ojo! La función as.Date utilizada aquí es la del paquete zoo, que enmascara a la del mismo nombre en R base. Para resolver cualquier ambigüedad podemos emplear el operador :: como en zoo::as.Date, base::as.Date.)
index(ParoReg) <- as.Date(index(ParoReg), frac=1)
index(ParoEPA) <- as.Date(index(ParoEPA), frac=1) frac=1 significa: “Alinea con el día que corresponde a dejar pasar la totalidad del periodo”. Esto alinea con el último dia del mes, del trimestre, etc. Observemos ahora:
head(ParoReg)## 1939-09-30 1939-10-31 1939-11-30 1939-12-31 1940-01-31 1940-02-29
## 632510 608040 464140 441388 479338 486834head(ParoEPA)## 2002-03-31 2002-06-30 2002-09-30 2002-12-31 2003-03-31 2003-06-30
## 2152.8 2103.3 2196.0 2232.4 2328.5 2216.0Ahora va a haber fechas identicas en ambas series que van a permitir el empalme (aunque no todas las fechas de una serie aparezcan en la otra). ¡SOLO el poder hacer uso de la línea a continuación justifica el trabajo que nos hemos tomado en fechar!
paro <- cbind(ParoReg/1000,ParoEPA)
paro <- window(paro,start="2002-03-31",end="2013-07-31")
plot(na.approx(paro),screens=1,col=1:2)
paro <- log(paro)
plot(na.approx(paro),screens=1,col=1:2)
Creemos un modelo bivariante:
library(dlm)
crearMod <- function(p) {
ep <- exp(p)
mod1 <- dlmModPoly(1)
mod2 <- dlmModPoly(1)
mod <- mod1 %+% mod2
mod$FF <- matrix(c(1,1),2,1)
mod$GG <- matrix(1,1,1)
mod$W <- matrix(ep[1],1,1)
mod$V <- diag(rep(ep[2],2))
return(mod)
}
modelo <- crearMod(rep(-5,2))
modelo## $FF
## [,1]
## [1,] 1
## [2,] 1
##
## $V
## [,1] [,2]
## [1,] 0.006737947 0.000000000
## [2,] 0.000000000 0.006737947
##
## $GG
## [,1]
## [1,] 1
##
## $W
## [,1]
## [1,] 0.006737947
##
## $m0
## [1] 0 0
##
## $C0
## [,1] [,2]
## [1,] 1e+07 0e+00
## [2,] 0e+00 1e+07¡ %+% es diferente de + !
estim <- dlmMLE(paro,c(0,0),build=crearMod)
modelo <- crearMod(estim$par)
#
dlmFilter(paro,modelo)$m[1:6]## 2002-03-31 2002-04-30 2002-05-31 2002-06-30 2002-07-31
## 0.000000 7.658069 7.645967 7.627698 7.621299 7.606837paro[1:6,]## ParoReg/1000 ParoEPA
## 2002-03-31 7.641614 7.674525
## 2002-04-30 7.630495 NA
## 2002-05-31 7.602363 NA
## 2002-06-30 7.582210 7.651263
## 2002-07-31 7.581644 NA
## 2002-08-31 7.592861 NAsuavizado.dlm <- dlmSmooth(paro,modelo)$s[-1]require(KFAS)## Loading required package: KFAS## Please cite KFAS in publications by using:
##
## Jouni Helske (2017). KFAS: Exponential Family State Space Models in R. Journal of Statistical Software, 78(10), 1-39. doi:10.18637/jss.v078.i10.mod1 <- SSModel(paro ~
-1 # Sin el nivel "de oficio"
+ SSMcustom(
Z=matrix(c(1,1), 2, 1), # Matriz de observación
T=diag(1), # Matriz de transición
Q=diag(1), # Matriz covarianzas del estado
P1inf=matrix(1,1,1) # Priori difusa
),
H=diag(2) # Matriz de covarianzas obs
)Examinemos la estructura del objeto resultante:
mod1## Call:
## SSModel(formula = paro ~ -1 + SSMcustom(Z = matrix(c(1, 1), 2,
## 1), T = diag(1), Q = diag(1), P1inf = matrix(1, 1, 1)), H = diag(2))
##
## State space model object of class SSModel
##
## Dimensions:
## [1] Number of time points: 137
## [1] Number of time series: 2
## [1] Number of disturbances: 1
## [1] Number of states: 1
## Names of the states:
## [1] custom1
## Distributions of the time series:
## [1] gaussian
##
## Object is a valid object of class SSModel.mod1$Z## , , 1
##
## custom1
## ParoReg/1000 1
## ParoEPA 1mod1$T## , , 1
##
## custom1
## custom1 1mod1$Q## , , 1
##
## [,1]
## [1,] 1mod1$H## , , 1
##
## [,1] [,2]
## [1,] 1 0
## [2,] 0 1Para ajustar el modelo, necesitamos la updatefn más abajo, cuyo cometido es reemplazar valores de los parametros en lugares específicos de las matrices que hay que estimar. Nuestro modelo tiene dos parámetros, una varianza en Q y dos varianzas en H que supondremos iguales.
Obsérvese que estamos imponiendo que las dos varianzas de la ecuación de observación sean iguales, lo que tiene sentido si les atribuimos la misma precisión: no tiene porqué ser el caso.
updatefn <- function(pars, model) {
ep <- exp(pars)
model["Q"] <- matrix(ep[1],1,1)
model["H"] <- diag(c(ep[2],ep[2]))
return(model)
}Una vez definida esta función, proporcionamos valores iniciales para los dos parámetros en init e invocamos fitSSM.
init <- rep(0.1, 2)
fit <- fitSSM(mod1, inits=init, updatefn=updatefn)Los resultados se obtienen en fit cuyo contenido es:
fit## $optim.out
## $optim.out$par
## [1] -6.935134 -5.670112
##
## $optim.out$value
## [1] -216.5231
##
## $optim.out$counts
## function gradient
## 115 NA
##
## $optim.out$convergence
## [1] 0
##
## $optim.out$message
## NULL
##
##
## $model
## Call:
## SSModel(formula = paro ~ -1 + SSMcustom(Z = matrix(c(1, 1), 2,
## 1), T = diag(1), Q = diag(1), P1inf = matrix(1, 1, 1)), H = diag(2))
##
## State space model object of class SSModel
##
## Dimensions:
## [1] Number of time points: 137
## [1] Number of time series: 2
## [1] Number of disturbances: 1
## [1] Number of states: 1
## Names of the states:
## [1] custom1
## Distributions of the time series:
## [1] gaussian
##
## Object is a valid object of class SSModel.En $optim.out$value vemos el valor de la verosimilitud maximizada y en $optim.out$par los valores “brutos” de los parámetros (sin hacer las transformaciones indicadas más arriba). Podemos reemplazarlos en el modelo mediante,
mod1 <- updatefn(pars = fit$optim.out$par, model = mod1)
mod1$H## , , 1
##
## [,1] [,2]
## [1,] 0.003447479 0.000000000
## [2,] 0.000000000 0.003447479mod1$Q## , , 1
##
## [,1]
## [1,] 0.0009729929Con los valores de los parámetros ya estimados podemos correr el filtro de Kalman para obtener los valores suavizados del estado:
mod1.out <- KFS(mod1)En mod1.out tenemos ahora bastante información, que incluye los valores suavizados del estado (en el componente alphahat).
str(mod1.out$V)## num [1, 1, 1:137] 0.00099 0.0009 0.000832 0.000702 0.000803 ...suavizado.kfas <- mod1.out$alphahatPodemos comparar con el resultado de dlm. La inicialización aproximada de dlm y la exacta de KFAS proporcionan resultados prácticamente indistinguibles.
cbind(suavizado.kfas,suavizado.dlm)## Time Series:
## Start = 1
## End = 137
## Frequency = 1
## suavizado.kfas suavizado.dlm
## 1 7.641083 7.641083
## 2 7.631495 7.631495
## 3 7.622190 7.622190
## 4 7.618479 7.618480
## 5 7.615753 7.615754
## 6 7.622653 7.622654
## 7 7.637961 7.637962
## 8 7.646846 7.646846
## 9 7.659738 7.659738
## 10 7.673170 7.673169
## 11 7.678990 7.678989
## 12 7.681860 7.681859
## 13 7.683229 7.683228
## 14 7.665970 7.665970
## 15 7.652705 7.652706
## 16 7.649080 7.649081
## 17 7.640845 7.640846
## 18 7.644453 7.644455
## 19 7.658010 7.658010
## 20 7.665335 7.665335
## 21 7.677532 7.677532
## 22 7.691838 7.691837
## 23 7.696418 7.696417
## 24 7.696911 7.696910
## 25 7.695135 7.695134
## 26 7.681382 7.681382
## 27 7.668308 7.668308
## 28 7.661740 7.661741
## 29 7.649038 7.649039
## 30 7.647921 7.647922
## 31 7.653156 7.653156
## 32 7.653972 7.653972
## 33 7.659265 7.659265
## 34 7.664441 7.664440
## 35 7.666086 7.666085
## 36 7.662248 7.662247
## 37 7.653295 7.653294
## 38 7.637564 7.637564
## 39 7.618956 7.618957
## 40 7.604402 7.604403
## 41 7.599789 7.599790
## 42 7.596358 7.596359
## 43 7.588961 7.588963
## 44 7.605291 7.605291
## 45 7.615496 7.615495
## 46 7.616644 7.616644
## 47 7.632950 7.632948
## 48 7.635081 7.635079
## 49 7.623927 7.623926
## 50 7.613741 7.613740
## 51 7.596696 7.596696
## 52 7.577827 7.577828
## 53 7.576062 7.576062
## 54 7.573711 7.573712
## 55 7.565999 7.566001
## 56 7.578376 7.578376
## 57 7.585407 7.585407
## 58 7.584816 7.584816
## 59 7.598645 7.598644
## 60 7.600427 7.600426
## 61 7.591648 7.591647
## 62 7.589374 7.589373
## 63 7.580601 7.580602
## 64 7.569903 7.569904
## 65 7.580533 7.580534
## 66 7.589633 7.589633
## 67 7.591586 7.591587
## 68 7.614606 7.614607
## 69 7.634721 7.634722
## 70 7.651355 7.651356
## 71 7.687060 7.687059
## 72 7.712347 7.712346
## 73 7.727766 7.727766
## 74 7.749553 7.749553
## 75 7.769161 7.769161
## 76 7.790312 7.790313
## 77 7.818280 7.818281
## 78 7.852993 7.852995
## 79 7.892510 7.892512
## 80 7.945718 7.945719
## 81 7.999469 7.999469
## 82 8.052285 8.052284
## 83 8.100327 8.100325
## 84 8.145619 8.145617
## 85 8.188173 8.188170
## 86 8.198995 8.198993
## 87 8.209229 8.209228
## 88 8.223685 8.223684
## 89 8.221244 8.221244
## 90 8.232407 8.232408
## 91 8.253639 8.253639
## 92 8.264905 8.264905
## 93 8.281803 8.281803
## 94 8.304646 8.304645
## 95 8.316213 8.316213
## 96 8.330635 8.330634
## 97 8.346313 8.346312
## 98 8.339446 8.339445
## 99 8.335516 8.335516
## 100 8.338656 8.338656
## 101 8.325404 8.325405
## 102 8.327527 8.327528
## 103 8.341247 8.341247
## 104 8.341810 8.341810
## 105 8.349850 8.349850
## 106 8.365962 8.365961
## 107 8.370044 8.370044
## 108 8.379729 8.379729
## 109 8.393233 8.393232
## 110 8.381614 8.381615
## 111 8.376317 8.376317
## 112 8.381163 8.381163
## 113 8.372672 8.372673
## 114 8.380799 8.380801
## 115 8.404320 8.404320
## 116 8.411657 8.411658
## 117 8.427805 8.427806
## 118 8.453495 8.453494
## 119 8.462108 8.462108
## 120 8.478718 8.478718
## 121 8.501207 8.501205
## 122 8.493710 8.493710
## 123 8.494406 8.494406
## 124 8.505287 8.505286
## 125 8.493533 8.493534
## 126 8.499405 8.499407
## 127 8.522222 8.522222
## 128 8.521959 8.521960
## 129 8.532599 8.532600
## 130 8.552839 8.552839
## 131 8.549428 8.549428
## 132 8.556202 8.556202
## 133 8.571724 8.571722
## 134 8.551795 8.551796
## 135 8.542243 8.542244
## 136 8.545986 8.545986
## 137 8.525974 8.525975Podemos representar las observaciones y el valor suavizado del estado junto con su intervalo de confianza 95%.
plot(paro[,1], xlab="Fecha", ylab="Paro (log)", type="p",
cex=0.5, col="green",
main="Paro registrado y EPA")
legend(x="bottomright",legend=c("Registrado","EPA",expression(alpha[t])),
col=c("green","blue","black"), pch=c(1,1,NA), pt.cex=0.5,
seg.len=1, lty=c(NA,NA,1), cex=0.9, lwd=c(NA,NA,2))
points(paro[,2], col="blue", cex=0.5)
lines(zoo(mod1.out$alphahat, order.by=index(paro)), lwd=2)
ciw <- 1.96*sqrt(mod1.out$V[1,1,])
lines(zoo(mod1.out$alphahat + ciw, order.by=index(paro)), lty=3)
lines(zoo(mod1.out$alphahat - ciw, order.by=index(paro)), lty=3)
PIB y consumo total de energia
Lectura datos PIB (del Banco de España):
temp <- read.table(file="PIB.txt",skip=2,header=FALSE,dec=",")
head(temp)## V1 V2
## 1 1995:1 157.257
## 2 1995:2 158.022
## 3 1995:3 158.572
## 4 1995:4 159.538
## 5 1996:1 160.417
## 6 1996:2 161.615tail(temp)## V1 V2
## 71 2012:3 232.995
## 72 2012:4 231.201
## 73 2013:1 230.480
## 74 2013:2 230.169
## 75 2013:3 230.345
## 76 2013:4 230.747PIB <- zoo(temp[,2],order.by=as.yearqtr(1995 + seq(0,75)/4))
head(PIB)## 1995 Q1 1995 Q2 1995 Q3 1995 Q4 1996 Q1 1996 Q2
## 157.257 158.022 158.572 159.538 160.417 161.615tail(PIB)## 2012 Q3 2012 Q4 2013 Q1 2013 Q2 2013 Q3 2013 Q4
## 232.995 231.201 230.480 230.169 230.345 230.747A continuación, datos de producción de energía eléctrica mensuales (del Banco de España):
temp <- read.table(file="elec.txt",skip=2,header=FALSE,dec=",")
head(temp)## V1 V2
## 1 1977:01 7863
## 2 1977:02 7147
## 3 1977:03 7452
## 4 1977:04 6726
## 5 1977:05 6857
## 6 1977:06 6541tail(temp)## V1 V2
## 440 2013:08 22046.6
## 441 2013:09 20937.2
## 442 2013:10 21120.9
## 443 2013:11 21651.0
## 444 2013:12 22886.6
## 445 2014:01 23118.4ELEC <- zoo(temp[,2],order.by=as.yearmon(1977 + seq(0,444)/12))
head(ELEC)## ene 1977 feb 1977 mar 1977 abr 1977 may 1977 jun 1977
## 7863 7147 7452 6726 6857 6541tail(ELEC)## ago 2013 sep 2013 oct 2013 nov 2013 dic 2013 ene 2014
## 22046.6 20937.2 21120.9 21651.0 22886.6 23118.4Tenemos ahora que alinear las series. Podemos fecharlas el último día de cada periodo:
index(PIB) <- as.Date(index(PIB), frac=1)
index(ELEC) <- as.Date(index(ELEC), frac=1)
head(PIB)## 1995-03-31 1995-06-30 1995-09-30 1995-12-31 1996-03-31 1996-06-30
## 157.257 158.022 158.572 159.538 160.417 161.615head(ELEC)## 1977-01-31 1977-02-28 1977-03-31 1977-04-30 1977-05-31 1977-06-30
## 7863 7147 7452 6726 6857 6541Al tener un fechado homogéneo, podemos emparejarlas y se alinearán automáticamente:
datos <- cbind(100*PIB,ELEC)
tail(datos)## 100 * PIB ELEC
## 2013-08-31 NA 22046.6
## 2013-09-30 23034.5 20937.2
## 2013-10-31 NA 21120.9
## 2013-11-30 NA 21651.0
## 2013-12-31 23074.7 22886.6
## 2014-01-31 NA 23118.4plot(datos, log="y")
plot(na.approx(datos), col=1:2, main="Interpolando datos")
plot(na.approx(datos),screen=1,col=1:2, main="PIB y energía bruta",ylab="Escala log",xlab="Fecha")
datos <- window(datos, start=as.Date("1995-01-31"), end=as.Date("2014-01-31"))
plot(na.approx(datos),screen=1,col=1:2, main="PIB y energía bruta",xlab="Fecha")
Podríamos ahora proceder como en el caso del paro registrado y EPA. ¡Observad! El nivel de las series no es (ni tiene porque ser) homogéneo. El multiplicar por 100 una permite representarlas en una escala numérica común, pero todavía tenemos que dar cuenta de su diferente nivel (en media, la serie “negra” por encima de la “roja”). Una posibilidad podría ser incluir un parámetro en la ecuación de observación que ajustara los niveles respectivos:
Modelos en espacio de estado generalizados
Leamos primero los datos ya preparados:
load("Pseudomona.rda") # Contiene una dataframe de nombre 'datos'.
datos## N IPM Tiempo DDDe DDS
## 1998 Q1 27 25 1 2.5035765 0
## 1998 Q2 27 25 2 14.9647887 0
## 1998 Q3 27 25 3 20.5081001 0
## 1998 Q4 27 25 4 19.0707351 0
## 1999 Q1 47 42 5 13.9973082 0
## 1999 Q2 47 42 6 13.5743520 0
## 1999 Q3 47 42 7 3.9556962 0
## 1999 Q4 47 42 8 10.2005533 0
## 2000 Q1 7 6 9 14.8351648 0
## 2000 Q2 5 4 10 1.6430412 0
## 2000 Q3 2 1 11 28.3512064 0
## 2000 Q4 7 7 12 9.5737483 0
## 2001 Q1 7 6 13 23.0658199 0
## 2001 Q2 12 9 14 7.6243455 0
## 2001 Q3 11 9 15 10.7898010 0
## 2001 Q4 15 9 16 9.7634763 0
## 2002 Q1 6 5 17 5.8291771 0
## 2002 Q2 4 2 18 6.2586926 1
## 2002 Q3 8 5 19 3.7323944 1
## 2002 Q4 6 5 20 2.8236915 1
## 2003 Q1 8 6 21 1.8839713 1
## 2003 Q2 7 7 22 0.0000000 1
## 2003 Q3 3 3 23 2.1084337 1
## 2003 Q4 6 6 24 7.0175439 1
## 2004 Q1 6 5 25 6.4862104 1
## 2004 Q2 7 5 26 7.6923077 1
## 2004 Q3 6 4 27 9.3608597 1
## 2004 Q4 7 5 28 6.6192560 1
## 2005 Q1 5 4 29 6.8983402 1
## 2005 Q2 10 8 30 2.2522523 1
## 2005 Q3 5 5 31 18.6969697 1
## 2005 Q4 3 3 32 5.2356021 1
## 2006 Q1 5 3 33 42.7612994 1
## 2006 Q2 6 5 34 7.4417010 1
## 2006 Q3 7 1 35 12.5282167 1
## 2006 Q4 8 3 36 5.2492047 1
## 2007 Q1 5 1 37 3.9043584 1
## 2007 Q2 4 1 38 3.9874082 1
## 2007 Q3 6 4 39 12.1268657 1
## 2007 Q4 10 3 40 4.1832669 1
## 2008 Q1 8 4 41 7.0106908 1
## 2008 Q2 10 3 42 5.7588076 1
## 2008 Q3 16 5 43 1.2012012 1
## 2008 Q4 9 4 44 1.4065817 1
## 2009 Q1 4 2 45 1.4317181 1
## 2009 Q2 4 2 46 1.8936170 1
## 2009 Q3 16 7 47 1.2486127 1
## 2009 Q4 11 6 48 1.8671193 1
## 2010 Q1 4 3 49 0.8970727 1
## 2010 Q2 6 4 50 11.0146252 1
## 2010 Q3 10 5 51 2.5945471 1
## 2010 Q4 7 5 52 1.5903308 1
## 2011 Q1 7 3 53 0.8354756 1
## 2011 Q2 8 5 54 0.0000000 1
## 2011 Q3 13 8 55 0.1062925 1
## 2011 Q4 12 8 56 1.3605442 1
## 2012 Q1 4 2 57 0.7207578 1
## 2012 Q2 11 2 58 0.9647651 1
## 2012 Q3 18 11 59 1.8041237 1
## 2012 Q4 9 4 60 0.5859375 1
## 2013 Q1 4 2 61 0.0000000 1
## 2013 Q2 6 3 62 1.1585994 1
## 2013 Q3 1 1 63 0.0000000 1
## 2013 Q4 6 4 64 2.0360825 1A continuación especificamos el modelo deseado, con ayuda de funciones constructoras disponibles en el paquete KFAS. En los lugares en que han de aparecer algunos parámetros ponemos valores arbitrarios (ceros); luego serán rellenados por los valores estimados.
En u tenemos sujetos en riesgo. El modelo que sigue especifica que la sensibilidad evoluciona como un nivel aleatorio mas partes explicadas por la DDS y DDDe –que supondremos fijas–. La ecuación de observación ya no es la conocida función lineal a que estamos habituados, sino que lo que tenemos es que la observación se distribuye como Poisson con \(\lambda\) igual a la sensibilidad multiplicada por \(N\) (número de sujetos en riesgo).
f <- formula(~ -1 + DDDe + DDS)
mod2<- SSModel(IPM ~ 1
+ SSMregression(rformula = f,
Q = matrix(0, 2, 2),
P1inf=diag(2),
data=datos)
+ SSMtrend(degree = 1,
Q = list(matrix(0, 1, 1)),
P1inf = matrix(1, 1, 1)) ,
distribution="poisson", u = datos[,"N"], data=datos
)
str(mod2)## List of 14
## $ y : Time-Series [1:64, 1] from 1 to 64: 25 25 25 25 42 42 42 42 6 4 ...
## $ Z : num [1, 1:3, 1:64] 2.5 0 1 15 0 ...
## ..- attr(*, "dimnames")=List of 3
## .. ..$ : NULL
## .. ..$ : chr [1:3] "DDDe" "DDS" "level"
## .. ..$ : NULL
## $ H : chr "Omitted"
## $ T : num [1:3, 1:3, 1] 1 0 0 0 1 0 0 0 1
## ..- attr(*, "dimnames")=List of 3
## .. ..$ : chr [1:3] "DDDe" "DDS" "level"
## .. ..$ : chr [1:3] "DDDe" "DDS" "level"
## .. ..$ : NULL
## $ R : num [1:3, 1:3, 1] 1 0 0 0 1 0 0 0 1
## ..- attr(*, "dimnames")=List of 3
## .. ..$ : chr [1:3] "DDDe" "DDS" "level"
## .. ..$ : NULL
## .. ..$ : NULL
## $ Q : num [1:3, 1:3, 1] 0 0 0 0 0 0 0 0 0
## $ a1 : num [1:3, 1] 0 0 0
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:3] "DDDe" "DDS" "level"
## .. ..$ : NULL
## $ P1 : num [1:3, 1:3] 0 0 0 0 0 0 0 0 0
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:3] "DDDe" "DDS" "level"
## .. ..$ : chr [1:3] "DDDe" "DDS" "level"
## $ P1inf : num [1:3, 1:3] 1 0 0 0 1 0 0 0 1
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:3] "DDDe" "DDS" "level"
## .. ..$ : chr [1:3] "DDDe" "DDS" "level"
## $ u : Time-Series [1:64, 1] from 1 to 64: 27 27 27 27 47 47 47 47 7 5 ...
## $ distribution: chr "poisson"
## $ tol : num 1.49e-08
## $ call : language SSModel(formula = IPM ~ 1 + SSMregression(rformula = f, Q = matrix(0, 2, 2), P1inf = diag(2), data = datos) | __truncated__ ...
## $ terms :Classes 'terms', 'formula' language IPM ~ 1 + SSMregression(rformula = f, Q = matrix(0, 2, 2), P1inf = diag(2), data = datos) + SSMtrend(degree | __truncated__ ...
## .. ..- attr(*, "variables")= language list(IPM, SSMregression(rformula = f, Q = matrix(0, 2, 2), P1inf = diag(2), data = datos), SSMtrend(degree =| __truncated__ ...
## .. ..- attr(*, "factors")= int [1:3, 1:2] 0 1 0 0 0 1
## .. .. ..- attr(*, "dimnames")=List of 2
## .. .. .. ..$ : chr [1:3] "IPM" "SSMregression(rformula = f, Q = matrix(0, 2, 2), P1inf = diag(2), data = datos)" "SSMtrend(degree = 1, Q = list(matrix(0, 1, 1)), P1inf = matrix(1, 1, 1))"
## .. .. .. ..$ : chr [1:2] "SSMregression(rformula = f, Q = matrix(0, 2, 2), P1inf = diag(2), data = datos)" "SSMtrend(degree = 1, Q = list(matrix(0, 1, 1)), P1inf = matrix(1, 1, 1))"
## .. ..- attr(*, "term.labels")= chr [1:2] "SSMregression(rformula = f, Q = matrix(0, 2, 2), P1inf = diag(2), data = datos)" "SSMtrend(degree = 1, Q = list(matrix(0, 1, 1)), P1inf = matrix(1, 1, 1))"
## .. ..- attr(*, "specials")=Dotted pair list of 6
## .. .. ..$ SSMregression: int 2
## .. .. ..$ SSMtrend : int 3
## .. .. ..$ SSMseasonal : NULL
## .. .. ..$ SSMcycle : NULL
## .. .. ..$ SSMarima : NULL
## .. .. ..$ SSMcustom : NULL
## .. ..- attr(*, "order")= int [1:2] 1 1
## .. ..- attr(*, "intercept")= int 1
## .. ..- attr(*, "response")= int 1
## .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
## - attr(*, "class")= chr "SSModel"
## - attr(*, "p")= int 1
## - attr(*, "m")= int 3
## - attr(*, "k")= int 3
## - attr(*, "n")= int 64
## - attr(*, "tv")= int [1:5] 1 1 0 0 0
## - attr(*, "state_types")= chr [1:3] "regression" "regression" "level"
## - attr(*, "eta_types")= chr [1:3] "regression" "regression" "level"Esta es la matriz de observación (temporalmente variable):
mod2$Z[,,1:5]## [,1] [,2] [,3] [,4] [,5]
## DDDe 2.503577 14.96479 20.5081 19.07074 13.99731
## DDS 0.000000 0.00000 0.0000 0.00000 0.00000
## level 1.000000 1.00000 1.0000 1.00000 1.00000Matriz de covarianzas de las observaciones, no hay.
mod2$H## [1] "Omitted"Matriz de covarianzas de los elementos del vector de estado:
mod2$Q## , , 1
##
## [,1] [,2] [,3]
## [1,] 0 0 0
## [2,] 0 0 0
## [3,] 0 0 0Matriz de transición:
mod2$T## , , 1
##
## DDDe DDS level
## DDDe 1 0 0
## DDS 0 1 0
## level 0 0 1Definimos una función updatefn para reemplazar los valores estimados de los parámetros en sus ubicaciones.
updatefn <- function(pars, model, ...) {
Q <- diag(c(0,0,exp(pars[1])))
model["Q"] <- Q
return(model)
}Proporcionamos valores iniciales al único parámetro del modelo y llevamos a cabo la estimación. Es mucho más onerosa que en el caso lineal gaussiano, que el filtro de Kalman maneja directamente; ahora la estimación se hace por simulación de modelos que aproximan el que deseamos (detalles en Durbin-Koopman, cap. 11).
init <- c(-6)
fit <- fitSSM(model=mod2, inits=init, updatefn=updatefn, method="BFGS", tol=10^-4, nsim=100)El resultado del ajuste puede verse en fit; contiene los valores de los parámetros (en $optim.out$par), el valor optimizado de la verosimilitud (en $optim.out$value), un diagnóstico sobre la convergencia (optim.out$convergence: valor cero significa que todo fue bien) e información acerca del modelo.
fit## $optim.out
## $optim.out$par
## [1] -5.913961
##
## $optim.out$value
## [1] 134.4912
##
## $optim.out$counts
## function gradient
## 13 3
##
## $optim.out$convergence
## [1] 0
##
## $optim.out$message
## NULL
##
##
## $model
## Call:
## SSModel(formula = IPM ~ 1 + SSMregression(rformula = f, Q = matrix(0,
## 2, 2), P1inf = diag(2), data = datos) + SSMtrend(degree = 1,
## Q = list(matrix(0, 1, 1)), P1inf = matrix(1, 1, 1)), data = datos,
## u = datos[, "N"], distribution = "poisson")
##
## State space model object of class SSModel
##
## Dimensions:
## [1] Number of time points: 64
## [1] Number of time series: 1
## [1] Number of disturbances: 3
## [1] Number of states: 3
## Names of the states:
## [1] DDDe DDS level
## Distributions of the time series:
## [1] poisson
##
## Object is a valid object of class SSModel.Podemos ahora reemplazar los valores estimados de los parámetros en el modelo y correr el filtro de Kalman:
mod2 <- updatefn(pars=fit$optim.out$par, model=mod2)
mod2.out <- KFS(mod2)En mod2.out$a tenemos los valores filtrados, en mod2.out$alphahat valores suavizados y en mod2.out$muhat valores esperados de la variable respuesta (\(N_t\lambda_t\), valor esperado de sensibles en el periodo \(t\)). Se imprimen los primeros valores de las respectivas series.
head(mod2.out$a)## NULLhead(mod2.out$alphahat)## DDDe DDS level
## [1,] 0.001250509 -0.1343691 -0.1117110
## [2,] 0.001250509 -0.1343691 -0.1138130
## [3,] 0.001250509 -0.1343691 -0.1171291
## [4,] 0.001250509 -0.1343691 -0.1214188
## [5,] 0.001250509 -0.1343691 -0.1270862
## [6,] 0.001250509 -0.1343691 -0.1324245head(cbind(mod2.out$model$y,mod2.out$muhat))## mod2.out$model$y mod2.out$muhat
## [1,] 25 24.22189
## [2,] 25 24.55063
## [3,] 25 24.63956
## [4,] 25 24.49003
## [5,] 42 42.12179
## [6,] 42 41.87537Para decidir sobre el posible efecto de la DDS, examinaremos su valor suavizado (estimado con todas las observaciones, pasadas, presentes y futuras).
mod2.out$alphahat[,2][1]## [1] -0.1343691y compararemos con su desviación típica
sqrt(mod2.out$V[2,2,])[1]## [1] 0.2117715El efecto sería aproximadamente significativo al nivel 5% si tuviera una magnitud de al menos dos desviaciones típicas, lo que aquí no sucede.
Mostramos a continuación dos gráficos, comparando los sensibles observados y esperados de acuerdo con el modelo y la estimación de la sensibilidad \(\lambda_t\) a lo largo del tiempo, con un intervalo de confianza 95% para la misma.
o <- zoo(mod2.out$model$y, order.by=index(datos)) # Observados
e <- zoo(mod2.out$muhat, order.by=index(datos)) # Esperados
plot(o, xlab="Fecha", ylab="Sensibles", type="p",
cex=0.6, col="green",
main="Sensibles esperados y observados (PseudoAer vs. IPM)")
points(e, cex=0.6, col="black")
legend(x="topright",
legend=c("Sensibles obs.", expression(N[t]*lambda[t])),
col=c("green","black"), pch=c(1,1), pt.cex=0.5)
a <- zoo(datos[,"N"], order.by=index(datos))
sensibilidad <- 100 * (e / a)
plot(sensibilidad, xlab="Fecha", ylab="Sensibles", type="l",
cex=0.6, col="black", lwd=2, ylim=c(0,100),
main="Sensibilidad ajustada y observada (PseudoAer vs. IPM)")
ciw <- 1.96 * sqrt(mod2.out$V_mu[1,1,]) * (100 / a)
lines(sensibilidad + ciw, lty=2)
lines(sensibilidad - ciw, lty=2)
points(100*o/a,col="green",cex=0.5)
abline(v=2002.25,col="red")
text(2002.25,0,"DDS",pos=4,col="red")