(“Estimación de parámetros” = “Identificación de modelos” en la literatura por y para ingenieros. Entre económetras, “identificación” tiene un significado radicalmente diferente: un modelo se dice “identificable” cuando se pueden estimar sus parámetros sin ambigüedad con la información disponible, y no identificable en caso contrario.)

A partir de primeros principios

Modifiquemos la función kalman2 en uno de los ejemplos precedente de manera que calcule la verosimilitud de un modelo univariante de nivel local:

kalman3 <- function(y,x0,P0,Q,H) {
  n <- length(y)
  x <- P <- rep(0,n)
  x[1] <- x0
  P[1] <- P0
  vero <- 0
  for (i in 1:n) {
    if (!is.na(y[i])) {
    inno <- y[i] - x[i]
    F    <- P[i] + H
    K    <- P[i] / F
    x[i+1] <- x[i] + K * inno
    P[i+1] <- P[i]*(1-K) + Q
    vero <- vero + log(F) + (inno^2)  / F
    } else {
      x[i+1] <- x[i]
      P[i+1] <- P[i] + Q
    }
  }
  vero <- -0.5*vero
  return(list(x=x,P=P,vero=vero))
}

Podemos invocarla y obtener, además de los valores del estado y matriz de covarianzas, el valor de la función de verosimilitud (supuesta normalidad de los ruidos).

result   <- kalman3(y=Nile,0,10^7,Q=150,H=15000)
result
## $x
##   [1]    0.0000 1118.3225 1139.2495 1079.5581 1113.2815 1123.1712 1129.8541
##   [8] 1078.9378 1100.9634 1137.2302 1137.5806 1120.4573 1099.1077 1100.3189
##  [15] 1088.8256 1081.5557 1068.9572 1080.2888 1051.9469 1042.5792 1052.2108
##  [22] 1056.9027 1071.8477 1079.4411 1095.9494 1111.7779 1122.1930 1113.3392
##  [29] 1112.0604 1079.6964 1056.7758 1039.3146 1006.3513 1000.0215  984.0962
##  [36]  957.1148  953.1976  928.3191  937.0495  947.8032  949.8209  938.5116
##  [43]  918.2868  874.2947  869.5089  853.5705  878.9200  899.9539  893.4888
##  [50]  881.1698  875.4456  865.2240  863.3001  863.3667  863.2367  847.5180
##  [57]  847.2785  837.4539  833.5105  853.1531  844.1967  838.1850  840.7358
##  [64]  841.1415  850.9259  863.5847  866.7633  862.5052  876.5356  866.4966
##  [71]  848.3756  829.4099  830.9881  829.1818  820.8887  818.9968  840.0197
##  [78]  841.9203  844.9719  845.2599  849.5158  839.4786  830.8719  831.5499
##  [85]  852.3300  858.5768  870.6980  863.6874  869.3295  879.3814  873.2572
##  [92]  887.2161  889.0029  890.1441  916.7654  916.3121  900.1111  901.9079
##  [99]  884.4137  868.2031  856.0078
## 
## $P
##   [1] 10000000.000    15127.534     7681.748     5230.130     4027.975
##   [6]     3325.305     2871.896     2560.401     2337.081     2172.037
##  [11]     2047.303     1951.431     1876.784     1818.076     1771.538
##  [16]     1734.414     1704.654     1680.699     1661.357     1645.698
##  [21]     1632.994     1622.670     1614.268     1607.423     1601.842
##  [26]     1597.287     1593.567     1590.529     1588.045     1586.015
##  [31]     1584.354     1582.996     1581.885     1580.976     1580.232
##  [36]     1579.623     1579.124     1578.716     1578.382     1578.109
##  [41]     1577.885     1577.702     1577.552     1577.429     1577.328
##  [46]     1577.246     1577.178     1577.123     1577.078     1577.041
##  [51]     1577.011     1576.986     1576.966     1576.949     1576.935
##  [56]     1576.924     1576.915     1576.908     1576.901     1576.896
##  [61]     1576.892     1576.889     1576.886     1576.884     1576.882
##  [66]     1576.881     1576.879     1576.878     1576.878     1576.877
##  [71]     1576.876     1576.876     1576.876     1576.875     1576.875
##  [76]     1576.875     1576.875     1576.874     1576.874     1576.874
##  [81]     1576.874     1576.874     1576.874     1576.874     1576.874
##  [86]     1576.874     1576.874     1576.874     1576.874     1576.874
##  [91]     1576.874     1576.874     1576.874     1576.874     1576.874
##  [96]     1576.874     1576.874     1576.874     1576.874     1576.874
## [101]     1576.874
## 
## $vero
## [1] -554.2235

Si en lugar de utilizar valores arbitrarios de P y Q queremos estimarlos, escribimos una función objetivo y la pasamos a (por ejemplo) optim.

objetivo <- function(p) {
  Q <- exp(p[1])                 # varianza no negativa
  H <- exp(p[2])                 # varianza no negativa
  return(kalman3(y=Nile,0,10^10,
                 Q=Q,
                 H=H)$vero)
       }

La optimización numérica pude presentar dificultades. Observad los diferentes resultados. Conviene probar diferentes valores iniciales y/o hacer más rigurosa la condición de término de la búsqueda (mediante reltol)

parest <- optim(par=c(4,5),
                fn=objetivo, control=list(fnscale=-1))
parest
## $par
## [1] 7.297411 9.622523
## 
## $value
## [1] -553.0837
## 
## $counts
## function gradient 
##       63       NA 
## 
## $convergence
## [1] 0
## 
## $message
## NULL
parest <- optim(par=c(4,5), method="L-BFGS",
                fn=objetivo, control=list(fnscale=-1))
parest
## $par
## [1] 7.292454 9.622353
## 
## $value
## [1] -553.0837
## 
## $counts
## function gradient 
##       18       18 
## 
## $convergence
## [1] 0
## 
## $message
## [1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"
parest <- optim(par=c(4,5), fn=objetivo,
                control=list(fnscale=-1, trace=2, reltol=10^-10))
##   Nelder-Mead direct search function minimizer
## function value for initial parameters = 4165.944027
##   Scaled convergence tolerance is 4.16594e-07
## Stepsize computed as 0.500000
## BUILD              3 4165.944027 2944.069378
## EXTENSION          5 3668.062799 2157.344621
## EXTENSION          7 2944.069378 1183.656021
## EXTENSION          9 2157.344621 696.500626
## EXTENSION         11 1183.656021 570.524366
## LO-REDUCTION      13 696.500626 565.290275
## LO-REDUCTION      15 615.440402 565.290275
## REFLECTION        17 570.524366 563.434497
## LO-REDUCTION      19 565.290275 556.534338
## HI-REDUCTION      21 563.434497 555.947636
## LO-REDUCTION      23 556.534338 555.464702
## HI-REDUCTION      25 555.947636 554.940713
## EXTENSION         27 555.464702 554.264972
## HI-REDUCTION      29 554.940713 554.264972
## EXTENSION         31 554.600565 553.391760
## LO-REDUCTION      33 554.264972 553.391760
## REFLECTION        35 553.595021 553.122361
## LO-REDUCTION      37 553.391760 553.122361
## LO-REDUCTION      39 553.275811 553.122361
## LO-REDUCTION      41 553.215429 553.122361
## REFLECTION        43 553.122385 553.088911
## HI-REDUCTION      45 553.122361 553.088911
## LO-REDUCTION      47 553.093317 553.088911
## HI-REDUCTION      49 553.090196 553.084454
## LO-REDUCTION      51 553.088911 553.084454
## HI-REDUCTION      53 553.085031 553.084454
## HI-REDUCTION      55 553.084987 553.084048
## HI-REDUCTION      57 553.084454 553.083936
## HI-REDUCTION      59 553.084048 553.083729
## LO-REDUCTION      61 553.083936 553.083729
## HI-REDUCTION      63 553.083745 553.083729
## REFLECTION        65 553.083745 553.083705
## HI-REDUCTION      67 553.083729 553.083704
## HI-REDUCTION      69 553.083705 553.083701
## HI-REDUCTION      71 553.083704 553.083699
## HI-REDUCTION      73 553.083701 553.083699
## HI-REDUCTION      75 553.083699 553.083698
## LO-REDUCTION      77 553.083699 553.083698
## Exiting from Nelder Mead minimizer
##     79 function evaluations used
parest
## $par
## [1] 7.291885 9.622393
## 
## $value
## [1] -553.0837
## 
## $counts
## function gradient 
##       79       NA 
## 
## $convergence
## [1] 0
## 
## $message
## NULL

