Licenciado en Matemáticas (UCM). Doctorado en estadística (UGR).
Encargado de la visualización y análisis de datos covid del Principado de Asturias (2021-2022).
Miembro de la Sociedad Española de Estadística e IO y la Real Sociedad Matemática Española.
Actualmente, investigador y docente en la Facultad de Estadística de la UCM. Divulgando por Twitter e Instagram
Objetivos
Entender el concepto de serie temporal y sus diferencias con la regresión → lo que te equivocaste ayer influye en lo que te equivocarás hoy
Entender conceptos teóricos básicos de procesos estocásticos
Aprender a manejar paquetes estadísticos de R de series temporales → la aplicabilidad de la teoría será tu valor en el futuro
Introducirnos en la metodología Box-Jenkins → los datos deben ser estacionarios
Evaluación
Evaluación continua: 3 entregas individuales a ordenador en clase (20%-25%-35%), y una entrega individual teórica a papel en clase (20%). Asistencia no obligatoria pero se valorará positivamente la participación.
Examen final: la nota ponderará en función de tu evaluación continua.
Más de un 7 -> podrás decidir peso del final entre un 0% y un 100% de la nota (es decir, no será obligatorio el final).
Entre 6 y 7 -> podrás decidir peso del final entre un 35% y un 100%.
Entre 5 y 6 -> podrás decidir peso del final entre un 60% y un 100%.
Entre 3.5 y 5 -> podrás decidir peso del final entre un 80% y un 100%.
Por debajo del 3.5 -> el peso del final será del 100%
Si tienes que hacer el examen final, será obligatorio presentarse y sacar más de un 3 para aprobar.
Planificación entregas
Entrega I (20%): 8 de octubre (120 minutos).
Entrega II (25%): … (120 minutos).
Entrega III (35%): … (120 minutos).
Entrega teórica (20%): … (120 minutos).
Examen final: 14 de enero (10:00-13:30)
Se podrán modificar las fechas por saturación con otras asignaturas siempre y cuando el/la delegado/a lo solicite con más de 7 días de antelación.
Diapositivas: diapositivas en Quarto disponibles y actualizadas en https://javieralvarezliebana.es/docencia/time-series. En el menú de las diapositivas (abajo a la izquierda) tienes una opción para descargarlas en pdf en Tools
airquality del paquete {datasets} (ya instalado por defecto): medidas diarias (153 observaciones) de la calidad del aire en Nueva York, de mayo a septiembre de 1973. Se midieron 6 variables: ozono, radiación solar, viento, temperatura, mes y día.
Durante la carrera es probable que hayas tratado con multitud de datos pero hay uno muy especial que trataremos en esta asignatura de manera diferente: las series temporales.
Vamos a cargar el fichero retiro_temp.csv donde tenemos los datos de temperaturas diarios (AEMET) desde 1980 hasta 2024 de la estación instalada en El Retiro (Madrid).
Código
library(readr) # de tidyverse# en tidyverse, read_ en lugar de read.# tendremos datos en formato tibble en lugar de data.frameretiro <-read_csv(file ="./datos/retiro_temp.csv")retiro
# A tibble: 16,314 × 8
fecha id_station nombre provincia altitud tmed tmin tmax
<date> <dbl> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
1 2000-01-01 3195 MADRID, RETIRO MADRID 667 5.4 0.3 10.4
2 2000-01-02 3195 MADRID, RETIRO MADRID 667 5 0.3 9.6
3 2000-01-03 3195 MADRID, RETIRO MADRID 667 3.5 0.1 6.9
4 2000-01-04 3195 MADRID, RETIRO MADRID 667 4.3 1.4 7.2
5 2000-01-05 3195 MADRID, RETIRO MADRID 667 0.6 -0.4 1.6
6 2000-01-06 3195 MADRID, RETIRO MADRID 667 3.8 -1.1 8.8
7 2000-01-07 3195 MADRID, RETIRO MADRID 667 6.2 0.6 11.7
8 2000-01-08 3195 MADRID, RETIRO MADRID 667 5.4 -0.1 11
9 2000-01-09 3195 MADRID, RETIRO MADRID 667 5.5 3 8
10 2000-01-10 3195 MADRID, RETIRO MADRID 667 4.8 1.8 7.8
# ℹ 16,304 more rows
¿Qué es una serie temporal?
¿Qué analizar de estos datos?
Podemos por ejemplo visualizar un boxplot de las temperaturas medias de cada día durante estos últimos 44 años…
Código
library(tidyverse)ggplot(retiro) +geom_boxplot(aes(y = tmed)) +scale_y_continuous(labels = scales::label_number(suffix ="ºC")) +theme_minimal() +labs(title ="Temperatura desde 1980 hasta 2024",x ="Cuatrimestre", y ="Temperatura media diaria")
¿Qué es una serie temporal?
… la densidad de la temperatura durante todo ese tiempo…
Código
ggplot(retiro) +geom_density(aes(x = tmed)) +scale_x_continuous(labels = scales::label_number(suffix ="ºC")) +theme_minimal() +labs(title ="Temperatura desde 1980 hasta 2024",x ="Temperatura media diaria")
¿Qué es una serie temporal?
… pero también podríamos querer relacionar la temperatura media con el mes (por ejemplo con una regresión)…
Código
ggplot(retiro |>mutate(mes =as_factor(lubridate::month(fecha))) |>summarise(mean_temp =mean(tmed, na.rm =TRUE),.by ="mes")) +geom_col(aes(x = mes, y = mean_temp)) +theme_minimal() +labs(title ="Temperatura media por mes",x ="Mes", y ="ºC (media)")
¿Qué es una serie temporal?
… o analizar cómo la temperatura media va incrementándose en cada década…
Código
ggplot(retiro |>mutate(periodo =if_else(fecha <as_date("1990-01-01"),"1980-1990",if_else(fecha <as_date("2000-01-01"),"1990-2000",if_else(fecha <as_date("2010-01-01"),"2000-2010",if_else(fecha <as_date("2020-01-01"),"2010-2020", "después de 2020")))))) +geom_boxplot(aes(x = periodo, y = tmed)) +theme_minimal() +labs(title ="Temperatura media según periodo",x ="periodo", y ="ºC (media)")
¿Qué es una serie temporal?
En todos ejemplos anteriores hemos analizado una variable continua (temperatura) en función de una variable discreta o de grupo (periodo, década, etc).
¿Pero y si queremos relacionarla con una variable temporal “continua” como es la propia fecha?
¿Qué es una serie temporal?
Código
ggplot(retiro) +geom_line(aes(x = fecha, y = tmed), linewidth =0.3, alpha =0.7) +theme_minimal() +labs(title ="Temperatura media como SERIE TEMPORAL",x ="t (fecha)", y ="ºC (media)")
¿Qué es una serie temporal?
Fíjate bien…¿qué elementos detectas?
Código
ggplot(retiro) +geom_line(aes(x = fecha, y = tmed), linewidth =0.3, alpha =0.7) +theme_minimal() +labs(title ="Temperatura media como SERIE TEMPORAL",x ="t (fecha)", y ="ºC (media)")
¿Qué es una serie temporal?
Tendencia: lo que ajustarías con un modelo clásico (por ejemplo, una regresión lineal) y representa el comportamiento global de la serie, algo así como un nivel base respecto al que la serie oscila.
(en nuestro caso: la temperatura global aumenta con el paso de los años)
Código
ggplot(retiro, aes(x = fecha, y = tmed)) +geom_line(linewidth =0.3, alpha =0.7) +geom_smooth(method ="lm", se =FALSE) +theme_minimal() +labs(title ="Temperatura media como SERIE TEMPORAL",x ="t (fecha)", y ="ºC (media)")
¿Qué es una serie temporal?
Estacionalidad: al margen de esa tendencia general, si hacemos zoom, en muchas series podemos observar un patrón que se repite cada x unidades temporales. En el caso de la temperatura, hay un patrón anual: diciembre hace más frío que en agosto.
Código
ggplot(retiro |>filter(between(fecha, as_date("2020-01-01"), as_date("2023-12-31"))),aes(x = fecha, y = tmed)) +geom_line(linewidth =0.3, alpha =0.7) +geom_smooth(method ="loess") +theme_minimal() +labs(title ="Temperatura media diaria de 2020 a 2023",x ="t (fecha)", y ="ºC (media)")
¿Qué es una serie temporal?
Atípicos: como sucede siempre en estadística será importantísimo analizar y tratar los datos atípicos muy alejados de lo esperado. Por ejemplo, en nuestro caso, Filomena.
Código
ggplot(retiro |>filter(between(fecha, as_date("2020-01-01"), as_date("2023-12-31"))) |>mutate(filomena =between(fecha, as_date("2020-12-25"), as_date("2021-01-22"))),aes(x = fecha, y = tmed)) +geom_line(linewidth =0.3, alpha =0.7) +geom_point(aes(alpha = filomena), color ="#991545") +scale_alpha_manual(values =c(0, 1)) +guides(alpha ="none") +theme_minimal() +labs(title ="Temperatura media diaria de 2020 a 2023",x ="t (fecha)", y ="ºC (media)")
¿Qué es una serie temporal?
Intervenciones: incluso podría suceder que la serie tuviese un corte o salto en su comportamiento. Por ejemplo, imagina que de repente el aparato de medición empieza a medir +25 grados de la temperatura real.
Código
ggplot(retiro |>filter(between(fecha, as_date("2020-01-01"), as_date("2023-12-31"))) |>mutate(tmed =if_else(fecha <="2021-12-31", tmed, tmed +25))) +geom_line(aes(x = fecha, y = tmed), linewidth =0.3, alpha =0.7) +guides(alpha ="none") +theme_minimal() +labs(title ="Temperatura media diaria de 2020 a 2023",subtitle ="Error de +25ºC a partir de 2022",x ="t (fecha)", y ="ºC (media)")
Ejemplos de series
En esta asignatura será fundamental un concepto: estacionariedad. Diremos que una serie es estacionaria si oscila de manera estable con una media y varianza constante.
Distintos objetivos
Análisis descriptivo
Visualización: ¿cómo son los datos? ¿Existe algún ausente o valor atípico?
¿Se puede descomponer la serie en series más sencillas?
Análisis probabilístico:
¿Existe un modelo teórico tal que lo que observamos sea simplemente una muestra dicho modelo probabilístico?
Aunque los datos sean aleatorios, ¿podemos modelizar de manera teórica alguna de sus característica?
Predicción
Conociendo su comportamiento pasado, ¿cuánto valdrá su valor mañana?
¿Cuánto me estoy equivocando? ¿Cómo medir ese error?
Bloques del curso
Bloque I: analisis exploratorio. Descomposición y suavizado
Paso 2: para Mac basta con que hacer click en el archivo .pkg, y abrirlo una vez descargado. Para sistemas Windows, debemos clickar en install R for the first time y después en Download R for Windows. Una vez descargado, abrirlo como cualquier archivo de instalación.
Paso 3: abrir el ejecutable de instalación.
Cuidado
Siempre que tengas que descargar algo de CRAN (ya sea el propio R o un paquete), asegúrate de tener conexión a internet.
Instalación de R Studio
RStudio será el Word que usaremos para escribir (lo que se conoce como un IDE: entorno integrado de desarrollo).
Paso 1: entra la web oficial de RStudio (ahora llamado Posit) y selecciona la descarga gratuita.
Paso 2: selecciona el ejecutable que te aparezca acorde a tu sistema operativo.
Paso 3: tras descargar el ejecutable, hay que abrirlo como otro cualquier otro y dejar que termine la instalación.
Scripts (documentos .R)
Un script será el documento en el que programamos, nuestro archivo .doc (aquí con extensión .R) donde escribiremos las órdenes. Para abrir nuestro primero script, haz click en el menú en File < New File < R Script.
Cuidado
Es importante no abusar de la consola: todo lo que no escribas en un script, cuando cierres, lo habrás perdido.
Cuidado
R es case-sensitive: es sensible a mayúsculas y minúsculas por lo que x y X representa variables distintas.
Ejecutando el primer script
Ahora tenemos una cuarta ventana: la ventana donde escribiremos nuestros códigos. ¿Cómo ejecutarlo?
Escribimos el código a ejecutar.
Guardamos el archivo .R haciendo click en Save current document.
El código no se ejecuta salvo que se lo indiquemos. Tenemos tres opciones de ejecutar un script:
Copiar y pegar en consola.
Seleccionar líneas y Ctrl+Enter
Activar Source on save a la derecha de guardar: no solo guarda sino que ejecuta el código completo.
Sé organizado: proyectos
De la misma manera que en el ordenador solemos trabajar de manera ordenada por carpetas, en RStudio podemos hacer lo mismo para trabajar de manera eficaz creando proyectos.
Un proyecto será una «carpeta» dentro de RStudio, de manera que nuestro directorio raíz automáticamente será la propia carpeta de proyecto (pudiendo pasar de un proyecto a otro con el menu superior derecho).
Podemos crear uno en una carpeta nueva o en una carpeta ya existente.
Buenas prácticas
Tip 1: asignar, evaluar y comparar no es lo mismo. Si te has fijado en R estamos usando <- para asignar valores a variables. Usaremos = para evaluar argumentos en funciones y == para saber si dos elementos son iguales.
x <-1# asignarx =1# evaluarx ==1# comparar
Tip 2: programa como escribes. Al igual que cuando redactas en castellano, acostúmbrate a incorporar espacios y saltos de línea paranoquedarteciego (es una buena práctica y no un requisito porque R no procesa los espacios)
x <-1# óptimox<-1# regux<-1# peor (decídete)
Buenas prácticas
Tip 3: no seas caótico, estandariza nombres, acostúmbrate siempre a hacerlo igual. El único requisito es que debe empezar siempre por una letra (y sin tildes). La forma más recomendable es la conocida como snake_case
Tip 4: facilita la lectura y escritura, pon márgenes. En Tools < Global Options puedes personalizar algunas opciones de RStudio. En Code < Display podemos indicarle en Show margin (no interacciona con el código).
Buenas prácticas
Tip 5: el tabulador es tu mejor amigo. En RStudio tenemos una herramienta maravillosa: si escribes parte del nombre de una variable o función y tabulas, RStudio te autocompleta
Buenas prácticas
Tip 6: ni un paréntesis soltero. Siempre que abras un paréntesis deberás cerrarlo. Para facilitar esta tarea entra en Tools < Global Options < Code < Display y activa la opción Rainbow parentheses
Buenas prácticas
Tip 7: fíjate en el lateral izquierdo. No solo podrás ver la línea de código por la que vas sino que, en caso de estar cometiendo un error de sintaxis, el propio RStudio te avisará.
Tip 8: intenta trabajar siempre por proyectos (para esta clase, crea un script clase2.R en el proyecto que creamos en la anterior clase)
¿Existen variables más allá de los números en la ciencia de datos? Piensa por ejemplo en los datos que podrías guardar de una persona:
La edad o el peso será un número.
edad <-33
Su nombre será una cadena de texto (conocida como string o char).
nombre <-"javi"
A la pregunta «¿estás matriculado en la Facultad?» la respuesta será lo que llamamos una variable lógica (TRUE si está matriculado o FALSE en otro caso).
matriculado <-TRUE
Su fecha de nacimiento será precisamente eso, una fecha, un tipo de variable crucial en esta asignatura
Variables de fecha
Un tipo de datos muy especial: los datos de tipo fecha.
fecha_char <-"2021-04-21"
Parece una simple cadena de texto pero debería representar un instante en el tiempo. ¿Qué debería suceder si sumamos un 1 a una fecha?
fecha_char +1
Error in fecha_char + 1: non-numeric argument to binary operator
Las fechas NO pueden ser texto: debemos convertir la cadena de texto a fecha.
Para trabajar con fechas usaremos el paquete {lubridate}, que deberemos instalar antes de poder usarlo.
install.packages("lubridate")
Variables de fecha
Una vez instalado, de todos los paquetes (libros) que tenemos, le indicaremos que nos cargue ese concretamente.
library(lubridate) # instala si no lo has hecho
Para convertir a tipo fecha usaremos la función as_date() del paquete {lubridate} (por defecto en formato yyyy-mm-dd)
# ¡no es una fecha, es un texto!fecha_char +1
Error in fecha_char + 1: non-numeric argument to binary operator
class(fecha_char)
[1] "character"
fecha <-as_date("2023-03-28")fecha +1
[1] "2023-03-29"
class(fecha)
[1] "Date"
Variables de fecha
En as_date() el formato de fecha por defecto es yyyy-mm-dd así si la cadena de texto no se introduce de manera adecuada…
as_date("28-03-2023")
[1] NA
Para cualquier otro formato debemos especificarlo en el argumento opcional format = ... tal que %d representa días, %m meses, %Y en formato de 4 años y %y en formato de 2 años.
as_date("28-03-2023", format ="%d-%m-%Y")
[1] "2023-03-28"
as_date("28-03-23", format ="%d-%m-%y")
[1] "2023-03-28"
as_date("03-28-2023", format ="%m-%d-%Y")
[1] "2023-03-28"
as_date("28/03/2023", format ="%d/%m/%Y")
[1] "2023-03-28"
Variables de fecha
En dicho paquete tenemos funciones muy útiles para manejar fechas:
Con today() podemos obtener directamente la fecha actual.
today()
[1] "2024-12-09"
Con now() podemos obtener la fecha y hora actual
now()
[1] "2024-12-09 19:22:45 CET"
Con year(), month() o day() podemos extraer el año, mes y día
fecha <-today()year(fecha)
[1] 2024
month(fecha)
[1] 12
Resúmenes de paquetes
Amplia contenido
Tienes un resumen en pdf de los paquetes más importantes en la carpeta correspondiente en el campus
Vectores: concatenar
Cuando trabajamos con datos normalmente tendremos columnas que representan variables: llamaremos vectores a una concatenación de celdas (valores) del mismo tipo (lo que sería una columna de una tabla).
La forma más sencilla es con el comando c() (c de concatenar), y basta con introducir sus elementos entre paréntesis y separados por comas
edades <-c(32, 27, 60, 61)edades
[1] 32 27 60 61
Consejo
Un número individual x <- 1 (o bien x <- c(1)) es en realidad un vector de longitud uno –> todo lo que sepamos hacer con un número podemos hacerlo con un vector de ellos.
💻 Tu turno
Intenta realizar los siguientes ejercicios sin mirar las soluciones
📝 Define el vector x como la concatenación de los 5 primeros números impares. Calcula la longitud del vector
Código
# Dos formasx <-c(1, 3, 5, 7, 9)x <-seq(1, 9, by =2)length(x)
📝 Accede al tercer elemento de x. Accede al último elemento (sin importar la longitud, un código que pueda ejecutarse siempre). Elimina el primer elemento.
Código
x[3]x[length(x)]x[-1]
📝 Obtén los elementos de x mayores que 4. Calcula el vector 1/x y guárdalo en una variable.
Código
x[x >4]z <-1/xz
📝 Crea un vector que represente los nombres de 5 personas, de los cuales uno es desconocido.
📝 Encuentra del vector x de ejercicios anteriores los elementos mayores (estrictos) que 1 Y ADEMÁS menores (estrictos) que 7. Encuentra una forma de averiguar si todos los elementos son o no positivos.
Código
x[x >1& x <7]all(x >0)
📝 Dado el vector x <- c(1, -5, 8, NA, 10, -3, 9), ¿por qué su media no devuelve un número sino lo que se muestra en el código inferior?
x <-c(1, -5, 8, NA, 10, -3, 9)mean(x)
[1] NA
📝 Dado el vector x <- c(1, -5, 8, NA, 10, -3, 9), extrae los elementos que ocupan los lugares 1, 2, 5, 6.
📝 Dado el vector x del ejercicio anterior, ¿cuales tienen un dato ausente? Pista: las funciones is.algo() comprueban si el elemento es tipo algo (tabula)
Código
is.na(x)
📝 Define el vector x como la concatenación de los 4 primeros números pares. Calcula el número de elementos de x menores estrictamente que 5.
Código
x[x <5] sum(x <5)
📝 Calcula el vector 1/x y obtén la versión ordenada (de menor a mayor) de las dos formas posibles
Código
z <-1/xsort(z)z[order(z)]
📝 Encuentra del vector x los elementos mayores (estrictos) que 1 y menores (estrictos) que 6. Encuentra una forma de averiguar si todos los elementos son o no negativos.
Código
x[x >1& x <7]all(x >0)
Primera base de datos
Cuando analizamos datos solemos tener varias variables de cada individuo: necesitamos una «tabla» que las recopile. La opción más inmediata son las matrices: concatenación de variables del mismo tipo e igual longitud.
Imagina que tenemos estaturas y pesos de 4 personas. ¿Cómo crear un dataset con las dos variables?
La opción más habitual es usando cbind(): concatenamos (bind) vectores en forma de columnas (c)
También podemos construir la matriz por filas con la función rbind() (concatenar - bind - por filas - rows), aunque lo recomendable es tener cada variable en columna e individuo en fila como luego veremos.
rbind(estaturas, pesos) # Construimos la matriz por filas
Podemos comprobar las dimensiones con dim(), nrow() y ncol(): las matrices son un tipo de datos tabulados (organizados en filas y columnas)
dim(datos_matriz)
[1] 4 2
nrow(datos_matriz)
[1] 4
ncol(datos_matriz)
[1] 2
Segundo intento: data.frame
Las matrices tienen el mismo problema que los vectores: si juntamos datos de distinto tipo, se perturba la integridad del dato ya que los convierte (fíjate en el código inferior: las edades y los TRUE/FALSE los ha convertido a texto)
edades soltero nombres
[1,] "14" "TRUE" "javi"
[2,] "24" NA "laura"
[3,] NA "FALSE" "lucía"
De hecho al no ser números ya no podemos realizar operaciones aritméticas
matriz +1
Error in matriz + 1: non-numeric argument to binary operator
Segundo intento: data.frame
Para poder trabajar con variables de distinto tipo tenemos en R lo que se conoce como data.frame: concatenación de variables de igual longitud pero que pueden ser de tipo distinto.
tabla <-data.frame(edades, soltero, nombres)class(tabla)
[1] "data.frame"
tabla
edades soltero nombres
1 14 TRUE javi
2 24 NA laura
3 NA FALSE lucía
Segundo intento: data.frame
Dado que un data.frame es ya un intento de «base de datos» las variables no son meros vectores matemáticos: tienen un significado y podemos (debemos) ponerles nombres que describan su significado
edad estado nombre f_nacimiento
1 14 TRUE javi 1989-09-10
2 24 NA laura 1992-04-01
3 NA FALSE lucía 1980-11-27
Intento final: tibble
Las tablas en formato data.frame tienen algunas limitaciones. La principal es que no permite la recursividad: imagina que definimos una base de datos con estaturas y pesos, y queremos una tercera variable con el IMC
📝 Carga del paquete {datasets} el conjunto de datos airquality (variables de la calidad del aire de Nueva York desde mayo hasta septiembre de 1973). ¿Es el conjunto de datos airquality de tipo tibble? En caso negativo, conviértelo a tibble (busca en la documentación del paquete en https://tibble.tidyverse.org/index.html).
📝 Una vez convertido a tibble obtén el nombre de las variables y las dimensiones del conjunto de datos. ¿Cuántas variables hay? ¿Cuántos días se han medido?
📝 Filtra solo los datos del mes de agosto. ¿Cómo indicarle que queremos solo las filas que cumplan una condición concreta? (pista: en realidad todo son vectores “formateados”)
Código
airquality_tb[Month ==8, ]
📝 Selecciona aquellos datos que no sean ni de julio ni de agosto.
Del paquete {Biostatistics} usaremos el conunto de datos pinniped, que guarda los datos de peso de cuerpo y cerebro (desagregado por sexo y mono/poligamia) de 33 especies de mamíferos marinos.
Intenta responder a las preguntas planteadas en el workbook
Comunicar: rmd y Quarto
Una de las principales fortalezas de R es la facilidad para generar informes, libros, webs, apuntes y hasta diapositivas (este mismo material por ejemplo). Para ello instalaremos antes
el paquete {rmarkdown} (para generar archivos .rmd)
install.packages("rmarkdown")
instalar Quarto (si ya conocías R, el «nuevo» .rmd ahora como .qmd)
Comunicar: rmd y Quarto
Hasta ahora solo hemos programado en scripts (archivos .R) dentro de proyectos, pero en muchas ocasiones no trabajaremos solos y necesitaremos comunicar los resultados en diferentes formatos:
Los archivos de extensión .qmd (o .rmd antes) nos permitirán fácilmente combinar:
Markdown: lenguaje tipado que nos permite crear contenido simple (tipo wordpress, con texto, negritas, cursivas, etc) con un diseño legible.
Matemáticas (latex): lenguaje para escribir notación matemática como \(x^2\) o \(\sqrt{y}\) o \(\int_{a}^{b} f(x) dx\)
Código y salidas: podremos no solo mostrar el paso final sino el código que has ido realizando (en R, Python, C++, Julia, …), con cajitas de código llamadas CHUNKS.
Imágenes, gráficas, tablas, estilos (css, js), etc.
Comunicar: rmd y Quarto
La principal ventaja de realizar este tipo de material en Quarto/Rmarkdown es que, al hacerlo desde RStudio, puedes generar un informe o una presentación sin salirte del entorno de programación en el que estás trabajando
De esta forma podrás analizar los datos, resumirlos y a la vez comunicarlos con la misma herramienta.
Recientemente el equipo de RStudio desarrolló Quarto, una versión mejorada de Rmarkdown (archivos .qmd), con un formato un poco más estético y simple. Tienes toda la documentación y ejemplos en https://quarto.org/
Vamos a crear el primer fichero rmarkdown con Quarto con extensión .qmd. Para ello solo necesitaremos hacer click en
File << New File << Quarto Document
Nuestro primer informe
Tras hacerlo nos aparecerán varias opciones de formatos de salida:
archivo .pdf
archivo .html (recomendable): documento dinámico, permite la interacción con el usuario, como una «página web».
archivo .doc (nada recomendable)
De momento dejaremos marcado el formato HTML que viene por defecto, y escribiremos el título de nuestro documento. Tras ello tendremos nuestro archivo .qmd (ya no es un script .R como los que hemos abierto hasta ahora).
Nuestro primer informe
Deberías tener algo similar a la captura de la imagen con dos modos de edición: Source (con código, la opción recomendada hasta que lo domines) y Visual (más parecido a un blog)
Para ejecutar TODO el documento debes clickar Render on Save y darle a guardar.
Salida de Quarto
Deberías haber obtenido una salida en html similar a esta (y se te ha generado en tu ordenador un archivo html)
Editor: source vs visual
Como se indicaba, tienes dos formas de trabajar: con código puro y algo parecido a un Notion (blog)
Un fichero .qmd se divide básicamente en tres partes:
Cabecera: la parte que tienes al inicio entre ---.
Texto: que podremos formatear y mejorar con negritas (escrito como negritas, con doble astérisco al inicio y final), cursivas (cursivas, con barra baja al inicio y final) o destacar nombres de funciones o variables de R. Puedes añadir ecuaciones como \(x^2\) (he escrito $x^2$, entre dólares).
Código R
Cabecera de un qmd
La cabecera están en formato YAML y contiene los metadatos del documento
title y subtitle: el título/subtítulo del documento
author: autor del mismo
format: formato de salida (podremos personalizar)
theme: si tienes algún archivo de estilos
toc: si quieres índice o no
toc-location: posición del índice
toc-title: título del índice
editor: si estás en modo visual o source.
---title:"prueba"format:html:editor: visual---
Cabecera de un qmd
La cabecera están en formato YAML y contiene los metadatos del documento
title y subtitle: el título/subtítulo del documento
Respecto a la escritura solo hay una cosa importante: salvo que indiquemos lo contrario, TODO lo que vamos a escribir es texto (normal). No código R.
Vamos a empezar escribiendo una sección al inicio (# Intro y detrás por ej. la frase
Este material ha sido diseñado por el profesor Javier Álvarez Liébana, docente en la Universidad Complutense de Madrid
Además al Running Code le añadiremos una almohadilla #: las almohadillas FUERA DE CHUNKS nos servirán para crear epígrafes (secciones) en el documento
Índice de un qmd
Para que el índice capture dichas secciones modificaremos la cabecera del archivo como se observa en la imagen (puedes cambiar la localización del índice y el título si quieres para probar).
Texto en un qmd
Vamos a personalizar un poco el texto haciendo lo siguiente:
Vamos a añadir negrita al nombre (poniendo ** al inicio y al final).
Vamos añadir cursiva a la palabra material (poniendo _ al inicio y al final).
Vamos añadir un enlacehttps://www.ucm.es, asociándolo al nombre de la Universidad. Para ello el título lo ponemos entre corchetes y justo detrás el enlace entre paréntesis [«Universidad Complutense de Madrid»](https://www.ucm.es)
Código en un qmd
Para añadir código R debemos crear nuestras cajas de código llamadas chunks: altos en el camino en nuestro texto markdown donde podremos incluir código de casi cualquier lenguaje (y sus salidas).
Para incluir uno deberá de ir encabezado de la siguiente forma tienes un atajo Command + Option + I (Mac) o Ctrl + Shift + I (Windows)
Código en un qmd
Dentro de dicha cajita (que tiene ahora otro color en el documento) escribiremos código R como lo veníamos haciendo hasta ahora en los scripts.
Vamos por ejemplo a definir dos variables y su suma de la siguiente manera, escribiendo dicho código en nuestro .qmd (dentro de ese chunk)
# Código Rx <-1y <-2x + y
[1] 3
Etiquetando chunks
Los chunks pueden tener un nombre o etiqueta, de forma que podamos referenciarlos de nuevo para no repetir código.
Ejecutando chunks
En cada chunk aparecen dos botones:
botón de play: activa la ejecución y salida de ese chunk particular (lo puedes visualizar dentro de tu propio RStudio)
botón de rebobinar: activa la ejecución y salida de todos los chunk hasta ese (sin llegar a él)
Además podemos incluir código R dentro de la línea de texto (en lugar de mostrar el texto x ejecuta el código R mostrando la variable).
Personalización de chunks
Los chunks podemos personalizarlos con opciones al inicio del chunk precedido de #|:
#| echo: false: ejecuta código y se muestra resultado pero no visualiza código en la salida.
#| include: false: ejecuta código pero no muestra resultado y no visualiza código en la salida.
#| eval: false: no ejecuta código, no muestra resultado pero sí visualiza código en la salida.
#| message: false: ejecuta código pero no muestra mensajes de salida.
#| warning: false: ejecuta código pero no muestra mensajes de warning.
#| error: true: ejecuta código y permite que haya errores mostrando el mensaje de error en la salida.
Estas opciones podemos aplicarlas chunk a chunk o fijar los parámetros de forma global con knitr::opts_chunk$set() al inicio del documento (dentro de un chunk).
Personalizando chunks
Si queremos que aplique la opción a todos los chunks por defecto debemos incluirlo al final de la cabecera, como opciones de ejecución
Además de texto y código podemos introducir lo siguiente:
Ecuaciones: puedes añadir además ecuaciones como \(x^2\) (he escrito $x^2$, la ecuación entre dólares).
Listas: puedes itemizar elementos poniendo *
* Paso 1: ...
* Paso 2: ...
Cross-references: puedes etiquetar partes del documento (la etiqueta se construye con {#nombre-seccion}) y llamarlas luego con [Sección](@nombre-seccion)
Gráficas/imágenes en qmd
Por último, también podemos añadir pies de gráficas o imágenes añadiendo #| fig-cap: "..."
Fíjate que el caption está en el margen (por ejemplo). Puedes cambiarlo introduciendo ajustes en la cabecera (todo lo relativo a figuras empieza por fig-, y puedes ver las opciones tabulando). Tienes más información en https://quarto.org/
El archivo de estilos debe estar en la misma carpeta que el archivo .qmd
Añadir estilos
También puedes hacerlo de manera sencilla añadiendo a los textos un poco de HTML. Por ejemplo, para personalizar el color de un texto va entre corchetes y justo tras el texto, entre llaves, las opciones de estilo
Esta palabra es [roja]{style="color:red;"} ...
... y esta [verde y en negrita]{style="color:green; font-weight: bold;"}
Esta palabra es roja …
… y esta verde y en negrita
Revealjs
Puedes añadir algunas «animaciones» usando lo que se conoce como Revealjs (javascript), especifcándolo en la cabecera y usando bloques de dicho lenguaje delimitados por ::: al inicio y final, y la palabra de la «herramienta» a usar. Por ejemplo {.incremental} hace una transición de los elementos.
format: revealjs
::: {.incremental}- Me- llamo- Javi:::
Me
llamo
Javi
Bloques de llamada
También puedes usar los bloques de llamada que por defecto son note, tip, warning, caution e important (aunque los puedes crear y personalizar). Para ello basta con usar :::{.callout-tipo} y el tipo que quieras
:::{.callout-tip}Note that there are five types of callouts, including: `note`, `tip`, `warning`, `caution`, and `important`.:::
Consejo
Recuerda que los 5 tipos son note, tip, warning, caution e important.
Peligro
Úsalos con cabeza, a veces mucho recursos estético puede marear.
# install.packages("reticulate")library(reticulate)install_python("3.9.12") # Instalar python en PC sino lo tienes# Instalar paquetes de Pythonreticulate::py_install("numpy")reticulate::py_install("matplotlib")
Vamos a realizar un pequeño simulacro antes de la entrega usando el dataset starwars del paquete {dplyr}
Ejemplo de entrega
library(dplyr)starwars
# A tibble: 87 × 14
name height mass hair_color skin_color eye_color birth_year sex gender
<chr> <int> <dbl> <chr> <chr> <chr> <dbl> <chr> <chr>
1 Luke Sk… 172 77 blond fair blue 19 male mascu…
2 C-3PO 167 75 <NA> gold yellow 112 none mascu…
3 R2-D2 96 32 <NA> white, bl… red 33 none mascu…
4 Darth V… 202 136 none white yellow 41.9 male mascu…
5 Leia Or… 150 49 brown light brown 19 fema… femin…
6 Owen La… 178 120 brown, gr… light blue 52 male mascu…
7 Beru Wh… 165 75 brown light blue 47 fema… femin…
8 R5-D4 97 32 <NA> white, red red NA none mascu…
9 Biggs D… 183 84 black light brown 24 male mascu…
10 Obi-Wan… 182 77 auburn, w… fair blue-gray 57 male mascu…
# ℹ 77 more rows
# ℹ 5 more variables: homeworld <chr>, species <chr>, films <list>,
# vehicles <list>, starships <list>
En él tenemos diferentes variables de los personajes de Star Wars, con características de su pelo, piel, altura, nombre, etc.
Ejemplo de entrega
Crea un documento .qmd con nombre, título, formato e índice. Cada ejercicio posterior será una subsección del documento. Ejecuta los chunks que consideres y comenta las salidas para responder a cada pregunta
Ejercicio 1. ¿Cuántos personajes hay guardados en la base de datos? ¿Cuántas características se han medido de cada uno?
Ejercicio 2. Extrae en dos variables distintas nombres y edades las variables correspondientes de la tabla. ¿De qué tipo es la variable nombre? ¿Y la variable birth_year?
Ejercicio 3. Obtén el vector de nombres de los personajes ordenados de mayores a jóvenes.
Ejemplo de entrega
Ejercicio 4. Busca ayuda de la función unique(). Úsala para saber que modalidades tiene la variable cualitativa correspondiente al color de ojos. ¿Cuántos distintos hay?
Ejercicio 5. ¿Existe ALGÚN valor ausente en la variable de color ojos?
Ejercicio 6. Calcula la media y desviación típica de las variables de estatura y peso (cuidado con los ausentes). Define un nuevo tibble con esas dos variables e incorpora una tercera variable que se llame “IMC” que calcule el índice de masa corporal. Incorpora con $ $ la fórmula usada para el IMC.
Estructuras de control
Una estructura de control se compone de una serie de comandos orientados a decidir el camino que tu código debe recorrer
Si se cumple la condición A, ¿qué sucede?
¿Y si sucede B?
¿Cómo puedo repetir una misma expresión (dependiendo de una variable)?
Si has programado antes, quizás te sea familiar las conocidas como estructuras condicionales tales como if (blabla) {...} else {...} o buclesfor/while (a evitar siempre que podamos).
Estructura If
Una de las estructuras de control más famosas son las conocidas como estructuras condicionalesif.
SI (IF) un conjunto de condiciones se cumple (TRUE), entonces ejecuta lo que haya dentro de las llaves
Por ejemplo, la estructura if (x == 1) { código A } lo que hará será ejecutar el código A entre llaves pero SOLO SI la condición entre paréntesis es cierta (solo si x es 1). En cualquier otro caso, no hará nada.
Por ejemplo, definamos un vector de edades de 8 personas
edad <-c(14, 17, 24, 56, 31, 20, 87, 73)edad <18
[1] TRUE TRUE FALSE FALSE FALSE FALSE FALSE FALSE
Estructura If
Nuestra estructura condicional hará lo siguiente: si existe algún menor de edad, imprimirá por pantalla un mensaje.
if (any(edad <18)) { print("Existe alguna persona menor de edad")}
[1] "Existe alguna persona menor de edad"
Estructura If
if (any(edad <18)) { print("Existe alguna persona menor de edad")}
En caso de que las condiciones no sean ciertas dentro de if() (FALSE), no sucede nada
if (all(edad >=18)) { print("Todos son mayores de edad")}
No obtenemos ningún mensaje porque la condición all(edad >= 18) no es TRUE, así que no ejecuta nada.
Estructura If-else
La estructura if (condicion) { código A } puede combinarse con un else { código B }: cuando la condición no está verificada, se ejecutará el código alternativo B dentro de else { }, permitiéndonos decidir que sucede cuando se cumple y cuando no.
Por ejemplo, if (x == 1) { código A } else { código B } ejecutará A si x es igual a 1 y B en cualquier otro caso.
if (all(edad >=18)) { print("Todos son mayores de edad")} else {print("Existe alguna persona menor de edad")}
[1] "Existe alguna persona menor de edad"
Estructura If-else
Esta estructura if - else puede ser anidada: imagina que queremos ejecutar un código si todos son menores; si no sucede, pero todos son mayores de 16, hacer otra cosa; en cualquier otra cosa, otra acción.
if (all(edad >=18)) { print("Todos son mayores de edad")} elseif (all(edad >=16)) {print("Hay algún menor de edad pero todos con 16 años o más")} else { print("Hay alguna persona con menos de 16 años") }
[1] "Hay alguna persona con menos de 16 años"
Truco
Puedes colapsar las estructuras haciendo click en la flecha a la izquierda que aparece en tu script.
If-else vectorizado
Esta estructura condicional se puede vectorizar (en una sola línea) con if_else() (del paquete {dplyr}), cuyos argumentos son
la condición a evaluar
lo que sucede cuando se cumple y cuando no
un argumento opcional para cuando la condición a evaluar es NA
Vamos a etiquetar sin son mayores/menores y un “desconocido” cuando no conocemos
En R base existe ifelse(): no deja especificar que hacer con los ausentes pero permite especificar distintos tipos de datos en TRUE y en FALSE.
Bucles
Aunque en la mayoría de ocasiones se pueden reemplazar por otras estructuras más eficientes y legibles, es importante conocer una de las expresiones de control más famosas: los bucles.
for { }: permite repetir el mismo código en un número prefijado y conocido de veces.
while { }: permite repetir el mismo código pero en un número indeterminado de veces (hasta que una condición deje de cumplirse).
Bucles for
Un bucle for es una estructura que permite repetir un conjunto de órdenes un número finito, prefijado y conocido de veces dado un conjunto de índices.
Vamos a definir un vector x <- c(0, -7, 1, 4) y otra variable vacía y. Tras ello definiremos un bucle for con for () { }: dentro de los paréntesis indicaremos un índice y unos valores a recorrer, dentro de las llaves el código a ejecutar en cada iteración (en este caso, rellenar y como x + 1)
x <-c(0, -7, 1, 4)y <-c()
Bucles for
Un bucle for es una estructura que permite repetir un conjunto de órdenes un número finito, prefijado y conocido de veces dado un conjunto de índices.
Vamos a definir un vector x <- c(0, -7, 1, 4) y otra variable vacía y. Tras ello definiremos un bucle for con for () { }: dentro de los paréntesis indicaremos un índice y unos valores a recorrer, dentro de las llaves el código a ejecutar en cada iteración (en este caso, rellenar y como x + 1)
x <-c(0, -7, 1, 4)y <-c()for (i in1:4) {}
Bucles for
Un bucle for es una estructura que permite repetir un conjunto de órdenes un número finito, prefijado y conocido de veces dado un conjunto de índices.
Vamos a definir un vector x <- c(0, -7, 1, 4) y otra variable vacía y. Tras ello definiremos un bucle for con for () { }: dentro de los paréntesis indicaremos un índice y unos valores a recorrer, dentro de las llaves el código a ejecutar en cada iteración (en este caso, rellenar y como x + 1)
x <-c(0, -7, 1, 4)y <-c()for (i in1:4) { y[i] <- x[i] +1}
Bucles for
Fíjate que debido a que R funciona de manera vectorial por defecto, el bucle es lo mismo que hacer x + 1 directamente.
x <-c(0, -7, 1, 4)y <-c()for (i in1:4) { y[i] <- x[i] +1}y
[1] 1 -6 2 5
y2 <- x +1y2
[1] 1 -6 2 5
Bucles for
Otra opción habitual es indicar los índices de manera «automática»: desde el primero 1 hasta el último (que corresponde con la longitud de x length(x))
x <-c(0, -7, 1, 4)y <-c()for (i in1:length(x)) { y[i] <- x[i] +1}y
[1] 1 -6 2 5
Bucles for
Así la estructura general de un bucle for será siempre la siguiente
for (índice in conjunto) { código (dependiente de i)}
SIEMPRE sabemos cuántas iteraciones tenemos (tantas como elementos haya en el conjunto a indexar)
Evitando bucles
Como ya hemos aprendido con el paquete{microbenchmark} podemos chequear como los bucles suelen ser muy ineficientes (de ahí que debamos evitarlos en la mayoría de ocasiones
library(microbenchmark)x <-1:1000microbenchmark(y <- x^2, for (i in1:100) { y[i] <- x[i]^2 },times =500)
Unit: microseconds
expr min lq mean median
y <- x^2 1.681 1.927 2.041144 1.968
for (i in 1:100) { y[i] <- x[i]^2 } 1358.371 1373.172 1447.217426 1383.176
uq max neval
2.091 5.125 500
1401.113 8067.529 500
Bucles for
Podemos ver otro ejemplo de bucle combinando números y textos: definimos un vector de edades y de nombres, e imprimimos el nombre y edad i-ésima.
nombres <-c("Javi", "Sandra", "Carlos", "Marcos", "Marta")edades <-c(33, 27, 18, 43, 29)library(glue)for (i in1:5) { print(glue("{nombres[i]} tiene {edades[i]} años")) }
Javi tiene 33 años
Sandra tiene 27 años
Carlos tiene 18 años
Marcos tiene 43 años
Marta tiene 29 años
Bucles for
Aunque normalmente se suelen indexar con vectors numéricos, los bucles pueden ser indexados sobre cualquier estructura vectorial, da igual de que tipo sea el conjunto
Vamos a combinar las estructuras condicionales y los bucles: usando el conjunto swiss del paquete {datasets}, vamos a asignar NA si los valores de fertilidad son mayores de 80.
for (i in1:nrow(swiss)) {if (swiss$Fertility[i] >80) { swiss$Fertility[i] <-NA }}
Esto es exactamente igual a un if_else() vectorizado
Otra forma de crear un bucle es con la estructura while { }, que nos ejecutará un bucle un número desconocido de veces, hasta que una condición deje de cumplirse (de hecho puede que nunca termine). Por ejemplo, vamos a inializar una variable ciclos <- 1, que incrementaremos en cada paso, y no saldremos del bucle hasta que ciclos > 4.
ciclos <-1while(ciclos <=4) {print(glue("No todavía, vamos por el ciclo {ciclos}")) ciclos <- ciclos +1}
No todavía, vamos por el ciclo 1
No todavía, vamos por el ciclo 2
No todavía, vamos por el ciclo 3
No todavía, vamos por el ciclo 4
Bucles while
Un bucle while será siempre como sigue
while(condición) { código a hacer mientras la condición sea TRUE# normalmente aquí se actualiza alguna variable}
Bucles while
¿Qué sucede cuando la condición nunca es FALSE? Pruébalo tu mismo
while (1>0) {print("Presiona ESC para salir del bucle")}
Cuidado
Un bucle while { } puede ser bastante «peligroso» sino controlamos bien cómo pararlo.
Bucles while
Contamos con dos palabras reservadas para abortar un bucle o forzar su avance:
break: permite abortar un bucle incluso si no se ha llegado a su final
for(i in1:10) {if (i ==3) {break# si i = 3, abortamos bucle }print(i)}
[1] 1
[1] 2
Bucles while
Contamos con dos palabras reservadas para abortar un bucle o forzar su avance:
next: fuerza un bucle a avanzar a la siguiente iteración
for(i in1:5) {if (i ==3) {next# si i = 3, la obvia y continua al siguiente }print(i)}
[1] 1
[1] 2
[1] 4
[1] 5
💻 Tu turno
Intenta realizar los siguientes ejercicios sin mirar las soluciones
La salida es 0 ya que sqrt(9) es igual 3, y dado que no es menor que 2, devuelve el segundo argumento que es 0
📝 ¿Cuál es la salida del siguiente código?
x <-c(1, NA, -1, 9)if_else(sqrt(x) <2, 0, 1)
Código
La salida es el vector c(0, NA, NA, 1) ya que sqrt(1) sí es menor que 2, sqrt(9) no lo es, y tanto en el caso de sqrt(NA) (raíz de ausente) como sqrt(-1) (devuelve NaN, not a number), su raíz cuadrada no puede verificarse si es menor que 2 o no, así que la salida es NA.
📝 Modifica el código inferior para que, cuando no se pueda verificar si la raíz cuadrada de un número es menor que 2, devuelva -1
x <-c(1, NA, -1, 9)if_else(sqrt(x) <2, 0, 1)
Código
x <-c(1, NA, -1, 9)if_else(sqrt(x) <2, 0, 1, missing =-1)
📝 ¿Cuál es son los valores de x e y del código inferior para z <- 1, z <- -1 y z <- -5?
z <--1if (z >0) { x <- z^3 y <--sqrt(z)} elseif (abs(z) <2) { x <- z^4 y <-sqrt(-z)} else { x <- z/2 y <-abs(z)}
Código
En primero caso x =1 e y =-1. En el segundo caso x =1 e y =1. En el tercer caso -1 y 2
📝 ¿Qué sucederá si ejecutamos el código inferior?
z <-"a"if (z >0) { x <- z^3 y <--sqrt(z)} elseif (abs(z) <2) { x <- z^4 y <-sqrt(-z)} else { x <- z/2 y <-abs(z)}
Código
# dará error ya que no es un argumento numéricoError in z^3: non-numeric argument to binary operator
📝 Del paquete {lubridate}, la función hour() nos devuelve la hora de una fecha dada, y la función now() nos devuelve fecha y hora del momento actual. Con ambas funciones haz que se imprima por pantalla (cat()) “buenas noches” solo a partir de las 21 horas.
Código
# Cargamos libreríalibrary(lubridate)# Fecha-hora actualfecha_actual <-now()# Estructura ifif (hour(fecha_actual) >21) {cat("Buenas noches") # print/cat dos formas de imprimir por pantalla}
📝 Modifica el código inferior para que se imprima un mensaje por pantalla si y solo si todos los datos de airquality son con mes distinto a enero
library(datasets)months <- airquality$Monthif (months ==2) {print("No hay datos de enero")}
Código
library(datasets)months <- airquality$Monthif (all(months !=1)) {print("No hay datos de enero")}
📝 Modifica el código inferior para guardar en una variable llamada temp_alta un TRUE si alguno de los registros tiene una temperatura superior a 90 grados Farenheit y FALSE en cualquier otro caso
temp <- airquality$Tempif (temp ==100) {print("Algunos de los registros tienen temperaturas superiores a 90 grados Farenheit")}
📝 Modifica el código inferior para diseñar un bucle for de 5 iteraciones que solo recorra los primeros 5 impares (y en cada paso del bucle los imprima)
for (i in1:5) {print(i)}
Código
for (i inc(1, 3, 5, 7, 9)) {print(i)}
📝 Modifica el código inferior para diseñar un bucle while que empiece con un contador count <- 1 y pare cuando llegue a 6
Intenta responder a las preguntas planteadas en el workbook donde tendrás que diseñar algunos estudios de simulación haciendo uso de bucles y estructuras condicionales
Creando funciones
No solo podemos usar funciones predeterminadas que vienen ya cargadas en paquetes, además podemos crear nuestras propias funciones para automatizar tareas. ¿Cómo crear nuestra propia función? Veamos su esquema básico:
Nombre: por ejemplo name_fun (sin espacios ni caracteres extraños). Al nombre le asignamos la palabra reservadafunction().
Definir argumentos de entrada (dentro de function()).
Cuerpo de la función dentro de { }.
Finalizamos la función con los argumentos de salida con return().
name_fun <-function() {}
Creando funciones
No solo podemos usar funciones predeterminadas que vienen ya cargadas en paquetes, además podemos crear nuestras propias funciones para automatizar tareas. ¿Cómo crear nuestra propia función? Veamos su esquema básico:
Nombre: por ejemplo name_fun (sin espacios ni caracteres extraños). Al nombre le asignamos la palabra reservadafunction().
Definir argumentos de entrada (dentro de function()).
Cuerpo de la función dentro de { }.
Finalizamos la función con los argumentos de salida con return().
name_fun <-function(arg1, arg2, ...) {}
Creando funciones
No solo podemos usar funciones predeterminadas que vienen ya cargadas en paquetes, además podemos crear nuestras propias funciones para automatizar tareas. ¿Cómo crear nuestra propia función? Veamos su esquema básico:
Nombre: por ejemplo name_fun (sin espacios ni caracteres extraños). Al nombre le asignamos la palabra reservadafunction().
Definir argumentos de entrada (dentro de function()).
Cuerpo de la función dentro de { }.
Finalizamos la función con los argumentos de salida con return().
name_fun <-function(arg1, arg2, ...) { código a ejecutar}
Creando funciones
No solo podemos usar funciones predeterminadas que vienen ya cargadas en paquetes, además podemos crear nuestras propias funciones para automatizar tareas. ¿Cómo crear nuestra propia función? Veamos su esquema básico:
Nombre: por ejemplo name_fun (sin espacios ni caracteres extraños). Al nombre le asignamos la palabra reservadafunction().
Definir argumentos de entrada (dentro de function()).
Cuerpo de la función dentro de { }.
Finalizamos la función con los argumentos de salida con return().
name_fun <-function(arg1, arg2, ...) { código a ejecutarreturn(var_salida)}
Creando funciones
arg1, arg2, ...: serán los argumentos de entrada, los argumentos que toma la función para ejecutar el código que tiene dentro
código: líneas de código que queramos que ejecute la función.
return(var_salida): se introducirán los argumentos de salida.
Todas las variables que definamos dentro de la función son variables LOCALES: solo existirán dentro de la función salvo que especifiquemos lo contrario.
Creando funciones
Veamos un ejemplo muy simple de función para calcular el área de un rectángulo.
Dado que el área de un rectángulo se calcula como el producto de sus lados, necesitaremos precisamente eso, sus lados: esos serán los argumentos de entrada y el valor a devolver será justo su área (\(lado_1 * lado_2\)).
# Definición del nombre de función y argumentos de entradacalcular_area <-function(lado_1, lado_2) {}
Creando funciones
Veamos un ejemplo muy simple de función para calcular el área de un rectángulo.
Dado que el área de un rectángulo se calcula como el producto de sus lados, necesitaremos precisamente eso, sus lados: esos serán los argumentos de entrada y el valor a devolver será justo su área (\(lado_1 * lado_2\)).
# Definición del nombre de función y argumentos de entradacalcular_area <-function(lado_1, lado_2) { area <- lado_1 * lado_2}
Creando funciones
Veamos un ejemplo muy simple de función para calcular el área de un rectángulo.
Dado que el área de un rectángulo se calcula como el producto de sus lados, necesitaremos precisamente eso, sus lados: esos serán los argumentos de entrada y el valor a devolver será justo su área (\(lado_1 * lado_2\)).
# Definición del nombre de función y argumentos de entradacalcular_area <-function(lado_1, lado_2) { area <- lado_1 * lado_2return(area)}
Uso de funciones
También podemos hacer una definición directa de las variables sin almacenar por el camino.
# Definición del nombre de función y argumentos de entradacalcular_area <-function(lado_1, lado_2) {return(lado_1 * lado_2)}
¿Cómo aplicar la función?
calcular_area(5, 3) # área de un rectángulo 5 x 3
[1] 15
calcular_area(1, 5) # área de un rectángulo 1 x 5
[1] 5
Uso de funciones
Consejo
Aunque no sea necesario, es recomendable hacer explícita la llamada de los argumentos, especificando en el código qué valor es para cada argumento para que no dependa de su orden, haciendo el código más legible
calcular_area(lado_1 =5, lado_2 =3) # área de un rectángulo 5 x 3
[1] 15
calcular_area(lado_2 =3, lado_1 =5) # área de un rectángulo 5 x 3
[1] 15
Argumentos por defecto
Imagina ahora que nos damos cuenta que el 90% de las veces usamos dicha función para calcular por defecto el área de un cuadrado (es decir, solo necesitamos un lado). Para ello, podemos definir argumentos por defecto en la función: tomarán dicho valor salvo que le asignemos otro.
¿Por qué no asignar lado_2 = lado_1por defecto, para ahorrar líneas de código y tiempo?
calcular_area <-function(lado_1, lado_2 = lado_1) {# Cuerpo de la función area <- lado_1 * lado_2# Resultado que devolvemosreturn(area)}
Argumentos por defecto
calcular_area <-function(lado_1, lado_2 = lado_1) {# Cuerpo de la función area <- lado_1 * lado_2# Resultado que devolvemosreturn(area)}
Ahora por defecto el segundo lado será igual al primero (si se lo añadimos usará ambos).
calcular_area(lado_1 =5) # cuadrado
[1] 25
calcular_area(lado_1 =5, lado_2 =7) # rectángulo
[1] 35
Salida múltiple
Compliquemos un poco la función y añadamos en la salida los valores de cada lado, etiquetados como lado_1 y lado_2, empaquetando la salida en una vector.
# Definición del nombre de función y argumentos de entradacalcular_area <-function(lado_1, lado_2 = lado_1) {# Cuerpo de la función area <- lado_1 * lado_2# Resultadoreturn(c("area"= area, "lado_1"= lado_1, "lado_2"= lado_2))}
Salida múltiple
Podemos complicar un poco más la salida añadiendo una cuarta variable que nos diga, en función de los argumentos, si rectángulo o cuadrado, teniendo que añadir en la salida una variable que de tipo caracter (o lógica).
# Definición del nombre de función y argumentos de entradacalcular_area <-function(lado_1, lado_2 = lado_1) {# Cuerpo de la función area <- lado_1 * lado_2# Resultadoreturn(c("area"= area, "lado_1"= lado_1, "lado_2"= lado_2,"tipo"=if_else(lado_1 == lado_2, "cuadrado", "rectángulo")))}calcular_area(5, 3)
area lado_1 lado_2 tipo
"15" "5" "3" "rectángulo"
Problema: al intentar juntar números y texto, lo convierte todo a números. Podríamos guardarlo todo en un tibble() como hemos aprendido o en un objeto conocido en R como listas
Orden de los argumentos
Antes nos daba igual el orden de los argumentos pero ahora el orden de los argumentos de entrada importa, ya que en la salida incluimos lado_1 y lado_2.
Recomendación
Como se comentaba, altamente recomendable hacer la llamada a la función indicando explícitamente los argumentos para mejorar legibilidad e interpretabilidad.
# Equivalente a calcular_area(5, 3)calcular_area(lado_1 =5, lado_2 =3)
area lado_1 lado_2 tipo
"15" "5" "3" "rectángulo"
Variables locales vs globales
Un aspecto importante sobre el que reflexionar con las funciones: ¿qué sucede si nombramos a una variable dentro de una función a la que se nos ha olvidado asignar un valor dentro de la misma?
Debemos ser cautos al usar funciones en R, ya que debido a la «regla lexicográfica», si una variable no se define dentro de la función, Rbuscará dicha variable en el entorno de variables.
x <-1funcion_ejemplo <-function() {print(x) # No devuelve nada, solo realiza la acción }funcion_ejemplo()
[1] 1
Variables locales vs globales
Si una variable ya está definida fuera de la función (entorno global), y además es usada dentro de cambiando su valor, el valor solo cambia dentro pero no en el entorno global.
x <-1funcion_ejemplo <-function() { x <-2print(x) # lo que vale dentro}
# lo que vale dentrofuncion_ejemplo() #<<
[1] 2
# lo que vale fueraprint(x) #<<
[1] 1
Variables locales vs globales
Si queremos que además de cambiar localmente lo haga globalmente deberemos usar la doble asignación (<<-).
x <-1y <-2funcion_ejemplo <-function() {# no cambia globalmente, solo localmente x <-3# cambia globalmente y <<-0#<<print(x)print(y)}funcion_ejemplo() # lo que vale dentro
[1] 3
[1] 0
x # lo que vale fuera
[1] 1
y # lo que vale fuera
[1] 0
💻 Tu turno
Intenta realizar los siguientes ejercicios sin mirar las soluciones
📝 Modifica el código inferior para definir una función llamada funcion_suma, de forma que dados dos elementos, devuelve su suma.
nombre <-function(x, y) { suma <-# código a ejecutarreturn()}# Aplicamos la funciónsuma(3, 7)
Código
funcion_suma <-function(x, y) { suma <- x + yreturn(suma)}funcion_suma(3, 7)
📝 Modifica el código inferior para definir una función llamada funcion_producto, de forma que dados dos elementos, devuelve su producto, pero que por defecto calcule el cuadrado
nombre <-function(x, y) { producto <-# código de la multiplicaciónreturn()}producto(3)producto(3, -7)
Código
funcion_producto <-function(x, y = x) { producto <- x * yreturn(producto)}funcion_producto(3)funcion_producto(3, -7)
📝 Define una función llamada igualdad_nombres que, dados dos nombres, nos diga si son iguales o no. Hazlo considerando importantes las mayúsculas, y sin que importen las mayúsculas. Usa el paquete {stringr}.
📝 Crea una función llamada calculo_IMC que, dados dos argumentos (peso y estatura en metros) y un nombre, devuelva una lista con el IMC (\(peso/(estatura_m^2)\)) y el nombre.
📝 Repite el ejercicio anterior pero con otro argumento opcional que se llame unidades (por defecto, unidades = "metros"). Desarrolla la función de forma que haga lo correcto si unidades = "metros" y si unidades = "centímetros".
📝 Crea un tibble ficticio de 7 personas, con tres variables (inventa nombre, y simula peso, estatura en centímetros), y aplica la función definida de forma que obtengamos una cuarta columna con su IMC.
Código
datos <-tibble("nombres"=c("javi", "sandra", "laura","ana", "carlos", "leo", NA),"peso"=rnorm(n =7, mean =70, sd =1),"estatura"=rnorm(n =7, mean =168, sd =5))datos |>mutate(IMC =calculo_IMC(nombres, peso, estatura, unidades ="centímetros")$IMC)
📝 Crea una función llamada atajo que tenga dos argumentos numéricos x e y. Si ambos son iguales, debes devolver "iguales" y hacer que la función acaba automáticamente (piensa cuándo una función sale). OJO: x e y podrían ser vectores. Si son distintos (de igual de longitud) calcula la proporción de elementos diferentes. Si son distintos (por ser distinta longitud), devuelve los elementos que no sean comunes.
Hasta ahora todo lo que hemos repasado en R lo hemos realizado en el paradigma de programación conocido como R base. Y es que cuando R nació como lenguaje, muchos de los que programaban en él imitaron formas y metodologías heredadas de otros lenguajes, basado en el uso de
Bucles for
Bucles while
Estructuras if-else
Y aunque conocer dichas estructuras puede sernos en algunos casos interesantes, en la mayoría de ocasiones han quedado caducas y vamos a poder evitarlas (en especial los bucles) ya que R está especialmente diseñado para trabajar de manera funcional (en lugar de elemento a elemento).
¿Qué es tidyverse?
En ese contexto de programación funcional, hace una década nacía {tidyverse}, un «universo» de paquetes para garantizar un flujo de trabajo eficiente, coherente y lexicográficamente sencillo de entender, basado en la idea de que nuestros datos están limpios y ordenados (tidy)
¿Qué es tidyverse?
{lubridate} manejo de fechas
{rvest}: web scraping
{tidymodels}: modelización/predicción
{tibble}: optimizando data.frame
{tidyr}: limpieza de datos
{readr}: carga datos rectangulares (.csv), {readxl} para importar archivos .xls y .xlsx
{dplyr}: gramática para depurar
{stringr}: manejo de textos
{purrr}: manejo de listas
{forcats}: manejo de cualitativas
{ggplot2}: visualización de datos
¿Qué es tidyverse?
{lubridate} manejo de fechas
{rvest}: web scraping
{tidymodels}: modelización/predicción
{tibble}: optimizando data.frame
{tidyr}: limpieza de datos
{readr}: carga datos rectangulares (.csv), {readxl} para importar archivos .xls y .xlsx
{dplyr}: gramática para depurar
{stringr}: manejo de textos
{purrr}: manejo de listas
{forcats}: manejo de cualitativas
{ggplot2}: visualización de datos
Filosofía base: tidy data
Tidy datasets are all alike, but every messy dataset is messy in its own way (Hadley Wickham, Chief Scientist en RStudio)
TIDYVERSE
El universo de paquetes {tidyverse} se basa en la idea introducida por Hadley Wickham (el Dios al que rezamos) de estandarizar el formato de los datos para
sistematizar la depuración
hacer más sencillo su manipulación.
código legible
Reglas del tidy data
Lo primero por tanto será entender qué son los conjuntos tidydata ya que todo {tidyverse} se basa en que los datos están estandarizados.
Cada variable en una única columna
Cada individuo en una fila diferente
Cada celda con un único valor
Cada dataset en un tibble
Si queremos cruzar múltiples tablas debemos tener una columna común
Tubería (pipe)
En {tidyverse} será clave el operador pipe (tubería) definido como |> (ctrl+shift+M): será una tubería que recorre los datos y los transforma.
En R base, si queremos aplicar tres funciones first(), second() y third() en orden, sería
third(second(first(datos)))
En {tidyverse} podremos leer de izquierda a derecha y separar los datos de las acciones
datos |>first() |>second() |>third()
Apunte importante
Desde la versión 4.1.0 de R disponemos de |>, un pipe nativo disponible fuera de tidyverse, sustituyendo al antiguo pipe%>% que dependía del paquete {magrittr} (bastante problemático).
Tubería (pipe)
La principal ventaja es que el código sea muy legible (casi literal) pudiendo hacer grandes operaciones con los datos con apenas código.
¿Pero qué aspecto tienen los datos no tidy? Vamos a cargar la tabla table4a del paquete {tidyr} (ya lo tenemos cargado del entorno tidyverse).
library(tidyr)table4a
# A tibble: 3 × 3
country `1999` `2000`
<chr> <dbl> <dbl>
1 Afghanistan 745 2666
2 Brazil 37737 80488
3 China 212258 213766
¿Qué puede estar fallando?
Pivotar: pivot_longer()
table4a
# A tibble: 3 × 3
country `1999` `2000`
<chr> <dbl> <dbl>
1 Afghanistan 745 2666
2 Brazil 37737 80488
3 China 212258 213766
❎ Cada fila representa dos observaciones (1999 y 2000) → las columnas 1999 y 2000 en realidad deberían ser en sí valores de una variable y no nombres de columnas.
Incluiremos una nueva columna que nos guarde el año y otra que guarde el valor de la variable de interés en cada uno de esos años. Y lo haremos con la función pivot_longer(): pivotaremos la tabla a formato long:
# A tibble: 6 × 3
country year cases
<chr> <chr> <dbl>
1 Afghanistan 1999 745
2 Afghanistan 2000 2666
3 Brazil 1999 37737
4 Brazil 2000 80488
5 China 1999 212258
6 China 2000 213766
cols: nombre de las variables a pivotar
names_to: nombre de la nueva variable a la quemandamos la cabecera de la tabla (los nombres).
values_to: nombre de la nueva variable a la que vamos a mandar los datos.
Datos SUCIOS: messy data
Veamos otro ejemplo con la tabla table2
table2
# A tibble: 12 × 4
country year type count
<chr> <dbl> <chr> <dbl>
1 Afghanistan 1999 cases 745
2 Afghanistan 1999 population 19987071
3 Afghanistan 2000 cases 2666
4 Afghanistan 2000 population 20595360
5 Brazil 1999 cases 37737
6 Brazil 1999 population 172006362
7 Brazil 2000 cases 80488
8 Brazil 2000 population 174504898
9 China 1999 cases 212258
10 China 1999 population 1272915272
11 China 2000 cases 213766
12 China 2000 population 1280428583
¿Qué puede estar fallando?
Pivotar: pivot_wider()
# A tibble: 12 × 4
country year type count
<chr> <dbl> <chr> <dbl>
1 Afghanistan 1999 cases 745
2 Afghanistan 1999 population 19987071
3 Afghanistan 2000 cases 2666
4 Afghanistan 2000 population 20595360
5 Brazil 1999 cases 37737
6 Brazil 1999 population 172006362
7 Brazil 2000 cases 80488
8 Brazil 2000 population 174504898
9 China 1999 cases 212258
10 China 1999 population 1272915272
11 China 2000 cases 213766
12 China 2000 population 1280428583
❎ Cada observación está dividido en dos filas → los registros con el mismo año deberían ser el mismo
Lo que haremos será lo opuesto: con pivot_wider()ensancharemos la tabla
# A tibble: 6 × 4
country year cases population
<chr> <dbl> <dbl> <dbl>
1 Afghanistan 1999 745 19987071
2 Afghanistan 2000 2666 20595360
3 Brazil 1999 37737 172006362
4 Brazil 2000 80488 174504898
5 China 1999 212258 1272915272
6 China 2000 213766 1280428583
Datos SUCIOS: messy data
Veamos otro ejemplo con la tabla table3
table3
# A tibble: 6 × 3
country year rate
<chr> <dbl> <chr>
1 Afghanistan 1999 745/19987071
2 Afghanistan 2000 2666/20595360
3 Brazil 1999 37737/172006362
4 Brazil 2000 80488/174504898
5 China 1999 212258/1272915272
6 China 2000 213766/1280428583
¿Qué puede estar fallando?
Separar: separate()
table3
# A tibble: 6 × 3
country year rate
<chr> <dbl> <chr>
1 Afghanistan 1999 745/19987071
2 Afghanistan 2000 2666/20595360
3 Brazil 1999 37737/172006362
4 Brazil 2000 80488/174504898
5 China 1999 212258/1272915272
6 China 2000 213766/1280428583
❎ Cada celda contiene varios valores
Lo que haremos será hacer uso de la función separate() para mandar separar cada valor a una columna diferente.
table3 |>separate(rate, into =c("cases", "pop"))
# A tibble: 6 × 4
country year cases pop
<chr> <dbl> <chr> <chr>
1 Afghanistan 1999 745 19987071
2 Afghanistan 2000 2666 20595360
3 Brazil 1999 37737 172006362
4 Brazil 2000 80488 174504898
5 China 1999 212258 1272915272
6 China 2000 213766 1280428583
Separar: separate()
table3 |>separate(rate, into =c("cases", "pop"))
# A tibble: 6 × 4
country year cases pop
<chr> <dbl> <chr> <chr>
1 Afghanistan 1999 745 19987071
2 Afghanistan 2000 2666 20595360
3 Brazil 1999 37737 172006362
4 Brazil 2000 80488 174504898
5 China 1999 212258 1272915272
6 China 2000 213766 1280428583
Fíjate que los datos, aunque los ha separado, los ha mantenido como texto cuando en realidad deberían ser variables numéricas. Para ello podemos añadir el argumento opcional convert = TRUE
table3 |>separate(rate, into =c("cases", "pop"), convert =TRUE)
# A tibble: 6 × 4
country year cases pop
<chr> <dbl> <int> <int>
1 Afghanistan 1999 745 19987071
2 Afghanistan 2000 2666 20595360
3 Brazil 1999 37737 172006362
4 Brazil 2000 80488 174504898
5 China 1999 212258 1272915272
6 China 2000 213766 1280428583
Datos SUCIOS: messy data
Veamos el último ejemplo con la tabla table5
table5
# A tibble: 6 × 4
country century year rate
<chr> <chr> <chr> <chr>
1 Afghanistan 19 99 745/19987071
2 Afghanistan 20 00 2666/20595360
3 Brazil 19 99 37737/172006362
4 Brazil 20 00 80488/174504898
5 China 19 99 212258/1272915272
6 China 20 00 213766/1280428583
¿Qué puede estar fallando?
Unir unite()
table5
# A tibble: 6 × 4
country century year rate
<chr> <chr> <chr> <chr>
1 Afghanistan 19 99 745/19987071
2 Afghanistan 20 00 2666/20595360
3 Brazil 19 99 37737/172006362
4 Brazil 20 00 80488/174504898
5 China 19 99 212258/1272915272
6 China 20 00 213766/1280428583
❎ Tenemos mismos valores divididos en dos columnas
Usaremos unite() para unir los valores de siglo y año en una misma columna
table5 |>unite(col = year_completo, century, year, sep ="")
# A tibble: 6 × 3
country year_completo rate
<chr> <chr> <chr>
1 Afghanistan 1999 745/19987071
2 Afghanistan 2000 2666/20595360
3 Brazil 1999 37737/172006362
4 Brazil 2000 80488/174504898
5 China 1999 212258/1272915272
6 China 2000 213766/1280428583
Ejemplo: relig_income
Vamos a realizar un ejemplo juntos con la tabla relig_income del paquete {tidyr}. Como se indica en la ayuda ? relig_income, la tabla representa la cantidad de personas que hay en cada tramo de ingresos anuales (20k = 20 000$) y en cada religión.
No lo es ya que en realidad solo deberíamos tener una variable de ingresos y la tenemos dividida en 11: todas ellas es la misma variable solo que adopta un valor diferente. ¿Cómo convertirla a tidy data?
Ejemplo: relig_income
La idea es pivotar todas las columnas de ingresos para que acaben en una sola columna llamada income, y los valores (el número de personas) en otra llamada people (por ejemplo). La tabla la haremos más larga y menos ancha así que…
Vamos a hilar más fino: ahora mismo en la variable income en realidad tenemos dos valores, el límite inferior y el superior de la renta. Vamos a separar dicha variable e ingresos en dos, llamadas income_inf y income_sup
Vamos a hilar más fino: ahora mismo en la variable income en realidad tenemos dos valores, el límite inferior y el superior de la renta. Vamos a separar dicha variable e ingresos en dos, llamadas income_inf y income_sup
relig_tidy |># Separamos por -separate(income, into =c("income_inf", "income_sup"), sep ="-")
Piensa ahora como podemos convertir los límites de ingresos a numéricas (eliminando símbolos, letras, etc)
Ejemplo: relig_income
Para ello usaremos el paquete {stringr}, en concreto la función str_remove_all(), a la que le podemos pasar los caracteres que queremos eliminar (fíjate que $ al ser un caracter reservado en R hay que indicárselo con \\$)
Fíjate que tenemos "Don't now/refused". ¿Qué deberíamos tener?
Debería ser un dato ausente así que usaremos if_else(): si contiene dicha frase, NA, en caso contrario su valor (consejo: str_detect() para detectar patrones en textos, y evitar tener que escribir toda la palabra sin errores)
¿Se te ocurre alguna forma de «cuantificar numéricamente» los valores ausentes que tenemos en este caso?
Si te fijas en realidad cuando hay ausente en el límite inferior en realidad podríamos poner un 0 (nadie puede ganar menos de eso) y cuando lo tenemos en el límite superior sería Inf
¿Por qué era importante tenerlo en tidydata? Lo veremos más adelante al visualizar los datos pero esto ya nos permite realizar filtros muy rápidos con muy poco código.
Por ejemplo: ¿cuántas personas agnósticas con ingresos superiores (o iguales) a 30 tenemos?
# una línea de códigosum(relig_tidy$people[relig_tidy$religion =="Agnostic"& relig_tidy$income_inf >=30])
[1] 609
💻 Tu turno
Intenta realizar los siguientes ejercicios sin mirar las soluciones
📝 Usa el dataset original relig_income y trata de responder a la última pregunta: ¿cuántas personas agnósticas con ingresos superiores (o iguales) a 30 tenemos? Compara el código a realizar cuando tenemos tidydata a cuando no. ¿Cuál es más legible si no supieses R? ¿Cuál tiene mayor probabilidad de error?
📝 Usando relig_tidy determina quién tiene más ingresos medios, ¿católicos (Catholic) o agnósticos (Agnostic)? Crea antes una variable avg_income (ingresos medios por intervalo): si hay 5 personas entre 20 y 30, y 3 personas entre 30 y 50, la media sería \((25*5 + 40*3)/8\) (si es Inf por arriba, NA)
📝 Echa un vistazo a la tabla table4b del paquete {tidyr}. ¿Es tidydata? En caso negativo, ¿qué falla? ¿Cómo convertirla a tidy data en caso de que no lo sea ya?
📝 Echa un vistazo a la tabla billboard del paquete {tidyr}. ¿Es tidydata? En caso negativo, ¿qué falla? ¿Cómo convertirla a tidy data en caso de que no lo sea ya?
Código
billboard |>pivot_longer(cols ="wk1":"wk76",names_to ="week",names_prefix ="wk",values_to ="position",values_drop_na =TRUE)
🐣 Caso práctico III
En el paquete {tidyr} contamos con el dataset who2 (dataset de la Organización Mundial de la Salud). Intenta responder a las preguntas planteadas en el workbook.
who2
# A tibble: 7,240 × 58
country year sp_m_014 sp_m_1524 sp_m_2534 sp_m_3544 sp_m_4554 sp_m_5564
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 Afghanistan 1980 NA NA NA NA NA NA
2 Afghanistan 1981 NA NA NA NA NA NA
3 Afghanistan 1982 NA NA NA NA NA NA
4 Afghanistan 1983 NA NA NA NA NA NA
5 Afghanistan 1984 NA NA NA NA NA NA
6 Afghanistan 1985 NA NA NA NA NA NA
7 Afghanistan 1986 NA NA NA NA NA NA
8 Afghanistan 1987 NA NA NA NA NA NA
9 Afghanistan 1988 NA NA NA NA NA NA
10 Afghanistan 1989 NA NA NA NA NA NA
# ℹ 7,230 more rows
# ℹ 50 more variables: sp_m_65 <dbl>, sp_f_014 <dbl>, sp_f_1524 <dbl>,
# sp_f_2534 <dbl>, sp_f_3544 <dbl>, sp_f_4554 <dbl>, sp_f_5564 <dbl>,
# sp_f_65 <dbl>, sn_m_014 <dbl>, sn_m_1524 <dbl>, sn_m_2534 <dbl>,
# sn_m_3544 <dbl>, sn_m_4554 <dbl>, sn_m_5564 <dbl>, sn_m_65 <dbl>,
# sn_f_014 <dbl>, sn_f_1524 <dbl>, sn_f_2534 <dbl>, sn_f_3544 <dbl>,
# sn_f_4554 <dbl>, sn_f_5564 <dbl>, sn_f_65 <dbl>, ep_m_014 <dbl>, …
Clase 3: introducción a series
Introducción al análisis descriptivo de series temporales
Como ya hemos visto, una serie temporal se puede definir de manera informal como una muestra de una variable (usualmente continua) recogida de manera secuencial en el tiempo
Código
ggplot(retiro) +geom_line(aes(x = fecha, y = tmed), linewidth =0.3, alpha =0.7) +theme_minimal() +labs(title ="Temperatura media como SERIE TEMPORAL",x ="t (fecha)", y ="ºC (media)")
Un poco de historia
Las primeras series temporales aparecieron en el siglo XIX, cuando el matemático Laplace (matemático, profesor de Napoleón y luego ministro) se dedicaba a estudiar la relación entre las fases de La Luna y la presión atmosférica (por el nivel del mar).
Aunque su primera aproximación fue un poco chapuza (intentó ajustar una función seno sin tener en cuenta el tiempo), Arthur Schuster decidió aplicar los trabajos de Fourier para que dichas sinusoidales dependiesen del tiempo (Fourier había demostrado que toda función periódica podía descomponerse como suma de senos y cosenos).
Un poco de historia
Treinta años más tarde, Yule y Slutsky aplicaron las ideas de la regresión desarrollados por Galton y Pearson al estudio de procesos cuya variable regresora es ella misma en otro instante temporal (procesos autoregresivos), aunque no fue hasta la llegada de Kolmogorov cuando se formalizó su definición matemática en el contexto de los procesos estocásticos.
Tras acabar la Segunda Guerra Mundial quedaron desclasificados algunos trabajos de Wiener, Kolmogorov, Bartlett y Tukey sobre la predicción de series temporales, así como su estudio en función del análisis de las correlaciones. Tras los trabajos de alisado de Holt y Winter en los años 60, en 1970 Box y Jenkins publican «La Biblia» de las series temporales, un manual donde se presenta una metodología para la identificación, estimación y predicción de series temporales (los conocidos como procesos SARIMA)
Métodos de análisis
Métodos descriptivos o clásicos: desarrollados entre 1940 y 1970, están enfocadas principalmente en la estimación de los valores de la serie, siendo bastante malos en la predicción futura.
Métodos de decomposición (tendencia-estacionalidad-error):
Suavizado clásico
STL (seasonal-trend decomposition procedure based on loess)
Métodos de alisado (dependencia del pasado disminuye con el tiempo):
Alisado exponencial simple
Alisado doble de Holt y triple de Holt-Winters . . .
Métodos probabilísticos: enmarcados dentro del análisis de los procesos estocásticos, se considera que la serie temporal observada es solo una muestra de un proceso estocástico. Necesitaremos hipótesis.
Métodos descriptivos
Los métodos descriptivos o clásicos se basan en entender, de manera casi empírica, el comportamiento de la serie, pareciéndose más una interpolación que a una metodología rigurosa de estimación y predicción.
Vamos a empezar denotando a la serie como \(X = \left\lbrace X_t \right\rbrace_{t}\), de la cual observamos una muestra\(\left(x_0, \ldots, x_{T} \right)\).
¿Cuál se te ocurre que sería el caso más sencillo de serie temporal?
Descomposición clásica: sin tendencia
El caso más sencillo es considerar que la serie es completamente aleatoria, es decir \(X_t = \varepsilon_t\).
Al término \(\left\lbrace \varepsilon_t \right\rbrace_{t \in T}\) le llamaremos error o innovación y, dado que se supone que no captura ningún patrón y que la serie debe ser finita, tendríamos
\[X_t = \varepsilon_t, \quad {\rm E} \left[X_t \right] = {\rm E} \left[\varepsilon_t \right] = 0, \quad {\rm Var} \left[\varepsilon_t \right] = \sigma_{\varepsilon}^{2} = cte < \infty\] donde normalmente el error sigue una distribución normal de varianza finita.
Además dicho término de error cumplirá que el pasado no proporciona ningún tipo de información sobre el futuro, es decir,
Paso 1: construye un tibble de 5 columnas, donde la primera columna contenga los valores \(t=1, 2, \ldots, 1000\); y donde la segunda columna contenga valores simulados según una normal \(N(0, \sigma = 0.5)\), la tercera con \(\sigma = 1\), la cuarta con \(\sigma = 2\) y la quinta con \(\sigma = 4\).
Código
n <-1000datos <-tibble("t"=1:n,"sd_05"=rnorm(n, mean =0, sd =0.5), "sd_1"=rnorm(n, mean =0, sd =1),"sd_2"=rnorm(n, mean =0, sd =2), "sd_4"=rnorm(n, mean =0, sd =4))
Sin tendencia
Paso 2: haz un gráfico (¿cuál harías?) solo considerando \(t\) y la primera serie sd_05
Código
ggplot(datos) +geom_line(aes(x = t, y = sd_05)) +theme_minimal() +labs(x ="t", y ="X_t",title ="Serie temporal X_t = eps_t con sd = 0.5")
Sin tendencia
Paso 3: ¿cómo deberíamos de transformar los datos para poder pintar todas las series a la vez? Hazte un borrador de cómo sería el código de ggplot para dibujarlo.
La idea es que si tenemos \(p\) series, en lugar de tener \(p\) variables distintas, tengamos una serie «debajo de» otra. Por ejemplo, vamos a definir la primera y pongamos debajo al segunda otra.
n <-1000datos <-tibble("t"=1:n, "X_t"=rnorm(n, mean =0, sd =0.5), "sd"="sd_0.5")datos_tidy <-rbind(datos, tibble("t"=1:n, "X_t"=rnorm(n, mean =0, sd =1), "sd"="sd_1"))datos_tidy
El resto las iremos concatenando de la misma manera, añadiendo filas al datos_tidy que ya tenemos.
datos_tidy <-rbind(datos_tidy, tibble("t"=1:n, "X_t"=rnorm(n, mean =0, sd =2), "sd"="sd_2"))datos_tidy <-rbind(datos_tidy, tibble("t"=1:n, "X_t"=rnorm(n, mean =0, sd =4), "sd"="sd_4"))
Lo anterior se pueda hacer más “conciso” con {tidyverse} haciendo uso de pivot_longer().
datos_tidy <- datos |>pivot_longer(cols ="sd_05":"sd_4", names_to ="sd", values_to ="X_t")
Sin tendencia
Código
ggplot(datos_tidy) +geom_line(aes(x = t, y = X_t, color = sd),alpha =0.7) + ggthemes::scale_color_colorblind() +facet_wrap(~sd) +theme_minimal() +labs(x ="t", y ="X_t", color ="Desv. típica",title ="Serie temporal X_t = eps_t con distintas varianzas")
Sin tendencia
Paso 4. Para automatizarlo, diseña una función tal que le introduzcas como argumento un tamaño muestral \(n\), un \(t\) y un vector de desviaciones típicas, y devuelva en formato tidy data los valores de las series temporales (tantas series como longitud tenga el vector de desviaciones)
Código
time_series_error <-function(n, t =1:n, sd_vec =1) { datos_tidy <-tibble()for (i in1:length(sd_vec)) { datos_tidy <- datos_tidy |>rbind(datos_tidy,tibble("t"= t, "sd"= glue::glue("sd_{sd_vec[i]}"),"X_t"=rnorm(n, mean =0, sd = sd_vec[i]))) }return(datos_tidy)}
Sin tendencia
Esta serie temporal \(X_t = \varepsilon_t\) es la más sencilla que podemos imaginar y no podemos predecirla ya que no hay ningún tipo de patrón determinístico que podamos capturar.
Código
datos <-time_series_error(n =1000, sd =c(0.5, 1, 2, 4))ggplot(datos) +geom_line(aes(x = t, y = X_t, color = sd), alpha =0.7) + ggthemes::scale_color_colorblind() +facet_wrap(~sd) +theme_minimal() +labs(x ="t", y ="X_t", color ="Desv. típica",title ="Serie temporal X_t = eps_t con distintas varianzas")
Clase 4: simulación error + tendencia
Introducción al análisis descriptivo de series temporales
Normalmente una serie temporal suele ser más complejo y lleva al menos incorporada una componente de tendencia o nivel \(\mu_t\) tal que
\[X_t = f \left(\mu_t, \varepsilon_t\right) =^{*} \mu_t + \varepsilon_t, \quad {\rm E} \left[X_t \right] = \mu_t\]\(*\) De momento estamos considerando una descomposición aditiva
Fíjate que ahora \({\rm E} \left[X_t \right] = \mu_t\) ya que la esperanza de la parte aleatoria (ruido) será asumida siempre nula: \(\mu_t\) es el nivel de la serie respecto a la que oscila en el infinito.
Dicha tendencia \(\mu_t\) puede ser a su vez modelada en función de \(t\) y de un vector de parámetros \(\beta\) tal que \(\mu_t := f \left(t, \beta \right)\). Esa función \(f \left( \cdot \right)\) puede ser cualquier función que se te ocurre pero algunas de las tendencias más habituales son:
Constante: \(\mu_t = \beta = \beta_0 = cte\)
Lineal: \(\mu_t = \beta_0 + \beta_1 t\)
Polinómica (no lineal): \(\mu_t = \beta_0 + \beta_1 t + \ldots + \beta_r t^{r}\)
No polinómica: \(\mu_t = \sin \left(\pi t \right)\)
\[X_t = f \left(\mu_t, \varepsilon_t\right) = \beta_0 + \beta_1 t + \varepsilon_t , \quad {\rm E} \left[X_t \right] = \beta_0 + \beta_1 t \to \pm \infty, \quad \widehat{X}_{t + k} = \widehat{\beta}_0 + \widehat{\beta}_1 \left(t + k \right)\]
Con tendencia
Polinómica (no lineal): \(\mu_t = \beta_0 + \beta_1 t + \ldots + \beta_r t^{r}\)
\[X_t = f \left(\mu_t, \varepsilon_t\right) = \beta_0 + \beta_1 t + \ldots + \beta_r t^{r} + \varepsilon_t , \quad \widehat{X}_{t + k} = \widehat{\beta}_0 + \widehat{\beta}_1 \left(t + k \right) + \ldots + \widehat{\beta}_r \left(t + k \right)^r\]
Con tendencia
Si nuestra tendencia está definida de forma paramétrica, para la predicción de los valores en un tiempo futuro \(t + k\) simplemente necesitamos realizar la estimación del vector de parámetros \(\widehat{\beta}\). Para ello recurriremos al método de los minimos cuadrados. Por ejemplo, en el caso de tendencia polinómica
💻 Si aún no lo has hecho, haz una función llamada time_series_error que simule una serie temporal solo con error. Los argumentos deben ser: tamaño muestral n, un vector temporal t y la desv típica sd (debes permitir que pueda ser un vector para simular varias a la vez, ya colocadas en tidydata)
Código
time_series_error <-function(n, t =1:n, sd_vec =1) { datos_tidy <-tibble()for (i in1:length(sd_vec)) { datos_tidy <- datos_tidy |>rbind(datos_tidy,tibble("t"= t, "sd"= glue::glue("sd_{sd_vec[i]}"),"X_t"=rnorm(n, mean =0, sd = sd_vec[i]))) }return(datos_tidy)}time_series_error(n =100, sd =c(0.5, 2))
💻 Tu turno
Ejercicio 2
💻 Usando la función anterior, define time_series_trend_error() que simule una serie temporal con tendencia cte y error, con solo 4 argumentos: n, t, desviación y la constante. Usa dicha función y dibuja.
Código
time_series_trend_error <-function(n =1000, t =1:n, beta_0 =0, sd =1) {# modo R base datos <-time_series_error(n = n, t = t, sd = sd) datos$X_t <- datos$X_t + beta_0# modo tidyverse# datos <- time_series_error(n = n, t = t, sd = sd) |> mutate(X_t = X_t + trend)return(datos)}datos <-time_series_trend_error(n =1000, beta_0 =3, sd =0.5)ggplot(datos, aes(x = t, y = X_t)) +geom_line(alpha =0.7) +geom_smooth(method ="lm", se =FALSE) + ggthemes::scale_color_colorblind() +theme_minimal() +labs(x ="t", y ="X_t", title ="Serie temporal X_t = mu_t + eps_t con mu_t = 3")
💻 Tu turno
Ejercicio 3
💻 Generaliza la función anterior para simular una serie temporal con error y tendencia lineal (donde antes definíamos solo una constante ahora será un vector de coeficientes). Fíjate que la línea de ajuste de ggplot es literal la estimación que haríamos si solo consideramos error + tendencia.
Código
time_series_trend_error <-function(n =1000, t =1:n, beta =c(1, -0.01), sd =1) { datos <-time_series_error(n = n, t = t, sd = sd) datos$X_t <- datos$X_t + (beta[1] + beta[2]*datos$t)return(datos)}datos <-time_series_trend_error(n =1000, beta =c(1, -0.01), sd =0.5)ggplot(datos, aes(x = t, y = X_t)) +geom_line(alpha =0.7) +geom_smooth(method ="lm", se =FALSE) + ggthemes::scale_color_colorblind() +theme_minimal() +labs(x ="t", y ="X_t", title ="Serie temporal X_t = mu_t + eps_t con mu_t = 1 - 0.01*t")
Clase 5: simulación error + tendencia
Introducción al análisis descriptivo de series temporales
💻 Generaliza la función anterior de manera que simule una serie temporal con error y tendencia polinómica (que acepte un vector de parámetros general).
Código
time_series_trend_error <-function(n =1000, t =1:n, beta =c(1, -0.01, 0.001, -0.0001), sd =1) { datos <-time_series_error(n = n, t = t, sd = sd)for (i in1:length(beta)) { datos$X_t <- datos$X_t + beta[i]*(datos$t^(i -1)) } return(datos)}datos <-time_series_trend_error(n =1000, beta =c(1, 0.01, 0.000001, -0.00000001), sd =0.5)ggplot(datos, aes(x = t, y = X_t)) +geom_line(alpha =0.7) +geom_smooth(formula = y ~poly(x, 3), se =FALSE) + ggthemes::scale_color_colorblind() +theme_minimal() +labs(x ="t", y ="X_t", title ="Serie temporal X_t = mu_t + eps_t con mu_t cúbica")
💻 Tu turno
Ejercicio 5
💻 Diseña una función estim_ts_trend_error() que, dada una serie (un tibble de dos columnas t y X_t), nos devuelva la misma tabla pero con una tercera columna con su estimación asumiendo una tendencia polinómica (necesitamos dos argumentos: la tabla y el grado del polinomio; haz uso dentro de poly(), chequea en la ayuda de la función ? poly())
Código
estim_ts_trend_error <-function(datos, degree =1) {if (degree ==0) { modelo <- datos |>lm(formula = X_t ~1) } else { modelo <- datos |>lm(formula = X_t ~poly(t, degree, raw =TRUE)) } datos$X_hat <-predict(modelo, tibble("t"= datos$t))return(datos)}datos <-time_series_trend_error(n =1000, beta =c(1, 0.01, 0.000001, -0.00000001), sd =0.5)# ajustamos un modelo polinómico de tendenciamodelo <- datos |>estim_ts_trend_error(degree =3)
Caso real: AEMET
Como ya te puedes estar imaginando, esta forma de estimar una serie temporal con un polinomio puede ser bastante imprecisa, máxime si aparece en nuestra serie una componente estacional (un patrón periódico).
Vamos a retomar por ejemplo nuestros datos de temperatura del AEMET
Código
ggplot(retiro) +geom_line(aes(x = fecha, y = tmed), linewidth =0.3, alpha =0.7) +theme_minimal() +labs(title ="Temperatura media como SERIE TEMPORAL",x ="t (fecha)", y ="ºC (media)")
Caso real: AEMET
💻 Aplica la función de estimación definida anteriormente a los datos reales del AEMET para incluir 3 nuevas columnas con los 3 métodos de estimación (tendencia constante, lineal y polinómica de grado 3)
retiro_estim <-tibble("fecha"= retiro$fecha, "t"=1:length(fecha), "X_t"= retiro$tmed) |># aplicamos función y renombramos variable de salida de la estimaciónestim_ts_trend_error(degree =0) |>rename(X_t_hat_0 = X_hat) |>estim_ts_trend_error(degree =1) |>rename(X_t_hat_1 = X_hat) |>estim_ts_trend_error(degree =3) |>rename(X_t_hat_3 = X_hat)
Caso real: AEMET
Como observas las predicciones no son precisas cuando hay una componente estacional ya que el ajuste realizado solo se fija en una tendencia con unos coef ctes.
¿Se te ocurre alguna idea para mejorar?
Código
ggplot(retiro_estim |>pivot_longer(-c(fecha, t), names_to ="type", values_to ="pred")) +geom_line(aes(x = fecha, y = pred, color = type),linewidth =0.4, alpha =0.75) + ggthemes::scale_color_colorblind() +theme_minimal() +labs(title ="Temperatura media como SERIE TEMPORAL",x ="t (fecha)", y ="ºC (media)")
Hasta ahora nos hemos centrado sobre todo en tres cosas
Repasar lo que supiéramos de R
Entender cómo simular y estimar series sencillas
Visualizar dichas series
Pero a partir de ahora pulsaremos un poco el acelerador (así quién tenga muchos problemas en la parte de programación, deberá empezar a usar del mail y las tutorías)
Recapitulando
¿Qué deberíamos saber hasta ahora?
Deberíamos tener una función ts_error() parecida a esta para simular un ruido de una varianza dada.
ts_error <-function(n, t =1:n, sd_vec =1) { datos_tidy <-tibble()for (i in1:length(sd_vec)) { datos_tidy <-rbind(datos_tidy,tibble("t"= t, "sd"= glue::glue("sd_{sd_vec[i]}"),"X_t"=rnorm(n, mean =0, sd = sd_vec[i]))) }return(datos_tidy)}ts_error(n =100, sd =c(0.5, 2))
Deberíamos tener una función ts_trend_error parecida a esta para simular una serie formada por tendencia polinómica más ruido, con una varianza dada y un vector de coeficientes dado.
ts_trend_error <-function(n =1000, t =1:n, beta =c(1, -0.01, 0.001, -0.0001), sd =1) { datos <-ts_error(n = n, t = t, sd = sd)for (i in1:length(beta)) { datos$X_t <- datos$X_t + beta[i]*(datos$t^(i -1)) } return(datos)}datos <-ts_trend_error(n =1000, beta =c(1, 0.01, 0.000001, -0.00000001), sd =0.5)
Recapitulando
¿Qué deberíamos saber hasta ahora?
Deberías ser capaz de entender cómo organizar los datos de manera que podamos visualizar de manera sencilla.
Y por último deberías tener una estim_ts_trend_error() similar a esta para que, dada una serie cualquiera, con dos columnas para \(t\) y \(X_t\) (con cualquier nombre) haga la estimación. Fíjate de los argumentos tag_estim y nuevos_valores que hacen para facilitarnos la vida a futuro
estim_ts_trend_error <-function(datos, degree =1, tag_estim =paste0("estim_poly_", degree),col_t ="t", col_X_t ="X_t", nuevos_valores =NULL) { datos <-# versión tidyverse datos |>select(all_of(col_t), all_of(col_X_t)) |>rename(t =all_of(col_t), X_t =all_of(col_X_t))# R base normal# datos <- datos[, c(col_t, col_X_t)]# names(datos)[names(datos) == col_t] <- "t"# names(datos)[names(datos) == col_X_t] <- "X_t"if (degree ==0) { modelo <- datos |>lm(formula = X_t ~1) } else { modelo <- datos |>lm(formula = X_t ~poly(t, degree, raw =TRUE)) } datos[, tag_estim] <-predict(modelo, tibble("t"= datos$t))if (!is.null(nuevos_valores)) { nuevos_datos <-tibble("t"= nuevos_valores, "X_t"=NA) nuevos_datos[, tag_estim] <-predict(modelo, tibble("t"= nuevos_datos$t)) datos <-rbind(datos, nuevos_datos) }return(datos)}
Caso real: AEMET
Con todo esto podemos aplicar nuestras funciones para estimar los datos reales del AEMET, estimando bajo 4 hipótesis: solo ruido, ruido + tendencia cte, ruido + tendencia lineal y ruido + tendencia polinómica de grado 3.
Fíjate que de cada estim_ts_trend_error() solo nos interesa la propia estimación ya que t y X_t es igual siempre, así que podríamos hacer un left_join() de las diferentes tablas.
ggplot(retiro_estim_tidy) +geom_line(aes(x = t, y = X_t, color = serie,linewidth = serie, alpha = serie, linetype = serie)) + ggthemes::scale_color_colorblind() +scale_alpha_manual(values =c(0.85, 0.85, 0.85, 0.85, 1)) +scale_linewidth_manual(values =c(1.1, 1.1, 1.1, 1.1, 0.1)) +scale_linetype_manual(values =c("dotted", "dotted", "dotted", "dotted", "solid")) +theme_minimal() +labs(x ="t", y ="Temperatura (ºC)", title ="Estimación serie AEMET")
Caso real: AEMET
La idea de los métodos de estimación es que podamos usarlos no solo para estimar sino también para predecir en instantes temporales futuros, haciendo uso de esos «nuevos valores» que podemos usar en estim_ts_trend_error().
Por ejemplo, en los datos tenemos solo hasta el 31 de agosto de 2024: ¿cuál es la predicción de los distintos métodos para todo el mes de septiembre, octubre y noviembre?
Caso real: AEMET
Por ejemplo, en los datos tenemos solo hasta el 31 de agosto de 2024: ¿cuál es la predicción de los distintos métodos para lo que queda de 2024 y 2025?
# filtro un poco para que se vea mejorggplot(retiro_predict_tidy |>filter(t >as_date("2015-01-01"))) +geom_line(aes(x = t, y = X_t, color = serie,linewidth = serie, alpha = serie, linetype = serie)) +geom_vline(xintercept =max(retiro$fecha), linetype ="twodash", color ="#a61d0f", alpha =0.5, linewidth =0.9) + ggthemes::scale_color_colorblind() +scale_alpha_manual(values =c(0.85, 0.85, 0.85, 0.85, 1)) +scale_linewidth_manual(values =c(1, 1, 1, 1, 0.2)) +scale_linetype_manual(values =c("dotted", "dotted", "dotted", "dotted", "solid")) +theme_minimal() +labs(x ="t", y ="Temperatura (ºC)", title ="Predicción serie AEMET")
Tendencia «dinámica»
El problema del ajuste anterior es que, amén de la parte puramente estocástica y la tendencia (más o menos compleja que pueda tener), existe una parte ESTACIONAL
Diremos que una serie tiene una componente estacional siempre que presente un patrón que se repite en periodos (aprox.) fijos en el tiempo tal que
Trataremos de manera general con los modelos aditivos ya que, en caso de ser multiplicativo, \(\log \left(X_t\right) =\log \left( \mu_t \right) + \log \left(S_t \right) + \log \left( \varepsilon_t \right)\)
Medias móviles
Existen diferentes estrategias para tener en cuenta la estacionalidad, muchas de ellas basadas en la idea de considerar que la tendencia no es algo estático
La más famosa (y sencilla) probablemente sea la idea de suavizado de medias móviles: en lugar de suavizar la serie considerando una media global, vamos mirar la serie por una pequeña ventana donde para cada \(t\) solo observemos un pequeño trozo de la serie.
La transformación es una media ponderada de \(q + s + 1\) valores donde, según avanzamos, se elimina el dato más antiguo y entra el más nuevo
Si los pesos son todos iguales\(a_j = \frac{1}{N}\), con \(N = q+s+1\), se conoce como media móvil de orden N.
Si \(q = s\) tal que \(a_{-j} = a_j\), para todo \(j=1,\ldots,q\) tal que \(k = 2*q + 1\), se conoce como k MA (moving average simétrica)
El problema de determinar los primeros/últimos valores se conoce como problema de los efectos terminales
Medias móviles
Para calcular una media móvil de orden \(k=2*q+1\) en R podemos hacerlo con filter() del paquete {stats} teniendo la variable ya ordenada (¡cuidado!: si tienes cargado {tidyverse} debes especificar que filter() es)
k <-3datos$x_linear <-predict(datos |>lm(formula = x ~ t), datos)datos <- datos[order(datos$t), ]datos$x_smooth_3ma <- stats::filter(datos$x, filter =rep(1/k, k))datos
Fíjate que cuando aumenta \(k\), la serie es más suavizada (más agresiva con las fluctuaciones) pero hay más datos ausentes (problema de efectos terminales)
💻 Realiza una estimación de la serie aplicando el suavizado de medias móviles con parámetro \(k=7, 14, 28, 365\) a los datos AEMET e incluye las 4 nuevas columnas con los 4 nuevos métodos de estimación al conjunto retiro_estim que teníamos de clases anteriores. Importante: la serie que le pases a stats::filter() no puede tener ausentes.
Código
# versión tidyverse (en R base tabla$variable_nueva <- valor)retiro_estim <- retiro_estim |>drop_na(X_t) |>arrange(t) |>mutate(x_smooth_7ma = stats::filter(X_t, filter =rep(1/7, 7)),x_smooth_14ma = stats::filter(X_t, filter =rep(1/14, 14)),x_smooth_28ma = stats::filter(X_t, filter =rep(1/28, 28)),x_smooth_365ma = stats::filter(X_t, filter =rep(1/365, 365)))ggplot(retiro_estim |>pivot_longer(-c(fecha, t), names_to ="type", values_to ="pred") |>filter(between(fecha, as_date("2016-01-01"), as_date("2022-01-01")))) +geom_line(aes(x = fecha, y = pred, color = type),linewidth =0.7, alpha =0.85) + ggthemes::scale_color_colorblind() +theme_minimal() +labs(title ="Temperatura media como SERIE TEMPORAL",x ="t (fecha)", y ="ºC (media)")
Clase 7: estacionalidad
Introduciendo y estimando la componente estacional
El problema todos los ajuste anteriores es que existe una parte ESTACIONAL que estamos ignorandome deliberadamente. Diremos que una serie tiene una componente estacional \(S_t\) siempre que presente un patrón que se repite en periodos (aprox.) fijos en el tiempo tal que
Trataremos de manera general con los modelos aditivos ya que, en caso de ser multiplicativo, \(\log \left(X_t\right) =\log \left( \mu_t \right) + \log \left(S_t \right) + \log \left( \varepsilon_t \right)\)
Fíjate que ahora \({\rm E} \left[X_t \right] = \mu_t + S_t\) y por tanto nuestra estimación muestral vendrá definida como
¿Qué significa que \(S_t\) estacional (o periódico)?
Componente estacional
En general, diremos que una función \(S_t := S(t)\) es estacional (o periódica) de periodo \(s\) siempre que \(S_t = S_{t + s} = S_{t + 2s}\): si es estacional de periodo \(s\) significa que, cada \(s\) valores, la función se repite.
Fíjate que si \(S_t = S_{t + s}\) también entonces sucede que \(S_t = S_{t - s}\).
Importante
Es importante entender que determinar el periodo no es solo determinar si es «anual» o «mensual»: el valor \(s\) es el número de valores de la componente \(S_t\) que pasan hasta que se repita. Por ejemplo, una serie puede tener una periodicidad anual pero si los datos son recogidos de manera trimestral, \(s = 4\); la misma periodicidad para unos datos recogidos de manera diaria será \(s = 365\).
En el caso de los datos del AEMET, \(s = 365\), no solo porque la temperatura tengan una periodicidad anual obviamente, sino porque los datos son diarios (la misma variable, recogida de manera mensual, tendría \(s = 12\)).
Componente estacional
¿Cómo estimar dicha componente estacional? Los métodos de descomposición clásica tienen la siguiente estructura:
Estimación de la tendencia. Dada la serie original \(X_t\) se realiza una estimación de su nivel o tendencia \(\mu_t\). Como hemos visto tenemos distintas alternativas:
medias móviles (será clave determinar la ventana \(k\))
tendencia polinómica
regresión local (conocida como regresión LOESS o LOWESS: ajusta a los datos una regresión polinómica pero de manera LOCAL, en cada punto solo se utiliza un % de los datos). Ver https://www.statology.org/lowess-smoothing-r/
Tras estimar la tendencia se construye la serie centrada
Fíjate que \(\hat{Y}_t\) es ya una serie sin tendencia, y cuyos valores ya no representan la serie original sino la anomalía que tendría cada \(t\) respecto al nivel global (en el caso de las temperaturas, por ejemplo la anomalía promedio de temperatura entre los distintos días del año y el nivel general de la serie).
Estimar los coeficientes estacionales. El objetivo será obtener un conjunto de \(s\) coeficientes \(\left(S_1, S_2, \ldots, S_{s-2}, S_{s-1}, S_s \right)\) que cumplirán por definición dos condiciones:
Se repiten cada \(s\) valores (por solo necesitamos estimar un tramo).
Su suma es cero (ya que representan las anomalías respecto al nivel general, por lo que si hay valores por encima tienen que existir por debajo)
Componente estacional
Estimar los coeficientes estacionales. El objetivo será obtener un conjunto de \(s\) coeficientes \(\left(S_1, S_2, \ldots, S_{s-2}, S_{s-1}, S_s \right)\).
Para su estimación lo que haremos será calcular, con la serie centrada, la diferencia entre la media de cada periodo estacional y la media general
En el caso del AEMET: \(\widehat{Y}_j\) representa la media (de la serie centrada) de todos los días \(j=1, \ldots, 365\), es decir, la media de los 1-enero, sin importar el año; de los 2-enero, …, de los 31-diciembre (la anomalía de temperatura que en promedio hace el 1-enero respecto a la tendencia general, y así para cada día del año). En este caso \(n = 45\) ya que tenemos datos de 45 años (salvo el final de 2024).
Componente estacional
Estimar las innovaciones. Una vez estimada la tendencia y estacionalidad
Fíjate que la componente estacional estimada cumple también que \(\widehat{S}_t = \widehat{S}_{t + s}\), es decir, solo necesitamos \(s\) coeficientes estimados (fíjate en el subíndice del \(\hat{S}\))
\(\widehat{Y}_t\) representa la estimación de la serie centrada: lo que falta por modelizar es componente estacional
\(X_t - \widehat{S}_t\) representa la estimación de la serie DESESTACIONALIZADA: lo que falta por modelizar es tendencia sin tener en cuenta el efecto estacional (ejemplo: tasa de paro sin el efecto que tienen periodos como verano o navidad)
Algunas observaciones
Cuidado
El método descrito se conoce método clásico de descomposición (o STL si se usa una regresión local para estimar la tendencia en lugar de medias móviles) pero la componente estacional podría ser también estimada por cualquier función periódica o armónica de periodo \(s\), por ejemplo, \(sin(\frac{2 \pi t}{s})\), se repite cada \(s\) valores).
Cuidado
Todo lo anterior también está diseñado bajo la hipótesis de que solo tenemos una periodicidad \(s\): nuestra serie podría tener distintas periodicidades superpuestas (Fourier: toda función periódica puede representarse como suma de funciones seno de distinta frecuencia y amplitud). De momento supondremos siempre un único ciclo.
💻 Tu turno
Ejercicio 1
💻 Renombra la función de estimación de la tendencia polinómica que ya tenemos hecha como estim_ts_trend_poly()
estim_ts_trend_poly <-function(datos, degree =1, tag_estim =paste0("estim_poly_", degree),col_t ="t", col_X_t ="X_t", nuevos_valores =NULL) {# R base normal datos <- datos[, c(col_t, col_X_t)]names(datos)[names(datos) == col_t] <-"t"names(datos)[names(datos) == col_X_t] <-"X_t"if (degree ==0) { modelo <- datos |>lm(formula = X_t ~1) } else { modelo <- datos |>lm(formula = X_t ~poly(t, degree, raw =TRUE)) } datos[, tag_estim] <-predict(modelo, tibble("t"= datos$t))if (!is.null(nuevos_valores)) { nuevos_datos <-tibble("t"= nuevos_valores, "X_t"=NA) nuevos_datos[, tag_estim] <-predict(modelo, tibble("t"= nuevos_datos$t)) datos <-rbind(datos, nuevos_datos) }return(datos)}
💻 Tu turno
Ejercicio 2
💻 Diseña una función general estim_ts_trend() que use estim_ts_trend_poly() que ya tenemos hecha. Dicha función debe tener los mismos parámetros que teníamos en esa función con dos argumentos extras: tipo_trend_estim (que haga la polinómica si "poly" y medias móviles si "MA") y k (la anchura de la ventana; la haremos siempre simétrica y con pesos uniformes de momento).
Código
estim_ts_trend <-function(datos, degree =1, k =NA, tipo_trend_estim ="poly",tag_estim =paste0("estim_", tipo_trend_estim, if_else(is.na(k), degree, k)),col_t ="t", col_X_t ="X_t", nuevos_valores =NULL) {if (tipo_trend_estim =="poly") { estim <-estim_ts_trend_poly(datos, degree, tag_estim, col_t, col_X_t, nuevos_valores) } elseif (tipo_trend_estim =="MA") { estim <- datos |>select(all_of(col_t), all_of(col_X_t)) |>rename(t = col_t, X_t = col_X_t) |># hay que ordenar!arrange(t) |># con !! + := le indicamos que el nombre sale de una variable# (es una técnica llamada lazyeval)mutate(!!tag_estim := stats::filter(X_t, filter =rep(1/k, k), sides =2)) }return(estim) }temp_retiro <-read_csv(file ="./datos/retiro_temp.csv")ejemplo1 <- temp_retiro |>estim_ts_trend(degree =1, tipo_trend_estim ="poly", col_t ="fecha", col_X_t ="tmed")ejemplo2 <- temp_retiro |>estim_ts_trend(k =7, tipo_trend_estim ="MA", col_t ="fecha", col_X_t ="tmed")
💻 Tu turno
Ejercicio 3
💻 Diseña una función ts_detrend() que use estim_ts_trend() que ya tenemos hecha (por lo que tendrá que tener, mínimo, los mismos parámetros que teníamos en esa función), de manera que dado un conjunto de datos, devuelva una nueva columna con la tendencia estimada (según el método que decidamos) y una nueva columna con la serie centrada
Código
ts_detrend <-function(datos, degree =1, k =NA, tipo_trend_estim ="poly",tag_estim =paste0("estim_", tipo_trend_estim, if_else(is.na(k), degree, k)),col_t ="t", col_X_t ="X_t", nuevos_valores =NULL) { estim <-# estimamos tendenciaestim_ts_trend(datos, degree, k, tipo_trend_estim, tag_estim, col_t, col_X_t, nuevos_valores) |># la llamamos siempre estim_trendrename(estim_trend = tag_estim) |># el tipo de estimación lo guardamos en otra variable apartemutate("tipo_trend_estim"= tag_estim,# calculamos serie centrada# si es NA, ponemos la media global en la estimación de la tendencia"estim_trend"=if_else(is.na(estim_trend), mean(X_t, na.rm =TRUE), estim_trend),"detrend"= X_t - estim_trend)return(estim) }ejemplo1 <- temp_retiro |>ts_detrend(degree =1, col_t ="fecha", col_X_t ="tmed")ejemplo2 <- temp_retiro |>ts_detrend(k =365, tipo_trend_estim ="MA", col_t ="fecha", col_X_t ="tmed")
💻 Diseña una función ts_deseason() que, dada una serie cualquiera y un periodo \(s\), nos devuelva una nueva columna con los coeficientes estacionales estimados (con el método visto anteriormente) y una nueva columna con la serie desestacionalizada. Vas a necesitar darle como argumento no solo el nombre de las columnas donde esté \(t\) y \(Y_t\) sino como se llamará la variable de grupo que debes crear dentro usando \(s\).
ts_deseason <-function(datos, s, col_group, col_t ="t", col_Y_t ="Y_t") {# primero renombramos como siempre y ordenamos estim <- datos |>rename(t = col_t, Y_t = col_Y_t) |>arrange(t)# Contruimos variable de grupo estim <- estim |>rowid_to_column(var ="id") |>mutate(!!col_group := ((id -1) %% s) +1) |>select(-id)# Estimamos la estacionacionalidad estim_season <- estim |>summarise("season_coef"=mean(Y_t, na.rm =TRUE), .by = col_group) |>mutate("estim_season"= season_coef -mean(season_coef, na.rm =TRUE))# Lo únimos a los datos estim <- estim |>left_join(estim_season, by = col_group) |>mutate("tipo_estim_season"="clasica","deseason"= Y_t - estim_season) return(estim)}
💻 Tu turno
Ejercicio 2
💻 Tomando los datos de retiro: a) aplica la función ts_detrend() de manera que la tendencia sea estimada con \(MA(k = 365)\); b) a esa tabla aplícale ts_deseason() para estimar \(S_t\), donde la variable de grupo se llama "mes_dia"; c) tras acabar renombra Y_t como detrend
💻 Con estim_retiro del ejercicio anterior, calcula una última columna que sea la estimación total de la serie estim_X_t (estimación de la tendencia + estimación de la estacionalidad). Tras ello pivota como consideres para poder dibujar en la misma gráfica la serie real y la estimación
💻 Tu turno
Ejercicio 3
estim_retiro |>mutate("estim_X_t"= estim_trend + estim_season) |># solo queremos dos de las curvasselect(t, X_t, estim_X_t) |>pivot_longer(cols =-t, names_to ="serie", values_to ="X_t") |>ggplot() +geom_line(aes(x = t, y = X_t, color = serie, alpha = serie)) + ggthemes::scale_color_colorblind() +scale_alpha_manual(values =c(0.9, 0.35)) +theme_minimal() +labs(x ="fecha", y ="Temperatura (ºC)",title ="Estimación decomposición clásica",subtitle ="Tendencia estimada con MA(k = 365)")
💻 Tu turno
Ejercicio 4
💻 Con estim_retiro del ejercicio 2, calcula una última columna que sea la estimación total de la serie estim_X_t (estimación de la tendencia + estimación de la estacionalidad) pero ahora selecciona y pivota como consideres para poder hacer luego una visualización de 6 gráficas (por separado pero en el mismo plot): i) la serie real, ii) la estimación de la tendencia, iii) la estimación de la estacionalidad, iv) la serie sin tendencia, v) la serie sin tendencia ni estacionalidad, vi) la estimación
💻 Tu turno
Ejercicio 4
estim_retiro |>mutate("estim_X_t"= estim_trend + estim_season) |># solo queremos dos de las curvasselect(t, X_t, estim_trend, detrend, estim_season, deseason, estim_X_t) |>pivot_longer(cols =-t, names_to ="serie", values_to ="X_t") |>ggplot() +geom_line(aes(x = t, y = X_t, color = serie)) + ggthemes::scale_color_colorblind() +facet_wrap(~serie, scales ="free_y") +theme_minimal() +labs(x ="fecha", y ="Temperatura (ºC)",title ="Estimación decomposición clásica",subtitle ="Tendencia estimada con MA(k = 365)")
💻 Tu turno
Ejercicio 6
💻 Repite los ejercicios anteriores pero haciendo una pequeña modificación en la función que estima la estacionalidad, permitiendo que pueda estimarla de manera clásica o de manera sinuosoidal \(\widehat{S}_t = \sin(2 \pi t/s)\).
Clase 9: métricas de error
Calibrando nuestra estimación: train/validación y errores
Hasta ahora todo lo que hemos hecho ha sido suponer que el comportamiento de nuestra serie temporal se podía explicar en términos de subpatrones o conductas, que en su forma aditiva, pueden tener la siguiente estructura
\(\widehat{T}_t\) la hemos estimado mediante un polinomio o medias móviles
\(\widehat{S}_t\) la hemos estimado asumiendo que solo hay un periodo
Siempre hemos ponderado las observaciones por igual pero…¿no tendríamos que dar más importancia a los datos más recientes?
Clasificación de Pegel
Imagen extraída de González Velasco, M., & Puerto García, I. M. del. (2009). Series temporales. Universidad de Extremadura.
Estacionalidad multiplicativa = serie heterocedástica (varianza no constante)
Alisado exponencial
Para superar esas limitaciones, en los años 50 y 60 se propusieron una serie de métodos que, más allá de asumir una estructura en los datos, su objetivo era describir la serie en términos de sus propios cambios.
El objetivo de los distintos métodos de alisado/suavizado será intuir la inercia de la serie, eliminando sus posibles fluctuaciones aleatorias, asumiendo que si la tendencia de la serie es ascendente, probablemente la serie siga subiendo, teniendo en consideración la pendiente con la que crece o decrece.
El nombre de alisado/suavizado EXPONENCIAL se debe a que vamos a ponderar el pasado de la serie con un conjunto de pesos que, normalmente, decrecerán de manera exponencial (ejemplo: si la observación inmediatamente anterior tiene un peso de \(0.5\), la siguiente tendrá un peso de \(0.5^2 = 0.25\))
Alisado simple
Imagina que tenemos el valor de la serie en \(X_t\) y su estimación \(\widehat{X}_t\)… ¿cómo poder predecir \(\widehat{X}_{t+1}\) solo usando ambos valores?
La idea que tuvo Holt (1956) fue la de hacer una media ponderada del valor real \(X_t\) y su estimación \(\widehat{X}_t\)
Si \(\theta \to 1\), entonces \(\widehat{X}_{t+1} \to \widehat{X}_t\) el modelo de alisado produce predicciones casi constantes sin muchas variaciones
Si \(\theta \to 0\), entonces \(\widehat{X}_{t+1} \to X_t\) el modelo de alisado produce predicciones muy variables dependientes del último valor observado (con mucha fluctuación causada por la propia aleatoriedad de la serie)
Alisado simple
Si hacemos lo mismo para la estimación en \(X_t\) tenemos
Fíjate en la serie geométrica\(\sum_{j=0}^{\infty} \theta^j\). ¿Recuerdas cuál era la suma de una serie geométrica cuya razón es menor que 1?
Alisado simple
Dado que \(\theta < 1\), tenemos que \(\sum_{j=0}^{\infty} \theta^j = \frac{1}{1 - \theta} \rightarrow \left( 1 - \theta \right) \sum_{j=0}^{\infty} \theta^j = 1\), así que \(\left( 1 - \theta \right) \sum_{j=0}^{t-1} \theta^j X_{t-j}\) representa una media ponderada del pasado (pesos decrecientes geométricamente con suma 1 en el infinito).
Si te fijas además todas las observaciones son influyentes pero su influencia va decreciendo.
Alisado simple
Fíjate que la fórmula anterior se puede generalizar para predicción a horizonte \(h\): si solo disponemos de información hasta el instante \(t\), tenemos que
El método de alisado simple produce «flat predictions»: siempre devuelve a horizonte \(h\) el último valor predicho \(\widehat{X}_{t+1|t}\) (asume que la tendencia es localmente constante).
Alisado simple
Fíjate que el método de alisado simple es una simple media ponderada, de una forma muy particular, pero que podríamos considerar otro tipo de métodos basados en promedios:
Si consideramos el conocido como método naive o trivial, tenemos que la predicción a horizonte \(h\) (con la información disponible hasta \(t\)) es
\[\widehat{X}_{t+h|t} = X_t\]
la única observación importante es la última. Esto es equivalente a un alisado simple con \(\theta = 0\)
Alisado simple
Fíjate que el método de alisado simple es una simple media ponderada, de una forma muy particular, pero que podríamos considerar otro tipo de métodos basados en promedios:
Si consideramos el conocido como mean method, lo que hemos llamado poly0 (hace simplemente la media cte de todos los valores), la predicción a horizonte \(h\) (con la información disponible hasta \(t\)) es
una tendencia que va variando muy lentamente con el tiempo (ya que el futuro se parece más a los valores pasados inmediatos que a los valores más lejanos).
Para el método de alisado simple necesitamos conocer el valor inicial y encontrar el mejor \(\theta\) (por ejemplo, mediante mínimos cuadrados de los errores observados para distintos valores).
Supondremos \(\theta\) constante pero existen métodos (ver aquí) para considerar que \(\theta\) se adapta a su vez a lo largo del tiempo
Ejemplo en R
Vamos a trabajar con el paquete {tsibbledata} que nos permite usar varios ejemplos de series temporales. En concreto vamos a usar el dataset global_economy que contiene distintas estadísticas económicas para distintos años países, y nos vamos a centrar solo en datos de exportaciones de Algeria
library(tsibbledata)library(tsibble)algeria_economy <- global_economy |>filter(Country =="Algeria") |># solo hay un país, eliminamos keyupdate_tsibble(key =NULL)algeria_economy
# A tsibble: 58 x 9 [1Y]
Country Code Year GDP Growth CPI Imports Exports Population
<fct> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 Algeria DZA 1960 2723648552. NA NA 67.1 39.0 11124888
2 Algeria DZA 1961 2434776646. -13.6 NA 67.5 46.2 11404859
3 Algeria DZA 1962 2001468868. -19.7 NA 20.8 19.8 11690153
4 Algeria DZA 1963 2703014867. 34.3 NA 36.8 24.7 11985136
5 Algeria DZA 1964 2909351793. 5.84 NA 29.4 25.1 12295970
6 Algeria DZA 1965 3136258897. 6.21 NA 25.8 22.6 12626952
7 Algeria DZA 1966 3039834559. -4.80 NA 24.7 26.0 12980267
8 Algeria DZA 1967 3370843066. 9.45 NA 21.6 23.4 13354197
9 Algeria DZA 1968 3852115817. 10.8 NA 24.2 23.1 13744387
10 Algeria DZA 1969 4257218772. 8.43 2.57 28.1 23.8 14144438
# ℹ 48 more rows
Ejemplo en R
Si observamos la serie no se aprecia una tendencia clara ni una estacionalidad evidente, perfecto para nuestro método de alisado simple.
Código
ggplot(algeria_economy) +geom_line(aes(x = Year, y = Exports)) +scale_x_continuous(n.breaks =20) +theme_minimal()
Ejemplo en R
¿Cómo realizar la estimación de alisado simple?
Vamos a usar el universo de paquetes {tidyverts} (una forma de trabajar en modo tidyverse con series temporales), de los cuales ya conocemos {tsibbledata} y {tsibble}. Vamos a instalar también el paquete {fable} que nos proporciona herramientas para la predicción de series temporales
library(fable)
Ejemplo en R
Dentro de este paquete existe una función llamada model() que nos permite ajustar distintos modelos a los datos.
ETS(var_objetivo ~ componentes del modelo): para ajustar los métodos de alisado/suavizado exponencial necesitamos especificarle dentro de model() que nuestro modelo es de tipo ETS() (exponential time-series smoothing model)
# NO EJECUTES QUE ESTÁ AÚN SIN COMPLETARfit_algeria <- algeria_economy |>model(ETS(Exports ~ ...))
Ejemplo en R
Una vez que hemos determinado nuestra variable objetivo, las componentes del modelo se pueden incorporar con las funciones error(), trend() y season()
Cada una de las funciones admite dentro los siguientes tipos:
"A": aditivo
"M": multiplicativo
"N": ninguno (sin esa componente)
Ejemplo en R
En nuestro caso hemos visto que el alisado simple es para predecir una serie asumiendo que no hay tendencia ni estacionalidad ("N" en ambos)
fit_algeria <- algeria_economy |># fíjate que le podemos poner nombre al modelomodel(alisado_simple =ETS(Exports ~error("A") +trend("N") +season("N")))fit_algeria
# A mable: 1 x 1
alisado_simple
<model>
1 <ETS(A,N,N)>
Ejemplo en R
Se guarda en un objeto de tipo modelo (mdl_df o mable): para poder ver la info en un formato tabulado basta hacer augment(): .fitted guarda las estimaciones, .resid residuales estimados (en este caso coincide con .innov)
Tras modelizar la serie podemos predecir el futuro con forecast(h = n_instantes_futuros) indicándole los valores futuros (para convertirlo a un formato tabla, dado que ya no es un modelo como antes, debemos convertirlo con as_tibble() o as_tsibble())
Fíjate que en el caso de alisado simple todas las predicciones son la misma (ya que asume que es localmente constante)
Podemos visualizar el ajuste (del pasado) haciendo augment() + pivot_longer() (seleccionándole antes las columnas a pivotar)
Código
fit_algeria |>augment() |>as_tibble() |>select(Year, Exports, .fitted) |>rename(export_real = Exports, export_fit = .fitted) |>pivot_longer(cols =-Year, names_to ="serie", values_to ="X_t") |>ggplot() +geom_line(aes(x = Year, y = X_t, color = serie)) + ggthemes::scale_color_colorblind() +scale_x_continuous(n.breaks =20) +theme_minimal() +labs(y ="% del PIB", title ="Exports: Algeria")
Fíjate como el alisado simple “llega tarde” a los cambios ya que lo que hace es asumir que todo va a seguir igual
Ejemplo en R
Podemos unir las predicciones futuras a los valores ajustados del pasado con un full_join() (conviertiendo lsa tablas a tibble, seleccionando las columnas que nos interesa y renombrándolas)
# A tibble: 65 × 4
Year export_real export_fit export_predict
<dbl> <dbl> <dbl> <dbl>
1 1960 39.0 39.5 NA
2 1961 46.2 39.1 NA
3 1962 19.8 45.1 NA
4 1963 24.7 23.8 NA
5 1964 25.1 24.6 NA
6 1965 22.6 25.0 NA
7 1966 26.0 23.0 NA
8 1967 23.4 25.5 NA
9 1968 23.1 23.8 NA
10 1969 23.8 23.2 NA
# ℹ 55 more rows
Ejemplo en R
Tras ello podemos visualizar (haciendo un pivot_longer() previo as usual)
Código
fit_predict_algeria |>pivot_longer(cols =-Year, names_to ="serie", values_to ="X_t") |>ggplot() +geom_line(aes(x = Year, y = X_t, color = serie)) + ggthemes::scale_color_colorblind() +scale_x_continuous(n.breaks =20) +theme_minimal() +labs(y ="% del PIB", title ="Exports: Algeria")
Ejemplo en R
También podemos simplificar el código de la visualización haciendo uso de autoplot(), que te incluye además unos intervalos de confianza.
Código
predict_algeria |>autoplot(algeria_economy) +geom_line(data = fit_algeria |>augment(),aes(y = .fitted), col ="#D55E00") +scale_x_continuous(n.breaks =20) +theme_minimal() +labs(y ="% del PIB", title ="Exports: Algeria")
Clase 11: modelos con fable
Modelos con fable: MEAN()
No solo vamos a poder aplicar el alisado en esta nueva lógica de programación, sino que vamos a poder programar de manera sencilla los modelos que hemos investigado hasta ahora
MEAN(var_objetivo): la predicción a horizonte \(h\) se define como la media de todos los valores conocidos de la serie \(\widehat{X}_{t+h|t} = \frac{1}{t} \sum_{j=1}^{t} X_t\) hasta el instante \(t\) (predicción constante, la que antes llamábamos poly0)
TSLM(var_objetivo ~ formula): la predicción a horizonte \(h\) se define como una regresión polinómica tal que \(\widehat{X}_{t+h|t} = \beta_0 + \sum_{j=1}^{p} \beta_j (t+h)^p\)
RW(var_objetivo ~ drift()): modelo conocido como random walk o paseo aleatorio con drift o tendencia tal que \(X_{t+1} = c + X_{t} + \varepsilon_t\). La predicción a horizonte \(h\) se define entonces como
\[\widehat{X}_{t+h|t} = c + X_{(t+h-1)|t} \simeq c + \widehat{X}_{(t+h-1)|t} \simeq 2*c + \widehat{X}_{(t+h-2)|t} \simeq ... \simeq c*h + X_{t}\]
NAIVE(var_objetivo): el modelo naive o trivial es una simplificación del random walk con \(c=0\) (sin tendencia, de hecho NAIVE() y RW() sin drift hacen lo mismo. La predicción a horizonte \(h\) se define entonces como simplemente el último valor conocido en tiempo \(t\).
SNAIVE(var_objetivo ~ lag()): similar al modelo anterior solo que en lugar de repetir el último valor repite los últimos valores estacionales (por ejemplo, el último mes, el último año, etc).
La predicción a horizonte \(h\) se define entonces como
\(k\) es la parte entera de \((h-1)/s\) (el número de periodos completos que han pasado hasta \(t+h\)). Por ejemplo, si \(h = 26\) y \(s = 12\), \((h-1)/s = 25/12 = 2.08333\), cuya parte entera es \(k = 2\) (han pasado 2 años completos hasta \(h = 26\)). Si \(s = 12\), la predicción de cualquier de los enero futuros será igual a la predicción del último enero conocido.
Modelos con fable: SNAIVE()
SNAIVE(var_objetivo ~ lag()): similar al modelo anterior solo que en lugar de repetir el último valor repite los últimos valores estacionales (por ejemplo, el último mes, el último año, etc).
También podemos con el paquete {feasts} hacer uso de las decomposiciones que hemos aprendido
classical_decomposition(var_objetivo ~ season(s), type = ...): descomposición clásica aprendida de tipo multiplicativa o aditiva. La predicción a horizonte \(h\) se define tal que
library(feasts)fit_airpass <- airpass |>model(cts =classical_decomposition(value ~season(12), type ="additive"))
Métodos de decomposición
Para este tipo particular de métodos (métodos de decomposición) existe una función components() que automáticamente nos devuelve la tendencia estimada (trend), la estacionalidad (seasonal), el residuo estimado (random) y la serie desestacionalizada (season_adjust, serie - seasonal)
fit_airpass <- airpass |>model(cts =classical_decomposition(value ~season(12), type ="additive"))fit_airpass |>components()
# A dable: 144 x 7 [1M]
# Key: .model [1]
# : value = trend + seasonal + random
.model index value trend seasonal random season_adjust
<chr> <mth> <dbl> <dbl> <dbl> <dbl> <dbl>
1 cts 1949 Jan 112 NA -24.7 NA 137.
2 cts 1949 Feb 118 NA -36.2 NA 154.
3 cts 1949 Mar 132 NA -2.24 NA 134.
4 cts 1949 Apr 129 NA -8.04 NA 137.
5 cts 1949 May 121 NA -4.51 NA 126.
6 cts 1949 Jun 135 NA 35.4 NA 99.6
7 cts 1949 Jul 148 127. 63.8 -42.6 84.2
8 cts 1949 Aug 148 127. 62.8 -42.1 85.2
9 cts 1949 Sep 136 128. 16.5 -8.48 119.
10 cts 1949 Oct 119 129. -20.6 11.1 140.
# ℹ 134 more rows
Métodos de decomposición
Con autoplot() nos visualiza cada componente
fit_airpass <- airpass |>model(cts =classical_decomposition(value ~season(12), type ="additive"))fit_airpass |>components() |>autoplot()
Métodos de decomposición
Así quedaría en modo multiplicativo
fit_airpass <- airpass |>model(cts =classical_decomposition(value ~season(12), type ="mult"))fit_airpass |>components() |>autoplot()
Métodos de decomposición
Fíjate que el modo multiplicativo es equivalente a tomar el logaritmo del valor objetivo
fit_airpass <- airpass |>model(cts =classical_decomposition(log(value) ~season(12), type ="additive"))fit_airpass |>components() |>autoplot()
Métodos de decomposición
STL(var_objetivo ~ trend(window = ..., degre = ...) + season(period = ..., window = "periodic")): la predicción a horizonte \(h\) se define también como una descomposición clásica en tendencia y estacionalidad pero ahora la tendencia es estimada mediante una regresión local conocida como regresión LOESS o LOWESS. Dicha regresión ajusta a los datos una regresión polinómica (degree) pero de manera LOCAL, en cada punto solo se utiliza un número window de observaciones).
Para este tipo particular de métodos (métodos de decomposición STL) components() nos devuelve la tendencia estimada (trend), la estacionalidad (season_s), el residuo estimado (remainder) y la serie desestacionalizada (season_adjust, serie - season_s)
La ecuación anterior la podemos expresar de manera iterativa como
\[\begin{eqnarray}\widehat{X}_{t+h|t} &=& \ell_t \quad \text{predicción en base a componentes} \nonumber \\ \ell_t &=& \theta \ell_{t-1}+ \left( 1 - \theta \right)X_t \quad \text{suavizado de componentes} \end{eqnarray}\]
Alisado doble (lineal) de Holt
\[\begin{eqnarray}\widehat{X}_{t+h|t} &=& \ell_t \quad \text{predicción (en base a componentes)} \nonumber \\ \ell_t &=& \theta \ell_{t-1}+ \left( 1 - \theta \right)X_t \quad \text{suavizado de componentes} \end{eqnarray}\]
En 1957 Holt propuso extender este método a datos con tendencia, incluyendo ahora dos componentes: nivel \(\ell_t\) y tendencia \(\mu_t\) (realizando un doble suavizado iterativo)
Predicción a horizonte \(h\): último nivel conocido más \(h\) veces la última tendencia conocida (en \(h\) instantes avanza \(h\) veces la tendencia)
Estimación nivel: media ponderada entre el último valor de la serie y \(\widehat{X}_{t|t-1}\) (predicción a horizonte \(h = 1\) en tiempo \(t\))
Estimación tendencia: media ponderada entre el último valor de la tendencia y la estimación de la tendencia a tiempo \(t\) (diferencia de nivel \(\ell_t\) y \(\ell_{t-1}\))
# A tsibble: 288 x 6 [1M]
# Key: .model [2]
.model index value .fitted .resid .innov
<chr> <mth> <dbl> <dbl> <dbl> <dbl>
1 alisado_simple 1949 Jan 112 112. -0.0948 -0.0948
2 alisado_simple 1949 Feb 118 112. 6.00 6.00
3 alisado_simple 1949 Mar 132 118. 14.0 14.0
4 alisado_simple 1949 Apr 129 132. -3.00 -3.00
5 alisado_simple 1949 May 121 129. -8.00 -8.00
6 alisado_simple 1949 Jun 135 121. 14.0 14.0
7 alisado_simple 1949 Jul 148 135. 13.0 13.0
8 alisado_simple 1949 Aug 148 148. 0.00130 0.00130
9 alisado_simple 1949 Sep 136 148. -12.0 -12.0
10 alisado_simple 1949 Oct 119 136. -17.0 -17.0
# ℹ 278 more rows
Alisado doble (lineal) de Holt
Código
predicciones |>autoplot(airpass, level =NULL) +geom_line(data = estimaciones, aes(x = index, y = .fitted, color = .model)) +geom_line(data = predicciones, aes(x = index, y = .mean, color = .model)) +scale_color_manual(values = ggthemes::colorblind_pal()(3)[-1]) +theme_minimal()
Alisado doble de Gardner-McKenzie
Las predicciones generadas por el alisado doble de Holt nos dan una tendencia constante (creciente o decreciente) indefinidamente hacia el infinito. Sin embargo, la evidencia empírica indica que estos métodos tienden a sobrepronosticar a horizontes de previsión más largos.
Para solventar esto Gardner y McKenzie (1985) introdujeron un parámetro de amortiguación \(\phi\): la predicción a futuro empieza siendo una línea recta que termina doblándose hasta acabar en una recta.
# A tsibble: 432 x 6 [1M]
# Key: .model [3]
.model index value .fitted .resid .innov
<chr> <mth> <dbl> <dbl> <dbl> <dbl>
1 alisado_simple 1949 Jan 112 112. -0.0948 -0.0948
2 alisado_simple 1949 Feb 118 112. 6.00 6.00
3 alisado_simple 1949 Mar 132 118. 14.0 14.0
4 alisado_simple 1949 Apr 129 132. -3.00 -3.00
5 alisado_simple 1949 May 121 129. -8.00 -8.00
6 alisado_simple 1949 Jun 135 121. 14.0 14.0
7 alisado_simple 1949 Jul 148 135. 13.0 13.0
8 alisado_simple 1949 Aug 148 148. 0.00130 0.00130
9 alisado_simple 1949 Sep 136 148. -12.0 -12.0
10 alisado_simple 1949 Oct 119 136. -17.0 -17.0
# ℹ 422 more rows
Alisado doble de Gardner-McKenzie
Código
predicciones |>autoplot(airpass, level =NULL) +geom_line(data = estimaciones, aes(x = index, y = .fitted, color = .model)) +geom_line(data = predicciones, aes(x = index, y = .mean, color = .model)) +scale_color_manual(values = ggthemes::colorblind_pal()(4)[-1]) +theme_minimal()
Clase 13: métodos de alisado
Alisado triple de Holt-Winters
Veamos un repaso de los métodos de alisado explicados.
Alisado simple: la predicción de la serie a horizonte \(h\) es el último nivel de la serie estimado (donde se estimó que la serie la última vez). Dicho nivel se obtiene de manera iterativa ponderando el nivel previo (la estimación anterior) y la última observación real
\[\begin{eqnarray}\color{red}{\widehat{X}_{t+h|t}} &=& \color{purple}{\ell_t} \quad \text{predicción (en base a componentes)} \nonumber \\ \color{purple}{\ell_t} &=& \theta \color{purple}{\ell_{t-1}}+ \left( 1 - \theta \right)X_t \quad \text{suavizado de componentes} \end{eqnarray}\]
Alisado triple de Holt-Winters
Alisado doble: la predicción de la serie a horizonte \(h\) es el último nivel de la serie estimado (donde se estimó la serie la última vez) y la última pendiente estimada: para predecir donde estará en la carretera a las 5h uso donde estaba a las 4h y la pendiente de la carretera en ese momento.
Como antes, el nivel se obtiene iterativamente ponderando el nivel previo (estimación anterior) y la última observación. La pendiente se estima de la misma forma, usando la predicción de la estimación y la última estimación disponible.
Unos años más tarde, Holt y Winters (1960) ampliaron los métodos anteriores para poder capturar la estacionalidad. Se conoce como alisado triple ya que tendremos ahora tres ecuaciones de suavizado: suavizado del nivel, suavizado de la tendencia y suavizado de la componente estacional.
Ahora tendremos 4 parámetros:
parámetros de suavizado \(\left(\theta_1, \theta_2, \theta_3 \right)\)
parámetro de estacionalidad \(s\)
Y tendremos dos tipos de formas de incluir la estacionalidad: aditivo (varianza constante a lo largo del tiempo, la componente estacional se expresa en términos absolutos) y multiplicativo (la componente estacional se expresa en términos relativos, en porcentajes).
Alisado triple de Holt-Winters
La idea es que ahora, al nivel y la tendencia, para predecir la serie a horizonte \(h\) le sumaremos una componente estacional. Imagina que tenemos una serie mensual con \(s= 12\) (periodicidad anual). ¿Cuánto valdrá la predicción de \(t+1\) con los datos hasta \(t\)?
Fíjate que en el **nivel y la tendencia el valor más reciente que podemos usar es \(t\)* pero \(s_{t+2}\) sí está disponible ya que es periódica, es decir, \(s_{t+2} = s_{(t+2) - s}\)
¿Y si \(h = 13\)?
Alisado triple de Holt-Winters
\[\color{red}{\widehat{X}_{t+13|t}} = \color{purple}{\ell_t} + 13*\color{green}{\mu_t} + \color{blue}{s_{t+13}}\] Pero como tenemos una periodicidad \(s = 12\), entonces \(s_{t+13} = s_{t+13-12} = s_{t+1}\): le hemos restado la cantidad de años completos que había en \(h = 13\).
De manera coloquial podemos expresarlo de manera general como
\[\color{red}{\widehat{X}_{t+h|t}} = \color{purple}{\ell_t} + h*\color{green}{\mu_t} + \color{blue}{s_{(t+h) - s*per.comp}}\] Se conoce como alisado triple ya que tendremos que suavizar de manera iterativa tres componentes: nivel, tendencia y estacionalidad.
Si entendemos el nivel como «en que punto está la serie», y dado que la componente estacional es cíclica («por dónde va la tasa de paro» debería ser ajeno a si es navidad o verano), el nivel será suavizado ponderando la última predicción y la serie desestacionalizada por el periodo previo
Por último, el suavizado de la estacionalidad será similar: una ponderación entre la estimación de la misma y la última estacionalidad (hace s periodos)
# A tsibble: 432 x 6 [1M]
# Key: .model [3]
.model index value .fitted .resid .innov
<chr> <mth> <dbl> <dbl> <dbl> <dbl>
1 alisado_simple 1949 Jan 112 112. -0.0948 -0.0948
2 alisado_simple 1949 Feb 118 112. 6.00 6.00
3 alisado_simple 1949 Mar 132 118. 14.0 14.0
4 alisado_simple 1949 Apr 129 132. -3.00 -3.00
5 alisado_simple 1949 May 121 129. -8.00 -8.00
6 alisado_simple 1949 Jun 135 121. 14.0 14.0
7 alisado_simple 1949 Jul 148 135. 13.0 13.0
8 alisado_simple 1949 Aug 148 148. 0.00130 0.00130
9 alisado_simple 1949 Sep 136 148. -12.0 -12.0
10 alisado_simple 1949 Oct 119 136. -17.0 -17.0
# ℹ 422 more rows
Alisado triple de Holt-Winters
Fíjate que las «montañitas» ahora suben pero no incrementan su altura, simplemente acaban un poco más arriba ya que tenemos una componente de tendencia.
Código
predicciones |>autoplot(airpass, level =NULL) +geom_line(data = estimaciones, aes(x = index, y = .fitted, color = .model)) +geom_line(data = predicciones, aes(x = index, y = .mean, color = .model)) +scale_color_manual(values = ggthemes::colorblind_pal()(4)[-1]) +theme_minimal()
# A tsibble: 576 x 6 [1M]
# Key: .model [4]
.model index value .fitted .resid .innov
<chr> <mth> <dbl> <dbl> <dbl> <dbl>
1 alisado_simple 1949 Jan 112 112. -0.0948 -0.0948
2 alisado_simple 1949 Feb 118 112. 6.00 6.00
3 alisado_simple 1949 Mar 132 118. 14.0 14.0
4 alisado_simple 1949 Apr 129 132. -3.00 -3.00
5 alisado_simple 1949 May 121 129. -8.00 -8.00
6 alisado_simple 1949 Jun 135 121. 14.0 14.0
7 alisado_simple 1949 Jul 148 135. 13.0 13.0
8 alisado_simple 1949 Aug 148 148. 0.00130 0.00130
9 alisado_simple 1949 Sep 136 148. -12.0 -12.0
10 alisado_simple 1949 Oct 119 136. -17.0 -17.0
# ℹ 566 more rows
Alisado triple de Holt-Winters
Código
predicciones |>autoplot(airpass, level =NULL) +geom_line(data = estimaciones, aes(x = index, y = .fitted, color = .model)) +geom_line(data = predicciones, aes(x = index, y = .mean, color = .model)) +scale_color_manual(values = ggthemes::colorblind_pal()(5)[-1]) +theme_minimal()
Clase 14: evaluación
Diagnosis errores
Como ya hicimos una vez en clase, es importante realizar una diagnosis correcta de los residuales. Para ello tenemos la función gg_tsresiduals() del paquete {feasts} que nos permite visualizar su evolución temporal, sus autocorrelaciones y su distribución
Por ejemplo en el caso del alisado simple observamos como los errores van aumentando según avanza el tiempo debido a la heterocedasticidad de la serie, con magnitudes muy elevadas, y teniendo además unas autocorrelaciones muy altas en algunos retardos.
En el caso del alisado triple multiplicativo observamos como los errores ya no aumentan según avanza el tiempo (magnitudes pequeñas), y teniendo además unas autocorrelaciones casi todas dentro de la banda.
Como sucede en otro tipo de modelos, es importante darse cuenta de que a la hora de evaluar un modelo tendremos que considerar dos aspectos diferentes:
¿Cómo funciona el modelo con los datos que conoce? Es lo que hemos llamado hasta ahora estimaciones.
¿Cómo funcionaría el modelo con unos datos que no conoce? Es lo que hemos llamado hasta ahora predicciones
El problema es que las predicciones hasta ahora no podíamos evaluarlas ya que el dato real del futuro no lo tenemos…¿Y si partimos nuestras series (train y test) de forma que solo le dejamos usar una parte de la información para diseñar el modelo, y así poder evaluarlo en el otro subconjunto?
Train vs test
Vamos a dividir nuestra serie temporal de pasajeros aéreos en train y test. Normalmente el % de datos en test es igual al horizonte al que queremos evaluar cómo funciona (por ejemplo, si \(h = 24\), usaremos los 2 últimos años como test).
# A tibble: 4 × 10
.model .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
<chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 alisado_simple Training 2.22 33.6 25.7 0.375 8.96 0.802 0.925 0.303
2 alisado_doble Training 0.568 33.5 25.6 -0.362 9.00 0.800 0.923 0.302
3 alisado_triple_ad Training 0.964 17.0 12.8 0.354 5.22 0.400 0.469 0.188
4 alisado_triple_mult Training 1.66 11.3 7.91 0.410 2.80 0.247 0.310 0.256
accuracy(predicciones, airpass_test)
# A tibble: 5 × 10
.model .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
<chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 alisado_doble Test 93.2 117. 93.2 18.8 18.8 NaN NaN 0.708
2 alisado_simple Test 115. 137. 115. 23.6 23.6 NaN NaN 0.728
3 alisado_triple_ad Test 79.3 91.2 79.3 16.5 16.5 NaN NaN 0.720
4 alisado_triple_mult Test 34.6 39.9 34.8 7.54 7.58 NaN NaN 0.566
5 mean_cte Test 206. 219. 206. 44.2 44.2 NaN NaN 0.728
Error medio (ME): definido como \(\frac{1}{n} \sum_{t=1}^{n} \widehat{\varepsilon}_{t}\)
Error medio absoluto (MAE): definido como \(\frac{1}{n} \sum_{t=1}^{n} \left| \widehat{\varepsilon}_{t} \right|\)
Error cuadrático medio (RSME): definido como \(\frac{1}{n} \sum_{t=1}^{n} \widehat{\varepsilon}_{t}^2\)
Ambos dependen de la escala de los datos
accuracy(predicciones, airpass_test)
# A tibble: 5 × 10
.model .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
<chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 alisado_doble Test 93.2 117. 93.2 18.8 18.8 NaN NaN 0.708
2 alisado_simple Test 115. 137. 115. 23.6 23.6 NaN NaN 0.728
3 alisado_triple_ad Test 79.3 91.2 79.3 16.5 16.5 NaN NaN 0.720
4 alisado_triple_mult Test 34.6 39.9 34.8 7.54 7.58 NaN NaN 0.566
5 mean_cte Test 206. 219. 206. 44.2 44.2 NaN NaN 0.728
Ambos son adimensionales pero pueden tender a infinito si \(X_t \to 0\)
accuracy(predicciones, airpass_test)
# A tibble: 5 × 10
.model .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
<chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 alisado_doble Test 93.2 117. 93.2 18.8 18.8 NaN NaN 0.708
2 alisado_simple Test 115. 137. 115. 23.6 23.6 NaN NaN 0.728
3 alisado_triple_ad Test 79.3 91.2 79.3 16.5 16.5 NaN NaN 0.720
4 alisado_triple_mult Test 34.6 39.9 34.8 7.54 7.58 NaN NaN 0.566
5 mean_cte Test 206. 219. 206. 44.2 44.2 NaN NaN 0.728
Validación cruzada
Como suele ser habitual en el campo de la calibración de modelos, una opción muy habitual es la de la validación:
Construir distintos modelos con la información de train
Usar los conjuntos de la validación para evaluar los modelos (o qué configuración de hiperparámetros) y decidir cuál de ellos es mejor
Una vez elegido el modelo, volver a lanzarlo y evaluarlo en test
Validación cruzada
Una de las formas de validación más habitual es la validación cruzada: las observaciones del conjunto de train van rotando su rol.
Por ejemplo, si tenemos 100 observaciones en train, podemos hacer 100 iteraciones de validación, de manera que en cada una entrenamos el modelo con 99 de ellas y otra queda reservada solo para evaluar los modelos.
Validación cruzada
En el caso de las series temporales una estrategia habitual suele ser la siguiente:
Descartar las primeras \(n\) observaciones para validación: habrá un conjunto mínimo que siempre formará parte de train
Iteración i: entrenamos con las primeras \(n+i\) observaciones, evaluamos con una única observación \(n+i+1\).
Realizamos el promedio de las métricas de evaluación obtenidas de los conjuntos de validación.
Validación cruzada
Fíjate que lo anterior está basado en una one-step forecast (predicción a horizonte \(h = 1\)), pero quizás nuestro interés esté en ver cómo funciona nuestro método a horizontes de predicción mayores
Descartar las primeras \(n\) observaciones para validación: habrá un conjunto mínimo que siempre formará parte de train
Iteración i: entrenamos con las primeras \(n+i\) observaciones, evaluamos con una única observación \(n+i+h\).
Realizamos el promedio de las métricas de evaluación obtenidas de los conjuntos de validación.
Validación cruzada
El ejemplo inferior es para \(h = 4\).
Validación cruzada
Para generar los subconjuntos vamos primero como antes a dividir nuestro dataset en train y test.
Tras ello vamos a generar los subconjuntos de validación usando trian con stretch_tsibble(), indicándole el número de valores iniciales que siempre estarán en train, el tamaño que queremos incrementar los sucesivos conjuntos y un identificador de cada slot
Por ejemplo, vamos a reservar los 2 primeros años y vamos a avanzar a horizonte 1.
¿Cómo visualizar las métricas de validación cruzada
Código
airpass_fit |>accuracy() |>ggplot(aes(x = .model, y = RMSE, fill = .model, color = .model)) +geom_boxplot(alpha =0.5) +geom_jitter(width =0.25, alpha =0.7) + ggthemes::scale_color_colorblind() + ggthemes::scale_fill_colorblind() +theme_minimal()
Ejemplo real: AEMET again
Cargamos los datos
Código
library(readr) # de tidyverseretiro <-read_csv(file ="./datos/retiro_temp.csv")retiro
# A tibble: 16,314 × 8
fecha id_station nombre provincia altitud tmed tmin tmax
<date> <dbl> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
1 2000-01-01 3195 MADRID, RETIRO MADRID 667 5.4 0.3 10.4
2 2000-01-02 3195 MADRID, RETIRO MADRID 667 5 0.3 9.6
3 2000-01-03 3195 MADRID, RETIRO MADRID 667 3.5 0.1 6.9
4 2000-01-04 3195 MADRID, RETIRO MADRID 667 4.3 1.4 7.2
5 2000-01-05 3195 MADRID, RETIRO MADRID 667 0.6 -0.4 1.6
6 2000-01-06 3195 MADRID, RETIRO MADRID 667 3.8 -1.1 8.8
7 2000-01-07 3195 MADRID, RETIRO MADRID 667 6.2 0.6 11.7
8 2000-01-08 3195 MADRID, RETIRO MADRID 667 5.4 -0.1 11
9 2000-01-09 3195 MADRID, RETIRO MADRID 667 5.5 3 8
10 2000-01-10 3195 MADRID, RETIRO MADRID 667 4.8 1.8 7.8
# ℹ 16,304 more rows
Ejemplo real: AEMET again
Convertimos a tsibble los datos
Código
retiro_ts <- retiro |># index: variable temporal# key: si tuviéramos varias series a la vez (varias estaciones)# regular = TRUE: regular time interval (para que detecte# adecuadamente la periodicidad de la serie, en este caso [1D])as_tsibble(index = fecha, key =NULL, regular =TRUE)retiro_ts
# A tsibble: 16,314 x 8 [1D]
fecha id_station nombre provincia altitud tmed tmin tmax
<date> <dbl> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
1 1980-01-01 3195 MADRID, RETIRO MADRID 667 12.5 9.7 15.3
2 1980-01-02 3195 MADRID, RETIRO MADRID 667 8.4 5 11.8
3 1980-01-03 3195 MADRID, RETIRO MADRID 667 5.5 1 10
4 1980-01-04 3195 MADRID, RETIRO MADRID 667 5.5 1.2 9.8
5 1980-01-05 3195 MADRID, RETIRO MADRID 667 5.6 0.8 10.4
6 1980-01-06 3195 MADRID, RETIRO MADRID 667 2.7 0.6 4.8
7 1980-01-07 3195 MADRID, RETIRO MADRID 667 3 0.6 5.4
8 1980-01-08 3195 MADRID, RETIRO MADRID 667 2.3 -2 6.6
9 1980-01-09 3195 MADRID, RETIRO MADRID 667 3.8 0.8 6.8
10 1980-01-10 3195 MADRID, RETIRO MADRID 667 5.6 1.8 9.4
# ℹ 16,304 more rows
Ejemplo real: AEMET again
Preprocesamos los datos. ¿Hay huecos?
Código
retiro_ts |>count(year(fecha)) |>arrange(n)
Vemos que, amén de los datos que faltan (lógicamente) en 2024, hay datos que faltan en 2022, así que antes debemos rellenar los huecos: primero creando la fila (con valores vacíos) y luego rellenando la variable objetivo
Al margen de posibles outliers: no hay un patrón ni tendencia en la varianza
Ejemplo real: AEMET again
De hecho si realizamos el ajuste de un recta a la propia varianza agrupada vemos que no es significativo: no hay evidencias de que exista una tendencia en la varianza.
lm(data = resumen_var, formula = dispersion ~ year_month) |>summary()
Call:
lm(formula = dispersion ~ year_month, data = resumen_var)
Residuals:
Min 1Q Median 3Q Max
-0.16829 -0.08274 -0.02334 0.05620 0.87227
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.448e-01 1.414e-02 17.315 <2e-16 ***
year_month -1.808e-06 1.113e-06 -1.624 0.105
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.1214 on 534 degrees of freedom
Multiple R-squared: 0.004916, Adjusted R-squared: 0.003053
F-statistic: 2.638 on 1 and 534 DF, p-value: 0.1049
Ejemplo real: AEMET again
Recuerda que si tuviéramos heterocedasticidad tendríamos que aplicar algún tipo de transformación a los datos, conocidas como transformaciones Box-Cox (adaptadas por Bickel y Doksum (1981) para valores negativos):
¿Qué \(\lambda\) sería si no tuviéramos que hacer ninguna transformación (es decir, homocedástica)?
Ejemplo real: AEMET again
Para calcularlo vamos a usar features(variable, features = ...) del paquete {fabletools} (ya cargado en {fable}). En features = ... vamos a especificar la opción guerrero (un método para elegir el \(\lambda\) óptimo del paquete {feasts})
library(fable)library(feasts)retiro_ts |>features(tmed, features = guerrero)
# A tibble: 1 × 1
lambda_guerrero
<dbl>
1 0.922
Ese valor está muy cerca de \(1\) lo que apoya lo anterior: salvo outliers, no tenemos heterocedasticidad.
Ejemplo real: AEMET again
Determinar periodicidad\(s\)
pacf(retiro_ts$tmed)
No vemos a priori estacionalidad…porque recuerda que aquí tenemos datos diarios. ¿Y si aumentamos el número de retardos?
Ejemplo real: AEMET again
pacf(retiro_ts$tmed, lag.max =365*2)
Cuesta verlo pero…¿ves esa pequeña montañita en el 365?
Ejemplo real: AEMET again
Esto también lo podemos visualizar de manera conjunta con gg_tsdisplay()
retiro_ts |>gg_tsdisplay(y = tmed, lag =365)
Fíjate que además de mostrarnos un patrón en las autocorrelaciones cada 365 valores nos pinta la serie cada año (cada 365 valores corta una curva) y vemos que efectivamente comparten un patrón
Ejemplo real: AEMET again
Otra forma de comprobar que efectivamente es estacional \(s = 365\) es haciendo una diferenciación estacional de la serie: a cada valor le vamos a restar su valor \(s\) periodos previos.
La función que nos lo permite se llama difference() (del paquete {tsibble})
Ejemplo real: AEMET again
Si \(s\) está mal elegido seguirá presentando ese patrón (o uno aún más extraño)
retiro_ts |>gg_tsdisplay(difference(tmed, 123), lag =365)
Ejemplo real: AEMET again
Si \(s\) está bien elegido las autocorrelaciones debería de desplomarse
retiro_ts |>gg_tsdisplay(difference(tmed, 365), lag =365)
Ejemplo real: AEMET again
Si \(s\) está bien elegido las autocorrelaciones debería de desplomarse
retiro_ts |>gg_tsdisplay(difference(tmed, 365), lag =365,plot_type ="partial")
Ejemplo real: AEMET again
Una vez que tenemos que \(s = 365\) vamos a separar la muestra de momento solo en train-test, por ejemplo para evaluar cómo funcionan los modelos a 3 años vista (es decir, usando los últimos 3 años como test, en los que tenemos \(h = 974\) valores)
# code-fold: trueggplot(retiro_ts_monthly_train) +geom_line(aes(x = year_month, y = tmed)) +theme_minimal()
Ejemplo real: AEMET again
fit_retiro <-# entrenamos SOLO con train retiro_ts_monthly_train |>model("alisado_simple"=ETS(tmed ~error("A") +trend("N") +season("N")),"alisado_doble"=ETS(tmed ~error("A") +trend("A") +season("N")),"alisado_triple_ad"=ETS(tmed ~error("A") +trend("A") +season("A", period ="1 year")),"alisado_triple_m"=ETS(tmed ~error("M") +trend("A") +season("M", period ="1 year")))estimaciones <- fit_retiro |>augment()predicciones <- fit_retiro |>forecast(h =12*2+8)
Ejemplo real: AEMET again
Código
predicciones |>autoplot(retiro_ts |>index_by("year_month"=yearmonth(fecha)) |>summarise("tmed"=mean(tmed)), level =NULL) +geom_line(data = estimaciones, aes(x = year_month, y = .fitted, color = .model)) +geom_line(data = predicciones, aes(x = year_month, y = .mean, color = .model)) +scale_color_manual(values = ggthemes::colorblind_pal()(5)[-1]) +theme_minimal()
Ejemplo real: AEMET again
La evaluación numérica la haremos con accuracy(): el triple aditivo funciona ligeramente mejor en train pero ligeramente peor en test. Dado que funcionan similar y el aditivo es más simple, por principio de parsimonia, nos quedamos con el triple aditivo.
accuracy(fit_retiro)
# A tibble: 4 × 10
.model .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
<chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 alisado_simple Training 0.00644 4.08 3.43 -4.61 25.9 2.13 2.07 0.583
2 alisado_doble Training -0.0322 4.08 3.44 -5.01 26.1 2.14 2.08 0.581
3 alisado_triple_ad Training -0.0190 1.42 1.13 -1.33 9.24 0.704 0.722 0.182
4 alisado_triple_m Training -0.0264 1.46 1.17 -1.73 9.39 0.726 0.742 0.201
accuracy(predicciones, retiro_ts_monthly_test)
# A tibble: 4 × 10
.model .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
<chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 alisado_doble Test 7.46 10.5 8.16 30.3 40.1 NaN NaN 0.788
2 alisado_simple Test 7.93 10.9 8.48 33.6 41.5 NaN NaN 0.784
3 alisado_triple_ad Test 0.872 1.58 1.33 4.46 8.70 NaN NaN -0.0169
4 alisado_triple_m Test 0.654 1.40 1.17 3.88 8.12 NaN NaN -0.0341
Ejemplo real: AEMET again
El diagnóstico de los residuos (de momento solo visual, ya veremos con lso ARIMA como testar que los residuos son como los queremos, normales, etc) lo haremos con gg_tsresiduals() sobre el modelo “ganador”
retiro_ts_monthly_train |>model("alisado_triple_ad"=ETS(tmed ~error("A") +trend("A") +season("A", period ="1 year"))) |>gg_tsresiduals()
Casos reales
Practica para a ejecutar todo el proceso con estas series
Realizar todos los métodos de alisado conocidos así como otros modelos que sepas de fable. Pensar ANTES de ver su resultado cómo crees que van a funcionar y por qué
Evaluarlos en train-test
Usa el mejor de los métodos posibles en cada datasets y ahora ajustarlo pero con distintos valores de los parámetros \(\alpha\), \(\beta\) y \(\gamma\) y decidir los 3 mejores modelos haciendo uso de la validación cruzada
Casos reales: acciones de google
Vamos a ilustrar el inicio de cómo trabajar con la serie temporal que captura el precio de cierre de las acciones de Google (variable Close del paquete {tsibbledata::gafa_stock})
datos <- tsibbledata::gafa_stock |>filter(Symbol =="GOOG")datos
Si te fijas ningún año tiene 365 valores ya que los mercados bursátiles cierran los fines de semana.
Acciones de google
Para rellenar los huecos vamos a usar primero fill_gaps() rellenando las fechas que no tenemos (con valores vacíos) y luego fill(..., .direction = "down") para rellenar la variable objetivo con los valores previos (el sábadon y domingo será el valor de cierre del viernes).
datos <- datos |>fill_gaps() |>fill(Close, .direction ="down")datos |>count(year(Date)) |>arrange(n)
Si te fijas lo primero que estamos pidiendo es fill_gaps(): rellena huecos entre fechas. Por ejemplo, si los datos son diarios, debe buscar huecos entre días, pero para eso lo primero que tiene que saber es que la serie es diaria y saber cómo debe rellenar huecos.
Si te fijas en la cabecera del tsibble aparece esto:
# A tsibble: 1,258 x 8 [!]# Key: Symbol [1]
En Key: no debería figurar nada (ya que no tenemos distintas series al haber filtrado una sola) pero sobre todo…en [!]debería figurar[1D]` ¡pero no figura nada!
Acciones de google
Para arreglarlo vamos a redefinir el objeto de serie temporal indicándole dentro de as_tsibble() que
index variable temporal
key: si tuviéramos varias series a la vez (en este caso NULL)
regular = TRUE: regular time interval (para que detecte adecuadamente la periodicidad de la serie, en este caso [1D])
datos <- datos |>as_tsibble(index = Date, key =NULL, regular =TRUE)datos
ggplot(datos, aes(x = Date, y = Close)) +geom_line() +geom_smooth(method ="lm", se =FALSE) +theme_minimal()
Parece que sí, que existe una tendencia positiva
Acciones de google
ajuste_lineal <-lm(data = datos, formula = Close ~ Date)ajuste_lineal |>summary()
Call:
lm(formula = Close ~ Date, data = datos)
Residuals:
Min 1Q Median 3Q Max
-162.04 -53.29 -14.68 48.75 188.99
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -5.796e+03 5.134e+01 -112.9 <2e-16 ***
Date 3.876e-01 3.021e-03 128.3 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 68 on 1823 degrees of freedom
Multiple R-squared: 0.9003, Adjusted R-squared: 0.9002
F-statistic: 1.646e+04 on 1 and 1823 DF, p-value: < 2.2e-16
Existe una muy fuerte tendencia positiva. Si te fijas \(R^2 = 0.9\) y los contrastes son todos apoyando rechazar la hipótesis nula así parece que tiene sentido que la tendencia sea lineal (es decir, ADITIVA)
Acciones de google
Una forma de chequear si es aditiva o no es ver los errores tras ajustar tendencia: no se observa ningún patrón evidente que falte por modelizar (si fuese multiplicativa se vería una tendencia todavía)
Código
ggplot(tibble("fecha"= datos$Date, "res"= ajuste_lineal$residuals),aes(x= fecha, y = res)) +geom_line() +theme_minimal()
Acciones de google
Determinar homocedasticidad/heterocedasticidad y cómo evoluciona su varianza. En caso necesario aplicar transformación.
Una primera opción es visualización de la dispersión por periodos (por ejemplo, por meses)
Código
datos <- datos |>mutate("year_month"=yearmonth(Date))resumen_var <- datos |>index_by(year_month) |>summarise("dispersion"=sd(Close) /mean(Close))ggplot(resumen_var, aes(x = year_month, y = dispersion)) +geom_line() +geom_smooth(method ="lm", se =FALSE) +theme_minimal()
Al margen de posibles outliers: no hay un patrón ni tendencia en la varianza
Acciones de google
De hecho si realizamos el ajuste de un recta a la propia varianza agrupada vemos que no es significativo: no hay evidencias de que exista una varianza que cambien con el tiempo.
lm(data = resumen_var, formula = dispersion ~ year_month) |>summary()
Call:
lm(formula = dispersion ~ year_month, data = resumen_var)
Residuals:
Min 1Q Median 3Q Max
-0.015438 -0.008556 -0.003433 0.006435 0.078257
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.982e-02 5.878e-02 0.337 0.737
year_month 1.864e-07 3.462e-06 0.054 0.957
Residual standard error: 0.01414 on 58 degrees of freedom
Multiple R-squared: 4.995e-05, Adjusted R-squared: -0.01719
F-statistic: 0.002897 on 1 and 58 DF, p-value: 0.9573
Acciones de google
Determinar periodicidad\(s\)
pacf(datos$Close)
No vemos a priori estacionalidad…porque recuerda que aquí tenemos datos diarios. ¿Y si aumentamos el número de retardos?
Acciones de google
pacf(datos$Close, lag.max =365*3)
No se aprecia patrón estacional
Acciones de google
Esto también lo podemos visualizar de manera conjunta con gg_tsdisplay()
datos |>gg_tsdisplay(y = Close, lag =365)
Fíjate que además de no mostrarnos un patrón en las autocorrelaciones al pedirle que la serie cada año (cada 365 valores corta una curva) y vemos que efectivamente no comparten ningún patrón anual
Acciones de google
Otra forma de comprobar que efectivamente no es estacionales haciendo una diferenciación estacional de la serie: a cada valor le vamos a restar su valor \(s\) periodos previos.
La función que nos lo permite se llama difference() (del paquete {tsibble})
Acciones de google
Si \(s\) está mal elegido seguirá presentando ese patrón (o uno aún más extraño)
datos |>gg_tsdisplay(difference(Close, 30), lag_max =365*3)
Acciones de google
Si \(s\) está mal elegido seguirá presentando ese patrón (o uno aún más extraño)
datos |>gg_tsdisplay(difference(Close, 180), lag_max =365*3)
Acciones de google
Si \(s\) está mal elegido seguirá presentando ese patrón (o uno aún más extraño)
datos |>gg_tsdisplay(difference(Close, 365), lag_max =365*3)
Acciones de google
A partir de aquí todo tuyo. Piensa:
Dividir train/test
¿Qué modelos podrían tener sentido y por qué? ¿Algún alisado? ¿Alguno de los demás vistos en {fable}?
¿Cómo evaluar los modelos?
Clases 16: procesos estacionarios
Modelos ARIMA
En los años 50 y 60 las matemáticas y la probabilidad entraron de lleno (aún más) en el campo de las series temporales, introduciendo el concepto de proceso estocástico.
La idea es que la serie temporal \(X_t\) que observamos en un instante temporal \(t\) no es más que la realización de una variable aleatoria definida en dicho instante: el valor que observamos es un valor particular de los infinitos valores que podría haber tomado.
Veamos un ejemplo simulado
Proceso estocástico
Vamos a simular la siguiente serie temporal (lo que se conoce como paseo aleatorio: el valor anterior más una perturbación aleatoria)
La serie temporal que observamos no es más que una de las infinitas realizaciones que podría haber tomado
Proceso estocástico
Proceso estocástico: un conjunto de variables aleatorias \(\left\lbrace X_t \right\rbrace_{t \in T}\) definidas sobre el mismo espacio de probabilidades \(\left(\Omega, \mathcal{A}, P \right)\). El conjunto de índices \(T\) suele ser un espacio temporal continuo.
En el gráfico: todo en su conjunto
Proceso estocástico
Dado que \(\left\lbrace X_t \right\rbrace_{t \in T}\) es un conjunto de variables aleatorias, dependientes también de un espacio temporal, tendremos que
donde \(S\) es el espacio de estados (valores posibles que puede tomar).
Proceso estocástico
Realización o trayectoria de un proceso estocástico: si fijamos \(\omega \in \Omega\) (una aleatoriedad), tenemos
\[\begin{eqnarray}X\left(\cdot, \omega \right): ~T & \to & S \nonumber \\ t & \to & X_t (\omega) ~ \text{ una realización}\nonumber\end{eqnarray}\] En el gráfico: si interpretamos una sola curva y vemos su evolución en el tiempo
Proceso estocástico
Distribucion de probabilidad a tiempo t: si fijamos \(t \in T\)
Dado que para cada instante \(t\) tenemos una distribución de probabilidad, podremos caracterizar la serie con algunos de sus parámetros poblacionales:
Media del proceso: fijado un \(t\) tenemos que \(\mu_t = E \left[X_t \right]\) (si es constante en el tiempo –> sin tendencia)
Código
ggplot(serie_total) +geom_line(aes(x = t, y = X_t, color = id_serie),linewidth =0.5) +geom_line(data = serie_total |>summarise("mu_t"=mean(X_t), .by = t),aes(x = t, y = mu_t), linewidth =1.1) +guides(color ="none") +theme_minimal()
Dado que para cada instante \(t\) tenemos una distribución de probabilidad, podremos caracterizar la serie con algunos de sus parámetros poblacionales:
Varianza del proceso: fijado un \(t\) tenemos que \(\sigma_{t}^{2} = Var \left[X_t \right]\) (si es constante en el tiempo –> homocedástica)
Código
ggplot(serie_total) +geom_line(aes(x = t, y = X_t, color = id_serie),linewidth =0.5) +geom_line(data = serie_total |>summarise("var_t"=sd(X_t), .by = t),aes(x = t, y = var_t), linewidth =1.1) +guides(color ="none") +theme_minimal()
Proceso estocástico
Mientras que en los modelos de descomposición y en los alisados exponencial simplemente realizamos un ajuste de una curva observada (sin tener en cuenta su distribución estocástica subyacente), en el momento en el que consideramos modelos ARIMA (la serie como una realización de un conjunto de variables aleatorias correladas en el tiempo) será importante su estructura de autocorrelación: las diferentes distribuciones del proceso estocástico dependen entre sí para distintos \(t\).
Dados \(t_1\) y \(t_2\), llamaremos función de autocovarianzas a
Todo lo anterior es sobre la base de que yo pudiese generar una colección de curvas de manera que puedo observarlas de manera longitudinal (tiempo) o transversal (distribución probabilística a un tiempo dado).
Proceso estacionario
El problema es que normalmente es inviable obtener varias realizaciones: si tengo una serie temporal de las temperaturas a lo largo del año, esas temperaturas podrían haber sido otras pero no puedo volver al pasado y generar otras nuevas.
Así la única forma de estimar las características transversales (media, varianza, etc) de una serie haciendo uso de su evolución longitudinal es suponer que esas propiedades transversales (distribución) en cada instante temporal son estables (no varían)
Proceso estacionario
Por lo tanto el primer requisito que necesitaremos para la introducción de modelos ARIMA o modelos probablísticos será el concepto de estacionariedad: algo de la serie que no varíe con el tiempo
Diremos que \(\left\lbrace X_{t} \right\rbrace\) es un proceso estacionario (en el tiempo) cuando sus propiedades no dependen del instante \(t\) en el que las medimos. Dicho de manera informal, una serie será estacionaria cuando mire donde mire veo lo mismo, cuyo patrón no puede ser predicho
Desde un punto de vista matemático es estrictamente estacionario si para cualquier conjunto de instantes temporales y retardo \(h > 0\), tenemos que \(\left(X_{t_{1}}, \ldots, X_{t_{n}} \right)\) y \(\left(X_{t_{1}+h}, \ldots, X_{t_{n}+h} \right)\) tienen la misma distribución conjunta: son indistinguibles.
Proceso estacionario
Dado que está condición es demasiado estricta (nos obligaría a saber su distribución, si es una normal, una chi-cuadrado, etc), trabajaremos con procesos débilmente estacionarios
Media constante: \(\mu_{t} = \mu_{t+h} = cte\) para cualquier instante \(t\) y retardo \(h\).
Varianza constante: \(\sigma_{t}^{2} = \sigma_{t+h}^{2}\) para cualquier instante \(t\) y retardo \(h\).
Autocovarianzas independiente de t: \(\gamma_{s, t} = \gamma_{s+h, t+h} = \gamma_h\) la dependencia entre dos instantes solo depende de su distancia, al margen de en qué punto lo medimos.
Proceso estacionario
De lo anterior se derivan algunas propiedades sobre las autocorrelaciones (las barras que nos dibuja el acf)
Solo nos importan los \(\gamma_h = \gamma_{t, t+h}\) para cada \(h=0, 1, 2, ...\), tal que \(\gamma_h = \gamma_{-h}\)
De la misma forma \(\rho_h = \rho_{t, t+h} = \frac{\gamma_{t, t+h}}{\sigma_{t} \sigma_{t+h}} = \frac{\gamma_{t, t+h}}{\sigma^2} = \frac{\gamma_{h}}{\gamma_0}\) tal que \(\rho_{h} = \rho_{-h}\) y \(\rho_0 = 1\).
Dado que \(\left| \rho_{h} \right| \leq 1\) entonces \(\left|\gamma_h \right| = \left| \rho_h \right| * \gamma_0 \leq \gamma_0\) para todo \(h\) (por eso decrecen las barras del acf)
Si \(\left\lbrace X_t \right\rbrace\) es un proceso estacionario entonces su diferencia \(X_{t} - X_{t-1}\) también lo es.
Ruido blanco
Uno de los procesos estocásticos más importantes es el conocido como ruido blanco
Media 0: \(\mu_{t} = \mu_{t+h} = 0\) para cualquier instante \(t\) y retardo \(h\).
Varianza constante: \(\sigma_{t}^{2} = \sigma_{t+h}^{2} = \sigma \neq 0\) para cualquier instante \(t\) y retardo \(h\).
# A tsibble: 1,825 x 8 [1D]
# Key: Symbol [1]
Symbol Date Open High Low Close Adj_Close Volume
<chr> <date> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 GOOG 2014-01-02 554. 555. 551. 553. 553. 3666400
2 GOOG 2014-01-03 554. 555. 549. 549. 549. 3355000
3 GOOG 2014-01-04 NA NA NA 549. NA NA
4 GOOG 2014-01-05 NA NA NA 549. NA NA
5 GOOG 2014-01-06 553. 556. 550. 555. 555. 3561600
6 GOOG 2014-01-07 559. 566. 557. 566. 566. 5138400
7 GOOG 2014-01-08 569. 570. 563. 567. 567. 4514100
8 GOOG 2014-01-09 568. 568. 559. 561. 561. 4196000
9 GOOG 2014-01-10 566. 566. 557. 561. 561. 4314700
10 GOOG 2014-01-11 NA NA NA 561. NA NA
# ℹ 1,815 more rows
Diagnóstico de residuos
Ya vimos en su momento que tiene tendencia, por lo que no es un proceso estacionario ya que su distribución es distinta en distintos instantes temporales
ggplot(google) +geom_line(aes(x = Date, y = Close)) +theme_minimal()
Diagnóstico de residuos
Si pintamos las estimaciones muestrales (insesgadsa) de las autocorrelaciones tenemos que no decrecen, es decir, no es ruido blanco (queda algo por modelizar)
# alternativa a acf(google$Close) en forma gpglotgoogle |>ACF(Close, lag_max =100) |>autoplot() +labs(title ="Google closing stock price",y ="Estimación muestral de autocorrelaciones",x ="Retardos") +theme_minimal()
Diagnóstico de residuos
Para tener una cuantificación inferencial más rigurosa sobre si es o no ruido blanco podemos realizar el test de Ljung-Box que nos contrasta si un grupo cualquiera de autocorrelaciones son diferentes de cero (se conoce como contrastes Portmanteau a aquellos que en lugar de probar la aleatoriedad en cada retardo distinto, prueba la aleatoriedad “en general”)
Para ello usaremos features() con la opción ljung_box (ya la usamos para las transformaciones Box-Cox)
google |>features(Close, ljung_box, lag =10)
# A tibble: 1 × 3
Symbol lb_stat lb_pvalue
<chr> <dbl> <dbl>
1 GOOG 17980. 0
Según el contraste hay evidencia suficientes para rechazar que sean cero -> no es ruido blanco
Diagnóstico de residuos
¿Pero qué pasa si yo diferencio la serie (resto a cada valor su instante anterior)? Si ahora te diese dos fotos de la serie, ¿sabrías ahora distinguir su instante temporal?
ggplot(google) +geom_line(aes(x = Date, y =difference(Close))) +theme_minimal()
Diagnóstico de residuos
¡Todas se desploman!
# alternativa a acf(google$Close) en forma gpglotgoogle |>ACF(difference(Close), lag_max =100) |>autoplot() +labs(title ="Google closing stock price",subtitle ="Serie diferenciada",y ="Estimación muestral de autocorrelaciones",x ="Retardos") +theme_minimal()
Diagnóstico de residuos
google |>mutate(diff_close =difference(Close)) |>features(diff_close, ljung_box, lag =10)
# A tibble: 1 × 3
Symbol lb_stat lb_pvalue
<chr> <dbl> <dbl>
1 GOOG 31.7 0.000442
Simplemente aplicando una diferenciación nuestro proceso se ha convertido en ruido blanco
Clases 17: procesos MA
Correlogramas
Como hemos comentado los correlogramas van a ser ahora fundamentales ya que os da información de la estructura del proceso estocástico subyacente.
google |>ACF(Close, lag_max =100) |>autoplot() +labs(y ="Estimación muestral de autocorrelaciones",x ="Retardos") +theme_minimal()
google |>ACF(difference(Close), lag_max =100) |>autoplot() +labs(y ="Estimación muestral de autocorrelaciones",x ="Retardos") +theme_minimal()
Correlogramas
En lineas generales diremos que
Series con tendencia presentarán una caída muy lenta de las autocorrelaciones (la tendencia implica que el valor \(t+1\) depende del \(t, t-1, t-2, ...\))
Series con estacionalidad presentarán un patrón cíclico de caída a lo largo del tiempo.
Correlogramas
Código
datos <- timeSeriesDataSets::a10_ts |>as_tsibble(regular =TRUE)ggplot(datos) +geom_line(aes(x = index, y = value)) +theme_minimal()
Código
datos |>ACF(value, lag_max =100) |>autoplot() +labs(y ="Estimación muestral de autocorrelaciones",x ="Retardos") +theme_minimal()
Correlogramas
Código
datos <- timeSeriesDataSets::beer_ts |>as_tsibble(regular =TRUE)ggplot(datos) +geom_line(aes(x = index, y = value)) +theme_minimal()
Código
datos |>ACF(value, lag_max =100) |>autoplot() +labs(y ="Estimación muestral de autocorrelaciones",x ="Retardos") +theme_minimal()
Correlogramas
Código
datos <- timeSeriesDataSets::elec_ts |>as_tsibble(regular =TRUE)ggplot(datos) +geom_line(aes(x = index, y = value)) +theme_minimal()
Código
datos |>ACF(value, lag_max =100) |>autoplot() +labs(y ="Estimación muestral de autocorrelaciones",x ="Retardos") +theme_minimal()
Procesos lineales
En las próximas clases trabajaremos siempre suponiendo que nuestros procesos ya son estacionarios. Veremos más adelante qué hacer cuando no son estacionarios en media o en varianza
Un resultado importante (se conoce como Teorema de descomposición de Wald) es que cualquier proceso estacionario puede transformado en un proceso lineal definido como
\[X_t = \mu + \sum_{j = -\infty}^{j=\infty} \Psi_j \varepsilon_{t-j} \quad \sum_{j = -\infty}^{j=\infty} \left| \Psi_j \right| < \infty, \quad \left\lbrace \varepsilon_{t} \right\rbrace \text{ ruido blanco}\] Un proceso lineal es por tanto la combinación lineal de instantes de ruido blanco tal que los coeficientes sean absolutamente sumables.
Procesos lineales
\[X_t = \mu + \sum_{j = -\infty}^{j=\infty} \Psi_j \varepsilon_{t-j}, \quad \sum_{j = -\infty}^{j=\infty} \left| \Psi_j \right| < \infty, \quad \left\lbrace \varepsilon_{t} \right\rbrace \text{ ruido blanco}\] Si te fijas hay una constante inicial \(\mu\) que representa la media (a largo plazo) del proceso ya que, dado que la esperanza es lineal y \(\left\lbrace \varepsilon_{t} \right\rbrace\) es ruido blanco, entonces
donde \(q\) será el orden (lo llamaremos procesos de medias moviles de orden q o MA(q)). Si te fijas es un caso particular de proceso lineal donde todos los coeficientes son 0 salvo \(\left(\Psi_0 = 1, \Psi_1 = -\theta_1,\ldots, \Psi_q = -\theta_q \right)\).
Los procesos \(MA(q)\) son procesos cuyo futuro se describe promediando los errores del pasado y son tremendamente útiles para modelizar fenómenos influenciados por sucesos que producen un efecto inmediato de corta duracción (por ejemplo, variables económicas).
Fíjate que el signo de \(\rho_h\) será el contrario al signo de \(\theta_1\) y que, a partir de \(q\), las autocorrelaciones son 0 (se dice que es q-correlacionado)
MA(1)
💻 Diseña una función para simular \(n\) trayectorias de un \(MA(1)\) en función de \(\left(n,~\sigma,~\theta_1, \varepsilon_1\right)\)
💻 Aplica dicha función para \(n = 300\), \(\sigma = 3\), \(\varepsilon_1 = 0\) y \(\theta_1 = (-0.3, -1.5)\)
Código
set.seed(1234567)ggplot(MA_1_simul(300, 3, -0.3, 0)) +geom_line(aes(x = t, y = X_t)) +labs(title ="MA(1) con theta_1 = -0.3 y sigma = 3") +theme_minimal()
Código
set.seed(1234567)ggplot(MA_1_simul(300, 3, -1.5, 0)) +geom_line(aes(x = t, y = X_t)) +labs(title ="MA(1) con theta_1 = -1.5 y sigma = 3") +theme_minimal()
MA(1)
💻 Calcula las ACF de ambas series
Código
set.seed(1234567)MA_1_simul(300, 3, -0.3, 0) |>ACF(X_t) |>autoplot() +labs(title ="MA(1) con theta_1 = -0.3 y sigma = 3") +theme_minimal()
Código
set.seed(1234567)MA_1_simul(300, 3, -1.5, 0) |>ACF(X_t) |>autoplot() +labs(title ="MA(1) con theta_1 = -1.5 y sigma = 3") +theme_minimal()
¿Te has fijado en que las correlaciones caen a partir de \(q\) (q-correlacionado poblacionalmente aunque a nivel muestral puedan caer algunas todavía fuera de la banda de significación)?
MA(1)
Si te fijas, dado que \(X_t = \varepsilon_t - \theta_1 \varepsilon_{t-1}\), tenemos entonces que \(\varepsilon_t = X_t + \theta_1 \varepsilon_{t-1}\) por lo que sustituyendo tenemos que
El futuro de la serie acaba siendo modelizado con el último error y una ponderación del pasado, y dado que parece razonable asumir que los valores más cercanos influyan más que los valores más lejanos los pesos \(\left| \theta_{1}\right|^{j}\) deben disminuir con el retardo: necesitamos que \(\left| \theta_{1}\right| < 1\)
Condición de invertibilidad
Esta condición (veremos porque se llama de invertibilidad) será fundamental en los procesos MA ya que, aunque la función de autocorrelación será importante, tenemos un problema de identificación: dada una función de autocorrelación no existe un único proceso asociado
Hagamos un ejercicio rápido con dos \(MA(1)\) distintos
Usaremos la función de autocorrelaciones (ACF) para determinar si se trata de un proceso MA(1): si a partir de la primera las correlaciones caen drásticamente (q-correlacionado), es probable que el proceso sea un \(MA(1)\)
MA(1)
Recuerda que estamos visualizando las estimaciones muestrales por lo que cuanto más aumente \(\theta_1\) más lenta será la caída (si \(\theta_1 > 0\) las autocorrelaciones irán cambiando de signo).
por lo que los coeficientes de \(\Pi (B)\) se determinarán tal que \(\Pi \left(B \right) \Theta_q (B) = 1\).
MA(q): invertibilidad
Si \(\Pi \left(B \right) \Theta_q (B) = 1\), entonces \(\Pi (B) = \Theta_{q}^{-1} (B)\).
Es decir, dicho polinomio \(\Theta_{q}(B)\) tiene que ser por tanto invertible. Se puede demostrar como un proceso MA(q) es invertible si y solo si las raíces de
son (en módulo) mayores que la unidad (\(\left| z \right| > 1\)).
MA(q) invertibilidad
El paquete {pracma} nos permite determinar las raíces de un polinomio con la función roots() introduciendo los coeficientes del polinomio de más grado a menos.
Por ejemplo si queremos calcular las raíces de \(-z^3 +2 z^2 +1.5* z - 2\)
pracma::roots(c(-1, 2, 1.5, -2))
[1] 2.2728039 -1.0843343 0.8115305
Si queremos definir un \(MA(q)\), por ejemplo, \(X_t = \varepsilon_t - 0.2 \varepsilon_{t-1} + 0.5 \varepsilon_{t-2}\), y ver si cumple las condiciones de invertibilidad basta con meter los paramétros de antiguo a reciente y cambiando el signo
pracma::roots(c(-0.5, 0.2, 1))
[1] 1.628286 -1.228286
MA(q)
💻 Diseña una función para simular \(n\) trayectorias de un \(MA(q)\) en función de \(q\), \(n\), \(\sigma\), \(\left(\theta_1,~\theta_2, \ldots, ~\theta_q\right)\) y \(\left(\varepsilon_1, \ldots, \varepsilon_q\right)\)
Código
MA_q_simul <-function(q, n, sigma, theta, eps_q) { X_t <- eps_q eps <-c(eps_q, rnorm(n = n - q, mean =0, sd = sigma))for (i in (q +1):n) { X_t[i] <- eps[i] -sum(theta * eps[(i -1):(i - q)]) } ts <-tibble("t"=1:n, "X_t"= X_t) |>as_tsibble(index = t)return(ts)}
MA(q)
💻 Aplica dicha función para simular un \(MA(1)\) con \(n = 10000\), \(\sigma = 1\), \(\varepsilon_1 = 0\) y \(\theta_1 = -0.7\), y visualiza el proceso
Código
set.seed(1234567)ggplot(MA_q_simul(1, 10000, 1, theta =-0.7, eps =rep(0, 1))) +geom_line(aes(x = t, y = X_t)) +labs(title ="MA(1) con theta_1 = -0.7 y sigma = 1") +theme_minimal()
MA en fable
Parece impredecible pero…¡no!
Vamos a realizar el proceso completo con el MA(1) anteriormente generado
Rechazamos la hipótesis nula: no es ruido blanco –> seguimos
MA en fable
¿Es un proceso estacionario?
Para ello aplicaremos el conocido como Kwiatkowski-Phillips-Schmidt-Shin (KPSS) test (Kwiatkowski et al., 1992)] incluyendo unitroot_kpss en features(). En dicho test su hipótesis nula es que el proceso es estacionario (cualquier p-valor por encima de 0.1 te lo trunca a 0.1)
No rechazamos la hipótesis nula de estacionariedad –> seguimos (no hay que aplicar transformaciones).
MA en fable
¿Cómo son las autocorrelaciones (ACF) para determinar el \(q\) adecuado?
datos |>ACF(X_t, lag_max =30) |>autoplot()
La primera autocorrelación que sale de la banda (y además muy evidente) es la primera –> ¿será \(q=1\)? (recuerda que el número de \(\rho_h\) no nulas nos caracteriza un MA)
MA en fable
Realizamos el ajuste ARIMA. De momento solo estamos en el final, la parte \(MA(q)\) pero el proceso final se denota como \(ARIMA(p, d, q)\). Para indicarle los índices usamos dentro de ARIMA(var_objetivo ~ pdq() la función pdq() (de momento los demás 0)
fit <- datos |>model("MA_1"=ARIMA(X_t ~pdq(0, 0, 1)))estimaciones <- fit |>augment()estimaciones
No se rechaza la hipótesis nula de ruido blanco –> hemos terminado. Podemos ver la estimación de \(\theta_1\) haciendo report() del modelo (los parámetros tienen el signo cambiado a nuestra ecuación)
Importante: un proceso MA(1) podría ser ajustado también con un MA(q) CON \(q > 1\) ya que basta con que poner el resto de coeficientes a 0 (fíjate en la estimación que nos devuelve el report()) –> principio de parsimonia: el modelo más sencillo
fit <- datos |>model("MA_4"=ARIMA(X_t ~pdq(0, 0, 4)))estimaciones <- fit |>augment()estimaciones |>features(.resid, ljung_box)
💻 Aplica dicha función de simulación para simular un \(MA(4)\) con \(n = 10000\), \(\sigma = 1\), \(\left(\varepsilon_1, \varepsilon_2, \varepsilon_3, \varepsilon_4 \right) = 0\) y \(\left(\theta_1, ~ \theta_2,~ \theta_3,~ \theta_4 \right) = (-0.1, 0.2, 0.1, 0.3)\) y comprueba que se cumple las condiciones de invertibilidad
No se rechaza la hipótesis nula de ruido blanco en ningún caso –> nos quedamos con el modelo más sencillo (MA(4)). Podemos ver la estimación de \(\theta_1\) haciendo report() del modelo (de hecho si lo hacemos con fit_5 en realidad es un MA(5) con el último muy pequeño).
Clases 19: procesos AR
Modelos AR
Hasta ahora hemos visto solo procesos cuya dependencia se contruye promediando los errores pasados. ¿Y si hacemos «lo mismo» pero promediando el pasado de la serie?
Dado que la idea es hacer una regresión con su propio pasado, a estos modelos se conoce como procesos autorregresivos (AR)
💻 Diseña una función para simular \(n\) trayectorias de un \(AR(p)\) en función de \(p\), \(n\), \(\sigma\), \(\left(\phi_1,~\phi_2, \ldots, ~\phi_p\right)\) y \(\left(X_1, \ldots, X_p\right)\)
Código
AR_p_simul <-function(p, n, sigma, phi, X) { X_t <- Xfor (i in (p +1):n) { X_t[i] <-sum(phi * X_t[(i -1):(i-p)]) +rnorm(1, mean =0, sd = sigma) } ts <-tibble("t"=1:n, "X_t"= X_t) |>as_tsibble(index = t)return(ts)}
AR(p)
💻 Aplica dicha función para simular un \(AR(2)\) con \(n = 300\), \(\sigma = 3\), \(\left(X_1, X_2 \right) = 0\) y \(\left(\phi_1, ~ \phi_2 \right) = (-0.3, 0.5)\)
Código
set.seed(1234567)AR_p_simul(p =2, n =300, sigma =3, phi =c(-0.3, 0.5), X =c(0, 0)) |>ggplot() +geom_line(aes(x = t, y = X_t)) +labs(title ="AR(2) con phi_1 = -0.3, phi_2 = 0.5 y sigma = 3") +theme_minimal()
AR(p)
💻 Aplica la función MA_q para simular un \(MA(2)\) con \(n = 10000\), \(\sigma = 3\), \(\left(\varepsilon_1, \varepsilon_2 \right) = 0\) y \(\left(\theta_1, ~ \theta_2 \right) = (-0.3, 0.5)\). Compara ambos procesos (mismos parámetros)
Código
set.seed(1234567)MA_q_simul(2, 10000, 3, theta =c(-0.3, 0.5), eps =rep(0, 2)) |>ggplot() +geom_line(aes(x = t, y = X_t)) +labs(title ="MA(2) con theta_1 = -0.3, theta_2 = 0.5 y sigma = 3") +theme_minimal()
Código
set.seed(1234567)AR_p_simul(p =2, n =10000, sigma =3, phi =c(-0.3, 0.5), X =c(0, 0)) |>ggplot() +geom_line(aes(x = t, y = X_t)) +labs(title ="AR(2) con phi_1 = -0.3, phi_2 = 0.5 y sigma = 3") +theme_minimal()
Ambos parecen algo impredible y aleatorio y ¡no es cierto!, simplemente son estacionarios (pero una parte de su comportamiento podemos predecirlo).
MA/AR vs ruido blanco
De hecho en ambos casos al realizar el contraste de ruido blanco de Ljung-Box obtenemos p-valores menores que 0.05 ==> rechazamos la hipótesis nula de ruido blanco
lo cual implica que \(\left| \rho_1 \right| < 1\) (ya que las correlaciones deberían ser menores que 1 en valor absoluto y decrecientes).
Fíjate que ahora las autocorrelaciones \(\rho_h\) NO decaen a 0
AR(1)
set.seed(1234567)MA_q_simul(q =1, n =10000, sigma =1, theta =c(0.75), eps =c(0)) |>ACF(X_t, lag_max =30) |>autoplot() +labs(title ="MA(1) con theta_1 = 0.75 y sigma = 1") +theme_minimal()
set.seed(1234567)AR_p_simul(p =1, n =10000, sigma =1, phi =c(0.75), X =c(0)) |>ACF(X_t, lag_max =30) |>autoplot() +labs(title ="AR(1) con phi_1 = 0.75 y sigma = 1") +theme_minimal()
¡La función ACF nunca va a anularse por completo!, de hecho tardará mucho en caer dentro de la banda de incorrelación.
AR(1)
Fíjate además que dado que \(\rho_h = \phi_{1}^{h}\), decrecen de manera exponencial con razón \(\phi_1\) (si aumentamos \(\phi_1\) tarda mucho más en caer)
set.seed(1234567)AR_p_simul(p =1, n =10000, sigma =1, phi =c(0.9), X =c(0)) |>ACF(X_t, lag_max =30) |>autoplot()
AR(p)
En general \(X_t = \phi_1 X_{t-1} + \phi_2 X_{t-2} + \ldots + \phi_p X_{t-p} + \varepsilon_t\)
Diremos que un proceso \(AR(p)\), definido como \(\Phi_p (B)X_t = \varepsilon_t\), es un proceso estacionario (o causal) si admite la representación \(MA(\infty)\) anterior, tal que
Al igual que antes, dado que \(\Phi_p X_t = \varepsilon_t\), los coeficientes de \(\Psi (B)\) se determinarán tal que $ _p (B ) = 1$.
AR(p): invertibilidad
Si \(\Phi_p (B) \Psi \left(B \right) = 1\), entonces \(\Psi (B) = \Phi_{p}^{-1} (B)\).
Es decir, dicho polinomio \(\Phi_{p}(B)\)tiene que ser de nuevo invertible. Se puede demostrar como un proceso AR(p) es invertible si y solo si las raíces de
Autocovarianzas: multiplicando la ecuación por \(X_{t-h}\) tenemos que
\[\begin{eqnarray}\gamma_h &=& E \left[X_{t-h}X_t \right]= E \left[X_{t-h} \left(\phi_1 X_{t-1} + \ldots + \phi_p X_{t-p} + \varepsilon_{t} \right) \right] \nonumber \\ &=& \phi_1 E \left[X_{t-h} X_{t-1} \right] + \ldots + \phi_p E \left[X_{t-h} X_{t-p} \right] + E \left[X_{t-h} \varepsilon_{t} \right]\ \nonumber \\ &=& \phi_1 \gamma_{h-1} + \ldots + \phi_p \gamma_{h-p} \quad h \geq 0 \end{eqnarray}\]
AR(p)
Autocorrelaciones:
\[\begin{eqnarray}\rho_h &=& \frac{\gamma_h}{\gamma_0} = \phi_1 \frac{\gamma_{h-1}}{\gamma_0} + \ldots + \phi_p \frac{\gamma_{h-p}}{\gamma_0} = \phi_1 \rho_{h-1} + \ldots + \phi_p \rho_{h-p} \quad h \geq 1 \nonumber \\ \Phi_p (B) \rho_h &=& 0 \end{eqnarray}\] Como vemos la estructura es más compleja ahora, de hecho deberíamos hacer un proceso iterativo complicado, pero lo que no sucede es que caigan a 0
Ecuaciones Yule-Walker
Si nos fijamos en las primeras \(p\) autocorrelaciones tenemos un sistema de ecuaciones conocido como ecuaciones de Yule-Walker
💻 Aplica la función AR_p para simular un \(AR(2)\) con \(n = 1000\), \(\sigma = 3\), \(\left(X_1, X_2 \right) = 0\) y \(\left(\phi_1, ~ \phi_2 \right) = (0.3, 0.5)\)
Código
set.seed(1234567)ggplot(AR_p_simul(p =2, 1000, 3, phi =c(0.3, 0.5), X =rep(0, 2))) +geom_line(aes(x = t, y = X_t)) +labs(title ="AR(2) con phi_1 = 0.3, phi_2 = 0.5 y sigma = 3") +theme_minimal()
pracma::roots(c(0.5, 0.3, 1))
[1] -0.3+1.382027i -0.3-1.382027i
Mod(pracma::roots(c(0.5, 0.3, 1)))
[1] 1.414214 1.414214
AR(p)
set.seed(1234567)AR_p_simul(p =2, 1000, 3, phi =c(0.3, 0.5), X =rep(0, 2)) |>ACF(X_t, lag_max =50) |>autoplot() +labs(title ="AR(2) con phi_1 = 0.3, phi_2 = 0.5 y sigma = 3") +theme_minimal()
Autocorrelación parcial
Es obvio que en el caso de los procesos \(MA(q)\) la función ACF nos permite su identificación pero no así en el caso de los procesos \(AR(p)\)
¿Cómo caracterizarlos? ¿Qué diferencia a un AR(1) de un AR(2)?
AR(1): \(X_t\) y \(X_{t-2}\) están relacionados indirectamente ya que \(X_t = \phi_1X_{t-1} + \varepsilon_{t} = \phi_{1}^2 X_{t-2} + \phi_{1} \varepsilon_{t-1} + \varepsilon_t\) pero dado que \(\left| \phi_1 \right| < 1\) es una relación que se va diluyendo
AR(2): sin embargo aquí \(X_t\) y \(X_{t-2}\)sí están relacionados directamente ya que \(X_t = \phi_1X_{t-1} + \phi_2 X_{t-2} + \varepsilon_{t}\) pero dado que \(\left| \phi_1 \right| < 1\) es una relación que se va diluyendo
Dicho de otra forma: en el caso de los \(AR(1)\), si conocemos \(X_{t-1}\), el valor de \(X_{t-2}\) es irrelevante (algo que no sucede en los \(AR(2)\)). Si yo pudiese calcular la autocorrelación entre \(X_t\) y \(X_{t-2}\) ELIMINANDO el efecto de \(X_{t-1}\), tendría que observar como dicho valor es nulo.
Autocorrelación parcial
¿Cómo eliminar de \(X_{t+k}\) el efecto de \(X_{t+1},\ldots,X_{t+k-1}\)?
En lugar de calular \(\rho_h = Cor[X_t, X_{t+h}]\) calcularemos
Dicho efecto lo estimaremos calculando el mejor predictor lineal óptimo de \(X_t\) y \(X_{t+h}\) en función de los retardos entre ellas \(\left(X_{t+1}, \ldots, X_{t+h-1} \right)\)
y se demuestra que \(\alpha_p = \phi_p\) y \(\alpha_h = 0\) para \(h > p\). Este truncamiento a partir de \(h > p\) no se da en un \(MA(q)\) ya que se puede entender como un \(AR(\infty)\).
set.seed(12345)AR_p_simul(p =2, n =10000, sigma =1, phi =c(-0.3, 0.4), X =rnorm(n =2, sd =1)) |>PACF(X_t, lag_max =30) |>autoplot() +labs(title ="AR(1) con phi_1 = -0.3, phi_2 = 0.4 y sigma = 1") +theme_minimal()
procesos \(MA(q)\): las autocorrelaciones parciales caen exponencialmente
procesos \(AR(p)\): las autocorrelaciones son aprox 0 a partir de \(p\) y la correlación correspondiente a \(p\) se sale de manera evidente
MA/AR en fable
¿Cómo predecir un MA o AR en fable?
Lo primero que haremos será generar un proceso \(AR(3)\) y dividir en train y test
set.seed(1234567)datos <-AR_p_simul(p =3, n =1000, sigma =3,phi =c(-0.2, -0.3, 0.4), X =c(0, 0, 0)) train <- datos |>slice(1:950)test <- datos |>slice(951:1000)
MA/AR en fable
Tras ello vamos a aplicar en model() el alisado simple y doble, y los modelos AR y MA haciendo uso de la función ARIMA() (basta con indicar la media y los órdenes p y q en la función pdq())
# A tibble: 16 × 10
.model .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
<chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 alisado_sim… Trai… 9.03e-2 3.83 3.07 Inf Inf 0.626 0.629 -0.264
2 alisado_dob… Trai… -9.03e-4 3.87 3.09 -Inf Inf 0.631 0.635 -0.252
3 AR1 Trai… -8.37e-3 3.70 2.94 64.2 153. 0.600 0.607 -0.126
4 AR2 Trai… -7.15e-3 3.25 2.57 44.8 172. 0.525 0.532 0.198
5 AR3 Trai… -7.96e-3 2.95 2.36 26.4 199. 0.481 0.484 0.00967
6 MA1 Trai… -8.18e-3 3.58 2.83 52.5 160. 0.578 0.587 0.00985
7 MA2 Trai… -5.57e-3 3.48 2.77 107. 168. 0.566 0.570 -0.171
8 MA3 Trai… -8.04e-3 3.15 2.50 46.1 161. 0.511 0.517 -0.0408
9 AR1 Test 6.66e-2 5.13 4.22 98.0 98.0 NaN NaN -0.311
10 AR2 Test 1.33e-2 4.89 4.09 105. 111. NaN NaN -0.312
11 AR3 Test 7.91e-2 4.76 4.01 48.9 128. NaN NaN -0.297
12 MA1 Test 3.47e-2 5.08 4.20 99.6 99.6 NaN NaN -0.317
13 MA2 Test 1.02e-3 5.07 4.21 110. 110. NaN NaN -0.320
14 MA3 Test 8.99e-2 5.02 4.16 87.5 106. NaN NaN -0.301
15 alisado_dob… Test -1.03e-1 5.16 4.26 99.4 112. NaN NaN -0.308
16 alisado_sim… Test 1.70e-1 5.16 4.24 99.8 102. NaN NaN -0.308
MA/AR en fable
Tanto en train como en test el mejor modelo es el AR(3) (lógico ya que los datos los hemos generado así)
fit_arima |>accuracy() |>bind_rows(predicciones |>accuracy(test)) |>slice_min(RMSE,n =1, by = .type)
# A tibble: 2 × 10
.model .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
<chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 AR3 Training -0.00796 2.95 2.36 26.4 199. 0.481 0.484 0.00967
2 AR3 Test 0.0791 4.76 4.01 48.9 128. NaN NaN -0.297
MA/AR en fable
Vamos entonces a repetir el proceso solo con dicho modelo. Además con report() vamos a poder obtener las estimación de los coeficientes (\(0.1914\), \(0.3162\) y \(-0.4154\), bastante preciso)
Tras el proceso, para asegurarnos de que aunque sea el mejor modelo lo que nos queda sin explicar (residuo) sea ruido blanco vamos a aplicar el test de Ljung-Box.
# antes de modelizardatos |>features(X_t, ljung_box)
MA(q): \(\Psi_0 = 1\) y \(\Psi_j = 0\) para todo \(j \geq q\)
AR(p): \(\Psi_0 = 1\) y \(\Psi_j = 0\) para todo \(j <0\), tal que \(\Phi_p(B) \Psi(B) = 1\) (se imponen condiciones en el decrecimiento de los coeficientes)
Los procesos ARMA intentarán combinar ambas propiedades para poder representar nuestros procesos con órdenes \(p\) y \(q\) bajos.
Procesos ARMA
Diremos que \(\left\lbrace X_t \right\rbrace_{t}\) es un proceso (mixto) autorregresivo-media móvil de orden \((p,q)\), denotado como \(ARMA(p,q)\) si
donde \(p \geq 1\) será el orden autorregresivo y \(q \geq 1\) será el orden de medias móvil
Usando el polinomio de retardos, podemos redefinir un ARMA(p, q) como
\[\Phi_p (B)X_t = \Theta_q (B) \varepsilon_t, \quad \left\lbrace \varepsilon_{t} \right\rbrace \text{ ruido blanco}\] tal que si \(\Phi_p (B)=1\) es un proceso puro MA(q) y si \(\Theta_q (B)=1\) es un proceso puro AR(p).
Procesos ARMA
Dado que \(\Phi_p (B)\) y \(\Theta_q (B)\) son polinomios, ambos pueden ser expresados en función de sus raíces tal que
Deben de ser diferentes: los polinomios \(\Phi_p (B)\) y \(\Theta_q (B)\) no pueden tener raíces comunes ya que si tuviesen, por ejemplo, \(\lambda_1 = \delta_1\), entonces podríamos cancelar en ambos lados
\[\left( 1- \lambda_2 B\right) \ldots \left( 1- \lambda_p B\right)X_t = \left( 1- \delta_2 B\right) \ldots \left( 1- \delta_q B\right) \varepsilon_t\] teniendo una sobreparametrización ya que en realidad el proceso sería ARMA(p-1, q-1).
Dado que este proceso puede ser farragoso podemos con ARIMA() realizar un ajuste automático dejando libre los órdenes e indicándole el \(p\) inicial con el que buscar, el \(q\) final y la métrica para elegir el mejor modelo (BIC penaliza modelos complejos y es consistente)
fit <- ARMA_2_1 |>model("ARMA"=ARIMA(X_t ~pdq(p_init =1, q_init =1), ic ="bic"))report(fit)
Por defecto este ajuste automático se realiza haciendo uso del algoritmo Hyndman-Khandakar que lo realiza stepwise y greedy: empezando por la esquina superior de modelos, selecciona el mejor según la métrica dada; tras ello busca en su entorno (moviendo y restado 1 p y q); si encuentro alguno mejor, cambia automáticamente.
Ajuste automático
si ARIMA(..., greedy = FALSE): evalua todo el entorno, no se cambia al primero que encuentre mejor
si ARIMA(..., stepwise = FALSE): evalua un grid más amplio de modeos
Ajuste automático
También puedes indicarle un vector concreto de \(p\) y \(q\) a probar
fit <- ARMA_2_1 |>model("ARMA"=ARIMA(X_t ~pdq(p =0:4, q =0:4), ic ="bic"))report(fit)
Como suele ser habitual en el campo de la calibración de modelos, una opción muy habitual es la de la validación:
Construir distintos modelos con la información de train
Usar los conjuntos de la validación para evaluar los modelos (o qué configuración de hiperparámetros) y decidir cuál de ellos es mejor
Una vez elegido el modelo, volver a lanzarlo y evaluarlo en test
Validación cruzada
Una de las formas de validación más habitual es la validación cruzada: las observaciones del conjunto de train van rotando su rol.
Por ejemplo, si tenemos 100 observaciones en train, podemos hacer 100 iteraciones de validación, de manera que en cada una entrenamos el modelo con 99 de ellas y otra queda reservada solo para evaluar los modelos.
Validación cruzada
En el caso de las series temporales una estrategia habitual suele ser la siguiente:
Descartar las primeras \(n\) observaciones para validación: habrá un conjunto mínimo que siempre formará parte de train
Iteración i: entrenamos con las primeras \(n+i\) observaciones, evaluamos con una única observación \(n+i+1\).
Realizamos el promedio de las métricas de evaluación obtenidas de los conjuntos de validación.
Validación cruzada
Fíjate que lo anterior está basado en una one-step forecast (predicción a horizonte \(h = 1\)), pero quizás nuestro interés esté en ver cómo funciona nuestro método a horizontes de predicción mayores
Descartar las primeras \(n\) observaciones para validación: habrá un conjunto mínimo que siempre formará parte de train
Iteración i: entrenamos con las primeras \(n+i\) observaciones, evaluamos con una única observación \(n+i+h\).
Realizamos el promedio de las métricas de evaluación obtenidas de los conjuntos de validación.
Validación cruzada
El ejemplo inferior es para \(h = 4\).
Validación cruzada
Vamos a generar un ARMA(2, 3) con la función arima.sim() para simular procesos ARMA
Tras ello vamos a generar los subconjuntos de validación usando trian con stretch_tsibble(), indicándole el número de valores iniciales que siempre estarán en train, el tamaño que queremos incrementar los sucesivos conjuntos y un identificador de cada slot
Por ejemplo, vamos a reservar los 350 primeros valores y vamos a avanzar a horizonte 1 (es decir, 100 slots de validación)
Tras generar los slots de validación entrenamos los modelos con dichos datos. Vamos a probar los siguientes modelos:
\(ARMA(1, 1)\)
\(ARMA(2, 2)\)
\(ARMA(p, q)\) automático stepwise greedy
\(ARMA(p, q)\) automático stepwise no greedy
\(ARMA(p, q)\) automático no stepwise
Validación cruzada
\(ARMA(1, 1)\)
\(ARMA(2, 2)\)
\(ARMA(p, q)\) automático stepwise greedy
\(ARMA(p, q)\) automático stepwise no greedy
\(ARMA(p, q)\) automático no stepwise
fit <- data_cv |>model("ARMA_11"=ARIMA(X_t ~pdq(p =1, d =0, q =1), ic ="bic"),"ARMA_22"=ARIMA(X_t ~pdq(p =2, d =0, q =2), ic ="bic"),"ARMA_stepwise_greedy"=ARIMA(X_t ~pdq(p_init =1, d =0, q_init =1), ic ="bic"),"ARMA_stepwise"=ARIMA(X_t ~pdq(p_init =1, d =0, q_init =1), ic ="bic", greedy =FALSE),"ARMA"=ARIMA(X_t ~pdq(p_init =1, d =0, q_init =1), ic ="bic", greedy =FALSE, stepwise =FALSE))estimaciones <- fit |>augment()
Validación cruzada
Tendremos las métricas para cada modelo y cada slot de cv que podemos promediar
¿Cómo visualizar las métricas de validación cruzada
Código
fit |>accuracy() |>ggplot(aes(x = .model, y = RMSE, fill = .model, color = .model)) +geom_boxplot(alpha =0.5) +geom_jitter(width =0.25, alpha =0.7) + ggthemes::scale_color_colorblind() + ggthemes::scale_fill_colorblind() +theme_minimal()
Si te fijas los mejores son el ARMA (completo) y el ARMA stepwise (pero no greedy): aunque tardan más los resultados son mejores. Fíjate que en el caso de los greedy tenemos un gap en la calidad según hacía la dirección a la que haya decidido orientarse: a veces llega al mejor pero no siempre
Validación cruzada
Dado que los mejores son el ARMA (completo) y el ARMA stepwise (pero no greedy) de manera similar, optaremos por este último ya que es más rápido.
fit <- train |>model("ARMA_stepwise"=ARIMA(X_t ~pdq(p_init =1, d =0, q_init =1), ic ="bic",greedy =FALSE))estimaciones <- fit |>augment()report(fit)
Hasta ahora estamos usando siempre la misma metodología
¿Es mi proceso estocástico ruido blanco? Un proceso de ruido blanco tiene a) media constante 0; b) varianza constante; c) instantes incorrelados.
❌ no puede serlo si tiene tendencia
❌ no puede serlo si es heterocedástico
❌ no puede serlo si tiene estacionariedad o algún patrón temporal
datos |>features(X_t, ljung_box)
Metodología Box-Jenkins
Hasta ahora estamos usando siempre la misma metodología
Si no fuese ruido blanco, ¿es estacionario? Un proceso estacionario tiene a) media constante; b) varianza constante; c) correlaciones solo dependientes del lag \(h\).
❌ no puede serlo si tiene tendencia
❌ no puede serlo si es heterocedástico
❌ no puede serlo si tiene regiones muy diferentes a otras
datos |>features(X_t, unitroot_kpss)
Metodología Box-Jenkins
Hasta ahora estamos usando siempre la misma metodología
Si es estacionario, ¿qué estructura ARMA (p, q) tiene?
Si tuviese solo parte MA(q) -> autocorrelaciones simples caerían a 0 a partir del retardo \(q\).
datos |>ACF(X_t, lag_max = ...) |>autoplot()
Si tuviese solo parte AR(p) -> autocorrelaciones parciales caerían a 0 a partir del retardo \(p\).
datos |>PACF(X_t, lag_max = ...) |>autoplot()
Metodología Box-Jenkins
Hasta ahora estamos usando siempre la misma metodología
Si es estacionario, ¿qué estructura ARMA (p, q) tiene?
En esta identificación hay que tener en cuenta que
Si ninguna de las dos se cumple implica que tenemos una estructura mixta ARMA(p, q) con \(p,q > 1\).
La identificación se realiza enfrentando las autocorrelaciones muestrales respecto a las que tendríamos si tuviésemos un ruido blanco simulado (bandas de confianza)
En esta fase el objetivo no es obtener un modelo adecuado a la primera ya que estabmos usando una muestra, sino restringir el conjunto de todos los posibles modelos ARMA a un subconjunto más pequeño
Metodología Box-Jenkins
Hasta ahora estamos usando siempre la misma metodología
Si es estacionario, ¿qué estructura ARMA (p, q) tiene?
La idea es empezar por el modelo más sencillo de los plausibles para luego proceder a su diagnóstico
Fase de diagnosis: ¿es el modelo ajustado el más verosimil teniendo en cuenta que tenemos una muestra de dicha distribución?
Metodología Box-Jenkins
Hasta ahora estamos usando siempre la misma metodología
Fase de diagnosis: ¿es el modelo ajustado el más verosimil teniendo en cuenta que tenemos una muestra de dicha distribución?
El objetivo es chequear que las hipótesis planteadas se cumplen
fit |>gg_tsresiduals()
¿Tiene el residuo media 0?
t.test(estimaciones$.resid, mu =0)
Metodología Box-Jenkins
Hasta ahora estamos usando siempre la misma metodología
Fase de diagnosis: ¿es el modelo ajustado el más verosimil teniendo en cuenta que tenemos una muestra de dicha distribución?
El objetivo es chequear que las hipótesis planteadas se cumplen
¿Tiene el residuo varianza constante?
lm(data = estimaciones, .resid ~ t) |>summary()
Metodología Box-Jenkins
Hasta ahora estamos usando siempre la misma metodología
Fase de diagnosis: ¿es el modelo ajustado el más verosimil teniendo en cuenta que tenemos una muestra de dicha distribución?
El objetivo es chequear que las hipótesis planteadas se cumplen
Hasta ahora estamos usando siempre la misma metodología
Fase de diagnosis: ¿es el modelo ajustado el más verosimil teniendo en cuenta que tenemos una muestra de dicha distribución?
En definitiva, ¿son nuestros residuos ruido blanco?
estimaciones |>features(.resid, ljung_box)
Metodología Box-Jenkins
Hasta ahora estamos usando siempre la misma metodología
Fase de diagnosis: ¿es el modelo ajustado el más verosimil teniendo en cuenta que tenemos una muestra de dicha distribución?
También podemos comprobar si los residuos son normales: si lo fuesen, sabríamos que incorrelación (lineal) implicaría independencia: nada por modelizar, hemos acabado.
Hasta ahora estamos usando siempre la misma metodología
Fase de diagnosis: ¿es el modelo ajustado el más verosimil teniendo en cuenta que tenemos una muestra de dicha distribución?
Un chequeo que puede ser también interesante es realizar una sobreparametrización del modelo: si hemos ajustado un ARMA(p, q) y cumple las hipótesis, ¿funciona sustancialmente mejor un ARMA(p + 1, q) o un ARMA(p , q+ 1)?
Si la fase de diagnosis no fuese exitosa deberíamos comenzar de nuevo el proceso ajustando un proceso de mayor orden
Metodología Box-Jenkins
Hasta ahora estamos usando siempre la misma metodología
Fase de evaluación: ¿cómo funciona nuestro modelo? ¿cuánto se equivoca?
Es aquí donde deberemos aplicar lo que conocemos respecto a validación de modelos
separamos en train y test
si tenemos varios posibles modelos paramétricos candidatos y tenemos dudas de cuál funcionaría mejor (o de si un modelo más complejo es rentable), realizamos una validación cruzada para obtener un métrica promedio
El modelo ganador lo volvemos a lanzar en train y test para obtener las métricas finales.
Metodología Box-Jenkins
Hasta ahora estamos usando siempre la misma metodología
Si no fuese ruido blanco, ¿es estacionario? Un proceso estacionario tiene a) media constante; b) varianza constante; c) correlaciones solo dependientes del lag \(h\).
❌ ¿Y si no fuese estacionario?
Procesos no estacionarios
Hasta ahora hemos supuesto que nuestro procesos eran estacionarios ya quelos procesos ARMA(p,q) solo pueden ser aplicados a procesos estacionarios.
Sin embargo una de las etapa más importantes es, en caso de que no lo sea, detectar la estructura no estacionaria de mi proceso
¿Falla en la estacionariedad en varianza?
¿Falla en la estacionariedad en media?
Procesos no estacionarios
Si nuestro proceso no fuera estacionario en varianza deberíamos evaluar que tendencia tiene la variabilidad de la serie respecto a su nivel para transformar los datos. Para ello tenemos dos opciones:
Agrupar los datos y calcular para cada grupo la media \(\overline{x}\) y la desviación típica \(s_{x}\): si suponemos que \(s_x = k \overline{x}^{\alpha}\), tomando logaritmos tenemos que \(\log(s_x) = \log(k) + \alpha \log( \overline{x})\), ¿qué \(\lambda = 1 - \alpha\) usar? (si fuese homocedástica, \(\alpha = 0\) y \(\lambda = 1\); si hubiese que aplicar logaritmo a los datos, \(\alpha = 1\) y \(\lambda = 0\))
La misma idea solo que ya directamente medimos \(CV_{x} = s_{x} / \overline{x}\) y chequeamos si tiene una tendencia constante (equivalente a lo anterior)
datos |>features(X_t, guerrero)
Procesos no estacionarios
Si nuestro proceso no fuera estacionario en media deberíamos evaluar que tendencia tiene la propia serie.
¿Cómo hacer la serie estacionaria en media?
La forma más sencilla es realizar lo que se conoce como diferenciaciones d ela serie. Dado un proceso estocástico \(X_t\) definiremos el operador diferencia de orden 1 como
\[Y_t:= \nabla X_t = (1-B)X_t = X_t - X_{t-1}\]
Dado un proceso estocástico \(X_t\) definiremos el operador diferencia de orden d como
Con dicho operador diremos que \(X_t\) es un proceso autoregresivo integrado de medias moviles ARIMA(p, d, q) si tras aplicar \(d\) diferencias con \(\nabla^d\), el proceso resultante es un proceso ARMA(p, q) invertible y estacionario y se define como
Fíjate que \(Y_t = \nabla^d X_t\) tiene que cumplir que sea un ARMA(p, q) estacionario.
Procesos ARIMA
Para identificar el orden \(d\) suele funcionar la regla de que \(d=1\) elimina tendencias lineales y \(d=2\) tendencias cuadráticas. También podemos visualizar la serie tras aplicar difference() o hacer uso del contraste unitroot_ndiffs que nos indica el orden \(d\) que identifica en nuestro proceso (o dejar libre el parámetro en ARIMA() y que el ajuste automática decida)
datos |>features(X_t, unitroot_ndiffs)
Clase 22: procesos SARIMA
Procesos SARIMA
Como hemos dicho un proceso mixto ARIMA(p, d, q) es definido como
\[\left(1 - \phi_1 B - \ldots - \phi_p B^p \right) \left( 1 - B \right)^d X_t = \left(1 - \theta_1 B - \ldots - \theta_q B^q \right) \varepsilon_t\] que de manera compacta puede ser definido como
Este enfoque es suficientemente flexible para modelizar procesos no estacionarios (pero que lo pueden ser tras diferenciar y aplicar una transformación) pero sin estacionalidad. ¿Cómo incluirla?
Para entenderlo vamos a empezar con lo que se conoce como procesos estacionales puros, procesos donde solo hay una componente estacional y nada más. Por ejemplo, pensemos en la ecuación de un AR(1)
\[X_t - \phi_1 X_{t-1} = (1 - \phi_1 B) X_t = \Phi_1(B)X_t = \varepsilon_t\] En un contexto puramente estacional, el concepto «dato anterior» no es el dato que corresponde al instante previo \(t-1\) sino al dato que corresponde al periodo estacional previo. Si tenemos una estacionalidad de periodo \(s\) entonces retroceder estacionalmente 1 periodo es \(X_{t-s}\).
Por denotaremos como \(AR(1)_s\) un proceso autorregresivo puramente estacional de periodo s y orden 1 al proceso
¿Cómo simularlos? Para ello vamos a definir un \(AR(1)_{s = 12}\) con \(\widetilde{\phi}_1 = 0.8\). Para ello vamos a hacer uso del paquete {astsa}
# install.packages("astsa")library(astsa)library(tsibble)library(tibble)n <-10000SAR_1 <-tibble("t"=1:n,"X_t"=sarima.sim(n = n, sar =0.8, S =12, sd =1.5)) |>as_tsibble(index = t)
¿Qué forma tendrán las autocorrelaciones?
Procesos SARIMA
ACF: en el caso de un AR(1) las correlaciones descendían a 0 asintóticamente; en el caso de un \(AR(1)_s\) dicho comportamiento asintótico solo sucede en las correlaciones múltiplos de s
PACF: en el caso de un AR(1) las correlaciones a partir de la primera caían a 0 drásticamente; en el caso de un \(AR(1)_s\)decaen drásticamente a 0 pero a partir de la correlación \(\rho_s\).
Procesos SARIMA
Vamos ahora a simular un \(AR(2)_{s = 12}\) con \(\widetilde{\phi}_1 = 0.5\) y \(\widetilde{\phi}_2 = -0.8\).
# install.packages("astsa")library(astsa)library(tsibble)library(tibble)n <-10000SAR_2 <-tibble("t"=1:n,"X_t"=sarima.sim(n = n, sar =c(0.5, -0.8), S =12, sd =1.5)) |>as_tsibble(index = t)
¿Cómo comprobar ahora que los coeficientes son válidos?
Fíjate que ahora al ser estacional los coeficientes afectan a \(B^{24}\) y \(B^{12}\) solo así que debo meter 0’s en medio.
ACF: en el caso de un AR(2) las correlaciones descendían a 0 asintóticamente; en el caso de un \(AR(2)_s\) dicho comportamiento sucede PERO solo en las correlaciones múltiplos de s.
PACF: en el caso de un AR(2) las correlaciones a partir de la segunda caían a 0 drásticamente; en el caso de un \(AR(1)_s\) dicho comportamiento sucede PERO a partir de la correlación \(\rho_s\).
Procesos SARIMA
💻 Te toca: simula distintos procesos \(MA\) puramente estacionales y «adivina» la forma de sus ACF/PACF antes de comprobarlo
Procesos SARIMA
Como pasa en los procesos ordinarios, en el caso de procesos puramente estacionales también podremos juntar ambas parte y definir un proceso mixto autorregresivo de medias móviles puramente estacional de periodo s al proceso \(ARMA(P, Q)_s\) definido como
Vamos ahora a simular un \(ARMA(2, 3)_{s = 12}\) con \(\widetilde{\phi}_1 = 0.5\) y \(\widetilde{\phi}_2 = -0.8\), y \(\widetilde{\theta}_1 = -0.2\)\(\widetilde{\theta}_2 = 0.5\) y \(\widetilde{\theta}_3 = -0.6\)
Código
n <-10000SARMA_23 <-tibble("t"=1:n,"X_t"=sarima.sim(n = n, sar =c(0.5, -0.8),sma =c(-0.2, 0.5, -0.6), S =12, sd =1.5)) |>as_tsibble(index = t)
Procesos SARIMA
ACF: en el caso de un AR(2) las correlaciones descendían a 0 asintóticamente; en el caso de un \(AR(2)_s\) dicho comportamiento sucede PERO solo en las correlaciones múltiplos de s.
PACF: en el caso de un AR(2) las correlaciones a partir de la segunda caían a 0 drásticamente; en el caso de un \(AR(1)_s\) dicho comportamiento sucede PERO a partir de la correlación \(\rho_s\).
Diferencias estacionales
Como pasa en los procesos ordinarios, en el caso de procesos puramente estacionales también podremos tener que el proceso no sea estacionario en media por lo que podremos definir un proceso mixto autorregresivo de medias móviles puramente estacional de periodo s con D diferencias al proceso \(ARIMA(P, D, Q)_s\) definido como un proceso \(ARMA(P, Q)_s\) tras \(D\) diferencias estacionales
Para comprobar si es necesario aplicar diferencias ordinarias
datos |>features(X_t, unitroot_ndiffs)
Para comprobar si es necesario aplicar diferencias estacionales
datos |>features(X_t, unitroot_nsdiffs)
Diferencias estacionales
Vamos ahora a simular un \(SARIMA(0, 2, 0)_{s = 12}\)
Código
n <-10000SARIMA_020 <-tibble("t"=1:n,"X_t"=sarima.sim(n = n, D =2, S =12, sd =1.5)) |>as_tsibble(index = t)
Diferencias estacionales
ACF: en el caso de un ARIMA(0, 2, 0) las correlaciones tenían tendencia sin decaer a 0; en el caso de un \(SARIMA(0, 2, 0)_s\) dicho comportamiento sucede en general PERO sobre todo en las correlaciones múltiplos de s.
Procesos SARIMA
Para terminar, vamos a juntar ambas partes (ordinaria y estacional): definiremos un proceso mixto autorregresivo integrado de medias móviles con componentes estacionales de periodo s al proceso \(SARIMA(p, d, q) \times (P, D, Q)_s\) definido como
\[ \Phi_p(B) \widetilde{\Phi}_P(B^s) \nabla_{s}^D \nabla^d X_t= \Theta_q(B) \widetilde{\Theta}_Q(B^s) \varepsilon_t\] tal que aplicando \(d\) diferencias ordinarias y \(D\) estacionales tenemos un proceso \(SARMA(p, q) \times (P, Q)_s\).
Modelización en fable
Vamos ahora a simular un \(ARIMA(2, 1, 3)_{s = 12}\) con \(\widetilde{\phi}_1 = 0.5\) y \(\widetilde{\phi}_2 = -0.8\), y \(\widetilde{\theta}_1 = -0.2\)\(\widetilde{\theta}_2 = 0.5\) y \(\widetilde{\theta}_3 = -0.6\). Tras ello vamos a realizar el ajusta en fable.
Código
n <-5000set.seed(1234567)SARMA_213 <-tibble("t"=1:n,"X_t"=sarima.sim(n = n, ar =0.5, ma =-0.9, d =2,sar =-0.7, sma =0.8, D =1,S =12, sd =1)) |>as_tsibble(index = t)
Lo primero es comprobar ya no solo si hay que aplicar una diferencia ordinaria sino si hace falta aplicar diferencias estacionales
.png)
