Información Sismos

Colombia es considerado un país de alto riesgo sísmico debido a que se ubica en el Cinturón de Fuego del Pacífico, una zona que concentra el 75% de los volcanes que existen en el mundo y donde suceden aproximadamente el 80% de los temblores más fuertes a nivel mundial.

Este es llamado también Anillo de Fuego del Pacífico y está integrado por la zona montañosa del oeste de Argentina, Chile, Perú, Ecuador, Colombia, Panamá, Costa Rica, Nicaragua, El Salvador, Honduras, Guatemala, México, Estados Unidos y Canadá, para luego doblar a la altura de las Islas Aleutianas y bajar por las costas e islas de Rusia, Japón Taiwán, Filipinas, Indonesia, Malasia, Timor Oriental, Brunéi, Singapur, Papúa Nueva Guinea, Islas Salomón, Tonga, Tuvalu y Nueva Zelanda.

En el caso específico de Colombia, el país se encuentra en dos áreas de subducción importantes, pues, por un lado, tiene la placa de Nazca con la Sudamericana y esta última que también choca con la placa del Caribe, lo que da paso a que tiemble constantemente.

Ante esta situación, los departamentos de Nariño, Chocó, Caldas y Santander son los lugares en donde más tiembla; en este último se encuentra el Municipio de Los Santos, considerado como la segunda zona más sísmica del mundo.

Por lo tanto, analizar y predecir esta situación es crucial para implementar medidas de precaución y manejo de riesgos ante posibles eventos sísmicos en Colombia.

Recurso Técnico de información:

Los datos para efectuar esta investigación se obtuvieron del Servicio Geológico Colombiano (SGC) a través del Catálogo de Sismicidad en formato Excel correspondiente al período desde el 28 de febrero de 2011 hasta el 28 de febrero de 2018. Esta base de datos contiene 33.848 registros.

El SGC es una entidad pública adscrita al Ministerio de Minas y Energía de Colombia, la cual tiene como función principal la investigación geocientífica del territorio colombiano para generar conocimiento sobre los recursos naturales y los riesgos geológicos del país. Teniendo en cuenta lo anterior, se eligió como fuente de información, ya que, es una fuente confiable y verificada que garantiza la precisión y completitud de los datos.

Autores

Estudiantes de la Maestría en Ciencia de Datos de la Universidad Javeriana de Cali.

Contenido de la Data

A continuación se presenta la base de datos a utilizar en la investigación.

#sismo <- read_csv("C:/Users/victo/Documents/Maestría/1.2. Segundo semestre/Series de tiempo/data_sismo.csv")#Victoria
#sismo <- read_csv("../../data_sismo.csv")#wladimir
sismo <- read_csv("D:/Documentos/GitHub/Series_de_tiempo_2/Series_de_tiempo/data_sismo.csv")#Mario
## Rows: 33848 Columns: 16
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr   (3): DEPARTAMENTO, MUNICIPIO, ESTADO
## dbl  (11): LATITUD (grados), LONGITUD (grados), PROFUNDIDAD (Km), MAGNITUD M...
## date  (1): FECHA
## time  (1): HORA_UTC
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
sismo <- mutate_if(sismo, is.character, tolower) # pasamos a minúsculas
glimpse(sismo)
## Rows: 33,848
## Columns: 16
## $ FECHA                    <date> 2011-02-28, 2011-02-28, 2011-02-28, 2011-02-…
## $ HORA_UTC                 <time> 02:34:30, 02:41:08, 06:20:53, 10:59:51, 11:0…
## $ `LATITUD (grados)`       <dbl> 6.789, 6.812, 6.746, 6.821, 6.800, 6.769, 6.7…
## $ `LONGITUD (grados)`      <dbl> -73.144, -73.129, -73.096, -73.095, -73.112, …
## $ `PROFUNDIDAD (Km)`       <dbl> 133.6, 144.3, 142.9, 143.9, 141.8, 147.3, 149…
## $ `MAGNITUD Ml`            <dbl> 1.9, 3.5, 1.7, 1.9, 1.4, 2.5, 2.8, 2.0, 1.9, …
## $ `MAGNITUD Mw`            <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ DEPARTAMENTO             <chr> "santander", "santander", "santander", "santa…
## $ MUNICIPIO                <chr> "los_santos", "los_santos", "los_santos", "lo…
## $ `# FASES`                <dbl> 6, 15, 7, 7, 4, 10, 11, 5, 5, 12, 6, 8, 4, 10…
## $ `RMS (Seg)`              <dbl> 0.2, 0.5, 0.5, 0.4, 0.1, 0.3, 0.4, 0.3, 0.2, …
## $ `GAP (grados)`           <dbl> 159, 104, 165, 129, 177, 128, 134, 137, 189, …
## $ `ERROR LATITUD (Km)`     <dbl> 3.0, 3.7, 5.1, 4.5, 3.3, 2.8, 3.4, 3.8, 2.1, …
## $ `ERROR LONGITUD (Km)`    <dbl> 5.5, 7.9, 11.6, 7.6, 13.6, 4.3, 5.7, 9.1, 8.2…
## $ `ERROR PROFUNDIDAD (Km)` <dbl> 8.3, 7.7, 8.8, 7.4, 3.9, 4.7, 6.5, 6.4, 4.3, …
## $ ESTADO                   <chr> "revisado", "revisado", "revisado", "revisado…

Observación de los datos

El Tiempo Universal Coordinado (HORA_UTC):

Es una escala de tiempo que es mantenida por los laboratorios de tiempo de todo el mundo y es determinada por relojes atómicos de alta precisión. El tiempo UTC es preciso a aproximadamente un nanosegundo (millonésima parte de un segundo) por día.

Magnitud Local (Ml):

Donde A es la máxima amplitud de la traza registrada y Ao la amplitud máxima que sería producida por un sismo patrón, siendo éste aquel que produciría una deflexión de 0.001 mm en un sismógrafo ubicado a 100 km del epicentro.

Magnitud de momento (Mw):

Esta magnitud es la más robusta; a diferencia de ML, mB y MS, la escala Mw no se satura, por lo que hoy en día es la más confiable y la más usada por las agencias dedicadas a la detección de sismos. También es la magnitud más usada por científicos para comparar los tamaños entre sismos.

Root Mean Square (RMS):

Es la medida del error, comparando la diferencia promedio entre el tiempo de arribo teórico y el tiempo de arribo observado en segundos, utilizando las lecturas de los sismogramas.

Brecha sísmica (GAP):

Zona geológica en la que no ha ocurrido un sismo fuerte durante un periodo prolongado de tiempo.

Descripción Datos 

skim(sismo)
Data summary
Name sismo
Number of rows 33848
Number of columns 16
_______________________
Column type frequency:
character 3
Date 1
difftime 1
numeric 11
________________________
Group variables None

Variable type: character

skim_variable n_missing complete_rate min max empty n_unique whitespace
DEPARTAMENTO 0 1 9 9 0 1 0
MUNICIPIO 0 1 10 10 0 1 0
ESTADO 0 1 8 8 0 1 0

Variable type: Date

skim_variable n_missing complete_rate min max median n_unique
FECHA 0 1 2011-02-28 2018-02-28 2014-08-31 2527

Variable type: difftime