Los parámetros que hemos estimado son los logaritmos de los que realmente deseamos. Hacemos la transformación inversa y obtenemos las predicciones

p <- exp(parest$par)
p
## [1]  1468.336 15099.145
result   <- kalman3(y=Nile,0,10^10,Q=p[1],H=p[2])
result
## $x
##   [1]    0.0000 1119.9983 1140.9266 1072.8020 1117.3066 1129.9686 1138.4540
##   [8] 1048.8709 1098.0300 1171.2859 1162.8933 1117.9534 1069.0401 1079.9840
##  [15] 1057.0182 1047.1327 1023.8666 1065.5549  994.3856  984.6709 1026.1423
##  [22] 1045.8615 1089.6847 1105.7882 1144.2911 1175.1840 1187.1493 1145.1923
##  [29] 1133.1265 1037.2440  984.5822  955.0580  885.3586  899.9472  882.0731
##  [36]  833.7287  855.6942  811.9898  867.5260  916.2444  930.3295  903.8097
##  [43]  856.3366  749.4515  769.3551  751.3721  849.7914  916.5941  894.0085
##  [50]  859.2978  849.0727  827.4273  832.1190  840.6308  846.3361  806.7322
##  [57]  816.9492  797.4726  797.0795  861.9364  834.4536  820.1821  832.1480
##  [64]  835.5793  864.5263  896.4244  896.5781  876.6666  912.2650  874.5490
##  [71]  821.5387  775.4729  794.3028  799.0278  783.8020  788.3937  855.5696
##  [78]  856.7525  861.3573  857.7911  866.3905  833.7137  811.0961  818.2791
##  [85]  880.1459  890.2525  915.8159  884.0935  894.4811  915.9787  889.0186
##  [92]  923.9891  919.1862  914.3307  982.5914  963.7443  905.6092  909.1844
##  [99]  858.1404  819.6566  798.3892
## 
## $P
##   [1] 1.000000e+10 1.656746e+04 9.367966e+03 7.249495e+03 6.366226e+03
##   [6] 5.946458e+03 5.734616e+03 5.624464e+03 5.566300e+03 5.535336e+03
##  [11] 5.518781e+03 5.509910e+03 5.505150e+03 5.502594e+03 5.501222e+03
##  [16] 5.500485e+03 5.500088e+03 5.499876e+03 5.499761e+03 5.499700e+03
##  [21] 5.499667e+03 5.499649e+03 5.499640e+03 5.499634e+03 5.499632e+03
##  [26] 5.499630e+03 5.499629e+03 5.499629e+03 5.499629e+03 5.499629e+03
##  [31] 5.499629e+03 5.499628e+03 5.499628e+03 5.499628e+03 5.499628e+03
##  [36] 5.499628e+03 5.499628e+03 5.499628e+03 5.499628e+03 5.499628e+03
##  [41] 5.499628e+03 5.499628e+03 5.499628e+03 5.499628e+03 5.499628e+03
##  [46] 5.499628e+03 5.499628e+03 5.499628e+03 5.499628e+03 5.499628e+03
##  [51] 5.499628e+03 5.499628e+03 5.499628e+03 5.499628e+03 5.499628e+03
##  [56] 5.499628e+03 5.499628e+03 5.499628e+03 5.499628e+03 5.499628e+03
##  [61] 5.499628e+03 5.499628e+03 5.499628e+03 5.499628e+03 5.499628e+03
##  [66] 5.499628e+03 5.499628e+03 5.499628e+03 5.499628e+03 5.499628e+03
##  [71] 5.499628e+03 5.499628e+03 5.499628e+03 5.499628e+03 5.499628e+03
##  [76] 5.499628e+03 5.499628e+03 5.499628e+03 5.499628e+03 5.499628e+03
##  [81] 5.499628e+03 5.499628e+03 5.499628e+03 5.499628e+03 5.499628e+03
##  [86] 5.499628e+03 5.499628e+03 5.499628e+03 5.499628e+03 5.499628e+03
##  [91] 5.499628e+03 5.499628e+03 5.499628e+03 5.499628e+03 5.499628e+03
##  [96] 5.499628e+03 5.499628e+03 5.499628e+03 5.499628e+03 5.499628e+03
## [101] 5.499628e+03
## 
## $vero
## [1] -553.0837
Nile.ll <- ts(result$x,start=1871)
plot(Nile)
lines(Nile.ll,col="red")

Una observación interesante es que obtenemos exactamente los mismos valores filtrados si alteramos las varianzas pero dejando su cociente constante:

