(“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.)
Modifiquemos la función kalman2
en uno de los ejemplos precedente de manera que calcule la verosimilitud de un modelo univariante de nivel local:
<- function(y,x0,P0,Q,H) {
kalman3 <- length(y)
n <- P <- rep(0,n)
x 1] <- x0
x[1] <- P0
P[<- 0
vero for (i in 1:n) {
if (!is.na(y[i])) {
<- y[i] - x[i]
inno <- P[i] + H
F <- P[i] / F
K +1] <- x[i] + K * inno
x[i+1] <- P[i]*(1-K) + Q
P[i<- vero + log(F) + (inno^2) / F
vero else {
} +1] <- x[i]
x[i+1] <- P[i] + Q
P[i
}
}<- -0.5*vero
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).
<- kalman3(y=Nile,0,10^7,Q=150,H=15000)
result 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
.
<- function(p) {
objetivo <- exp(p[1]) # varianza no negativa
Q <- exp(p[2]) # varianza no negativa
H 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
)
<- optim(par=c(4,5),
parest 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
<- optim(par=c(4,5), method="L-BFGS",
parest 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"
<- optim(par=c(4,5), fn=objetivo,
parest 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
<- exp(parest$par)
p p
## [1] 1468.336 15099.145
<- kalman3(y=Nile,0,10^10,Q=p[1],H=p[2])
result 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
<- ts(result$x,start=1871)
Nile.ll 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:
<- 3.6
k <- k * p
p <- kalman3(y=Nile, 0, k*10^10, Q=p[1], H=p[2])
result 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\).
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
).
<- StructTS(x=Nile,type="level", init=c(1469,15100))
ajuste.ll 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
;
$coef ajuste.ll
## 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
:
<- StructTS(x=Nile,type="trend", init=c(1469,1,15100))
ajuste.llt <- ts(fitted(ajuste.llt)[,1],start=1871)
Nile.llt 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):
$coef ajuste.llt
## level slope epsilon
## 1416.832 0.000 15080.406
No hay problema con valores perdidos:
<- Nile
Nile.bis 20:30] <- NA
Nile.bis[<- StructTS(x=Nile.bis,type="trend", init=c(1000, 0.5, 10000))
ajuste.llt 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))
<- StructTS(log10(UKgas), type = "BSM")
UKgas.bsm plot(cbind(fitted(UKgas.bsm), resids=resid(UKgas.bsm)), main = "Consumo de gas en el Reino Unido")
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)
<- dlm(FF=1, GG=1, V=15100, W=1469, m0=0, C0=10^7)
mod.Nilo.ll <- dlmFilter(Nile,mod.Nilo.ll)
res.Nilo.ll plot(res.Nilo.ll$y)
lines(res.Nilo.ll$f,col="red")
<- dlmSmooth(Nile,mod.Nilo.ll)
smo.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
<- dlmModPoly(2, dV=15047, dW=c(1427,0), m0=c(0,0),
mod.Nilo.llt C0=10^7*diag(2))
<- dlmFilter(Nile,mod.Nilo.llt)
res.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:
<- function(x) {
crearMod 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
.
<- dlmMLE(Nile, parm=c(0,0), build=crearMod)
ajuste.Nilo.ll <- dlmMLE(Nile, parm=c(0,0), build=crearMod,
ajuste.Nilo.ll method="Nelder-Mead",
control=list(reltol=10^-10))
<- crearMod(ajuste.Nilo.ll$par)
mod.Nilo.ll 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.
<- dget(file="Cordero.dge")
Cordero 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
<- function(x) {
crearMod <- dlmModPoly(1,dV=exp(x[1]),dW=exp(x[2])) +
mod dlmModSeas(12,dV=0,dW=c(exp(x[3]),rep(0,10)))
return(mod)
}<- dlmMLE(Cordero,parm=c(0,0,0), build=crearMod)
ajuste.Cordero.lls 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"
<- crearMod(ajuste.Cordero.lls$par)
mod.Cordero.lls <- dlmFilter(Cordero,mod.Cordero.lls)
res.Cordero.lls plot(Cordero)
lines(res.Cordero.lls$f,col="red")
<- dlmSmooth(Cordero,mod.Cordero.lls)
smo.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).
<- function(x) {
crearMod <- matrix(ifelse(time(Nile)<1899,0,1),nc=1)
X <- dlmModPoly(1,dV=exp(x[1]),dW=exp(x[2]))
mod1 <- dlmModReg(X,dV=0,addInt=FALSE)
mod2 return(mod1+mod2)
}<- crearMod(c(0,0))
mod.Nilo.salto <- dlmMLE(Nile, parm=c(0,0), build=crearMod)
ajuste.Nilo.salto 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"
<- crearMod(ajuste.Nilo.salto$par)
mod.Nilo.salto 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
<- dlmFilter(Nile,mod.Nilo.salto) res.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.
<- ts(res.Nilo.salto$a[,2], start=c(1871,1))
efecto <- rep(NA,length(efecto))
Var.efecto for (i in 1:length(efecto)) {
<- ( res.Nilo.salto$U.C[[i]] %*% diag(res.Nilo.salto$D.C[i,]^2)
Var.efecto[i] %*% t(res.Nilo.salto$U.C[[i]]) )[2,2]
}plot(efecto)
<- 1.96 * sqrt(Var.efecto)
ancho.ic 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.
<- dlmSmooth(res.Nilo.salto)
res.Nilo.smooth <- ts(res.Nilo.smooth$s[,2], start=c(1871,1))
smooth <- dlmSvd2var(res.Nilo.smooth$U.S, res.Nilo.smooth$D.S)
temp <- rep(NA,length(smooth))
Var.smooth for (i in 1:length(smooth)) {
<- (temp[[i]])[2,2]
Var.smooth[i] }
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.
1] smooth[
## [1] -248.1385
sqrt(Var.smooth[1])
## [1] 28.75747