skim_variable n_missing complete_rate min max median n_unique
HORA_UTC 0 1 1 secs 86376 secs 35560.5 secs 27764

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
LATITUD (grados) 0 1.00 6.79 0.03 6.73 6.77 6.79 6.81 6.88 ▂▇▇▂▁
LONGITUD (grados) 0 1.00 -73.14 0.03 -73.19 -73.16 -73.14 -73.12 -73.04 ▃▇▅▁▁
PROFUNDIDAD (Km) 0 1.00 145.06 5.16 1.80 142.00 145.30 148.10 188.20 ▁▁▁▇▁
MAGNITUD Ml 0 1.00 1.89 0.55 0.60 1.50 1.80 2.20 6.30 ▇▇▁▁▁
MAGNITUD Mw 31826 0.06 3.43 0.47 1.50 3.10 3.40 3.70 6.40 ▁▇▇▁▁
# FASES 0 1.00 9.89 6.34 3.00 6.00 8.00 12.00 72.00 ▇▁▁▁▁
RMS (Seg) 0 1.00 0.39 0.12 0.00 0.30 0.40 0.50 1.60 ▅▇▁▁▁
GAP (grados) 0 1.00 144.34 45.45 47.00 133.00 138.00 155.00 319.00 ▂▇▂▁▁
ERROR LATITUD (Km) 0 1.00 4.12 1.63 1.30 3.00 3.80 4.80 33.80 ▇▁▁▁▁
ERROR LONGITUD (Km) 0 1.00 4.85 1.89 1.20 3.60 4.50 5.60 41.70 ▇▁▁▁▁
ERROR PROFUNDIDAD (Km) 0 1.00 6.43 1.92 2.10 5.10 6.10 7.40 89.90 ▇▁▁▁▁

La profundidad de sismos en kilómetros muestra una media de 145.0610908, con una desviación estándar de 5.16222274.Por su aprte, las magnitudes Ml y Mw, tienen medias de 1.8909980 y 3.4296736 respectivamente. La magnitud Mw muestra una tasa de completitud más baja, con un 5.97% de datos faltantes.

El número de fases registradas por evento tiene una media de 9.8856949, el RMS (Root Mean Square) y GAP (Brecha) presentan una media de 0.3912964 y 144.3370952 respectivamente.

Imputar Datos

  • Existen 31826 datos faltantes NA para la variable Magnitud Mw por tanto se decide remover de la base de datos
  • Para las demás variables de tipo categórico, fecha, tiempo, y numérico no hay datos faltantes
  • Las variables categóricas corresponden a la descripción del departamento de Santander, en su municipio de Los_Santos y para la columna estado se encuentran todos los datos en revisado, por tanto se procede a removerlas de la data original sismo
  • Para el ejercicio académico se ha decidido remover las variables de error de latitud, longitud, y profundidad
  • También se renombra las variables que contienen comillas en su estructura de columnas
sismo_2 <- sismo
# las columnas retiradas fueron 7,8,9,13,14,15,16
sismo_2 <- sismo_2 [, c(1,2,3,4,5,6,10,11,12)] # data con 9 variables
colnames(sismo_2)[3] <-"LATITUD"  # renombrar variable
colnames(sismo_2)[4] <-"LONGITUD"  # renombrar variable
colnames(sismo_2)[5] <-"PROFUNDIDAD"  # renombrar variable
colnames(sismo_2)[6] <-"MAGNITUD_ML"  # renombrar variable
colnames(sismo_2)[7] <-"NUM_FASES"  # renombrar variable
colnames(sismo_2)[8] <-"RMS"  # renombrar variable
colnames(sismo_2)[9] <-"GAP"  # renombrar variable
head(sismo_2)

Variable Numéricas

Los siguientes gráficos muestran la distribución de la latitud, longitud, profundidad, magnitud, número de fases, Root Mean Square ( RMS) y la brecha sismica (GAP) de los sismos en el municipio de Los Santos, Santander. Cada gráfico representa la frecuencia de sismos en el eje Y un dato específico relacionado con cada una de las variables en el eje x.

La profundidad presenta una media de 145 km, con relación a la magnitud es una forma asimetrica a la derecha con una media de 1.8, el número de fases es una forma asimetrica a la izquierda y tiene una media de 9.8, por su parte el RMS y el GAP presentan una forma multimodal ya que aparecen varios picos y sus medias son de 0.3 y de 144 respectivamente.

g1 = ggplot(sismo_2, aes(x = LATITUD)) + geom_histogram(fill="royalblue") + theme_gray()
g2 = ggplot(sismo_2, aes(x = LONGITUD)) + geom_histogram(fill="royalblue") + theme_gray()
g3 = ggplot(sismo_2, aes(x = PROFUNDIDAD)) + geom_histogram(fill="royalblue") + theme_gray()
g4 = ggplot(sismo_2, aes(x = MAGNITUD_ML)) + geom_histogram(fill="royalblue") + theme_gray()
g5 = ggplot(sismo_2, aes(x = NUM_FASES)) + geom_histogram(fill="royalblue") + theme_gray()
g6 = ggplot(sismo_2, aes(x = RMS)) + geom_histogram(fill="royalblue") + theme_gray()
g7 = ggplot(sismo_2, aes(x = GAP)) + geom_histogram(fill="royalblue") + theme_gray()


ggarrange(g1, g2, g3, g4, g5, g6, g7, labels = c("A", "B", "C","D", "E", "F","G"), ncol = 2, nrow = 2)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## $`1`

## 
## $`2`

## 
## attr(,"class")
## [1] "list"      "ggarrange"

Variable Serie Tiempo

Análisis de la Data

A continuación se se seleccionan las variables de análisis para el presente ejercicio:

Profundidad: esta variable se selcciono dado que representa la profundidad del foco sísmico y este influye en cómo se siente un terremoto en la superficie, los terremotos superficiales, son aquellos que tienen un foco a menos de 70 km de profundidad, estos generalmente causan mas daño ya que la energía que liberan no se dispera tanto antes de llegar a la superficie.

Magnitud: la magnitud en Colombia se reorta utilizando la escala de Richter, esta medida permite conocer cuanta es la energía liberada por un sismo.

sismo_3 <- sismo_2 # copiamos data2
sismo_3 <- sismo_3 [, c(1,5,6)] # data con 3 variables 1:fecha, 5:profundidad, 6:magnitud
#sismo_3$fecha_hora <- paste(sismo_3$FECHA," ",sismo_3$HORA_UTC)# unir columnas 1 y 2
#sismo_3$fecha_hora <- as.POSIXct(sismo_3$fecha_hora) # convertir a tipo date time
#sismo_3$FECHA_UNIX <- as.POSIXct(sismo_3$FECHA)
sismo_3 <- sismo_3[order(sismo_3$FECHA), ] #ordenamo sismo_3 de menor a mayor
head(sismo_3,2)

Se agrupan las fechas para obtener un dato único por día, con relación a la magnitud se decidió calcular el maximo de cada día para conocer el sismo que libera la mayor cantidad de energía y con relación a la profundidad se decidió calcular el minimo con el fin de determianr el sismo que esta mas cerca a la superficie y puede causar mas daño.

sismo_4_max <- sismo_3 %>% # Agrupamos por fecha y obtenemos el máximo de la magnitud para cada fecha
  group_by(FECHA) %>% 
  summarise(MAX_MAGNITUD_ML = max(MAGNITUD_ML),
            MIN_PROFUNDIDAD = min(PROFUNDIDAD))
head(sismo_4_max, 2)

Dado que no hay información para todos los días del año entre el 28 de febrero de 2011 y el 28 de febrero de 2018 primero se creo un data frame con todas las fechas en ese rango y posteriormente los datos faltantes en las variables MAGNITUD Y PROFUNCIDAD se llenaron de la siguiente forma:

Profundidad: la profundidad se decio llenar con la media que es 145.0610908. Esta decición esta fundamentada en que la revisar la gráfica de series de tiempo, si se llenaran con profundidad 0 generaria mucho ruido y no permitiria visualizar adecuadamente los simso que realmente presentan una profundida baja y cercana a la superficie.

Magnitud: la magnitud se decidio llenar con 0.

# Crear un data frame con todas las fechas en el rango deseado
todas_las_fechas <- data.frame(FECHA = seq(as.Date("2011-02-28"), as.Date("2018-02-28"), by = "day"))

sismo_4_max_completo <- todas_las_fechas %>%
  left_join(sismo_4_max, by = "FECHA") %>%
  mutate(MAX_MAGNITUD_ML = ifelse(is.na(MAX_MAGNITUD_ML), 0, MAX_MAGNITUD_ML),
         MIN_PROFUNDIDAD = ifelse(is.na(MIN_PROFUNDIDAD), 145.0610908, MIN_PROFUNDIDAD))

