Trabajo Fin de Curso
Javier Gómez Dorado
2022-04-19
Instrucciones Generales:
En el fichero SeriesMarzo2022CursoINE.txt tenéis 5 series temporales. En la primera linea del fichero está el nombre de la serie, que podéis cambiar en los programas a y1, y2, y3, y4, y5, Cada serie está en una columna del fichero, por ejemplo, la segunda serie empieza con 11.7203056607791 y 13.6656699978044 y termina con 93.3096271455892.
De las 5, tenéis que analizar las 2 que queráis.
Las series tienen las siguientes características:
Todas son simuladas, mensuales, de longitud 252 meses, es decir, 21 años completos.
Se supone, para que podamos comparar todos las fechas de las intervenciones, que todas empiezan en enero de 1980. De todas formas, yo prefiero que nos refiramos a las intervenciones por el número de observación (del 1 al 252). Por ejemplo, a78, e165 o t128 para AO, LS y TC, respectivamente.
En ninguna hay que tomar logaritmos.
El modelo de todas es un \(ARIMA(p,d,q)\times(P,D,Q)_{12}\), con cada uno de los p, d, q, P, D y Q que vale o 0 o 1. De todas formas, en una de las series funciona mejor un modelo donde uno de esos parámetros vale 2 (a veces pasa eso, que un modelo que no es el del que se ha simulado funciona mejor que el verdadero)
Cada serie tiene exactamente una intervención, que puede ser, o un outlier aditivo, o un cambio de nivel, o un cambio transitorio con parámetro 0.7 y que dura 7 meses, el cual podemos modelizar en R ayudándonos de un vector que valga
c(1, 0.7, 0.7^2, 0.7^3, 0.7^4, 0.7^5, 0.7^6)En el caso de un cambio transitorio, también nos puede ayudar a verlo en nuestra serie el ajustar varios outlier aditivos en puntos sucesivos y ver que el parámetro estimado de cada uno de ellos va decreciendo.
Al ajustar modelos a la serie podéis emplear estos consejos que acabo de dar junto con lo que hemos hecho en el curso con otros ejemplos.
Tenéis que entregar los modelos e intervención de cada una de las 4 series que hayáis escogido , así como los programas R y su salida (abreviada, sin gráficos, solo la parte numérica) con los que hayáis llegado a esa conclusión, que incluya comentarios entre las lineas de R donde se diga lo que vais haciendo y por qué. Tiene que incluir:
La lectura de las series.
Las pruebas de diferenciación regular y/o estacional.
Los intentos de ajuste de modelos.
La búsqueda de intervenciones, usando los residuos de algún(os) modelo(s) que parezca que pueda(n) ser adecuado(s).
Las directrices son las mismas que hemos usado en clase, pasar los
cuatro test del programa TesRes (hay una serie en que solo
se pasan tres test, pero el cuarto no se pasa por poco), que la
correlación máxima entre parámetros estimados no sea demasiado grande
(pongamos que no supere 0.7, por ejemplo) y, si se puede escoger entre
varios modelos que cumplan todo esto, minimizar el AIC.
Análisis de la Serie 1:
Lectura de la serie.
Mediante la siguiente expresión denominamos por y1 a la
serie Serie1Mar2022, la convertimos en un objeto de series
temporales indicando que su periodicidad es anual y que comienza en
enero de 1980.
Series <- read.csv("SeriesMarzo2022CursoINE.txt", header = TRUE, sep = ",",dec =".")
head(Series); tail(Series)## Serie1Mar2022 Serie2Mar2022 Serie3Mar2022 Serie4Mar2022 Serie5Mar2022
## 1 5.624840 11.72031 -1.4547821 -0.22666393 -17.42202
## 2 7.180966 13.66567 0.5404285 -2.34302419 -14.07551
## 3 10.174619 17.49035 2.6619062 0.76358500 -14.70060
## 4 14.748377 12.68971 1.9476337 -0.30329438 -18.67759
## 5 15.978658 21.62860 3.2746635 -1.19379181 -10.70948
## 6 9.665338 19.33614 3.8791944 0.01854653 -10.79782## Serie1Mar2022 Serie2Mar2022 Serie3Mar2022 Serie4Mar2022 Serie5Mar2022
## 247 220.7427 96.48861 11.070914 -52.28748 -39.30166
## 248 220.4529 90.23218 12.052707 -49.90917 -34.48425
## 249 232.3560 90.68834 4.183355 -53.48277 -35.47603
## 250 239.8574 98.07596 3.415131 -54.33832 -26.96899
## 251 241.0920 96.93653 -3.365590 -53.92812 -31.42053
## 252 247.5119 93.30963 -1.993826 -53.27820 -41.02509y1 <- ts(Series$Serie1Mar2022, frequency = 12, start = c(1980, 1))A continuación exponemos en cuatro gráficos la representación de la
serie y1, representamos también el pseudoespectro, para el
cual hemos utilizado un núcleo sencillo. También representamos el
autocrrelograma y el autocorrelagrama parcial.
par(mfrow=c(2,2))
plot(y1)
k2 = kernel("modified.daniell", c(1))
spec.pgram(y1, k2, taper=0,log='yes', detrend=FALSE, demean=TRUE)
acf(y1);pacf(y1)
Pruebas de diferenciación regular y/o estacional.
A la vista de representación de la serie y1, esta parece
homocedástica por lo que tal y como nos dice el enunciado no es
necesario tomar logaritmos. En el gráfico del pseudoespectro podemos
observar picos espectrales en la tendencia y en las frecuencias
estacionales. Tenemos autocrrelaciones bastante altas y con un
decrecimiento lento lo cual suguiere que necesitamos diferenciar la
serie pues no parece que sea estacionaria.
Tomamos una diferenciación regular de la serie original, a la que
llamamos dyy volvemos a generar los cuatro gráficos
anteriores esta vez para la serie diferenciada dy
par(mfrow=c(2,2))
dy<-diff(y1); plot(dy)
spec.pgram(dy, k2, taper=0,log='yes', detrend=FALSE, demean=TRUE, fast=FALSE)
acf(dy);pacf(dy)
En el pseudoespectro se ha reducido el pico espectral de la tendencia aunque se han acentuado los picos de las frecuencias estacionales. El autocorrelograma contiene autocrrelaciones grandes y una estructura sinusoidal que nos sugiere que hemos de hacer una diferenciación estacional.
Realizamos una diferenciación estacional de la serie original a la
que denominamos d12yy volvemos a generar para ella los
cuatro gráficos anteriores.
par(mfrow=c(2,2))
d12y<-diff(y1,lag=12);plot(d12y)
spec.pgram(d12y, k2, taper=0,log='yes', detrend=FALSE, demean=TRUE, fast=FALSE)
acf(d12y);pacf(d12y)
En el pseudoespectro volvemos a tener alto el pico espectral de la tendencia y los valores altos y de decrecimiento lento del autocorrelograma nos sugieren que la serie no es estacionaria.
Tomamos la decisión de diferenciar tanto estacional como de forma
regular la serie original y la denominamos d12dy. Y
volvemos a generar los cuatro gráficos.
par(mfrow=c(2,2))
d12dy<-diff(dy,lag=12);plot(d12dy)
spec.pgram(d12dy, k2, taper=0,log='yes', detrend=FALSE, demean=TRUE, fast=FALSE)
acf(d12dy);pacf(d12dy)
Ahora podemos observar cómo en el pseudoespectro se han suavizado los
picos de la tendencia y de las frecuencias espectrales y cómo no son
muchos las autocorrelaciones y autocrrelaciones parciales que se salen
de las bandas de confianza.
Ajuste de modelos.
En el apartado anterior llegamos a la conclusión de que en el modelo
\(SARIMA(p,d,q)\times(P,D,Q)_{12}\),
los valores para los parámetros \(d =
1\) y \(D = 1\) puesto que
decidimos tomar una diferenciación regular y otra estacional. Buscamos
modelos sencillos para ello construimos la siguiente tabla donde se
muestran las distintas combinaciones asignando a los parámetros \(p, P, q, Q\) los valores \(0\) ó \(1\). Se muestra en la tabla los p-valores
para los disntintos contrastes de la función TesRes.
Una
vez expuesta la tabla vamos a detenernos en el modelo elegido \(SARIMA(1, 1, 0)\times(0, 1, 1)_{12}\). En
primer lugar hemos descartados aquellos modelos que tenían una
correlación alta entre los parámetros estimados
max(abs(Corr)) < -0.7 y de los restantes hemos elegido
aquel que pasaba el mayor número de test y con p-valores más altos.
par(mfrow=c(3,2))
mod110011<-arima(y1,order=c(1,1,0),seasonal=c(0,1,1))
Temod110011<-TesRes(mod110011,2); Temod110011[6:8]## [[1]]
## Resumen VM
## 1 6.523975e-01 pvalLjungBox
## 2 2.835745e-07 pvalMcLeodLi
## 3 3.395805e-15 pvalShapiroWilk
## 4 2.097124e-01 pvalLjungEsta
## 5 -2.120360e-01 max(abs(Corr))
## 6 2.000000e+00 fila max
## 7 1.000000e+00 col max
##
## [[2]]
## Parametros Desv. tipicas Cocientes T p-valores
## ar1 0.4795322 0.05788048 8.284869 1.182392e-16
## sma1 -0.2294485 0.06825357 -3.361706 7.746247e-04
##
## [[3]]
## Var. residual log-vero aic
## 1 2.480535 -448.1441 902.2882plot(mod110011$residuals)
spec.pgram(mod110011$residuals, k2, taper=0,log='yes', detrend=FALSE, demean=TRUE, fast=FALSE)
acf(mod110011$residuals);pacf(mod110011$residuals)
acf(mod110011$residuals^2);pacf(mod110011$residuals^2)
No obstante el modelo elegido no pasa el test de normalidad
pvalShapiroWilk ni el test pvalMcLeodLi, cuya
hipótesis nula dice es que no existe heteroscedasticidad condicional
autorregresiva (ARCH) entre los retardos considerados. Aún así los
coeficientes del modelo son significativos pues presentan p-valores muy
bajos.
Búsqueda de intervenciones.
Analizamos los residuos y comprobamos cuales de ellos son mayores de 3 veces su desviación típica, luego probamos con 4 veces y con más de 5 veces su desviación típica.
ResEstan<-mod110011$residuals/sqrt(mod110011$sigma2)
par(mfrow=c(1,1))
plot(ResEstan)
which(abs(ResEstan)>5)## [1] 132 144which(abs(ResEstan)>4)## [1] 132 144which(abs(ResEstan)>3)## [1] 132 144 145Comprobamos que los residuos que exceden más de 5 veces su desviación típica son los que ocupan las posiciones 132, 144 y 145 respectivamente. El hecho de que las pasociones 144 y 145 sean consecutivas nos hace sospechar que la intervención de variable escalón podría ser adecuada.
Vamos a probar con una variabale escalón para ello definimos una matriz con todos ceros salvo de la posición 144 en adelante que serán todos unos
# A la vista de lo anterior metemos dos intervenciones
z1<-matrix(1,nrow=length(y1),ncol=1)
icual<-144 #pondremos 1 de 1992:1 en adelante
z1[1:icual,1]<-0 #de 1991:12 para atrás ponemos ceroUna vez construida la matriz la pasamos como regresor al modelo.
par(mfrow=c(3,2))
mod110011i1<-arima(y1,order=c(1,1,0),seasonal=c(0,1,1), xreg=z1)
Temod110011i1<-TesRes(mod110011i1,2); Temod110011i1[6:8]## [[1]]
## Resumen VM
## 1 6.335323e-01 pvalLjungBox
## 2 1.063291e-06 pvalMcLeodLi
## 3 4.739607e-15 pvalShapiroWilk
## 4 1.439422e-01 pvalLjungEsta
## 5 -1.820707e-01 max(abs(Corr))
## 6 2.000000e+00 fila max
## 7 1.000000e+00 col max
##
## [[2]]
## Parametros Desv. tipicas Cocientes T p-valores
## ar1 0.4853732 0.05740012 8.455961 2.767974e-17
## sma1 -0.2222557 0.06837665 -3.250461 1.152179e-03
## z1 -1.5820676 1.10584601 -1.430640 1.525334e-01
##
## [[3]]
## Var. residual log-vero aic
## 1 2.459992 -447.1334 902.2669plot(mod110011i1$residuals)
spec.pgram(mod110011i1$residuals, k2, taper=0,log='yes', detrend=FALSE, demean=TRUE, fast=FALSE)
acf(mod110011i1$residuals);pacf(mod110011i1$residuals)
acf(mod110011i1$residuals^2);pacf(mod110011i1$residuals^2)
No parece mejorar mucho el modelo y de hecho el coeficiente que para la intervención no parece significativo.
Vamos a optar por considerar los tres valores como outliers y para ello definimos una matriz con todos cero salvo en las posiciones donde los residuos exceden 5 veces sus varianza donde completaremos con unos.
z2<-matrix(0,nrow=length(y1),ncol=1)
z2[132,1]<-1
z2[144,1]<-1
z2[145,1]<-1par(mfrow=c(3,2))
mod110011i2<-arima(y1,order=c(1,1,0),seasonal=c(0,1,1), xreg=z2)
Temod110011i2<-TesRes(mod110011i2,2); Temod110011i2[6:8]## [[1]]
## Resumen VM
## 1 2.595672e-02 pvalLjungBox
## 2 4.839237e-06 pvalMcLeodLi
## 3 8.153841e-14 pvalShapiroWilk
## 4 3.218501e-01 pvalLjungEsta
## 5 -1.695552e-01 max(abs(Corr))
## 6 2.000000e+00 fila max
## 7 1.000000e+00 col max
##
## [[2]]
## Parametros Desv. tipicas Cocientes T p-valores
## ar1 0.5372060 0.05522276 9.727981 2.290943e-22
## sma1 -0.2342128 0.06323087 -3.704089 2.121519e-04
## z2 -2.6621813 0.51545093 -5.164762 2.407451e-07
##
## [[3]]
## Var. residual log-vero aic
## 1 2.238059 -435.9052 879.8104plot(mod110011i2$residuals)
spec.pgram(mod110011i2$residuals, k2, taper=0,log='yes', detrend=FALSE, demean=TRUE, fast=FALSE)
acf(mod110011i2$residuals);pacf(mod110011i2$residuals)
acf(mod110011i2$residuals^2);pacf(mod110011i2$residuals^2)
En este último modelo sí que parece que mejoran los p-valores de los contrastes y todos los coeficients son significativos por lo que lo elegimos como mejor modelo. Con todo lo expuesto el modelo elegido es:
\[ SARIMA(p, d, q) \times(P, D, Q)_{s} \] donde \(p=1\), \(d=1\), \(q=0\), \(P=0\), \(D=1\), \(Q=1\) y \(s=12\). Cuya representación es:
\[ \phi(B) \Phi\left(B^{s}\right) \Delta^{d} \Delta_{s}^{D} y_{t}=\theta(B) \Theta\left(B^{s}\right) \varepsilon_{t} \] donde \(\phi=-0.5372060\), \(\Phi=0\), \(\theta=0\), \(\Theta= -0.2342128\), se incorpora al modelo además un outlier aditivo cuyo coeficiente en el modelo es \(\alpha=-2.6621813\).
Análisis de la Serie 2:
Lectura de la serie.
Mediante la siguiente expresión denominamos por y2 a la
serie Serie2Mar2022, la convertimos en un objeto de series
temporales indicando que su periodicidad es anual y que comienza en
enero de 1980.
y2 <- ts(Series$Serie2Mar2022, frequency = 12, start = c(1980, 1))A continuación exponemos en cuatro gráficos la representación de la
serie y1, representamos también el pseudoespectro, para el
cual hemos utilizado un núcleo sencillo. También representamos el
autocrrelograma y el autocorrelagrama parcial.
par(mfrow=c(2,2))
plot(y2)
k2 = kernel("modified.daniell", c(1))
spec.pgram(y2, k2, taper=0,log='yes', detrend=FALSE, demean=TRUE)
acf(y2);pacf(y2)
Pruebas de diferenciación regular y/o estacional.
A la vista de representación de la serie y2, esta parece
homocedástica por lo que tal y como nos dice el enunciado no es
necesario tomar logaritmos. En el gráfico del pseudoespectro podemos
observar picos espectrales en la tendencia y en las frecuencias
estacionales. Tenemos autocrrelaciones bastante altas y con un
decrecimiento lento lo cual suguiere que necesitamos diferenciar la
serie pues no parece que sea estacionaria.
Tomamos una diferenciación regular de la serie original, a la que
llamamos dyy volvemos a generar los cuatro gráficos
anteriores esta vez para la serie diferenciada dy
par(mfrow=c(2,2))
dy2<-diff(y2); plot(dy2)
spec.pgram(dy2, k2, taper=0,log='yes', detrend=FALSE, demean=TRUE, fast=FALSE)
acf(dy2);pacf(dy2)
En el pseudoespectro se ha reducido el pico espectral de la tendencia aunque se han acentuado los picos de las frecuencias estacionales. El autocorrelograma contiene autocrrelaciones grandes y una estructura sinusoidal que nos sugiere que hemos de hacer una diferenciación estacional.
Realizamos una diferenciación estacional de la serie original a la
que denominamos d12y2y volvemos a generar para ella los
cuatro gráficos anteriores.
par(mfrow=c(2,2))
d12y2<-diff(y2,lag=12);plot(d12y2)
spec.pgram(d12y2, k2, taper=0,log='yes', detrend=FALSE, demean=TRUE, fast=FALSE)
acf(d12y2);pacf(d12y2)
En el pseudoespectro volvemos a tener alto el pico espectral de la tendencia y los valores altos y de decrecimiento lento del autocorrelograma nos sugieren que la serie no es estacionaria.
Tomamos la decisión de diferenciar tanto estacional como de forma
regular la serie original y la denominamos d12dy2. Y
volvemos a generar los cuatro gráficos.
par(mfrow=c(2,2))
d12dy2<-diff(dy2,lag=12);plot(d12dy2)
spec.pgram(d12dy2, k2, taper=0,log='yes', detrend=FALSE, demean=TRUE, fast=FALSE)
acf(d12dy2);pacf(d12dy2)
Ahora podemos observar cómo en el pseudoespectro se han suavizado los
picos de la tendencia y de las frecuencias espectrales y cómo no son
muchos las autocorrelaciones y autocrrelaciones parciales que se salen
de las bandas de confianza.
Ajuste de modelos.
En el apartado anterior llegamos a la conclusión de que en el modelo
\(SARIMA(p,d,q)\times(P,D,Q)_{12}\),
los valores para los parámetros \(d =
1\) y \(D = 1\) puesto que
decidimos tomar una diferenciación regular y otra estacional. Buscamos
modelos sencillos para ello construimos la siguiente tabla donde se
muestran las distintas combinaciones asignando a los parámetros \(p, P, q, Q\) los valores \(0\) ó \(1\). Se muestra en la tabla los p-valores
para los disntintos contrastes de la función TesRes.
Una
vez expuesta la tabla vamos a detenernos en el modelo elegido \(SARIMA(1, 1, 0)\times(0, 1, 1)_{12}\). En
primer lugar hemos descartados aquellos modelos que tenían una
correlación alta entre los parámetros estimados
max(abs(Corr)) < -0.7 y de los restantes hemos elegido
aquel que pasaba el mayor número de test y con p-valores más altos.
par(mfrow=c(3,2))
mod110011a<-arima(y2,order=c(1,1,0),seasonal=c(0,1,1))
Temod110011a<-TesRes(mod110011a,2); Temod110011a[6:8]## [[1]]
## Resumen VM
## 1 6.076551e-01 pvalLjungBox
## 2 8.545989e-01 pvalMcLeodLi
## 3 6.229303e-16 pvalShapiroWilk
## 4 7.340771e-01 pvalLjungEsta
## 5 1.016174e-01 max(abs(Corr))
## 6 2.000000e+00 fila max
## 7 1.000000e+00 col max
##
## [[2]]
## Parametros Desv. tipicas Cocientes T p-valores
## ar1 -0.3712163 0.06029342 -6.156828 7.421616e-10
## sma1 -0.4421692 0.05986385 -7.386247 1.510314e-13
##
## [[3]]
## Var. residual log-vero aic
## 1 2.069658 -427.4281 860.8561plot(mod110011a$residuals)
spec.pgram(mod110011a$residuals, k2, taper=0,log='yes', detrend=FALSE, demean=TRUE, fast=FALSE)
acf(mod110011a$residuals);pacf(mod110011a$residuals)
acf(mod110011a$residuals^2);pacf(mod110011a$residuals^2)
No obstante el modelo elegido no pasa el test de normalidad
pvalShapiroWilk, aún así los coeficientes del modelo son
significativos pues presentan p-valores muy bajos.
Búsqueda de intervenciones.
Analizamos los residuos y comprobamos cuales de ellos son mayores de 3 veces su desviación típica, luego probamos con 4 veces y con más de 5 veces su desviación típica.
ResEstana<-mod110011a$residuals/sqrt(mod110011a$sigma2)
par(mfrow=c(1,1))
plot(ResEstana)
which(abs(ResEstana)>4)## [1] 167which(abs(ResEstana)>3)## [1] 167 179which(abs(ResEstana)>2)## [1] 104 167 169 170 173 179 191Comprobamos que el residuos que exceden más de 4 veces su desviación típica es el que ocupa la posición 167, vamos a probar si existe un outlier aditivo en esa posición. Para ello definimos una matriz con todo ceros salvo en la posición 167 que pondremos un uno.
z1a<-matrix(0,nrow=length(y2),ncol=1)
z1a[167,1]<-1par(mfrow=c(3,2))
mod110011ai1<-arima(y2,order=c(1,1,0),seasonal=c(0,1,1), xreg=z1a)
Temod110011ai1<-TesRes(mod110011ai1,2); Temod110011ai1[6:8]## [[1]]
## Resumen VM
## 1 5.180981e-01 pvalLjungBox
## 2 1.139428e-05 pvalMcLeodLi
## 3 1.565377e-05 pvalShapiroWilk
## 4 7.115647e-01 pvalLjungEsta
## 5 1.609525e-01 max(abs(Corr))
## 6 3.000000e+00 fila max
## 7 1.000000e+00 col max
##
## [[2]]
## Parametros Desv. tipicas Cocientes T p-valores
## ar1 -0.3983845 0.06017197 -6.620765 3.573450e-11
## sma1 -0.3608190 0.05912539 -6.102608 1.043517e-09
## z1a -8.3586015 0.81448974 -10.262378 1.041109e-24
##
## [[3]]
## Var. residual log-vero aic
## 1 1.442082 -383.7971 775.5942plot(mod110011ai1$residuals)
spec.pgram(mod110011ai1$residuals, k2, taper=0,log='yes', detrend=FALSE, demean=TRUE, fast=FALSE)
acf(mod110011ai1$residuals);pacf(mod110011ai1$residuals)
acf(mod110011ai1$residuals^2);pacf(mod110011ai1$residuals^2)
No parece que hayamos mejorado puesto que ahora tampoco se pasa el
test pvalMcLeodLi, cuya hipótesis nula dice es que no
existe heteroscedasticidad condicional autorregresiva (ARCH) entre los
retardos considerados. Por lo que nos quedaremos con el modelo sin
intervenciones
\[ SARIMA(p, d, q) \times(P, D, Q)_{s} \] donde \(p=1\), \(d=1\), \(q=0\), \(P=0\), \(D=1\), \(Q=1\) y \(s=12\). Cuya representación es:
\[ \phi(B) \Phi\left(B^{s}\right) \Delta^{d} \Delta_{s}^{D} y_{t}=\theta(B) \Theta\left(B^{s}\right) \varepsilon_{t} \] donde \(\phi=0.3712163\), \(\Phi=0\), \(\theta=0\), \(\Theta= -0.4421692\).