k <- 3.6
p <- k * p
result   <- kalman3(y=Nile, 0, k*10^10, Q=p[1], H=p[2])
result
## $x
##   [1]    0.0000 1119.9983 1140.9266 1072.8020 1117.3066 1129.9686 1138.4540
##   [8] 1048.8709 1098.0300 1171.2859 1162.8933 1117.9534 1069.0401 1079.9840
##  [15] 1057.0182 1047.1327 1023.8666 1065.5549  994.3856  984.6709 1026.1423
##  [22] 1045.8615 1089.6847 1105.7882 1144.2911 1175.1840 1187.1493 1145.1923
##  [29] 1133.1265 1037.2440  984.5822  955.0580  885.3586  899.9472  882.0731
##  [36]  833.7287  855.6942  811.9898  867.5260  916.2444  930.3295  903.8097
##  [43]  856.3366  749.4515  769.3551  751.3721  849.7914  916.5941  894.0085
##  [50]  859.2978  849.0727  827.4273  832.1190  840.6308  846.3361  806.7322
##  [57]  816.9492  797.4726  797.0795  861.9364  834.4536  820.1821  832.1480
##  [64]  835.5793  864.5263  896.4244  896.5781  876.6666  912.2650  874.5490
##  [71]  821.5387  775.4729  794.3028  799.0278  783.8020  788.3937  855.5696
##  [78]  856.7525  861.3573  857.7911  866.3905  833.7137  811.0961  818.2791
##  [85]  880.1459  890.2525  915.8159  884.0935  894.4811  915.9787  889.0186
##  [92]  923.9891  919.1862  914.3307  982.5914  963.7443  905.6092  909.1844
##  [99]  858.1404  819.6566  798.3892
## 
## $P
##   [1] 3.600000e+10 5.964285e+04 3.372468e+04 2.609818e+04 2.291841e+04
##   [6] 2.140725e+04 2.064462e+04 2.024807e+04 2.003868e+04 1.992721e+04
##  [11] 1.986761e+04 1.983567e+04 1.981854e+04 1.980934e+04 1.980440e+04
##  [16] 1.980174e+04 1.980032e+04 1.979955e+04 1.979914e+04 1.979892e+04
##  [21] 1.979880e+04 1.979874e+04 1.979870e+04 1.979868e+04 1.979867e+04
##  [26] 1.979867e+04 1.979867e+04 1.979866e+04 1.979866e+04 1.979866e+04
##  [31] 1.979866e+04 1.979866e+04 1.979866e+04 1.979866e+04 1.979866e+04
##  [36] 1.979866e+04 1.979866e+04 1.979866e+04 1.979866e+04 1.979866e+04
##  [41] 1.979866e+04 1.979866e+04 1.979866e+04 1.979866e+04 1.979866e+04
##  [46] 1.979866e+04 1.979866e+04 1.979866e+04 1.979866e+04 1.979866e+04
##  [51] 1.979866e+04 1.979866e+04 1.979866e+04 1.979866e+04 1.979866e+04
##  [56] 1.979866e+04 1.979866e+04 1.979866e+04 1.979866e+04 1.979866e+04
##  [61] 1.979866e+04 1.979866e+04 1.979866e+04 1.979866e+04 1.979866e+04
##  [66] 1.979866e+04 1.979866e+04 1.979866e+04 1.979866e+04 1.979866e+04
##  [71] 1.979866e+04 1.979866e+04 1.979866e+04 1.979866e+04 1.979866e+04
##  [76] 1.979866e+04 1.979866e+04 1.979866e+04 1.979866e+04 1.979866e+04
##  [81] 1.979866e+04 1.979866e+04 1.979866e+04 1.979866e+04 1.979866e+04
##  [86] 1.979866e+04 1.979866e+04 1.979866e+04 1.979866e+04 1.979866e+04
##  [91] 1.979866e+04 1.979866e+04 1.979866e+04 1.979866e+04 1.979866e+04
##  [96] 1.979866e+04 1.979866e+04 1.979866e+04 1.979866e+04 1.979866e+04
## [101] 1.979866e+04
## 
## $vero
## [1] -581.3785

Los valores filtrados son idénticos; la varianza del estado, lógicamente, se ha multiplicado por \(k\).

Filtrado y suavizado de Kalman mediante funciones en R

Función StructTS (sólo modelos univariantes estructurales)

Naturalmente, hay funciones de librería que hacen todo con mínima intervención por nuestra parte. Para modelos estructurales básicos (nivel local, tendencia local lineal, y estacionales) de series univariantes, la más cómoda y rápida opción es StructTS (cargada por defecto en R).

ajuste.ll <- StructTS(x=Nile,type="level", init=c(1469,15100))
plot(Nile)
lines(fitted(ajuste.ll),col="red")

Por defecto, nos devuelve los valores filtrados. Podemos obtener los valores estimados de los parámetros directamente del objeto retornado por StructTS;

ajuste.ll$coef
##   level epsilon 
##    1469   15100

como vemos coindiden con los estimados anteriormente.

La función auxiliar tsSmooth proporciona los valores suavizados, y tsdiag proporciona algunos diagnósticos.

plot(Nile)
lines(fitted(ajuste.ll),col="red")
lines(tsSmooth(ajuste.ll), col="green")

tsdiag(ajuste.ll)

Para ajustar un modelo de tendencia local sólo tenemos que alterar el argumento type= de StructTS:

ajuste.llt <- StructTS(x=Nile,type="trend", init=c(1469,1,15100))
Nile.llt <- ts(fitted(ajuste.llt)[,1],start=1871)
plot(Nile)
lines(Nile.llt,col="red")

No parece haber razón para ajustar una tendencia lineal local (soluciones “frontera” son frecuentes en este tipo de modelos):

ajuste.llt$coef
##     level     slope   epsilon 
##  1416.832     0.000 15080.406

No hay problema con valores perdidos:

Nile.bis <- Nile
Nile.bis[20:30] <- NA
ajuste.llt <- StructTS(x=Nile.bis,type="trend", init=c(1000, 0.5, 10000))
fitted(ajuste.llt)
## Time Series:
## Start = 1871 
## End = 1970 
## Frequency = 1 
##          level        slope
## 1871 1120.0000   0.00000000
## 1872 1147.6924   9.23125671
## 1873 1049.9268 -17.52429249
## 1874 1114.6348  -1.06688445
## 1875 1132.1468   2.03414159
## 1876 1143.2886   3.33888918
## 1877 1041.3218  -9.88874221
## 1878 1088.1760  -3.53398622
## 1879 1159.1781   4.01029478
## 1880 1157.6074   3.49386287
## 1881 1124.0087   0.32451551
## 1882 1084.8651  -2.81731725
## 1883 1087.5166  -2.40842375
## 1884 1068.2958  -3.59797885
## 1885 1056.8832  -4.12492203
## 1886 1037.3290  -5.12303869
## 1887 1055.6880  -3.65659510
## 1888 1013.4872  -5.99403144
## 1889 1000.2392  -6.42357016
## 1890  993.8156  -6.42357016
## 1891  987.3921  -6.42357016
## 1892  980.9685  -6.42357016
## 1893  974.5449  -6.42357016
## 1894  968.1214  -6.42357016
## 1895  961.6978  -6.42357016
## 1896  955.2742  -6.42357016
## 1897  948.8506  -6.42357016
## 1898  942.4271  -6.42357016
## 1899  936.0035  -6.42357016
## 1900  929.5799  -6.42357016
## 1901  906.4885  -7.17757627
## 1902  843.7551  -9.65772184
## 1903  858.2745  -8.58230815
## 1904  846.3579  -8.73133698
## 1905  813.1211  -9.83986071
## 1906  821.7948  -8.98743140
## 1907  794.3536  -9.85597134
## 1908  818.3801  -8.22118735
## 1909  843.0040  -6.59404604
## 1910  853.8481  -5.70662585
## 1911  845.9599  -5.82062323
## 1912  825.9901  -6.57889720
## 1913  775.2885  -8.99996898
## 1914  773.1834  -8.61338681
## 1915  757.1852  -9.03552404
## 1916  791.6480  -6.50639077
## 1917  821.7437  -4.34669922
## 1918  819.0877  -4.24569128
## 1919  808.9701  -4.60017959
## 1920  806.2884  -4.48337323
## 1921  797.9060  -4.72238131
## 1922  799.1629  -4.35392384
## 1923  802.8007  -3.85952849
## 1924  806.2331  -3.40715779
## 1925  790.6884  -4.16151124
## 1926  793.3062  -3.73973940
## 1927  784.2771  -4.06894745
## 1928  782.0433  -3.95474239
## 1929  808.5533  -2.05982974
## 1930  800.9650  -2.40342357
## 1931  796.5162  -2.53041673
## 1932  802.2599  -2.01728534
## 1933  805.4586  -1.69420856
## 1934  820.1079  -0.68314745
## 1935  838.6030   0.50181756
## 1936  845.8500   0.91808372
## 1937  843.8833   0.74025897
## 1938  863.8771   1.92587875
## 1939  854.7709   1.24717890
## 1940  835.0803  -0.03981897
## 1941  813.4138  -1.36814626
## 1942  815.9906  -1.12600742
## 1943  814.5319  -1.14641258
## 1944  805.1009  -1.65438804
## 1945  803.1627  -1.67178100
## 1946  829.1423   0.02243638
## 1947  832.7378   0.24130850
## 1948  837.7304   0.53228928
## 1949  839.3901   0.60132152
## 1950  845.7791   0.95568393
## 1951  834.8490   0.22800347
## 1952  825.1216  -0.38147715
## 1953  826.2733  -0.28761609
## 1954  851.8809   1.29771769
## 1955  860.6701   1.75636710
## 1956  876.7053   2.63059141
## 1957  869.8236   2.04817113
## 1958  877.7778   2.40980590
## 1959  891.1386   3.08038016
## 1960  885.0700   2.52011641
## 1961  902.8809   3.45653887
## 1962  906.2985   3.45415218
## 1963  908.7420   3.39225325
## 1964  941.9082   5.21588683
## 1965  943.0687   4.96748634
## 1966  924.7098   3.53865354
## 1967  927.1807   3.47324613
## 1968  906.1024   1.96928051
## 1969  885.6657   0.59671487
## 1970  869.3762  -0.43773630
plot(Nile.bis)
lines(fitted(ajuste.llt)[,1],col="red")