# Filtrar desde el 28 de febrero de 2017 hasta el 28 de febrero de 2018
sismo_4_max_completo_filtrado <- sismo_4_max_completo %>%
  filter(FECHA >= as.Date("2017-02-28") & FECHA <= as.Date("2018-02-28"))

head(sismo_4_max_completo_filtrado, 2)

Se revisa que el número de filas concurde para comprobar que se encuentran todos los días entre el 28 de febrero de 2011 y el 28 de febrero de 2018.

nrow(todas_las_fechas)
## [1] 2558
nrow(sismo_4_max_completo_filtrado)
## [1] 366

Serie de tiempo Magnitud

indice.ts <- ts(sismo_4_max_completo_filtrado$MAX_MAGNITUD_ML,start = c(2017,2),end = c(2018,2), frequency = 365)
#head(indice.ts,10) 
head(indice.ts)
## Time Series:
## Start = c(2017, 2) 
## End = c(2017, 7) 
## Frequency = 365 
## [1] 3.9 2.6 3.7 3.7 2.7 1.7

Datos de la serie

"Tipo ts"
## [1] "Tipo ts"
class(indice.ts)
## [1] "ts"
"Inicio serie"
## [1] "Inicio serie"
start(indice.ts)
## [1] 2017    2
"fin serie"
## [1] "fin serie"
end(indice.ts)
## [1] 2018    2

Gráfica de la serie

hist(indice.ts,main="",ylab="valor", col='orange', xlab ="MAGNITUD_ML")
title(main="MAGNITUD DEL SISMO EN EL TIEMPO")

La distribución de la magnitud de los sismos en el municipio de Los Santos Santander evidencia que la mayoria de los sismos presentados son de magnitud baja y no presnetan un mayor riezgo para la población. Hay presencia de sismos de magnitud fuerte sin embargo hay que revisar la profundidad de estos para determinar si causan daños a la población.

plot(indice.ts,main="",ylab="valor", col='orange', xlab ="FECHA")
title(main="MAGNITUD DEL SISMO EN EL TIEMPO")

Serie de tiempo Profundidad

P_indice.ts <- ts(sismo_4_max_completo_filtrado$MIN_PROFUNDIDAD,start = c(2017,2),end = c(2018,2), frequency = 365)
#head(indice.ts,10) 
head(P_indice.ts)
## Time Series:
## Start = c(2017, 2) 
## End = c(2017, 7) 
## Frequency = 365 
## [1] 136.8 139.6 138.5 133.9 136.1 136.8

Datos de la serie

"Tipo ts"
## [1] "Tipo ts"
class(P_indice.ts)
## [1] "ts"
"Inicio serie"
## [1] "Inicio serie"
start(P_indice.ts)
## [1] 2017    2
"fin serie"
## [1] "fin serie"
end(P_indice.ts)
## [1] 2018    2

Gráfica de la serie

hist(P_indice.ts,main="",ylab="valor", col='saddlebrown', xlab ="PROFUNDIDA")
title(main="PROFUNDIDAD DEL SISMO (KM) EN EL TIEMPO")

La distribución de la profundiad de los sismos en el municipio de Los Santos Santander evidencia que la mayoria de los sismos presentados son de alta profunidad y no presnetan un mayor riezgo para la población.

plot(P_indice.ts,main="",ylab="valor", col='saddlebrown', xlab ="FECHA")
title(main="PROFUNDIDAD DEL SISMO (KM) EN EL TIEMPO")

Gráfica de sismos

Profundidad

library(ggplot2)
library(TSstudio)


data_profundidad_from_2012 <- sismo_4_max_completo[-(1:308), ]

data_serie <- ts(data_profundidad_from_2012$MIN_PROFUNDIDAD, frequency=365, start=2012)
plot(data_serie)

autoplot(data_serie)+
        labs(title = "Serie de tiempo",       
             x = "Tiempo",
             y = "Valor",
             colour = "#00a0dc")+
        theme_bw() 

fit <- decompose(data_serie, type='additive')

autoplot(data_serie, series="Serie tiempo") + 
        autolayer(trendcycle(fit), series="Tendencia") +
        labs(title = "Serie de tiempo",      
             x = "Tiempo",
             y = "Valor"
        ) + 
        theme_bw()
## Warning: Removed 364 rows containing missing values or values outside the scale range
## (`geom_line()`).

ggseasonplot(data_serie)

ts_seasonal(data_serie, type = "all")
## Warning: No trace type specified and no positional attributes specified
## No trace type specified:
##   Based on info supplied, a 'scatter' trace seems appropriate.
##   Read more about this trace type -> https://plotly.com/r/reference/#scatter
## No scatter mode specifed:
##   Setting the mode to markers
##   Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode

Magnitud

library(ggplot2)
library(TSstudio)

data_magnitud_from_2012 <- sismo_4_max_completo[-(1:308), ]


data_serie <- ts(data_magnitud_from_2012$MAX_MAGNITUD_ML, frequency=365, start=2012)
plot(data_serie)

autoplot(data_serie)+
        labs(title = "Serie de tiempo",       
             x = "Tiempo",
             y = "Valor",
             colour = "#00a0dc")+
        theme_bw() 

fit <- decompose(data_serie, type='additive')

autoplot(data_serie, series="Serie tiempo") + 
        autolayer(trendcycle(fit), series="Tendencia") +
        labs(title = "Serie de tiempo",      
             x = "Tiempo",
             y = "Valor"
        ) + 
        theme_bw()
## Warning: Removed 364 rows containing missing values or values outside the scale range
## (`geom_line()`).

ggseasonplot(data_serie)

ts_seasonal(data_serie, type = "all")
## Warning: No trace type specified and no positional attributes specified
## No trace type specified:
##   Based on info supplied, a 'scatter' trace seems appropriate.
##   Read more about this trace type -> https://plotly.com/r/reference/#scatter
## No scatter mode specifed:
##   Setting the mode to markers
##   Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode

Holt-Winter

A continuación, se aplica la metolodgía Holt-Winter y de suavizamiento a la variable tiempo. Esta técnica de análisis de series temporales se utiliza para predecir datos con componentes de tendencia y estacioanlidad.

Magnitud

library(aTSA)
## 
## Adjuntando el paquete: 'aTSA'
## The following objects are masked from 'package:tseries':
## 
##     adf.test, kpss.test, pp.test
## The following object is masked from 'package:forecast':
## 
##     forecast
## The following object is masked from 'package:graphics':
## 
##     identify
# Transformar la serie (si es necesario)
adf.test(diff(indice.ts))
## Augmented Dickey-Fuller Test 
## alternative: stationary 
##  
## Type 1: no drift no trend 
##      lag   ADF p.value
## [1,]   0 -34.4    0.01
## [2,]   1 -24.0    0.01
## [3,]   2 -18.2    0.01
## [4,]   3 -16.6    0.01
## [5,]   4 -14.7    0.01
## [6,]   5 -13.5    0.01
## Type 2: with drift no trend 
##      lag   ADF p.value
## [1,]   0 -34.3    0.01
## [2,]   1 -23.9    0.01
## [3,]   2 -18.1    0.01
## [4,]   3 -16.5    0.01
## [5,]   4 -14.7    0.01
## [6,]   5 -13.5    0.01
## Type 3: with drift and trend 
##      lag   ADF p.value
## [1,]   0 -34.3    0.01
## [2,]   1 -23.9    0.01
## [3,]   2 -18.1    0.01
## [4,]   3 -16.5    0.01
## [5,]   4 -14.7    0.01
## [6,]   5 -13.5    0.01
## ---- 
## Note: in fact, p.value = 0.01 means p.value <= 0.01

Profundidad

library(aTSA)

