Ejemplos varios (II)
F. Tusell
Esta versión generada el: 2021-11-18
Paquetes utilizados
library(dlm)
library(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.library(lubridate)##
## Attaching package: 'lubridate'## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, unionlibrary(xts)## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlibrary(zoo)Reconstrucción de desgloses
Sucede en ocasiones que un agregado se observa con precisión, pero sus partes constituyentes no se observan todo el tiempo, o se observan con precisión menor.
Lectura
Emplearemos como ejemplo el desglose de la contratación en acciones del SIB (Sistema de Interconexión Bursátil Español) en cinco categorías:
cont <- read.csv(file="ContSIB.txt", header=TRUE, sep=";")
head(cont)## SIB Bancos Elec Const Comun Resto
## 1991:01 1246610 279735 260048 99669 89666 517501
## 1991:02 2363560 584725 547492 151214 154826 925304
## 1991:03 2509330 547520 511727 189694 163196 1097199
## 1991:04 1949460 395019 421632 155602 109580 867618
## 1991:05 2211420 763187 406980 163504 134659 743090
## 1991:06 1910650 458482 502456 148459 94741 706515cont <- zoo(cont, order.by=as.yearmon(rownames(cont), "%Y:%m"))
head(cont)## SIB Bancos Elec Const Comun Resto
## ene 1991 1246610 279735 260048 99669 89666 517501
## feb 1991 2363560 584725 547492 151214 154826 925304
## mar 1991 2509330 547520 511727 189694 163196 1097199
## abr 1991 1949460 395019 421632 155602 109580 867618
## may 1991 2211420 763187 406980 163504 134659 743090
## jun 1991 1910650 458482 502456 148459 94741 706515La primera columna es el total; las siguientes cinco, el desglose. Hay minúsculos descuadres de los que no nos preocupamos. La apariencia de las series es la que sigue:
plot(cont)
plot(as.xts(cont))
Introducción de valores perdidos
Vamos a introducir valores perdidos en diferentes intervalos en las series distintas del total. Solaparemos observaciones perdidas, de modo que no sea posible recuperar una componente por mera diferencia.
Podemos seleccionar observaciones por su número de fila o emplear fechas. Hay que recordar que cuando se fecha con año+mes (as.yearmon), los meses son fracciones de año:
as.yearmon("ene 2000") + 1/12## [1] "feb 2000"as.yearmon("ene 2000") + (1:4) / 12## [1] "feb 2000" "mar 2000" "abr 2000" "may 2000"También se puede emplear la función window para seleccionar rangos de fechas. En el código siguiente se emplean los tres métodos para insertar valores perdidos. Previamente, se almacenan las series completas en cont.orig:
cont.orig <- cont
cont[30:40, 2] <- NA
cont[35:55, 3] <- NA
cont[as.yearmon("ene 2000")+ (1:24)/12, 2] <- NA
cont[as.yearmon("ene 2001")+ (1:24)/12, 3] <- NA
cont[as.yearmon("ene 2010")+ (1:24)/12, 4] <- NA
cont[as.yearmon("ene 2011")+ (1:24)/12, 5] <- NA
cont[as.yearmon("ene 2013")+ (1:24)/12, 6] <- NA
cont[as.yearmon("ene 2014")+ (1:24)/12, 2] <- NA
window(cont[,"Comun"],
start=as.yearmon("ene 2006"),
end=as.yearmon("dic 2007")) <- NAplot(cont)
plot(as.xts(cont))
Modelización
mod1 <- SSModel(cont ~
-1 # Sin el nivel "de oficio"
+ SSMcustom(
Z=rbind(rep(1,5), diag(5)), # Matriz de observación
T=diag(5), # Matriz transición
Q=diag(5), # Covarianzas del estado
P1inf=diag(5) # Priori difusa
),
H=diag(6) # Covarianzas obs
)Examinemos la estructura del objeto resultante:
mod1## Call:
## SSModel(formula = cont ~ -1 + SSMcustom(Z = rbind(rep(1, 5),
## diag(5)), T = diag(5), Q = diag(5), P1inf = diag(5)), H = diag(6))
##
## State space model object of class SSModel
##
## Dimensions:
## [1] Number of time points: 333
## [1] Number of time series: 6
## [1] Number of disturbances: 5
## [1] Number of states: 5
## Names of the states:
## [1] custom1 custom2 custom3 custom4 custom5
## Distributions of the time series:
## [1] gaussian
##
## Object is a valid object of class SSModel.mod1$Z## , , 1
##
## custom1 custom2 custom3 custom4 custom5
## SIB 1 1 1 1 1
## Bancos 1 0 0 0 0
## Elec 0 1 0 0 0
## Const 0 0 1 0 0
## Comun 0 0 0 1 0
## Resto 0 0 0 0 1mod1$T## , , 1
##
## custom1 custom2 custom3 custom4 custom5
## custom1 1 0 0 0 0
## custom2 0 1 0 0 0
## custom3 0 0 1 0 0
## custom4 0 0 0 1 0
## custom5 0 0 0 0 1mod1$Q## , , 1
##
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 0 0 0 0
## [2,] 0 1 0 0 0
## [3,] 0 0 1 0 0
## [4,] 0 0 0 1 0
## [5,] 0 0 0 0 1mod1$H## , , 1
##
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 1 0 0 0 0 0
## [2,] 0 1 0 0 0 0
## [3,] 0 0 1 0 0 0
## [4,] 0 0 0 1 0 0
## [5,] 0 0 0 0 1 0
## [6,] 0 0 0 0 0 1Las matrices de covarianzas contienen los únicos parámetros que hay en el modelo. Tiene sentido en este problema pensar en ruidos de observación de muy reducida varianza, ya que las observaciones son practicamente exactas. Aunque las observaciones sean prácticamente exactas, no conviene fijar los términos diagonales de \(H\) a cero, en evitación de problemas numéricos.
Para 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. Supondremos que las precisiones en la observación de las series es similar y la varianza de los respectivos estados, en cambio, específica de cada uno.
updatefn <- function(pars, model) {
model["Q"] <- diag(exp(pars[1:5]))
return(model)
}Una vez definida esta función, proporcionamos valores iniciales para los dos parámetros en init e invocamos fitSSM.
init <- rep(1, 5)
fit <- fitSSM(mod1, inits=init, updatefn=updatefn)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] [,3] [,4] [,5] [,6]
## [1,] 1 0 0 0 0 0
## [2,] 0 1 0 0 0 0
## [3,] 0 0 1 0 0 0
## [4,] 0 0 0 1 0 0
## [5,] 0 0 0 0 1 0
## [6,] 0 0 0 0 0 1mod1$Q## , , 1
##
## [,1] [,2] [,3] [,4] [,5]
## [1,] 9986294 0.000 0.0000 0.00 0
## [2,] 0 1447.889 0.0000 0.00 0
## [3,] 0 0.000 444.3605 0.00 0
## [4,] 0 0.000 0.0000 55747.13 0
## [5,] 0 0.000 0.0000 0.00 9994469mod1.out <- KFS(mod1)Reconstrucción
En mod1.out tenemos ahora bastante información, que incluye los valores suavizados del estado (en el componente alphahat).
head(mod1.out$alphahat)## custom1 custom2 custom3 custom4 custom5
## [1,] 279681.0 260192.2 99730.82 89613.10 517447.0
## [2,] 584767.1 547311.2 151226.48 154867.07 925346.1
## [3,] 547552.5 511721.9 189564.01 163227.44 1097231.5
## [4,] 394996.0 421661.1 155672.71 109558.35 867595.0
## [5,] 763183.0 407052.0 163448.86 134653.93 743086.1
## [6,] 458522.6 502293.9 148454.88 94782.42 706555.6cont.r <- cbind(as.numeric(cont[,1]),mod1.out$alphahat)
colnames(cont.r) <- colnames(cont)
cont.r <- zoo(cont.r, order.by=index(cont))
plot(cont.r)
plot(as.xts(cont.r))
Podemos examinar con más detalle la reconstrucción efectuada representando los valores observados de las series y añadiendo en rojo los valores reconstruidos:
for (i in 2:6) {
plot(cont.r[,i], col="red", ylab=colnames(cont)[i])
lines(cont[,i])
legend("topleft",c("Obs","Reconst"), col=1:2, lty=c(1,1))
}
Comparemos ahora los valores reales con los reconstruidos:
for (i in 2:6) {
plot(cont.r[,i], col="red", ylab=colnames(cont)[i])
lines(cont.orig[,i], col="green")
legend("topleft",c("Original","Reconst"),
col=c("green","red"), lty=c(1,1))
}
Homogeneización de series
(Extracto de Pérez-Castroviejo and Tusell (2007))
load("PreciosVizcaya.rda")
colnames(PreciosVizcaya)## [1] "Bilbao.Pan.de.trigo.(Ptas/Kg)"
## [2] "Bilbao.Carne.de.ganado.vacuno"
## [3] "Bilbao.Carne.de.ganado.lanar"
## [4] "Bilbao.Carne.de.ganado.de.cerda"
## [5] "Bilbao.Manteca.(Ptas/Kg)"
## [6] "Bilbao.Tocino"
## [7] "Bilbao.Bacalao"
## [8] "Bilbao.Sardinas.saladas"
## [9] "Bilbao.Sardinas.frescas(docena)"
## [10] "Bilbao.Pesca.fresca.ordinaria"
## [11] "Bilbao.Arroz"
## [12] "Bilbao.Garbanzos"
## [13] "Bilbao.Patatas"
## [14] "Bilbao.Judías"
## [15] "Bilbao.Jabón"
## [16] "Bilbao.Lentejas"
## [17] "Bilbao.Sal"
## [18] "Bilbao.Huevos.(docena)"
## [19] "Bilbao.Azúcar"
## [20] "Bilbao.Café"
## [21] "Bilbao.Vinos.de.Jerez.(l.)"
## [22] "Bilbao.Vino.común.(l.)"
## [23] "Bilbao.Vino.rancio.(l.)"
## [24] "Bilbao.Sidra.(l)"
## [25] "Bilbao.Chacolí.(l.)"
## [26] "Bilbao.Aceite.(l.)"
## [27] "Bilbao.Leche.(l.)"
## [28] "Bilbao.Leña.(Kg)"
## [29] "Bilbao.Carbón.vegetal.(Kg)"
## [30] "Bilbao.Carbón.mineral.(Kg)"
## [31] "Bilbao.Cok.(Kg)"
## [32] "Bilbao.Petróleo.(l.)"
## [33] "Bilbao.Fluido.eléctrico.(5.bujías.mes)"
## [34] "Bilbao.Gas.(m3)"
## [35] "Bilbao.Merluza"
## [36] "Bilbao.Gallina.(unidad)"
## [37] "Barak91.Pan.Trigo.2"
## [38] "Barak91.Pan.trigo.3"
## [39] "Barak91.Pan.maiz"
## [40] "Barak91.Carne.vacuno"
## [41] "Barak91.Carne.lanar"
## [42] "Barak91.Carne.cerdo"
## [43] "Barak91.Tocino"
## [44] "Barak91.Bacalao"
## [45] "Barak91.Arroz"
## [46] "Barak91.Garbanzos"
## [47] "Barak91.Patatas"
## [48] "Barak91.Judias"
## [49] "Barak91.Habas"
## [50] "Barak91.Lentejas"
## [51] "Barak91.Aceite"
## [52] "Barak91.Vino"
## [53] "Barak91.Chacoli"
## [54] "Barak91.Sidra"
## [55] "Barak.Pan.de.trigo.(Kg)"
## [56] "Barak.Pan.de.cebada"
## [57] "Barak.Pan.de.maiz"
## [58] "Barak.Pan.de.centeno"
## [59] "Barak.Harina.de.trigo.(kg)"
## [60] "Barak.Harina.de.maíz"
## [61] "Barak.Harina.de.centeno"
## [62] "Barak.Carne.de.vaca.(Kg)"
## [63] "Barak.Carne.de.carnero.u.oveja"
## [64] "Barak.Carne.de.cabra"
## [65] "Barak.Carne.de.cerdo"
## [66] "Barak.Carne.despojos.de.reses"
## [67] "Barak.Tocino(Kg)"
## [68] "Barak.Embutidos.(Kg)"
## [69] "Barak.Pesaco.fresco.(Kg)"
## [70] "Barak.Pescado.en.escabeche.(Kg)"
## [71] "Barak.Bacalao.(Kg)"
## [72] "Barak.Frutas.(Kg)"
## [73] "Barak.Hortalizas.(Kg)"
## [74] "Barak.Patatas.(Kg)"
## [75] "Barak.Habas.(Kg)"
## [76] "Barak.Garbanzos.(Kg)"
## [77] "Barak.Arroz.(Kg)"
## [78] "Barak.Judías.(Kg)"
## [79] "Barak.Almortas.(Kg)"
## [80] "Barak.Lentejas.(Kg)"
## [81] "Barak.Vino.(l.)"
## [82] "Barak.Sidra.(l.)"
## [83] "Barak.Leche.(l.)"
## [84] "Barak.Petróleo.(l.)"
## [85] "Barak.Luz.(mes)"
## [86] "Barak.Carbón.vegetal.(kg)"
## [87] "Barak.Leña.(Kg)"
## [88] "Barak.Café.(Kg)"
## [89] "Barak.Huevos.(docena)"
## [90] "Barak.Azúcar.(Kg)"
## [91] "Barak.Jabón.(Kg)"
## [92] "Barak.Sal.(Kg)"
## [93] "Barak.Aceite.(l.)"
## [94] "Barak.Manteca.de.vaca.(Kg)"
## [95] "Barak.Manteca.de.cerdo.(Kg)"
## [96] "Barak.Achicoria.(Kg)"
## [97] "Barak.Velas.(Kg)"
## [98] "Barak.Carbón.mineral.(Kg)"
## [99] "Barak.Merluza"
## [100] "Basurto.PAN.ptas.kg.."
## [101] "Basurto.CARNE.s.hues."
## [102] "Basurto.CARNE.c.hues."
## [103] "Basurto.CARNE.Terne."
## [104] "Basurto.Chuletas.buey"
## [105] "Basurto.VINO..Ptas.l.."
## [106] "Basurto.ACEITE.Pts.l"
## [107] "Basurto.ALUBIA.Pts.k"
## [108] "Basurto.GARBAN.Pts.k"
## [109] "Basurto.PATATA.Pts.k"
## [110] "Basurto.LECHE.Pts.l."
## [111] "Basurto.HUEVOS.Pts.d."
## [112] "Basurto.ARROZ.Pts.k.."
## [113] "Basurto.AZUCAR.PtsK."
## [114] "Basurto.CAFE.Ptas.Kg"
## [115] "Basurto.BACALAO.P.K."
## [116] "Basurto.MERLUZ.P.K.."
## [117] "Basurto.BESUGO.P.K."
## [118] "Basurto.BONITO.P.K.."
## [119] "Basurto.SARDINA.una."
## [120] "Basurto.ANCHOA.P.K.."
## [121] "Basurto.ANGULA.P.K.."
## [122] "Basurto.JABON.ptas.Kg."
## [123] "Basurto.CARB.MIN..Ptas.Kg."
## [124] "Basurto.LEÑA.Ptas.Kg."
## [125] "Valmaseda.Garbanzo.Ptas.Kg"
## [126] "Valmaseda.Arroz.Ptas.Kg"
## [127] "Valmaseda.Aceite.ptas.l"
## [128] "Valmaseda.Vino.Ptas.l"
## [129] "Valmaseda.Vaca.Ptas.Kg"
## [130] "Valmaseda.Tocino.Ptas.Kg"
## [131] "AHV.Pan"
## [132] "AHV.Carne"
## [133] "AHV.Tocino"
## [134] "AHV.Alubias"
## [135] "AHV.Garbanzos"
## [136] "AHV.Aceite"
## [137] "AHV.Vino"
## [138] "AHV.Bacalao"
## [139] "AHV.Arroz"
## [140] "AHV.Patatas"
## [141] "AHV.Azúcar"
## [142] "AHV.Café"
## [143] "AHV.Jabón"
## [144] "AHV.Alpargatas.par."
## [145] "AHV.Botas.par."
## [146] "AHV.Jubones.uno."
## [147] "AHV.Calcetines.par."
## [148] "AHV.Terliz.m.."
## [149] "AHV.Mahón.m."
## [150] "AHV.Malaño.m."
## [151] "AHV.Francesilla.m."
## [152] "Textil.Lana"
## [153] "Textil.Algodón"
## [154] "Textil.Zapatillas"
## [155] "Textil.Lienzo"
## [156] "Textil.Terliz"
## [157] "Textil.Cretona"
## [158] "Textil.Mahón"
## [159] "Textil.MahónBaracaldo"
## [160] "Textil.Yute"
## [161] "Textil.Francesilla"
## [162] "Textil.Muselina"
## [163] "Textil.Cutí"
## [164] "Textil.Percal"
## [165] "Textil.Dril"
## [166] "Compl.Registros.alquiler"
## [167] "Compl.Jabón"
## [168] "Compl.Leña1"
## [169] "Compl.Carbonvegetal1"
## [170] "Compl.Carbon.mineral"
## [171] "Compl.Electricidad"class(PreciosVizcaya)## [1] "zoo"tmp <- grep("[Vv]aca|[Vv]acuno|[Bb]uey|Terne", colnames(PreciosVizcaya))
colnames(PreciosVizcaya)[tmp]## [1] "Bilbao.Carne.de.ganado.vacuno" "Barak91.Carne.vacuno"
## [3] "Barak.Carne.de.vaca.(Kg)" "Barak.Manteca.de.vaca.(Kg)"
## [5] "Basurto.CARNE.Terne." "Basurto.Chuletas.buey"
## [7] "Valmaseda.Vaca.Ptas.Kg"Vacuno <- PreciosVizcaya[,tmp[c(1,2,3,5,6,7)]]
plot(Vacuno)
Vacuno <- zoo(coredata(Vacuno),
order.by=year(index(Vacuno)))
head(Vacuno)## Bilbao.Carne.de.ganado.vacuno Barak91.Carne.vacuno
## 1862 NA NA
## 1863 NA NA
## 1864 NA NA
## 1865 NA NA
## 1866 NA NA
## 1867 NA NA
## Barak.Carne.de.vaca.(Kg) Basurto.CARNE.Terne. Basurto.Chuletas.buey
## 1862 NA NA NA
## 1863 NA NA NA
## 1864 NA NA NA
## 1865 NA NA NA
## 1866 NA NA NA
## 1867 NA NA NA
## Valmaseda.Vaca.Ptas.Kg
## 1862 0.89
## 1863 0.89
## 1864 0.89
## 1865 0.89
## 1866 0.89
## 1867 0.90Vacuno <- window(Vacuno, start=1862, end=1938)
Vacuno## Bilbao.Carne.de.ganado.vacuno Barak91.Carne.vacuno
## 1862 NA NA
## 1863 NA NA
## 1864 NA NA
## 1865 NA NA
## 1866 NA NA
## 1867 NA NA
## 1868 NA NA
## 1869 NA NA
## 1870 NA NA
## 1871 NA NA
## 1872 NA NA
## 1873 NA NA
## 1874 NA NA
## 1875 NA NA
## 1876 NA NA
## 1877 NA NA
## 1878 NA NA
## 1879 NA NA
## 1880 NA NA
## 1881 NA NA
## 1882 NA NA
## 1883 NA NA
## 1884 NA NA
## 1885 NA NA
## 1886 NA NA
## 1887 NA NA
## 1888 NA NA
## 1889 NA NA
## 1890 NA NA
## 1891 NA 1.40
## 1892 NA 1.60
## 1893 NA 1.40
## 1894 NA 1.50
## 1895 NA 1.45
## 1896 NA 1.45
## 1897 NA 1.55
## 1898 NA 1.20
## 1899 NA 1.30
## 1900 NA NA
## 1901 NA NA
## 1902 NA NA
## 1903 NA NA
## 1904 NA NA
## 1905 NA NA
## 1906 NA NA
## 1907 NA NA
## 1908 NA NA
## 1909 NA NA
## 1910 NA NA
## 1911 NA NA
## 1912 NA NA
## 1913 1.55 NA
## 1914 1.93 NA
## 1915 1.98 NA
## 1916 2.25 NA
## 1917 2.35 NA
## 1918 2.60 NA
## 1919 3.30 NA
## 1920 3.30 NA
## 1921 3.25 NA
## 1922 3.60 NA
## 1923 3.60 NA
## 1924 3.20 NA
## 1925 3.25 NA
## 1926 2.45 NA
## 1927 2.45 NA
## 1928 2.30 NA
## 1929 2.85 NA
## 1930 2.20 NA
## 1931 3.40 NA
## 1932 3.30 NA
## 1933 2.60 NA
## 1934 2.60 NA
## 1935 2.60 NA
## 1936 3.60 NA
## 1937 4.50 NA
## 1938 4.35 NA
## Barak.Carne.de.vaca.(Kg) Basurto.CARNE.Terne. Basurto.Chuletas.buey
## 1862 NA NA NA
## 1863 NA NA NA
## 1864 NA NA NA
## 1865 NA NA NA
## 1866 NA NA NA
## 1867 NA NA NA
## 1868 NA NA NA
## 1869 NA NA NA
## 1870 NA NA NA
## 1871 NA NA NA
## 1872 NA NA NA
## 1873 NA NA NA
## 1874 NA NA NA
## 1875 NA NA NA
## 1876 NA NA NA
## 1877 NA NA NA
## 1878 NA NA NA
## 1879 NA NA NA
## 1880 NA NA NA
## 1881 NA NA NA
## 1882 NA NA NA
## 1883 NA NA NA
## 1884 NA 2.25 NA
## 1885 NA NA NA
## 1886 NA 2.50 NA
## 1887 NA 2.31 NA
## 1888 NA 2.25 NA
## 1889 NA 2.25 NA
## 1890 NA 2.24 NA
## 1891 NA 2.25 NA
## 1892 NA 2.25 NA
## 1893 NA 2.25 NA
## 1894 NA 2.25 NA
## 1895 NA 2.25 NA
## 1896 NA 2.25 NA
## 1897 NA 2.25 NA
## 1898 NA NA NA
## 1899 NA NA NA
## 1900 NA NA NA
## 1901 NA NA NA
## 1902 NA NA NA
## 1903 NA NA NA
## 1904 NA NA NA
## 1905 NA NA NA
## 1906 1.50 NA NA
## 1907 1.50 NA NA
## 1908 1.55 NA 2.35
## 1909 1.54 NA 2.35
## 1910 1.43 NA 2.35
## 1911 1.50 NA 2.35
## 1912 1.45 NA 2.35
## 1913 1.70 NA 2.35
## 1914 1.53 NA 2.35
## 1915 1.75 1.85 2.35
## 1916 1.90 1.94 2.35
## 1917 1.90 1.88 2.45
## 1918 3.60 2.20 2.76
## 1919 2.90 2.81 3.01
## 1920 2.90 NA 3.29
## 1921 2.80 NA 3.60
## 1922 2.40 NA 3.48
## 1923 2.60 2.91 3.30
## 1924 3.00 3.72 3.45
## 1925 2.80 3.38 3.50
## 1926 2.80 3.19 3.50
## 1927 2.80 3.24 3.50
## 1928 NA 2.83 3.18
## 1929 NA 4.38 3.28
## 1930 NA 6.86 NA
## 1931 NA 5.55 3.45
## 1932 NA 5.42 3.32
## 1933 NA 5.13 3.16
## 1934 NA 5.21 3.08
## 1935 NA 5.26 3.06
## 1936 NA NA NA
## 1937 NA NA NA
## 1938 NA NA NA
## Valmaseda.Vaca.Ptas.Kg
## 1862 0.89
## 1863 0.89
## 1864 0.89
## 1865 0.89
## 1866 0.89
## 1867 0.90
## 1868 0.90
## 1869 0.87
## 1870 0.84
## 1871 0.84
## 1872 0.84
## 1873 NA
## 1874 NA
## 1875 NA
## 1876 0.98
## 1877 1.07
## 1878 1.00
## 1879 1.05
## 1880 1.10
## 1881 1.21
## 1882 1.14
## 1883 1.10
## 1884 1.18
## 1885 1.12
## 1886 1.20
## 1887 1.12
## 1888 1.10
## 1889 1.10
## 1890 1.20
## 1891 NA
## 1892 NA
## 1893 NA
## 1894 NA
## 1895 NA
## 1896 NA
## 1897 NA
## 1898 NA
## 1899 NA
## 1900 NA
## 1901 NA
## 1902 NA
## 1903 NA
## 1904 NA
## 1905 NA
## 1906 NA
## 1907 NA
## 1908 NA
## 1909 NA
## 1910 NA
## 1911 NA
## 1912 NA
## 1913 NA
## 1914 NA
## 1915 NA
## 1916 NA
## 1917 NA
## 1918 NA
## 1919 NA
## 1920 NA
## 1921 NA
## 1922 NA
## 1923 NA
## 1924 NA
## 1925 NA
## 1926 NA
## 1927 NA
## 1928 NA
## 1929 NA
## 1930 NA
## 1931 NA
## 1932 NA
## 1933 NA
## 1934 NA
## 1935 NA
## 1936 NA
## 1937 NA
## 1938 NAmod2 <- SSModel(Vacuno ~
-1 +
SSMcustom(Z=matrix(1,6,1),
T=diag(1),
Q=diag(1),
P1inf=diag(1)
),
H=diag(6)
)updatefn <- function(pars, model) {
model["Q"] <- matrix(exp(pars[1]),1,1)
model["Z"] <- matrix(c(exp(pars[2:6]),1),6,1)
model["H"] <- diag(exp(pars[7:12]))
return(model)
}Una vez definida esta función, proporcionamos valores iniciales para los dos parámetros en init e invocamos fitSSM.
init <- c(-.5,rep(0,11))
fit <- fitSSM(mod2, inits=init, updatefn=updatefn)
fit## $optim.out
## $optim.out$par
## [1] -3.1658879 0.5839930 -1.2067838 0.3827767 0.8495147 0.6723116
## [7] -4.0025694 4.6865276 -0.7076722 0.6705105 -1.0752786 -1.5611455
##
## $optim.out$value
## [1] 144.495
##
## $optim.out$counts
## function gradient
## 501 NA
##
## $optim.out$convergence
## [1] 1
##
## $optim.out$message
## NULL
##
##
## $model
## Call:
## SSModel(formula = Vacuno ~ -1 + SSMcustom(Z = matrix(1, 6, 1),
## T = diag(1), Q = diag(1), P1inf = diag(1)), H = diag(6))
##
## State space model object of class SSModel
##
## Dimensions:
## [1] Number of time points: 77
## [1] Number of time series: 6
## [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.mod2 <- updatefn(fit$optim.out$par, mod2)
mod2## Call:
## SSModel(formula = Vacuno ~ -1 + SSMcustom(Z = matrix(1, 6, 1),
## T = diag(1), Q = diag(1), P1inf = diag(1)), H = diag(6))
##
## State space model object of class SSModel
##
## Dimensions:
## [1] Number of time points: 77
## [1] Number of time series: 6
## [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.mod2$Z## , , 1
##
## custom1
## Bilbao.Carne.de.ganado.vacuno 1.7931843
## Barak91.Carne.vacuno 0.2991579
## Barak.Carne.de.vaca.(Kg) 1.4663505
## Basurto.CARNE.Terne. 2.3385118
## Basurto.Chuletas.buey 1.9587599
## Valmaseda.Vaca.Ptas.Kg 1.0000000str(mod2)## List of 14
## $ y : Time-Series [1:77, 1:6] from 1 to 77: NA NA NA NA NA NA NA NA NA NA ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:77] "1862-07-01" "1863-07-01" "1864-07-01" "1865-07-01" ...
## .. ..$ : chr [1:6] "Bilbao.Carne.de.ganado.vacuno" "Barak91.Carne.vacuno" "Barak.Carne.de.vaca.(Kg)" "Basurto.CARNE.Terne." ...
## ..- attr(*, "index")= num [1:77] 1862 1863 1864 1865 1866 ...
## $ Z : num [1:6, 1, 1] 1.793 0.299 1.466 2.339 1.959 ...
## ..- attr(*, "dimnames")=List of 3
## .. ..$ : chr [1:6] "Bilbao.Carne.de.ganado.vacuno" "Barak91.Carne.vacuno" "Barak.Carne.de.vaca.(Kg)" "Basurto.CARNE.Terne." ...
## .. ..$ : chr "custom1"
## .. ..$ : NULL
## $ H : num [1:6, 1:6, 1] 0.0183 0 0 0 0 ...
## $ T : num [1, 1, 1] 1
## ..- attr(*, "dimnames")=List of 3
## .. ..$ : chr "custom1"
## .. ..$ : chr "custom1"
## .. ..$ : NULL
## $ R : num [1, 1, 1] 1
## ..- attr(*, "dimnames")=List of 3
## .. ..$ : chr "custom1"
## .. ..$ : NULL
## .. ..$ : NULL
## $ Q : num [1, 1, 1] 0.0422
## $ a1 : num [1, 1] 0
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr "custom1"
## .. ..$ : NULL
## $ P1 : num [1, 1] 0
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr "custom1"
## .. ..$ : chr "custom1"
## $ P1inf : num [1, 1] 1
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr "custom1"
## .. ..$ : chr "custom1"
## $ u : chr "Omitted"
## $ distribution: chr [1:6] "gaussian" "gaussian" "gaussian" "gaussian" ...
## $ tol : num 1.49e-08
## $ call : language SSModel(formula = Vacuno ~ -1 + SSMcustom(Z = matrix(1, 6, 1), T = diag(1), Q = diag(1), P1inf = diag(1)), H = diag(6))
## $ terms :Classes 'terms', 'formula' language Vacuno ~ -1 + SSMcustom(Z = matrix(1, 6, 1), T = diag(1), Q = diag(1), P1inf = diag(1))
## .. ..- attr(*, "variables")= language list(Vacuno, SSMcustom(Z = matrix(1, 6, 1), T = diag(1), Q = diag(1), P1inf = diag(1)))
## .. ..- attr(*, "factors")= int [1:2, 1] 0 1
## .. .. ..- attr(*, "dimnames")=List of 2
## .. .. .. ..$ : chr [1:2] "Vacuno" "SSMcustom(Z = matrix(1, 6, 1), T = diag(1), Q = diag(1), P1inf = diag(1))"
## .. .. .. ..$ : chr "SSMcustom(Z = matrix(1, 6, 1), T = diag(1), Q = diag(1), P1inf = diag(1))"
## .. ..- attr(*, "term.labels")= chr "SSMcustom(Z = matrix(1, 6, 1), T = diag(1), Q = diag(1), P1inf = diag(1))"
## .. ..- attr(*, "specials")=Dotted pair list of 6
## .. .. ..$ SSMregression: NULL
## .. .. ..$ SSMtrend : NULL
## .. .. ..$ SSMseasonal : NULL
## .. .. ..$ SSMcycle : NULL
## .. .. ..$ SSMarima : NULL
## .. .. ..$ SSMcustom : int 2
## .. ..- attr(*, "order")= int 1
## .. ..- attr(*, "intercept")= int 0
## .. ..- attr(*, "response")= int 1
## .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
## - attr(*, "class")= chr "SSModel"
## - attr(*, "p")= int 6
## - attr(*, "m")= int 1
## - attr(*, "k")= int 1
## - attr(*, "n")= int 77
## - attr(*, "tv")= int [1:5] 0 0 0 0 0
## - attr(*, "state_types")= chr "custom"
## - attr(*, "eta_types")= chr "custom"mod2.out <- KFS(mod2)
Vacuno.std <- zoo(mod2.out$alphahat,
order.by=index(Vacuno))
plot(Vacuno.std, main="Valores suavizados carne de vacuno", xlab="Año")
Regresión con parámetros cambiantes
series <- read.table("17_Industry_Portfolios_Daily.txt",skip=9,header=TRUE)
series <- zoo(series,as.Date(rownames(series),format="%Y%m%d"))
head(series)## Food Mines Oil Clths Durbl Chems Cnsum Cnstr Steel FabPr Machn
## 1963-07-01 -0.38 -0.96 -0.88 -0.28 -1.05 -0.61 -0.33 -0.51 -0.84 -0.92 -1.37
## 1963-07-02 0.52 0.22 0.99 0.46 1.17 0.65 0.55 0.89 0.72 0.45 1.62
## 1963-07-03 0.66 0.47 0.89 0.48 1.04 0.72 0.51 0.55 0.45 0.47 0.75
## 1963-07-05 0.14 0.39 0.66 0.34 0.33 0.21 0.80 0.26 0.23 0.44 0.52
## 1963-07-08 -0.21 -0.74 -0.54 -0.48 -0.72 -0.65 -0.49 -0.43 -0.90 -0.37 -1.17
## 1963-07-09 0.32 0.90 0.74 0.54 0.72 0.27 0.25 0.06 0.88 0.35 0.86
## Cars Trans Utils Rtail Finan Other
## 1963-07-01 -0.49 -0.92 -0.18 -0.32 -0.46 -0.65
## 1963-07-02 0.87 0.97 0.49 0.43 0.49 0.78
## 1963-07-03 0.15 0.63 0.12 0.82 0.51 0.88
## 1963-07-05 0.57 -0.01 0.13 0.06 0.50 0.60
## 1963-07-08 -1.24 -1.04 0.02 -0.04 -0.40 -0.75
## 1963-07-09 0.31 0.99 0.21 0.22 0.62 0.28FFfact <- read.table("F-F_Research_Data_Factors_daily.txt",skip=3,header=TRUE)
FFfact <- zoo(FFfact,as.Date(rownames(FFfact),format="%Y%m%d"))
head(FFfact)## Mkt.RF SMB HML RF
## 1963-07-01 -0.66 0.16 -0.32 0.012
## 1963-07-02 0.78 -0.28 0.26 0.012
## 1963-07-03 0.63 -0.18 -0.10 0.012
## 1963-07-05 0.40 0.09 -0.29 0.012
## 1963-07-08 -0.62 0.05 -0.20 0.012
## 1963-07-09 0.45 0.01 0.11 0.012Según un modelo de valoración, el exceso de rendimiento de un valor (o grupo de valores) estaría relacionado con el exceso de rendimiento del mercado así:
\[E_v = \beta E_m \]
Los excesos de de rendimiento de un activo son los rendimientos del mismo menos el tipo de interés sin riesgo (RF, “risk free” en FFfact). Los excesos de rendimiento del mercado aparecen ya calculados en FFfact como Mkt.RF. Los de lps diferentes sectores se calculan así:
for (i in 1:ncol(series)) {
series[,i] <- series[,i] - FFfact[,"RF"]
}
series <- cbind(series,Mkt.RF=FFfact[,"Mkt.RF"])
head(series)## Food Mines Oil Clths Durbl Chems Cnsum Cnstr Steel
## 1963-07-01 -0.392 -0.972 -0.892 -0.292 -1.062 -0.622 -0.342 -0.522 -0.852
## 1963-07-02 0.508 0.208 0.978 0.448 1.158 0.638 0.538 0.878 0.708
## 1963-07-03 0.648 0.458 0.878 0.468 1.028 0.708 0.498 0.538 0.438
## 1963-07-05 0.128 0.378 0.648 0.328 0.318 0.198 0.788 0.248 0.218
## 1963-07-08 -0.222 -0.752 -0.552 -0.492 -0.732 -0.662 -0.502 -0.442 -0.912
## 1963-07-09 0.308 0.888 0.728 0.528 0.708 0.258 0.238 0.048 0.868
## FabPr Machn Cars Trans Utils Rtail Finan Other Mkt.RF
## 1963-07-01 -0.932 -1.382 -0.502 -0.932 -0.192 -0.332 -0.472 -0.662 -0.66
## 1963-07-02 0.438 1.608 0.858 0.958 0.478 0.418 0.478 0.768 0.78
## 1963-07-03 0.458 0.738 0.138 0.618 0.108 0.808 0.498 0.868 0.63
## 1963-07-05 0.428 0.508 0.558 -0.022 0.118 0.048 0.488 0.588 0.40
## 1963-07-08 -0.382 -1.182 -1.252 -1.052 0.008 -0.052 -0.412 -0.762 -0.62
## 1963-07-09 0.338 0.848 0.298 0.978 0.198 0.208 0.608 0.268 0.45Los betas estáticos se calcularían así:
for (i in 1:17) {
mod <- lm(series[,i] ~ -1 + Mkt.RF, data=series)
cat("Beta sector",colnames(series)[i],coefficients(mod),"\n")
}## Beta sector Food 0.7109178
## Beta sector Mines 0.8944075
## Beta sector Oil 0.9165465
## Beta sector Clths 0.9065931
## Beta sector Durbl 1.010685
## Beta sector Chems 0.9861027
## Beta sector Cnsum 0.8612434
## Beta sector Cnstr 1.040344
## Beta sector Steel 1.18822
## Beta sector FabPr 0.856544
## Beta sector Machn 1.289164
## Beta sector Cars 1.061703
## Beta sector Trans 0.9860509
## Beta sector Utils 0.5781272
## Beta sector Rtail 0.9647996
## Beta sector Finan 1.025672
## Beta sector Other 1.014342Hay abundante evidencia de que los betas no son constantes. Podríamos entonces plantearnos una regresión con un beta variable para ver su evolución. Trandicionalmente se han empleado “rolling windows”: podemos hacer uso de un modelo en espacio de estado. Por ejemplo, para el sector 1 (Food):
Especificamos el modelo. Uno de los raros casos en que queremos una regresión sin térm ino constante.
require(dlm)
mod <- dlmModReg(series[,"Mkt.RF"], addInt=FALSE)
mod ## $FF
## [,1]
## [1,] 1
##
## $V
## [,1]
## [1,] 1
##
## $GG
## [,1]
## [1,] 1
##
## $W
## [,1]
## [1,] 0
##
## $JFF
## [,1]
## [1,] 1
##
## $X
## X
## 1963-07-01 -0.66
## 1963-07-02 0.78
## ...
##
## $m0
## [1] 0
##
## $C0
## [,1]
## [1,] 1e+07Especificamos la función que reemplaza los valores de los parámetros en el modelo precedente. Nótese que dejamos la varianza V a 1; empleamos la verosimilitud condensada.
crearMod <- function(x) {
mod$W[1,1] <- exp(x[1])
return(mod)
}Maximizamos la verosimilitud y reemplazamos el único valor estimado en el modelo:
res <- dlmMLE(series[,"Food"], par=c(-2), build=crearMod)
res## $par
## [1] -9.207284
##
## $value
## [1] 1609.698
##
## $counts
## function gradient
## 10 10
##
## $convergence
## [1] 0
##
## $message
## [1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"mod <- crearMod(res$par)
mod## $FF
## [,1]
## [1,] 1
##
## $V
## [,1]
## [1,] 1
##
## $GG
## [,1]
## [1,] 1
##
## $W
## [,1]
## [1,] 0.0001003061
##
## $JFF
## [,1]
## [1,] 1
##
## $X
## X
## 1963-07-01 -0.66
## 1963-07-02 0.78
## ...
##
## $m0
## [1] 0
##
## $C0
## [,1]
## [1,] 1e+07Filtramos y suavizamos. Los valores filtrados serían de utilidad a un gestor en tiempo real; los valores suavizados son estimaciones retrospectivas.
filtered <- dlmFilter(series[,"Food"], mod)
Food.beta <- filtered$m
head(Food.beta)## 1963-07-01 1963-07-02 1963-07-03 1963-07-05 1963-07-08
## 0.0000000 0.5939393 0.6273569 0.7378824 0.6961089 0.6306361tail(Food.beta)## 2009-05-21 2009-05-22 2009-05-26 2009-05-27 2009-05-28 2009-05-29
## 0.5836754 0.5836754 0.5836754 0.5836754 0.5836754 0.5836754Food.beta <- Food.beta[-1]
Food.beta <- zoo(Food.beta, order.by=as.Date(names(Food.beta)))
plot(Food.beta)
El valor estático 0.71 estimado antes por regresión oculta una apreciable variabilidad, con un extraño desplome algo después del año 2000.
Referencias
Pérez-Castroviejo, Pedro, and Fernando Tusell. 2007. “Using Redundant and Incomplete Time Series for the Estimation of Cost of Living Indices in Biscaye, 1862–1940.” Review of Income and Wealth 53: 673–91.