Veamos un ejemplo con estacionalidad anual y datos trimestrales (¿como la averigua StructTS? Examinando los atributos del objeto ts. Es imperativo que lo que se pasa como argumento sea un objeto ts).

UKgas
##        Qtr1   Qtr2   Qtr3   Qtr4
## 1960  160.1  129.7   84.8  120.1
## 1961  160.1  124.9   84.8  116.9
## 1962  169.7  140.9   89.7  123.3
## 1963  187.3  144.1   92.9  120.1
## 1964  176.1  147.3   89.7  123.3
## 1965  185.7  155.3   99.3  131.3
## 1966  200.1  161.7  102.5  136.1
## 1967  204.9  176.1  112.1  140.9
## 1968  227.3  195.3  115.3  142.5
## 1969  244.9  214.5  118.5  153.7
## 1970  244.9  216.1  188.9  142.5
## 1971  301.0  196.9  136.1  267.3
## 1972  317.0  230.5  152.1  336.2
## 1973  371.4  240.1  158.5  355.4
## 1974  449.9  286.6  179.3  403.4
## 1975  491.5  321.8  177.7  409.8
## 1976  593.9  329.8  176.1  483.5
## 1977  584.3  395.4  187.3  485.1
## 1978  669.2  421.0  216.1  509.1
## 1979  827.7  467.5  209.7  542.7
## 1980  840.5  414.6  217.7  670.8
## 1981  848.5  437.0  209.7  701.2
## 1982  925.3  443.4  214.5  683.6
## 1983  917.3  515.5  224.1  694.8
## 1984  989.4  477.1  233.7  730.0
## 1985 1087.0  534.7  281.8  787.6
## 1986 1163.9  613.1  347.4  782.8
start(UKgas)
## [1] 1960    1
end(UKgas)
## [1] 1986    4
frequency(UKgas)
## [1] 4
plot(UKgas)

plot(log(UKgas))

UKgas.bsm <- StructTS(log10(UKgas), type = "BSM")
plot(cbind(fitted(UKgas.bsm), resids=resid(UKgas.bsm)), main = "Consumo de gas en el Reino Unido")

Uso de las funciones en dlm

El uso de dlm permite reproducir el trabajo anterior y ampliarlo en diversas direcciones. Así ajustaríamos el modelo de nivel local. Podemos obtener predicciones una etapa hacia adelante de las observaciones $f, del vector de estado ($a) o valores filtrados del vector de estado ($m). Podemos también suavizar.

library(dlm)
mod.Nilo.ll  <- dlm(FF=1, GG=1, V=15100, W=1469, m0=0, C0=10^7)
res.Nilo.ll <- dlmFilter(Nile,mod.Nilo.ll)
plot(res.Nilo.ll$y)
lines(res.Nilo.ll$f,col="red")
smo.Nilo.ll <- dlmSmooth(Nile,mod.Nilo.ll)
lines(smo.Nilo.ll$s,col="green")

Podemos emplear funciones auxiliares para especificar el modelo de interés:

dlmModPoly(1, dV=15099, dW=1469, m0=0, C0=10^7)
## $FF
##      [,1]
## [1,]    1
## 
## $V
##       [,1]
## [1,] 15099
## 
## $GG
##      [,1]
## [1,]    1
## 
## $W
##      [,1]
## [1,] 1469
## 
## $m0
## [1] 0
## 
## $C0
##       [,1]
## [1,] 1e+07
dlmModPoly(2, dV=15047, dW=c(1427,0.001), m0=c(0,0),
           C0=10^7*diag(2))
## $FF
##      [,1] [,2]
## [1,]    1    0
## 
## $V
##       [,1]
## [1,] 15047
## 
## $GG
##      [,1] [,2]
## [1,]    1    1
## [2,]    0    1
## 
## $W
##      [,1]  [,2]
## [1,] 1427 0.000
## [2,]    0 0.001
## 
## $m0
## [1] 0 0
## 
## $C0
##       [,1]  [,2]
## [1,] 1e+07 0e+00
## [2,] 0e+00 1e+07
mod.Nilo.llt <- dlmModPoly(2, dV=15047, dW=c(1427,0), m0=c(0,0),
                           C0=10^7*diag(2))