# Transformar la serie (si es necesario)
adf.test(diff(P_indice.ts))
## Augmented Dickey-Fuller Test 
## alternative: stationary 
##  
## Type 1: no drift no trend 
##      lag   ADF p.value
## [1,]   0 -32.4    0.01
## [2,]   1 -23.2    0.01
## [3,]   2 -18.8    0.01
## [4,]   3 -16.8    0.01
## [5,]   4 -14.8    0.01
## [6,]   5 -13.8    0.01
## Type 2: with drift no trend 
##      lag   ADF p.value
## [1,]   0 -32.4    0.01
## [2,]   1 -23.2    0.01
## [3,]   2 -18.8    0.01
## [4,]   3 -16.7    0.01
## [5,]   4 -14.7    0.01
## [6,]   5 -13.8    0.01
## Type 3: with drift and trend 
##      lag   ADF p.value
## [1,]   0 -32.4    0.01
## [2,]   1 -23.1    0.01
## [3,]   2 -18.7    0.01
## [4,]   3 -16.7    0.01
## [5,]   4 -14.7    0.01
## [6,]   5 -13.8    0.01
## ---- 
## Note: in fact, p.value = 0.01 means p.value <= 0.01

Teniendo en cuenta los resultados del test de Dickey-Fuller aumentado (ADF), se evidencio que todos los valores del estadístico ADF son negativos y los valores p son 0.01 para todos tipos de prueba:

  • Tipo 1: sin deriva y sin tendencia
  • Tipo 2: con deriva y sin tendencia
  • Tipo 3: con deriva y tendencia

Con base en lo anterior hay suficiente evidencia para rechazar la hipótesis nula de una raíz unitaria, concluyendo que la serie temporal es estacionaria.

ARIMA

A continuación, se aplica la metodología Box-Jenkins también conocida como ARIMA, este es un enfoque sistemático para construir modelos de series temporales que pueden predecir valores futuros basados en valores pasados y se utiliza especialmente cuando los datos no tienen una tendencia clara o una estacionalidad simple como es el caso de nuestros datos.

Teniendo en cuenta lo mencionado anteriormente y dado que el modelo ARIMA es particularmente adecuado para series temporales estacionarias, se justifica su uso en lugar del modelo Holt-Winters.

Este último es más adecuado para datos con componentes de tendencia y estacionalidad claros, que no es el caso de la serie temporal en estudio de de magnitud y produndidad de sismos en el municipio de los Santos, Santander, según los resultados de la prueba ADF.

Magnitud

library(forecast)

indice.ts2 <- head(indice.ts, -30)

# Ajustar el modelo ARIMA
modelo_magnitud_arima <- auto.arima(indice.ts2)

# Resumen del modelo
summary(modelo_magnitud_arima)
## Series: indice.ts2 
## ARIMA(2,1,2) 
## 
## Coefficients:
##           ar1      ar2      ma1      ma2
##       -0.9609  -0.0635  -0.0565  -0.9072
## s.e.   0.1022   0.0567   0.0878   0.0865
## 
## sigma^2 = 0.3394:  log likelihood = -294.06
## AIC=598.11   AICc=598.3   BIC=617.18
## 
## Training set error measures:
##                       ME      RMSE      MAE       MPE     MAPE MASE
## Training set -0.04803261 0.5782438 0.458386 -5.474935 16.43939  NaN
##                       ACF1
## Training set -0.0008658055

El modelo ARIMA(2,1,2) ajustado a la serie temporal indice.ts2 de (magnitud de los sismos) presenta los siguientes coeficientes estimados: AR1 (-0.9609), AR2 (-0.0635), MA1 (-0.0565) y MA2 (-0.9072), todos con errores estándar que sugieren una precisión razonable en las estimaciones. Las medidas de error en el conjunto de entrenamiento incluyen un ME de -0.048, un RMSE de 0.578, y un MAE de 0.458, indicando un ajuste razonable del modelo con errores moderadamente bajos. El MAPE de 16.44% sugiere que el modelo tiene una precisión aceptable en términos relativos, aunque no excelente.

# Hacer predicciones para los próximos 4 meses utilizando la función forecast del paquete forecast
pronostico_magnitud_arima <- forecast::forecast(modelo_magnitud_arima, h=30)

# Graficar las predicciones
plot(pronostico_magnitud_arima, main="Pronóstico de Magnitud de Sismos", ylab="Magnitud", xlab="Tiempo")

Se presenta la serie temporal de la magnitud de los sismos desde principios de 2017 hasta principios de 2018, junto con el pronóstico de magnitud para los últimos 30 días del período analizado. La línea azul representa las predicciones del modelo Prophet, mientras que las áreas sombreadas en gris claro y oscuro representan los intervalos de confianza al 80% y 95%, respectivamente, las bandas de confianza más amplias hacia el final del período indican una mayor incertidumbre en las predicciones futuras.

library(forecast)
library(tseries)

# Ajustar el modelo ARIMA automáticamente
modelo_magnitud_arima <- auto.arima(sismo_4_max_completo_filtrado$MAX_MAGNITUD_ML)

# Mostrar el resumen del modelo
summary(modelo_magnitud_arima)
## Series: sismo_4_max_completo_filtrado$MAX_MAGNITUD_ML 
## ARIMA(2,1,2) 
## 
## Coefficients:
##           ar1      ar2      ma1      ma2
##       -0.9673  -0.0619  -0.0484  -0.9152
## s.e.   0.0989   0.0546   0.0848   0.0836
## 
## sigma^2 = 0.3366:  log likelihood = -318.87
## AIC=647.75   AICc=647.92   BIC=667.25
## 
## Training set error measures:
##                       ME      RMSE       MAE      MPE    MAPE      MASE
## Training set -0.04863675 0.5761563 0.4552714 -5.54113 16.4699 0.7116663
##                      ACF1
## Training set -0.001168494

El modelo ARIMA de magnitud muestra un ajuste razonable, aunque con ciertas limitaciones en la predicción debido a la variabilidad de los datos sísmicos.

# Gráficos de diagnóstico
tsdiag(modelo_magnitud_arima)

# Gráfico de los residuos
plot(residuals(modelo_magnitud_arima), main="Residuos del Modelo ARIMA")

# ACF de los residuos
Acf(residuals(modelo_magnitud_arima), main="ACF de los Residuos del Modelo ARIMA")

Las dos gráficas indican que el modelo ARIMA está bien especificado, ya que los residuos no presentan autocorrelación significativa y están distribuidos de manera aleatoria alrededor de cero. Sin embargo, la presencia de algunos picos sugiere que podría haber eventos atípicos o anomalías pero dado que en los sismos es normal, no se consideran datos atípicos.

# Diagnóstico del modelo ARIMA: revisión de residuos
checkresiduals(modelo_magnitud_arima)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(2,1,2)
## Q* = 3.2928, df = 6, p-value = 0.7713
## 
## Model df: 4.   Total lags used: 10

La gráfica compuesta presenta un análisis de los residuos del modelo ARIMA(2,1,2) ajustado para predecir la magnitud de los sismos. La serie temporal de los residuos (parte superior) muestra que los residuos están distribuidos alrededor de cero sin patrones evidentes, aunque con algunos picos ocasionales. La función de autocorrelación (ACF) de los residuos (parte inferior izquierda) revela que casi todos los valores de autocorrelación caen dentro de las bandas de confianza, indicando que no hay autocorrelación significativa en los residuos. Esto sugiere que el modelo ha capturado bien las dependencias temporales en los datos. En conjunto, estos análisis indican que el modelo ARIMA(2,1,2) es adecuado, ya que los residuos no presentan autocorrelación significativa y siguen una distribución aproximadamente normal, aunque hay algunos picos que podrían ser eventos atípicos.

# Calcular pronósticos dentro de la muestra
fitted_values <- fitted(modelo_magnitud_arima)
actual_values <- sismo_4_max_completo_filtrado$MAX_MAGNITUD_ML

# Calcular errores
errors <- actual_values - fitted_values

# Calcular métricas de error
rmse <- sqrt(mean(errors^2))
mae <- mean(abs(errors))
mape <- mean(abs(errors / actual_values)) * 100

# Mostrar las métricas de error
cat("RMSE:", rmse, "\n")
## RMSE: 0.5761563
cat("MAE:", mae, "\n")
## MAE: 0.4552714
cat("MAPE:", mape, "%\n")
## MAPE: 16.4699 %

En conjunto, estas métricas sugieren que el modelo ARIMA(2,1,2) tiene un rendimiento razonablemente bueno para predecir la magnitud de los sismos.

# Prueba de Ljung-Box para autocorrelaciones en los residuos
Box.test(residuals(modelo_magnitud_arima), lag=10, type="Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  residuals(modelo_magnitud_arima)
## X-squared = 3.2928, df = 10, p-value = 0.9737

Los residuos se consideran independientes, lo que sugiere que el modelo ARIMA ha capturado adecuadamente las dependencias temporales en los datos.

# Residuos del modelo
residuales_magnitud_arima <- residuals(modelo_magnitud_arima)

# Prueba de Shapiro-Wilk para normalidad
shapiro.test(residuales_magnitud_arima)
## 
##  Shapiro-Wilk normality test
## 
## data:  residuales_magnitud_arima
## W = 0.96371, p-value = 6.991e-08
# Prueba de Ljung-Box para autocorrelación
Box.test(residuales_magnitud_arima, type="Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  residuales_magnitud_arima
## X-squared = 0.00050384, df = 1, p-value = 0.9821

Los resultados de los dos tests muestran que los residuos son independientes, lo que indica que el modelo ha capturado adecuadamente las dependencias temporales en los datos.

Profundidad

P_indice.ts2 <- head(P_indice.ts, -30)

# Ajustar el modelo ARIMA
modelo_profundidad_arima <- auto.arima(P_indice.ts2)

# Resumen del modelo
summary(modelo_profundidad_arima)
## Series: P_indice.ts2 
## ARIMA(0,1,1) 
## 
## Coefficients:
##           ma1
##       -0.9695
## s.e.   0.0143
## 
## sigma^2 = 124.6:  log likelihood = -1284.47
## AIC=2572.94   AICc=2572.98   BIC=2580.57
## 
## Training set error measures:
##                       ME    RMSE      MAE       MPE     MAPE MASE         ACF1
## Training set -0.08478686 11.1297 4.121583 -31.59247 34.47417  NaN -0.007245018

El modelo ARIMA(0,1,1) ajustado para la serie temporal P_indice.ts2 (profundidad) muestra que el término de media móvil es significativo y que el modelo ha capturado adecuadamente las dependencias temporales, como lo indica la baja autocorrelación en los residuos.

# Hacer predicciones para los próximos 4 meses utilizando la función forecast del paquete forecast
pronostico_profundidad_arima <- forecast::forecast(modelo_profundidad_arima, h=30)

# Graficar las predicciones
plot(pronostico_profundidad_arima, main="Pronóstico de Profundidad de Sismos", ylab="Profundidad", xlab="Tiempo")

En la serie temporal, la mayoría de los sismos ocurren a profundidades cercanas a 145 km, con algunas anomalías de eventos sísmicos a profundidades significativamente menores, sin embargo dado el contexto no se consideran datos atipicos o anomalias. El pronóstico, representado por la línea azul y las bandas de confianza sombreadas en gris claro y oscuro (al 80% y 95%, respectivamente), indica que se espera que las profundidades de los sismos se mantengan en niveles similares a los observados anteriormente, alrededor de 145 km.

# Ajustar el modelo ARIMA automáticamente
modelo_profundidad_arima <- auto.arima(sismo_4_max_completo_filtrado$MIN_PROFUNDIDAD)

# Mostrar el resumen del modelo
summary(modelo_profundidad_arima)
## Series: sismo_4_max_completo_filtrado$MIN_PROFUNDIDAD 
## ARIMA(0,1,1) 
## 
## Coefficients:
##           ma1
##       -0.9667
## s.e.   0.0135
## 
## sigma^2 = 115.4:  log likelihood = -1385.44
## AIC=2774.87   AICc=2774.91   BIC=2782.67
## 
## Training set error measures:
##                      ME     RMSE      MAE       MPE    MAPE      MASE
## Training set 0.04955531 10.71536 4.016825 -28.90284 31.8042 0.7223793
##                      ACF1
## Training set -0.007641542

El modelo ARIMA para profundidad muestra que los residuos no presentan autocorrelación significativa y el coeficiente MA1 es altamente significativo. Sin embargo, las métricas de error indican que las predicciones tienen un margen de error considerable.

# Gráficos de diagnóstico
tsdiag(modelo_profundidad_arima)

# Gráfico de los residuos
plot(residuals(modelo_profundidad_arima), main="Residuos del Modelo ARIMA")

# ACF de los residuos
Acf(residuals(modelo_profundidad_arima), main="ACF de los Residuos del Modelo ARIMA")

El modelo ha capturado adecuadamente las dependencias temporales en los datos. A pesar de estos valores atípicos debido a los datos de análisis, la mayoría de los residuos se mantienen estables y cercanos a cero, confirmando que el modelo ARIMA es adecuado en su ajuste.

# Diagnóstico del modelo ARIMA: revisión de residuos
checkresiduals(modelo_profundidad_arima)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,1,1)
## Q* = 1.5562, df = 9, p-value = 0.9967
## 
## Model df: 1.   Total lags used: 10

En la parte superior, la serie temporal de los residuos muestra que la mayoría de los residuos están cercanos a cero. La gráfica ACF (abajo a la izquierda) muestra que las autocorrelaciones de los residuos están dentro de las bandas de confianza, lo que sugiere que no hay autocorrelación significativa en los residuos, indicando un buen ajuste del modelo. Finalmente, el histograma de los residuos (abajo a la derecha) muestra una distribución aproximadamente normal, aunque con una ligera asimetría y algunos valores extremos.

# Residuos del modelo
residuales_profundidad_arima <- residuals(modelo_profundidad_arima)

# Prueba de Shapiro-Wilk para normalidad
shapiro.test(residuales_profundidad_arima)
## 
##  Shapiro-Wilk normality test
## 
## data:  residuales_profundidad_arima
## W = 0.32756, p-value < 2.2e-16
# Prueba de Ljung-Box para autocorrelación
Box.test(residuales_profundidad_arima, type="Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  residuales_profundidad_arima
## X-squared = 0.021548, df = 1, p-value = 0.8833

El test de Shapiro-Wilk tiene un estadístico W de 0.32756 y un p-valor menor a 2.2e-16, lo que indica que los residuos no siguen una distribución normal. Por otro lado, el test de Box-Ljung tiene un valor de chi-cuadrado de 0.021548 con 1 grado de libertad y un p-valor de 0.8833, sugiriendo que no hay autocorrelación significativa en los residuos.

Algoritmo Facebook’s Prophet

A continuación se utiliza una herramienta que ayuda a predecir lo que va a pasar en el futuro usando datos del pasado. Esta herramienta fue creada por Facebook y es muy buena para trabajar con datos que cambian de manera regular como los datos de los sismos. A continuación se utiliza el algoritmo Prophet para la variable magnitud y profundidad de los sismos.

Magnitud

# Cargar el paquete
library(prophet)
## Cargando paquete requerido: Rcpp
## Cargando paquete requerido: rlang
library(ggplot2)
library(forecast)

# Asumiendo que los datos están en un data.frame llamado sismo
# Convertir las fechas y profundidades en un formato adecuado para Prophet
data_magnitud <- data.frame(
  ds = as.Date(sismo_4_max_completo_filtrado$FECHA), 
  y = sismo_4_max_completo_filtrado$MAX_MAGNITUD_ML
)

# Quitar los últimos 30 días de los datos
train_data_magnitud <- head(data_magnitud, -30)

# Ajustar el modelo Prophet con los datos de entrenamiento
modelo_magnitud_prophet <- prophet(train_data_magnitud)
## Disabling yearly seasonality. Run prophet with yearly.seasonality=TRUE to override this.
## Disabling daily seasonality. Run prophet with daily.seasonality=TRUE to override this.
# Hacer predicciones para los próximos 30 días
future <- make_future_dataframe(modelo_magnitud_prophet, periods = 30)
pronostico_magnitud_prophet <- predict(modelo_magnitud_prophet, future)

# Graficar las predicciones
plot(modelo_magnitud_prophet, pronostico_magnitud_prophet) +
  labs(title = "Pronóstico de Magnitud de Sismos", y = "Magnitud", x = "Tiempo")

La línea azul central representa la predicción del modelo, mientras que las áreas sombreadas en azul claro y oscuro indican los intervalos de confianza del 80% y 95%, respectivamente. Los puntos negros representan las observaciones reales de la magnitud de los sismos durante el período de tiempo mostrado. La mayoría de las observaciones se encuentran dentro de los intervalos de confianza, lo que sugiere que el modelo tiene un buen desempeño en capturar la variabilidad de las magnitudes de los sismos.

# Calcular los residuos
residuales_magnitud_prophet <- data_magnitud$y[1:nrow(train_data_magnitud)] - pronostico_magnitud_prophet$yhat[1:nrow(train_data_magnitud)]

# Test de Shapiro-Wilk para normalidad de los residuos
shapiro_test <- shapiro.test(residuales_magnitud_prophet)
print(shapiro_test)
## 
##  Shapiro-Wilk normality test
## 
## data:  residuales_magnitud_prophet
## W = 0.95033, p-value = 3.206e-09
# Test de Box-Ljung para independencia de los residuos
box_ljung_test <- Box.test(residuales_magnitud_prophet, lag = 20, type = "Ljung-Box")
print(box_ljung_test)
## 
##  Box-Ljung test
## 
## data:  residuales_magnitud_prophet
## X-squared = 7.869, df = 20, p-value = 0.9927
# Graficar los residuos
residuals_df <- data.frame(time = train_data_magnitud$ds, residuals = residuales_magnitud_prophet)
ggplot(residuals_df, aes(x = time, y = residuals)) +
  geom_line() +
  labs(title = "Residuos del Modelo Prophet", x = "Tiempo", y = "Residuos")

# Graficar la función de autocorrelación de los residuos
acf(residuales_magnitud_prophet, main = "ACF de los Residuos del Modelo Prophet")

La gráfica que muestra la función de autocorrelación (ACF) de los residuos, indica que casi todos los valores de autocorrelación están dentro de las bandas de confianza, sugiriendo que no hay autocorrelación significativa en los residuos y que el modelo ha capturado adecuadamente las dependencias temporales en los datos. La segunda gráfica, que muestra la serie temporal de los residuos, muestra que los residuos están distribuidos alrededor de cero con algunas fluctuaciones pero sin patrones evidentes de autocorrelación. Por su parte, el test de Box-Ljung muestra que los residuos son independientes, lo que indica que el modelo Prophet ha capturado adecuadamente las dependencias temporales en los datos.

Profundidad

# Cargar el paquete
library(prophet)

# Asumiendo que los datos están en un data.frame llamado sismo_4_max_completo_filtrado
# Convertir las fechas y profundidades en un formato adecuado para Prophet
data_profundidad <- data.frame(
  ds = as.Date(sismo_4_max_completo_filtrado$FECHA), 
  y = sismo_4_max_completo_filtrado$MIN_PROFUNDIDAD
)

# Quitar los últimos 30 días de los datos
train_data_profundidad <- head(data_profundidad, -30)

# Ajustar el modelo Prophet con los datos de entrenamiento
modelo_profundidad_prophet <- prophet(train_data_profundidad)
## Disabling yearly seasonality. Run prophet with yearly.seasonality=TRUE to override this.
## Disabling daily seasonality. Run prophet with daily.seasonality=TRUE to override this.
# Hacer predicciones para los próximos 30 días
future <- make_future_dataframe(modelo_profundidad_prophet, periods = 30)
pronostico_profundidad_prophet <- predict(modelo_profundidad_prophet, future)

# Graficar las predicciones
plot(modelo_profundidad_prophet, pronostico_profundidad_prophet) +
  labs(title = "Pronóstico de Profundidad Mínima de Sismos", y = "Profundidad (km)", x = "Tiempo")

La línea azul central representa la predicción del modelo, mientras que las áreas sombreadas en azul claro indican los intervalos de confianza del 80% y 95%, respectivamente. Los puntos negros representan las observaciones reales de la profundidad de los sismos durante el período de tiempo mostrado. La mayoría de las observaciones se encuentran dentro de los intervalos de confianza, lo que sugiere que el modelo tiene un buen desempeño en capturar la variabilidad de la profundidad de los sismos.

# Calcular los residuos
residuales_profundidad_prophet <- train_data_profundidad$y - predict(modelo_profundidad_prophet, train_data_profundidad)$yhat

# Test de Shapiro-Wilk para normalidad de los residuos
shapiro_test <- shapiro.test(residuales_profundidad_prophet)
print(shapiro_test)
## 
##  Shapiro-Wilk normality test
## 
## data:  residuales_profundidad_prophet
## W = 0.34076, p-value < 2.2e-16
# Test de Box-Ljung para independencia de los residuos
box_ljung_test <- Box.test(residuales_profundidad_prophet, lag = 20, type = "Ljung-Box")
print(box_ljung_test)
## 
##  Box-Ljung test
## 
## data:  residuales_profundidad_prophet
## X-squared = 3.2081, df = 20, p-value = 1
# Graficar los residuos
residuals_df <- data.frame(time = train_data_profundidad$ds, residuals = residuales_profundidad_prophet)
ggplot(residuals_df, aes(x = time, y = residuals)) +
  geom_line() +
  labs(title = "Residuos del Modelo Prophet para Profundidad", x = "Tiempo", y = "Residuos")

# Graficar la función de autocorrelación de los residuos
acf(residuales_profundidad_prophet, main = "ACF de los Residuos del Modelo Prophet para Profundidad")

La gráfica muestra que, aunque la mayoría de los residuos se distribuyen alrededor de cero, hay algunos valores atípicos negativos que el modelo no ha capturado adecuadamente. La gráfica de autocorrelación (ACF) de los residuos, indica que no hay autocorrelación significativa, ya que todos los valores están dentro de las bandas de confianza.

Los resultados del test de Shapiro-Wilk (W = 0.34076, p-valor < 2.2e-16) sugieren que los residuos no siguen una distribución normal. Por otro lado, el test de Box-Ljung (X-squared = 3.2081, df = 20, p-valor = 1) indica que no hay autocorrelación significativa en los residuos.

Gráfica combinada Magnitud

# Cargar las librerías necesarias
library(prophet)
library(forecast)
library(ggplot2)

# Suponiendo que los datos están en el data.frame llamado sismo_4_max_completo
data_magnitud <- data.frame(
  ds = as.Date(sismo_4_max_completo$FECHA), 
  y = sismo_4_max_completo$MAX_MAGNITUD_ML
)

# Ajustar el modelo Prophet con estacionalidad diaria
modelo_magnitud_prophet <- prophet(data_magnitud, daily.seasonality = TRUE)

# Crear el dataframe futuro para predicciones de Prophet
future <- make_future_dataframe(modelo_magnitud_prophet, periods = 120)
pronostico_magnitud_prophet <- predict(modelo_magnitud_prophet, future)

# Ajustar el modelo ARIMA
modelo_arima <- auto.arima(data_magnitud$y)
pronostico_arima <- forecast::forecast(modelo_arima, h = 120)
# Crear un dataframe para los resultados de Prophet
resultados_prophet <- data.frame(
  ds = future$ds,
  yhat = pronostico_magnitud_prophet$yhat
)

# Crear un dataframe para los resultados de ARIMA
fechas_arima <- seq(max(data_magnitud$ds) + 1, by = "day", length.out = 120)
resultados_arima <- data.frame(
  ds = fechas_arima,
  yhat = pronostico_arima$mean
)

# Combinar los resultados con los datos originales
resultados_comb <- rbind(
  data.frame(ds = data_magnitud$ds, y = data_magnitud$y, modelo = "Datos"),
  data.frame(ds = resultados_prophet$ds, y = resultados_prophet$yhat, modelo = "Prophet"),
  data.frame(ds = resultados_arima$ds, y = resultados_arima$yhat, modelo = "ARIMA")
)
# Graficar los resultados
ggplot() +
  geom_line(data = resultados_comb[resultados_comb$modelo == "Datos", ], 
            aes(x = ds, y = y, color = modelo), size = 0.5) +
  geom_line(data = resultados_comb[resultados_comb$modelo == "Prophet", ], 
            aes(x = ds, y = y, color = modelo), size = 1) +
  geom_line(data = resultados_comb[resultados_comb$modelo == "ARIMA", ], 
            aes(x = ds, y = y, color = modelo), size = 1) +
  labs(title = "Pronóstico de Magnitud de Sismos: Prophet vs ARIMA",
       y = "Magnitud", x = "Fecha") +
  theme_minimal() +
  scale_color_manual(values = c("Datos" = "black", "Prophet" = "blue", "ARIMA" = "red"))
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

En general, ambos modelos parecen ser adecuados para capturar la tendencia general de la magnitud de los sismos.

Gráfica combinada Pronóstico

library(prophet)
library(forecast)
library(ggplot2)

# Asumiendo que los datos están en un data.frame llamado sismo
# Convertir las fechas y profundidades en un formato adecuado para Prophet
data_profundidad <- data.frame(
  ds = as.Date(sismo_4_max_completo$FECHA), 
  y = sismo_4_max_completo$MIN_PROFUNDIDAD
)

# Ajustar el modelo Prophet
modelo_profundidad_prophet <- prophet(data_profundidad, daily.seasonality = TRUE)

# Hacer predicciones para los próximos 4 meses (120 días)
future <- make_future_dataframe(modelo_profundidad_prophet, periods = 120)
pronostico_profundidad_prophet <- predict(modelo_profundidad_prophet, future)

# Ajustar el modelo ARIMA
modelo_arima <- auto.arima(data_profundidad$y)
pronostico_arima <- forecast::forecast(modelo_arima, h = 120)
# Crear un dataframe para los resultados de Prophet
resultados_prophet <- data.frame(
  ds = future$ds,
  yhat = pronostico_profundidad_prophet$yhat
)

# Crear un dataframe para los resultados de ARIMA
fechas_arima <- seq(max(data_profundidad$ds) + 1, by = "day", length.out = 120)
resultados_arima <- data.frame(
  ds = fechas_arima,
  yhat = pronostico_arima$mean
)

# Combinar los resultados con los datos originales
resultados_comb <- rbind(
  data.frame(ds = data_profundidad$ds, y = data_profundidad$y, modelo = "Datos"),
  data.frame(ds = resultados_prophet$ds, y = resultados_prophet$yhat, modelo = "Prophet"),
  data.frame(ds = resultados_arima$ds, y = resultados_arima$yhat, modelo = "ARIMA")
)
# Graficar los resultados
ggplot() +
  geom_line(data = resultados_comb[resultados_comb$modelo == "Datos", ], 
            aes(x = ds, y = y, color = modelo), size = 0.5) +
  geom_line(data = resultados_comb[resultados_comb$modelo == "Prophet", ], 
            aes(x = ds, y = y, color = modelo), size = 1) +
  geom_line(data = resultados_comb[resultados_comb$modelo == "ARIMA", ], 
            aes(x = ds, y = y, color = modelo), size = 1) +
  labs(title = "Pronóstico de Profundidad de Sismos: Prophet vs ARIMA",
       y = "Pronostico", x = "Fecha") +
  theme_minimal() +
  scale_color_manual(values = c("Datos" = "black", "Prophet" = "blue", "ARIMA" = "red"))

En general, ambos modelos parecen ser adecuados para capturar la tendencia general de la magnitud de los sismos

Modelo Elman

Modelo Elman Profundidad

library(ggplot2)
library(RSNNS)
library(dplyr)
data = data_profundidad

train_size <- floor(0.80 * nrow(data))

train_data <- data[1:train_size, ]
test_data <- data[(train_size + 1):nrow(data), ]


max_val <- max(data$y)
min_val <- min(data$y)

normalize <- function(x) {
  return ((x - min_val) / (max_val - min_val))
}

train_data$y <- normalize(train_data$y)
test_data$y <- normalize(test_data$y)

create_lagged_matrix <- function(data, lag) {
  inputs <- data %>%
    dplyr::mutate(lagged_y = dplyr::lag(y, lag)) %>%
    na.omit()
  
  inputs_matrix <- as.matrix(inputs$lagged_y)
  outputs_matrix <- as.matrix(inputs$y)
  
  return(list(inputs = inputs_matrix, outputs = outputs_matrix))
}

lag <- 2  

train_matrices <- create_lagged_matrix(train_data, lag)
test_matrices <- create_lagged_matrix(test_data, lag)

elman_model <- elman(train_matrices$inputs, train_matrices$outputs, size = c(5), maxit = 500, learnFuncParams = c(0.1), initFunc = "JE_Weights")

predictions <- predict(elman_model, test_matrices$inputs)

# Desnormalizar las predicciones
denormalize <- function(x) {
  return (x * (max_val - min_val) + min_val)
}

predictions <- denormalize(predictions)

# Comparar las predicciones con los valores reales
results <- data.frame(
  Date = test_data$ds[(lag + 1):nrow(test_data)],
  Actual = denormalize(test_data$y[(lag + 1):nrow(test_data)]),
  Predicted = predictions
)

head(results,10)
ggplot(data = results, aes(x = Date)) +
  geom_line(aes(y = Actual, color = "Actual")) +
  geom_line(aes(y = Predicted, color = "Predicted")) +
  labs(title = "Predicciones vs Valores Reales Profundidad",
       x = "Fecha",
       y = "Valor",
       color = "Leyenda") +
  theme_minimal()

test_data_new <- head(test_data, -2)


plot(data$ds, data$y, type = 'l', col = 'blue', main = 'Predicciones del Modelo Elman Profundidad con todo el historico', xlab = 'Fecha', ylab = 'Profundidad')
lines(test_data_new$ds, predictions, col = 'red')
legend("bottomleft", legend = c("Real","Predicción Prueba"), col = c("blue", "red"), lty = 1)

La primera gráfica compara las predicciones del modelo con los valores reales de profundidad de los sismos. Las líneas rojas representan los valores reales, mientras que las líneas turquesas representan las predicciones. Se observa que las predicciones siguen la tendencia general de los valores reales, aunque hay una considerable variabilidad en los datos reales que el modelo no logra capturar completamente. La segunda gráfica muestra las predicciones del modelo Elman para la profundidad de los sismos usando todo el histórico de datos. Las líneas azules representan los datos reales y las líneas rojas representan las predicciones del modelo. Se puede observar que el modelo Elman también sigue la tendencia general de los datos históricos, pero las predicciones son menos volátiles en comparación con los valores reales.

Modelo Elman Magnitud

#prueba de concepto Elman
library(ggplot2)
library(RSNNS)
library(dplyr)
data = data_magnitud

train_size <- floor(0.80 * nrow(data))

train_data <- data[1:train_size, ]
test_data <- data[(train_size + 1):nrow(data), ]


max_val <- max(data$y)
min_val <- min(data$y)

normalize <- function(x) {
  return ((x - min_val) / (max_val - min_val))
}

train_data$y <- normalize(train_data$y)
test_data$y <- normalize(test_data$y)

create_lagged_matrix <- function(data, lag) {
  inputs <- data %>%
    dplyr::mutate(lagged_y = dplyr::lag(y, lag)) %>%
    na.omit()
  
  inputs_matrix <- as.matrix(inputs$lagged_y)
  outputs_matrix <- as.matrix(inputs$y)
  
  return(list(inputs = inputs_matrix, outputs = outputs_matrix))
}

lag <- 2  

train_matrices <- create_lagged_matrix(train_data, lag)
test_matrices <- create_lagged_matrix(test_data, lag)

elman_model <- elman(train_matrices$inputs, train_matrices$outputs, size = c(5), maxit = 500, learnFuncParams = c(0.1), initFunc = "JE_Weights")

predictions <- predict(elman_model, test_matrices$inputs)

# Desnormalizar las predicciones
denormalize <- function(x) {
  return (x * (max_val - min_val) + min_val)
}

predictions <- denormalize(predictions)

# Comparar las predicciones con los valores reales
results <- data.frame(
  Date = test_data$ds[(lag + 1):nrow(test_data)],
  Actual = denormalize(test_data$y[(lag + 1):nrow(test_data)]),
  Predicted = predictions
)

head(results,10)
ggplot(data = results, aes(x = Date)) +
  geom_line(aes(y = Actual, color = "Actual")) +
  geom_line(aes(y = Predicted, color = "Predicted")) +
  labs(title = "Predicciones vs Valores Reales Magnitud",
       x = "Fecha",
       y = "Valor",
       color = "Leyenda") +
  theme_minimal()

test_data_new <- head(test_data, -2)


plot(data$ds, data$y, type = 'l', col = 'blue', main = 'Predicciones del Modelo Elman Magnitud con todo el historico', xlab = 'Fecha', ylab = 'Magnitud')
lines(test_data_new$ds, predictions, col = 'red')
legend("bottomleft", legend = c("Real","Predicción Prueba"), col = c("blue", "red"), lty = 1)

La primera gráfica compara las predicciones de magnitud de los sismos (en turquesa) con los valores reales (en rojo) para el año 2017 y principios de 2018. Se observa que las predicciones siguen de cerca la línea central de los valores reales, pero no capturan completamente la variabilidad y los picos de los datos observados. La segunda gráfica muestra las predicciones del modelo Elman para la magnitud de los sismos utilizando todo el histórico de datos. Las líneas azules representan los datos reales y las líneas rojas las predicciones del modelo. El modelo Elman sigue la tendencia general de los datos históricos pero con menos variabilidad.

Modelo Jordan Profundidad

library(ggplot2)
library(RSNNS)
library(dplyr)
data = data_profundidad

train_size <- floor(0.80 * nrow(data))

train_data <- data[1:train_size, ]
test_data <- data[(train_size + 1):nrow(data), ]


max_val <- max(data$y)
min_val <- min(data$y)

normalize <- function(x) {
  return ((x - min_val) / (max_val - min_val))
}

train_data$y <- normalize(train_data$y)
test_data$y <- normalize(test_data$y)

create_lagged_matrix <- function(data, lag) {
  inputs <- data %>%
    dplyr::mutate(lagged_y = dplyr::lag(y, lag)) %>%
    na.omit()
  
  inputs_matrix <- as.matrix(inputs$lagged_y)
  outputs_matrix <- as.matrix(inputs$y)
  
  return(list(inputs = inputs_matrix, outputs = outputs_matrix))
}

lag <- 2  

train_matrices <- create_lagged_matrix(train_data, lag)
test_matrices <- create_lagged_matrix(test_data, lag)

jordan_model <- jordan(train_matrices$inputs, train_matrices$outputs, size = c(5), maxit = 500, learnFuncParams = c(0.1), initFunc = "JE_Weights")

predictions <- predict(jordan_model, test_matrices$inputs)

# Desnormalizar las predicciones
denormalize <- function(x) {
  return (x * (max_val - min_val) + min_val)
}

predictions <- denormalize(predictions)

# Comparar las predicciones con los valores reales
results <- data.frame(
  Date = test_data$ds[(lag + 1):nrow(test_data)],
  Actual = denormalize(test_data$y[(lag + 1):nrow(test_data)]),
  Predicted = predictions
)

head(results,10)
ggplot(data = results, aes(x = Date)) +
  geom_line(aes(y = Actual, color = "Actual")) +
  geom_line(aes(y = Predicted, color = "Predicted")) +
  labs(title = "Predicciones vs Valores Reales Profundidad",
       x = "Fecha",
       y = "Valor",
       color = "Leyenda") +
  theme_minimal()

test_data_new <- head(test_data, -2)


plot(data$ds, data$y, type = 'l', col = 'blue', main = 'Predicciones del Modelo Jordan Profundidad con todo el historico', xlab = 'Fecha', ylab = 'Profundidad')
lines(test_data_new$ds, predictions, col = 'red')
legend("bottomleft", legend = c("Real","Predicción Prueba"), col = c("blue", "red"), lty = 1)

Estas dos gráficas muestran predicciones de la profundidad de los sismos utilizando diferentes modelos. En la primera gráfica, se comparan los valores reales de la profundidad (línea roja) con las predicciones realizadas (línea azul claro) a lo largo del tiempo, indicando que las predicciones siguen de cerca los valores reales pero con cierta desviación. La segunda gráfica presenta las predicciones realizadas con el modelo Jordan (línea roja) contra los datos históricos (línea azul). Aquí, se observa que las predicciones del modelo Jordan tienden a estabilizarse en un rango más estrecho, lo cual puede indicar una mayor consistencia, pero también una posible subestimación de las fluctuaciones observadas en los datos reales.

Modelo Jordan Magnitud

library(ggplot2)
library(RSNNS)
library(dplyr)
data = data_magnitud

train_size <- floor(0.80 * nrow(data))

train_data <- data[1:train_size, ]
test_data <- data[(train_size + 1):nrow(data), ]


max_val <- max(data$y)
min_val <- min(data$y)

normalize <- function(x) {
  return ((x - min_val) / (max_val - min_val))
}

train_data$y <- normalize(train_data$y)
test_data$y <- normalize(test_data$y)

create_lagged_matrix <- function(data, lag) {
  inputs <- data %>%
    dplyr::mutate(lagged_y = dplyr::lag(y, lag)) %>%
    na.omit()
  
  inputs_matrix <- as.matrix(inputs$lagged_y)
  outputs_matrix <- as.matrix(inputs$y)
  
  return(list(inputs = inputs_matrix, outputs = outputs_matrix))
}

lag <- 2  

train_matrices <- create_lagged_matrix(train_data, lag)
test_matrices <- create_lagged_matrix(test_data, lag)

jordan_model <- jordan(train_matrices$inputs, train_matrices$outputs, size = c(5), maxit = 500, learnFuncParams = c(0.1), initFunc = "JE_Weights")

predictions <- predict(jordan_model, test_matrices$inputs)

# Desnormalizar las predicciones
denormalize <- function(x) {
  return (x * (max_val - min_val) + min_val)
}

predictions <- denormalize(predictions)

# Comparar las predicciones con los valores reales
results <- data.frame(
  Date = test_data$ds[(lag + 1):nrow(test_data)],
  Actual = denormalize(test_data$y[(lag + 1):nrow(test_data)]),
  Predicted = predictions
)

head(results,10)
ggplot(data = results, aes(x = Date)) +
  geom_line(aes(y = Actual, color = "Actual")) +
  geom_line(aes(y = Predicted, color = "Predicted")) +
  labs(title = "Predicciones vs Valores Reales Magnitud",
       x = "Fecha",
       y = "Valor",
       color = "Leyenda") +
  theme_minimal()

test_data_new <- head(test_data, -2)


plot(data$ds, data$y, type = 'l', col = 'blue', main = 'Predicciones del Modelo Jordan Magnitud con todo el historico', xlab = 'Fecha', ylab = 'Magnitud')
lines(test_data_new$ds, predictions, col = 'red')
legend("bottomleft", legend = c("Real","Predicción Prueba"), col = c("blue", "red"), lty = 1)

Las dos gráficas muestran las predicciones de la magnitud de sismos utilizando dos modelos diferentes. La primera gráfica presenta las predicciones versus los valores reales, donde las predicciones están representadas por una línea (azul claro) y los valores reales por una (línea roja). Se observa que las predicciones del modelo siguen de cerca los valores reales, aunque con algunas desviaciones. La segunda gráfica muestra las predicciones del modelo Jordan comparadas con el histórico completo de los datos. En esta gráfica se observa que las predicciones del modelo Jordan tienden a mantenerse cerca de los valores reales, aunque con menos variabilidad que los datos históricos.