res.Nilo.llt <- dlmFilter(Nile,mod.Nilo.llt)
str(res.Nilo.llt, max.depth=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 10
##   ..$ m0 : num [1:2] 0 0
##   ..$ C0 : num [1:2, 1:2] 1e+07 0e+00 0e+00 1e+07
##   ..$ FF : num [1, 1:2] 1 0
##   ..$ V  : num [1, 1] 15047
##   ..$ GG : num [1:2, 1:2] 1 0 1 1
##   ..$ W  : num [1:2, 1:2] 1427 0 0 0
##   ..$ JFF: NULL
##   ..$ JV : NULL
##   ..$ JGG: NULL
##   ..$ JW : NULL
##   ..- attr(*, "class")= chr "dlm"
##  $ m  : Time-Series [1:101, 1:2] from 1870 to 1970: 0 1119 1162 1003 1128 ...
##  $ U.C:List of 101
##   ..$ : num [1:2, 1:2] 1 0 0 1
##   ..$ : num [1:2, 1:2] -1 0.00151 0.00151 1
##   ..$ : num [1:2, 1:2] 0.859 -0.511 0.511 0.859
##   ..$ : num [1:2, 1:2] -0.6 0.8 0.8 0.6
##   ..$ : num [1:2, 1:2] -0.433 0.901 0.901 0.433
##   ..$ : num [1:2, 1:2] -0.334 0.942 0.942 0.334
##   ..$ : num [1:2, 1:2] -0.27 0.963 0.963 0.27
##   ..$ : num [1:2, 1:2] -0.225 0.974 0.974 0.225
##   ..$ : num [1:2, 1:2] -0.192 0.981 0.981 0.192
##   ..$ : num [1:2, 1:2] -0.167 0.986 0.986 0.167
##   ..$ : num [1:2, 1:2] -0.147 0.989 0.989 0.147
##   ..$ : num [1:2, 1:2] -0.13 0.991 0.991 0.13
##   ..$ : num [1:2, 1:2] -0.117 0.993 0.993 0.117
##   ..$ : num [1:2, 1:2] -0.106 0.994 0.994 0.106
##   ..$ : num [1:2, 1:2] -0.0968 0.9953 0.9953 0.0968
##   ..$ : num [1:2, 1:2] -0.0888 0.996 0.996 0.0888
##   ..$ : num [1:2, 1:2] -0.082 0.997 0.997 0.082
##   ..$ : num [1:2, 1:2] -0.0761 0.9971 0.9971 0.0761
##   ..$ : num [1:2, 1:2] -0.0709 0.9975 0.9975 0.0709
##   ..$ : num [1:2, 1:2] -0.0664 0.9978 0.9978 0.0664
##   ..$ : num [1:2, 1:2] -0.0624 0.9981 0.9981 0.0624
##   ..$ : num [1:2, 1:2] -0.0588 0.9983 0.9983 0.0588
##   ..$ : num [1:2, 1:2] -0.0556 0.9985 0.9985 0.0556
##   ..$ : num [1:2, 1:2] -0.0527 0.9986 0.9986 0.0527
##   ..$ : num [1:2, 1:2] -0.0501 0.9987 0.9987 0.0501
##   ..$ : num [1:2, 1:2] -0.0477 0.9989 0.9989 0.0477
##   ..$ : num [1:2, 1:2] -0.0456 0.999 0.999 0.0456
##   ..$ : num [1:2, 1:2] -0.0436 0.999 0.999 0.0436
##   ..$ : num [1:2, 1:2] -0.0418 0.9991 0.9991 0.0418
##   ..$ : num [1:2, 1:2] -0.0401 0.9992 0.9992 0.0401
##   ..$ : num [1:2, 1:2] -0.0386 0.9993 0.9993 0.0386
##   ..$ : num [1:2, 1:2] -0.0371 0.9993 0.9993 0.0371
##   ..$ : num [1:2, 1:2] -0.0358 0.9994 0.9994 0.0358
##   ..$ : num [1:2, 1:2] -0.0346 0.9994 0.9994 0.0346
##   ..$ : num [1:2, 1:2] -0.0334 0.9994 0.9994 0.0334
##   ..$ : num [1:2, 1:2] -0.0324 0.9995 0.9995 0.0324
##   ..$ : num [1:2, 1:2] -0.0313 0.9995 0.9995 0.0313
##   ..$ : num [1:2, 1:2] -0.0304 0.9995 0.9995 0.0304
##   ..$ : num [1:2, 1:2] -0.0295 0.9996 0.9996 0.0295
##   ..$ : num [1:2, 1:2] -0.0287 0.9996 0.9996 0.0287
##   ..$ : num [1:2, 1:2] -0.0279 0.9996 0.9996 0.0279
##   ..$ : num [1:2, 1:2] -0.0271 0.9996 0.9996 0.0271
##   ..$ : num [1:2, 1:2] -0.0264 0.9997 0.9997 0.0264
##   ..$ : num [1:2, 1:2] -0.0257 0.9997 0.9997 0.0257
##   ..$ : num [1:2, 1:2] -0.0251 0.9997 0.9997 0.0251
##   ..$ : num [1:2, 1:2] -0.0245 0.9997 0.9997 0.0245
##   ..$ : num [1:2, 1:2] -0.0239 0.9997 0.9997 0.0239
##   ..$ : num [1:2, 1:2] -0.0233 0.9997 0.9997 0.0233
##   ..$ : num [1:2, 1:2] -0.0228 0.9997 0.9997 0.0228
##   ..$ : num [1:2, 1:2] -0.0223 0.9998 0.9998 0.0223
##   ..$ : num [1:2, 1:2] -0.0218 0.9998 0.9998 0.0218
##   ..$ : num [1:2, 1:2] -0.0213 0.9998 0.9998 0.0213
##   ..$ : num [1:2, 1:2] -0.0209 0.9998 0.9998 0.0209
##   ..$ : num [1:2, 1:2] -0.0205 0.9998 0.9998 0.0205
##   ..$ : num [1:2, 1:2] -0.02 1 1 0.02
##   ..$ : num [1:2, 1:2] -0.0197 0.9998 0.9998 0.0197
##   ..$ : num [1:2, 1:2] -0.0193 0.9998 0.9998 0.0193
##   ..$ : num [1:2, 1:2] -0.0189 0.9998 0.9998 0.0189
##   ..$ : num [1:2, 1:2] -0.0186 0.9998 0.9998 0.0186
##   ..$ : num [1:2, 1:2] -0.0182 0.9998 0.9998 0.0182
##   ..$ : num [1:2, 1:2] -0.0179 0.9998 0.9998 0.0179
##   ..$ : num [1:2, 1:2] -0.0176 0.9998 0.9998 0.0176
##   ..$ : num [1:2, 1:2] -0.0173 0.9999 0.9999 0.0173
##   ..$ : num [1:2, 1:2] -0.017 1 1 0.017
##   ..$ : num [1:2, 1:2] -0.0167 0.9999 0.9999 0.0167
##   ..$ : num [1:2, 1:2] -0.0164 0.9999 0.9999 0.0164
##   ..$ : num [1:2, 1:2] -0.0162 0.9999 0.9999 0.0162
##   ..$ : num [1:2, 1:2] -0.0159 0.9999 0.9999 0.0159
##   ..$ : num [1:2, 1:2] -0.0157 0.9999 0.9999 0.0157
##   ..$ : num [1:2, 1:2] -0.0154 0.9999 0.9999 0.0154
##   ..$ : num [1:2, 1:2] -0.0152 0.9999 0.9999 0.0152
##   ..$ : num [1:2, 1:2] -0.015 1 1 0.015
##   ..$ : num [1:2, 1:2] -0.0147 0.9999 0.9999 0.0147
##   ..$ : num [1:2, 1:2] -0.0145 0.9999 0.9999 0.0145
##   ..$ : num [1:2, 1:2] -0.0143 0.9999 0.9999 0.0143
##   ..$ : num [1:2, 1:2] -0.0141 0.9999 0.9999 0.0141
##   ..$ : num [1:2, 1:2] -0.0139 0.9999 0.9999 0.0139
##   ..$ : num [1:2, 1:2] -0.0137 0.9999 0.9999 0.0137
##   ..$ : num [1:2, 1:2] -0.0135 0.9999 0.9999 0.0135
##   ..$ : num [1:2, 1:2] -0.0134 0.9999 0.9999 0.0134
##   ..$ : num [1:2, 1:2] -0.0132 0.9999 0.9999 0.0132
##   ..$ : num [1:2, 1:2] -0.013 1 1 0.013
##   ..$ : num [1:2, 1:2] -0.0128 0.9999 0.9999 0.0128
##   ..$ : num [1:2, 1:2] -0.0127 0.9999 0.9999 0.0127
##   ..$ : num [1:2, 1:2] -0.0125 0.9999 0.9999 0.0125
##   ..$ : num [1:2, 1:2] -0.0124 0.9999 0.9999 0.0124
##   ..$ : num [1:2, 1:2] -0.0122 0.9999 0.9999 0.0122
##   ..$ : num [1:2, 1:2] -0.0121 0.9999 0.9999 0.0121
##   ..$ : num [1:2, 1:2] -0.0119 0.9999 0.9999 0.0119
##   ..$ : num [1:2, 1:2] -0.0118 0.9999 0.9999 0.0118
##   ..$ : num [1:2, 1:2] -0.0116 0.9999 0.9999 0.0116
##   ..$ : num [1:2, 1:2] -0.0115 0.9999 0.9999 0.0115
##   ..$ : num [1:2, 1:2] -0.0114 0.9999 0.9999 0.0114
##   ..$ : num [1:2, 1:2] -0.0113 0.9999 0.9999 0.0113
##   ..$ : num [1:2, 1:2] -0.0111 0.9999 0.9999 0.0111
##   ..$ : num [1:2, 1:2] -0.011 1 1 0.011
##   ..$ : num [1:2, 1:2] -0.0109 0.9999 0.9999 0.0109
##   ..$ : num [1:2, 1:2] -0.0108 0.9999 0.9999 0.0108
##   ..$ : num [1:2, 1:2] -0.0107 0.9999 0.9999 0.0107
##   .. [list output truncated]
##  $ D.C: num [1:101, 1:2] 3162.3 122.6 78.2 50.9 36.2 ...
##  $ a  : Time-Series [1:100, 1:2] from 1871 to 1970: 0 1679 1206 926 1137 ...
##  $ U.R:List of 100
##   ..$ : num [1:2, 1:2] -0.851 -0.526 0.526 -0.851
##   ..$ : num [1:2, 1:2] 0.708 0.706 0.706 -0.708
##   ..$ : num [1:2, 1:2] -0.851 -0.525 0.525 -0.851
##   ..$ : num [1:2, 1:2] -0.916 -0.401 -0.401 0.916
##   ..$ : num [1:2, 1:2] -0.948 -0.318 -0.318 0.948
##   ..$ : num [1:2, 1:2] -0.965 -0.261 -0.261 0.965
##   ..$ : num [1:2, 1:2] -0.976 -0.219 -0.219 0.976
##   ..$ : num [1:2, 1:2] -0.982 -0.188 -0.188 0.982
##   ..$ : num [1:2, 1:2] -0.987 -0.164 -0.164 0.987
##   ..$ : num [1:2, 1:2] -0.99 -0.144 -0.144 0.99
##   ..$ : num [1:2, 1:2] -0.992 -0.129 -0.129 0.992
##   ..$ : num [1:2, 1:2] -0.993 -0.116 -0.116 0.993
##   ..$ : num [1:2, 1:2] -0.994 -0.105 -0.105 0.994
##   ..$ : num [1:2, 1:2] -0.9954 -0.0959 -0.0959 0.9954
##   ..$ : num [1:2, 1:2] -0.9961 -0.0881 -0.0881 0.9961
##   ..$ : num [1:2, 1:2] -0.9967 -0.0814 -0.0814 0.9967
##   ..$ : num [1:2, 1:2] -0.9971 -0.0756 -0.0756 0.9971
##   ..$ : num [1:2, 1:2] -0.9975 -0.0705 -0.0705 0.9975
##   ..$ : num [1:2, 1:2] -0.998 -0.066 -0.066 0.998
##   ..$ : num [1:2, 1:2] -0.998 -0.062 -0.062 0.998
##   ..$ : num [1:2, 1:2] -0.9983 -0.0585 -0.0585 0.9983
##   ..$ : num [1:2, 1:2] -0.9985 -0.0553 -0.0553 0.9985
##   ..$ : num [1:2, 1:2] -0.9986 -0.0524 -0.0524 0.9986
##   ..$ : num [1:2, 1:2] -0.9988 -0.0499 -0.0499 0.9988
##   ..$ : num [1:2, 1:2] -0.9989 -0.0475 -0.0475 0.9989
##   ..$ : num [1:2, 1:2] -0.999 -0.0454 -0.0454 0.999
##   ..$ : num [1:2, 1:2] -0.9991 -0.0434 -0.0434 0.9991
##   ..$ : num [1:2, 1:2] -0.9991 -0.0416 -0.0416 0.9991
##   ..$ : num [1:2, 1:2] -0.999 -0.04 -0.04 0.999
##   ..$ : num [1:2, 1:2] -0.9993 -0.0384 -0.0384 0.9993
##   ..$ : num [1:2, 1:2] -0.999 -0.037 -0.037 0.999
##   ..$ : num [1:2, 1:2] -0.9994 -0.0357 -0.0357 0.9994
##   ..$ : num [1:2, 1:2] -0.9994 -0.0345 -0.0345 0.9994
##   ..$ : num [1:2, 1:2] -0.9994 -0.0333 -0.0333 0.9994
##   ..$ : num [1:2, 1:2] -0.9995 -0.0323 -0.0323 0.9995
##   ..$ : num [1:2, 1:2] -0.9995 -0.0312 -0.0312 0.9995
##   ..$ : num [1:2, 1:2] -0.9995 -0.0303 -0.0303 0.9995
##   ..$ : num [1:2, 1:2] -0.9996 -0.0294 -0.0294 0.9996
##   ..$ : num [1:2, 1:2] -0.9996 -0.0286 -0.0286 0.9996
##   ..$ : num [1:2, 1:2] -0.9996 -0.0278 -0.0278 0.9996
##   ..$ : num [1:2, 1:2] -1 -0.027 -0.027 1
##   ..$ : num [1:2, 1:2] -0.9997 -0.0263 -0.0263 0.9997
##   ..$ : num [1:2, 1:2] -0.9997 -0.0256 -0.0256 0.9997
##   ..$ : num [1:2, 1:2] -1 -0.025 -0.025 1
##   ..$ : num [1:2, 1:2] -0.9997 -0.0244 -0.0244 0.9997
##   ..$ : num [1:2, 1:2] -0.9997 -0.0238 -0.0238 0.9997
##   ..$ : num [1:2, 1:2] -0.9997 -0.0233 -0.0233 0.9997
##   ..$ : num [1:2, 1:2] -0.9997 -0.0227 -0.0227 0.9997
##   ..$ : num [1:2, 1:2] -0.9998 -0.0222 -0.0222 0.9998
##   ..$ : num [1:2, 1:2] -0.9998 -0.0217 -0.0217 0.9998
##   ..$ : num [1:2, 1:2] -0.9998 -0.0213 -0.0213 0.9998
##   ..$ : num [1:2, 1:2] -0.9998 -0.0208 -0.0208 0.9998
##   ..$ : num [1:2, 1:2] -0.9998 -0.0204 -0.0204 0.9998
##   ..$ : num [1:2, 1:2] -1 -0.02 -0.02 1
##   ..$ : num [1:2, 1:2] -0.9998 -0.0196 -0.0196 0.9998
##   ..$ : num [1:2, 1:2] -0.9998 -0.0192 -0.0192 0.9998
##   ..$ : num [1:2, 1:2] -0.9998 -0.0189 -0.0189 0.9998
##   ..$ : num [1:2, 1:2] -0.9998 -0.0185 -0.0185 0.9998
##   ..$ : num [1:2, 1:2] -0.9998 -0.0182 -0.0182 0.9998
##   ..$ : num [1:2, 1:2] -0.9998 -0.0179 -0.0179 0.9998
##   ..$ : num [1:2, 1:2] -0.9998 -0.0176 -0.0176 0.9998
##   ..$ : num [1:2, 1:2] -0.9999 -0.0172 -0.0172 0.9999
##   ..$ : num [1:2, 1:2] -1 -0.017 -0.017 1
##   ..$ : num [1:2, 1:2] -0.9999 -0.0167 -0.0167 0.9999
##   ..$ : num [1:2, 1:2] -0.9999 -0.0164 -0.0164 0.9999
##   ..$ : num [1:2, 1:2] -0.9999 -0.0161 -0.0161 0.9999
##   ..$ : num [1:2, 1:2] -0.9999 -0.0159 -0.0159 0.9999
##   ..$ : num [1:2, 1:2] -0.9999 -0.0156 -0.0156 0.9999
##   ..$ : num [1:2, 1:2] -0.9999 -0.0154 -0.0154 0.9999
##   ..$ : num [1:2, 1:2] -0.9999 -0.0152 -0.0152 0.9999
##   ..$ : num [1:2, 1:2] -0.9999 -0.0149 -0.0149 0.9999
##   ..$ : num [1:2, 1:2] -0.9999 -0.0147 -0.0147 0.9999
##   ..$ : num [1:2, 1:2] -0.9999 -0.0145 -0.0145 0.9999
##   ..$ : num [1:2, 1:2] -0.9999 -0.0143 -0.0143 0.9999
##   ..$ : num [1:2, 1:2] -0.9999 -0.0141 -0.0141 0.9999
##   ..$ : num [1:2, 1:2] -0.9999 -0.0139 -0.0139 0.9999
##   ..$ : num [1:2, 1:2] -0.9999 -0.0137 -0.0137 0.9999
##   ..$ : num [1:2, 1:2] -0.9999 -0.0135 -0.0135 0.9999
##   ..$ : num [1:2, 1:2] -0.9999 -0.0133 -0.0133 0.9999
##   ..$ : num [1:2, 1:2] -0.9999 -0.0132 -0.0132 0.9999
##   ..$ : num [1:2, 1:2] -1 -0.013 -0.013 1
##   ..$ : num [1:2, 1:2] -0.9999 -0.0128 -0.0128 0.9999
##   ..$ : num [1:2, 1:2] -0.9999 -0.0127 -0.0127 0.9999
##   ..$ : num [1:2, 1:2] -0.9999 -0.0125 -0.0125 0.9999
##   ..$ : num [1:2, 1:2] -0.9999 -0.0124 -0.0124 0.9999
##   ..$ : num [1:2, 1:2] -0.9999 -0.0122 -0.0122 0.9999
##   ..$ : num [1:2, 1:2] -0.9999 -0.0121 -0.0121 0.9999
##   ..$ : num [1:2, 1:2] -0.9999 -0.0119 -0.0119 0.9999
##   ..$ : num [1:2, 1:2] -0.9999 -0.0118 -0.0118 0.9999
##   ..$ : num [1:2, 1:2] -0.9999 -0.0116 -0.0116 0.9999
##   ..$ : num [1:2, 1:2] -0.9999 -0.0115 -0.0115 0.9999
##   ..$ : num [1:2, 1:2] -0.9999 -0.0114 -0.0114 0.9999
##   ..$ : num [1:2, 1:2] -0.9999 -0.0112 -0.0112 0.9999
##   ..$ : num [1:2, 1:2] -0.9999 -0.0111 -0.0111 0.9999
##   ..$ : num [1:2, 1:2] -1 -0.011 -0.011 1
##   ..$ : num [1:2, 1:2] -0.9999 -0.0109 -0.0109 0.9999
##   ..$ : num [1:2, 1:2] -0.9999 -0.0108 -0.0108 0.9999
##   ..$ : num [1:2, 1:2] -0.9999 -0.0106 -0.0106 0.9999
##   ..$ : num [1:2, 1:2] -0.9999 -0.0105 -0.0105 0.9999
##   .. [list output truncated]
##  $ D.R: num [1:100, 1:2] 5117 3167 326 210 165 ...
##  $ f  : Time-Series [1:100] from 1871 to 1970: 0 1679 1206 926 1137 ...
##  - attr(*, "class")= chr "dlmFiltered"

Una tendencia cuadrática local se especifica de manera igualmente fácil:

dlmModPoly(3, dV=15047, dW=1:3, m0=rep(0,3), C0=10^7*diag(3))
## $FF
##      [,1] [,2] [,3]
## [1,]    1    0    0
## 
## $V
##       [,1]
## [1,] 15047
## 
## $GG
##      [,1] [,2] [,3]
## [1,]    1    1    0
## [2,]    0    1    1
## [3,]    0    0    1
## 
## $W
##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    0    2    0
## [3,]    0    0    3
## 
## $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

Para estimar un modelo, debemos definir una función que lo re-cree para diferentes valores de los parámetros. Por ejemplo, en el caso del modelo de nivel local para el Nilo:

crearMod <- function(x) {
  return(dlmModPoly(1, dV=exp(x[1]), dW=exp(x[2]), m0=0, C0=10^7))
}

La función dlmMLE es un envoltorio que se ocupa de llamar a optim:

dlmMLE
## function (y, parm, build, method = "L-BFGS-B", ..., debug = FALSE) 
## {
##     logLik <- function(parm, ...) {
##         mod <- build(parm, ...)
##         return(dlmLL(y = y, mod = mod, debug = debug))
##     }
##     out <- optim(parm, logLik, method = method, ...)
##     return(out)
## }
## <bytecode: 0x55c70158c938>
## <environment: namespace:dlm>

Veamos como utilizarla. Observemos que como argumentos opcionales podemos pasar argumentos propios de optim.

ajuste.Nilo.ll <- dlmMLE(Nile, parm=c(0,0), build=crearMod)
ajuste.Nilo.ll <- dlmMLE(Nile, parm=c(0,0), build=crearMod,
                         method="Nelder-Mead",
                         control=list(reltol=10^-10))
mod.Nilo.ll <- crearMod(ajuste.Nilo.ll$par)
mod.Nilo.ll
## $FF
##      [,1]
## [1,]    1
## 
## $V
##          [,1]
## [1,] 15116.68
## 
## $GG
##      [,1]
## [1,]    1
## 
## $W
##          [,1]
## [1,] 1460.262
## 
## $m0
## [1] 0
## 
## $C0
##       [,1]
## [1,] 1e+07

Es igualmente fácil con estacionalidad.

Cordero <- dget(file="Cordero.dge")
dlmModSeas(12)
## $FF
##      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
## [1,]    1    0    0    0    0    0    0    0    0     0     0
## 
## $V
##      [,1]
## [1,]    1
## 
## $GG
##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
##  [1,]   -1   -1   -1   -1   -1   -1   -1   -1   -1    -1    -1
##  [2,]    1    0    0    0    0    0    0    0    0     0     0
##  [3,]    0    1    0    0    0    0    0    0    0     0     0
##  [4,]    0    0    1    0    0    0    0    0    0     0     0
##  [5,]    0    0    0    1    0    0    0    0    0     0     0
##  [6,]    0    0    0    0    1    0    0    0    0     0     0
##  [7,]    0    0    0    0    0    1    0    0    0     0     0
##  [8,]    0    0    0    0    0    0    1    0    0     0     0
##  [9,]    0    0    0    0    0    0    0    1    0     0     0
## [10,]    0    0    0    0    0    0    0    0    1     0     0
## [11,]    0    0    0    0    0    0    0    0    0     1     0
## 
## $W
##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
##  [1,]    1    0    0    0    0    0    0    0    0     0     0
##  [2,]    0    0    0    0    0    0    0    0    0     0     0
##  [3,]    0    0    0    0    0    0    0    0    0     0     0
##  [4,]    0    0    0    0    0    0    0    0    0     0     0
##  [5,]    0    0    0    0    0    0    0    0    0     0     0
##  [6,]    0    0    0    0    0    0    0    0    0     0     0
##  [7,]    0    0    0    0    0    0    0    0    0     0     0
##  [8,]    0    0    0    0    0    0    0    0    0     0     0
##  [9,]    0    0    0    0    0    0    0    0    0     0     0
## [10,]    0    0    0    0    0    0    0    0    0     0     0
## [11,]    0    0    0    0    0    0    0    0    0     0     0
## 
## $m0
##  [1] 0 0 0 0 0 0 0 0 0 0 0
## 
## $C0
##        [,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8]  [,9] [,10] [,11]
##  [1,] 1e+07 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00
##  [2,] 0e+00 1e+07 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00
##  [3,] 0e+00 0e+00 1e+07 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00
##  [4,] 0e+00 0e+00 0e+00 1e+07 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00
##  [5,] 0e+00 0e+00 0e+00 0e+00 1e+07 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00
##  [6,] 0e+00 0e+00 0e+00 0e+00 0e+00 1e+07 0e+00 0e+00 0e+00 0e+00 0e+00
##  [7,] 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 1e+07 0e+00 0e+00 0e+00 0e+00
##  [8,] 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 1e+07 0e+00 0e+00 0e+00
##  [9,] 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 1e+07 0e+00 0e+00
## [10,] 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 1e+07 0e+00
## [11,] 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 1e+07
dlmModPoly(1)
## $FF
##      [,1]
## [1,]    1
## 
## $V
##      [,1]
## [1,]    1
## 
## $GG
##      [,1]
## [1,]    1
## 
## $W
##      [,1]
## [1,]    1
## 
## $m0
## [1] 0
## 
## $C0
##       [,1]
## [1,] 1e+07
crearMod <- function(x) {
  mod <- dlmModPoly(1,dV=exp(x[1]),dW=exp(x[2]))  +
         dlmModSeas(12,dV=0,dW=c(exp(x[3]),rep(0,10)))
  return(mod)
}
ajuste.Cordero.lls <- dlmMLE(Cordero,parm=c(0,0,0), build=crearMod)
ajuste.Cordero.lls
## $par
## [1] -1.932482 12.719464 13.521408
## 
## $value
## [1] 703.5872
## 
## $counts
## function gradient 
##       29       29 
## 
## $convergence
## [1] 0
## 
## $message
## [1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"
mod.Cordero.lls <- crearMod(ajuste.Cordero.lls$par)
res.Cordero.lls <- dlmFilter(Cordero,mod.Cordero.lls)
plot(Cordero)
lines(res.Cordero.lls$f,col="red")
smo.Cordero.lls <- dlmSmooth(Cordero,mod.Cordero.lls)
lines(smo.Cordero.lls$s[,1]+smo.Cordero.lls$s[,2],col="green")

Rehagamos el ejemplo de análisis de intervención en una de las sesiones anteriores (efecto del cierre en 1898 de la primera presa de Assuan).

crearMod <- function(x) {
  X <- matrix(ifelse(time(Nile)<1899,0,1),nc=1)
  mod1 <- dlmModPoly(1,dV=exp(x[1]),dW=exp(x[2]))
  mod2 <- dlmModReg(X,dV=0,addInt=FALSE)
  return(mod1+mod2)
}
mod.Nilo.salto <- crearMod(c(0,0))
ajuste.Nilo.salto <- dlmMLE(Nile, parm=c(0,0), build=crearMod)
ajuste.Nilo.salto
## $par
## [1]  9.6987123 -0.5675894
## 
## $value
## [1] 544.237
## 
## $counts
## function gradient 
##       18       18 
## 
## $convergence
## [1] 0
## 
## $message
## [1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"
mod.Nilo.salto <- crearMod(ajuste.Nilo.salto$par)
mod.Nilo.salto
## $FF
##      [,1] [,2]
## [1,]    1    1
## 
## $V
##          [,1]
## [1,] 16296.61
## 
## $GG
##      [,1] [,2]
## [1,]    1    0
## [2,]    0    1
## 
## $W
##           [,1] [,2]
## [1,] 0.5668903    0
## [2,] 0.0000000    0
## 
## $JFF
##      [,1] [,2]
## [1,]    0    1
## 
## $X
##      [,1]
## [1,] 0   
## [2,] 0   
## [3,] ... 
## 
## $m0
## [1] 0 0
## 
## $C0
##       [,1]  [,2]
## [1,] 1e+07 0e+00
## [2,] 0e+00 1e+07
res.Nilo.salto <- dlmFilter(Nile,mod.Nilo.salto)

Veamos los resultados:

plot(dropFirst(Nile))
lines(dropFirst(res.Nilo.salto$f), col="red")

Una forma de valorar la influencia del cambio en 1898 es ver si el elemento del vector de estado que recoge su influencia es significativo. Dicho elemento es el segundo. Lo representaremos junto con un intervalo de confianza.

efecto <- ts(res.Nilo.salto$a[,2], start=c(1871,1))
Var.efecto <- rep(NA,length(efecto))
for (i in 1:length(efecto)) {
Var.efecto[i] <- ( res.Nilo.salto$U.C[[i]] %*% diag(res.Nilo.salto$D.C[i,]^2) 
         %*% t(res.Nilo.salto$U.C[[i]]) )[2,2]
}
plot(efecto)
ancho.ic <- 1.96 * sqrt(Var.efecto)
lines(efecto + ancho.ic, col="blue")
lines(efecto - ancho.ic, col="blue")

La gráfica anterior muestra de forma bastante palmaria que la magnitud del salto estimado es bastante grande en relación a su varianza, lo que sugiere que el cambio operado en 1898 es significativo.

Quizá en este caso sería más adecuado examinar no el valor filtrado del elemento de interés, sino el valor suavizado.

res.Nilo.smooth <- dlmSmooth(res.Nilo.salto)
smooth <- ts(res.Nilo.smooth$s[,2], start=c(1871,1))
temp <- dlmSvd2var(res.Nilo.smooth$U.S, res.Nilo.smooth$D.S)
Var.smooth  <- rep(NA,length(smooth))
for (i in 1:length(smooth)) {
   Var.smooth[i] <- (temp[[i]])[2,2]
}

Siendo el efecto buscado constante, el valor suavizado tanto del mismo como de su varianza lo serán también. Basta mirar el primero de sus valores: en general habría que mirar toda la trayectoria, y por eso se ha mostrado como calcular todos los valores, aunque aquí no se utilicen.

smooth[1]
## [1] -248.1385
sqrt(Var.smooth[1])
## [1] 28.75747