# Dudas a @crissthiandi <albertocenaa@gmail.com>
#' Carga paquetes para un análisis de series de tiempo
#'
#' Ideal para usar en trabajos similar a tidyverse
#' escribir código para análisis de series de tiempo
#'
#' @return Mensaje si todo sale bien
#' @export
#' @encoding UTF-8
#'
#' @examples
#'
#' paquetes.tsrutina()
paquetes.tsrutina <- function(){
require('tseries',warn.conflicts = F,quietly = F)
require('lmtest',warn.conflicts = F,quietly = F)
require('pracma',warn.conflicts = F,quietly = F)
require('ggfortify',warn.conflicts = F,quietly = F)
require('forecast',warn.conflicts = F,quietly = F,)
require('tseries',warn.conflicts = F,quietly = F)
require('greybox',warn.conflicts = F,quietly = F)
require('readr',warn.conflicts = F,quietly = F)
message("Se han cargado los paquetes necesarios")
}
#' Repite un texto y lo imprime en la consola
#'
#' Repite un texto de forma seriada y se imprime en consola, ideal para hacer una separación entre elementos impresos en consola
#'
#'
#' @param repite String, que se repetirá, idealmente del tipo "-" "_" "*" etc.
#' @param num_repetidas Int, Cuantas veces se repite el string anterior.
#' @param color Color del texto impresa, soportamos red, blue, green, magenta, silver y yellow
#'
#' @return Impresión en consola
#' @export
#' @encoding UTF-8
#'
#' @encoding UTF-8
#'
#' separador("-",10,"blue")
#' separador()
#'
separador<-function(repite="-",num_repetidas=45,color="red"){
b <- paste0("cat(","crayon::",color,"(",
"paste0(rep(","'",repite,"',",num_repetidas,")",")","),'\n')")
a <- str2expression(b)
eval(a)
}
#' Pruebas estadísticas a una serie de tiempo
#'
#' Esta función incluye el calculo y la decisión de dos pruebas estadísticas.
#'
#' Realiza las siguientes pruebas:
#' \itemize{\item{Dickey-Fuller:}{ La prueba de estacionariedad}}
#' \itemize{\item{Durbin-Watson:}{ La prueba de autocorrelacion}}
#'
#' Recibe un dataframe de dos columnas, la primera el tiempo y
#' la segunda el valor de la serie.
#'
#' Tambien soporta el uso de objetos TimeSeries.
#'
#' @param datos Data.frame o objeto TS a analizar
#' @param frecuencia Frecuencia de los datos, en caso de TS sobrescribe los valores
#' @param validar_ Indicador para verificar los datos [True/False]
#'
#' @return Solo arroja la decisión a tomar, por defecto con respecto a p=0.05
#' @export
#' @encoding UTF-8
#'
#' @importFrom lmtest dwtest
#' @importFrom tseries adf.test
#' @importFrom crayon green red yellow
#'
#' @examples
#' base=data.frame(tiempo=seq(Sys.Date(),by="days",length=20),valores=(rexp(50)+1)*sin(1:50))
#' serie_tiempo_pruebas(datos=base,frecuencia=4)
#'
serie_tiempo_pruebas <-function(datos,frecuencia=NULL,validar_=FALSE,
msg=TRUE,pausa_off=1){
if(validar_){
paquetes.tsrutina()
pausa()
conditional.tsrutina(datos)
}
if(is.data.frame(datos)){
pausa()
separador()
if(is.data.frame(datos)){
base <- datos
names(base) <- c("x","y")
report_data <- list()
base2 <- base
base2$xx <- 1:length(base$x)
regresion <- lm(y~xx,base2)
cat(crayon::green("\n Prueba de presencia de autocorrelación Durbin-Watson Test \n"))
cat(crayon::green("\n (Prueba aplicada los residuos de una regresión lineal) \n"))
prueba <- lmtest::dwtest(regresion)
print(prueba)
report_data$prueba_dwtest <- prueba
if(msg){
p_valor <- readline('Inserte un p valor, (intro para p=0.05): \n')
}else{
p_valor <- 0.05
}
if(p_valor==""){
p_valor <- 0.05
}else{
cat(crayon::red(sprintf("El valor de p= %s",p_valor)),"\n")
p_valor <- as.numeric(p_valor)
}
report_data$prueba_dwtest_msg <-
if(prueba$p.value>0.05){
mensaje("No se puede rechazar \"H0 = No hay presencia de autocorrelacion en los residuos\" \n")
}else{
mensaje("Se rechaza H0, se obta por \"H1 = Hay presencia de autocorrelacion en los residuos\" \n")
}
separador()
sprintf("\n \n")
#Probar Estacionariedad para proceder a los modelos no probabilisticos
#Prueba Dickey-Fuller:
#H0: Hay presencia de una raiz unitaria en la serie de tiempo
#H1: La serie de tiempo es Estacionaria.
#No se rechaza H0 pues pvalor=0.467 es mayor que 0.05 (No se rechaza H0)
basets <- ts(base$y,frequency=frecuencia)
cat(crayon::green("\n Prueba de presencia de estacionariedad Dickey-Fuller Test \n"))
prueba <- tseries::adf.test(basets)
print(prueba)
report_data$prueba_adf.test <- prueba
if(msg){#si no hay mensaje entonces predeterminado
p_valor<-readline('Inserte un p valor, (intro para p=0.05) \n')
}else{
p_valor <- 0.05
}
if(p_valor==""){
p_valor <- 0.05
}else{
cat(crayon::red(sprintf("El valor de p= %s",p_valor)),"\n")
p_valor <- as.numeric(p_valor)
}
report_data$prueba_adf.test_msg <-
if(prueba$p.value>p_valor){
mensaje("No se puede rechazar H0 = Hay presencia de una raiz unitaria \n")
}else{
mensaje("Se rechaza H0, se obta por H1 = La serie de tiempo es Estacionaria \n")
}
separador()
return(invisible(report_data))
}else{
stop("El objeto debe ser data frame")
}
}else{
# TODO remorver recursividad en funciones (?)
if(is.ts(datos) & !is.mts(datos)){
elementos <- tratamiento.ts_set(datos)
frecuencia <- ifelse(is.null(frecuencia),elementos$frecu,frecuencia)
serie_tiempo_pruebas(elementos$data,frecuencia,
validar_ = validar_,msg = msg,
pausa_off = pausa_off)
}else{
stop("El objeto debe ser un data frame con dos elementos o una serie de tiempo univariada")
}
}
}
#' Primer filtro de las funciones en TSRutina
#'
#' @param datos objeto a verificar
#'
#' @description Verifica si un data.frame tiene dos columnas y si la segunda de estas
#' es de carácter numérico, si es serie de tiempo no hace nada, retorna NULL.
#'
#' @return NULL or stop() event
#' @encoding UTF-8
#'
#' @examples
#'
#' conditional.tsrutina(datos)
#'
conditional.tsrutina <- function(datos){
if(is.ts(datos)){
return(NULL)
} else if(is.mts(datos)) {
stop("Esta paqueteria no soporta series de tiempo multivariadas")
}
if(ncol(datos)!=2)
stop("Los datos deben tener solo dos columnas, tiempo y valor en ese orden")
if(!is.numeric(datos[,2])){
stop("La segunda columna deben ser los valores, la primera el tiempo")}
message("\n Primera columna => variable de Tiempo
\n Segunda columna => variable de valor \n")
}
#' Convert ts to data.frame
#'
#' @param datosts objeto ts a tranformar en data.frame class
#'
#' @description Convierte un objeto serie de tiempo a data.frame con dos columnas
#' x => variable fecha
#' y => variable valor
#' La variable fecha es inicializada con el valor start() de la serie de tiempo y termina con el numero de
#' frecuencias que se pueden hacer
#'
#' @details Si el objeto TS contiene una frecuencia no compatible con las programadas, se te pedirá ingresar el nombre del vector
#' (previamente creado en tu enviroment a.k.a entorno de desarrollo) con las fechas a usar. Mira el ejemplo 2.
#'
#' @return a data.frame object
#' @export
#' @encoding UTF-8
#'
#' @examples
#'
#' tratamiento.ts_set(sunspot.year)
#'
#' ## Ejemplo 2
#' \dontrun{
#' ### Se crea el vector con las fechas
#' fechas <- c("2021-05-01","2020-10-05") %>% as.Date()
#' ### al correr tratamiento.ts_set(datos) se ingresa el nombre : fechas
#' }
#'
#'
tratamiento.ts_set <- function(datosts){
datos_conver <- as.data.frame(datosts)
year <- start(datosts)[1L]
mes <- start(datosts)[2L]
dia <- start(datosts)[3L]
dia <- ifelse(test = is.na(dia),yes = 1,no = dia)
fecha_inicio <- as.Date(paste(dia,mes,year, sep = "-"),
format = "%d-%m-%Y")
frecuencia <- frequency(datosts)
avance <- switch(as.character(frecuencia),
"1" = 'year',
"12" = 'month',
"4" = '3 month',
"54" = '7 days',
"3" = '4 month',
"365" = 'days',
T = FALSE
)
try(if(avance){
cat(crayon::red("No se encontro frecuencia compatible con función Seq() \n
Ingrese el nombre del vector con las fechas de tu serie de tiempo. \n"))
fecha_secuencia <- readline("El vector debe existir en tu Enviroment: ")
},silent = 1)
fecha_secuencia <- tryCatch(
expr = seq(fecha_inicio, by=avance, length=nrow(datos_conver)),
error = function(e) {
message("\n Se intentara conectar con el objeto/vector llamado: \n",
fecha_secuencia,"\n")# De no exitir hay error
busqueda <- try(get(fecha_secuencia,envir = parent.frame()),
outFile = cat("Error, no se encontro el vector ",fecha_secuencia))
return(busqueda)
}
)
datos_conver <- data.frame(x = fecha_secuencia,y = datos_conver$x)
elementos <- list()
elementos$data <- datos_conver
elementos$frecu <- frequency(datosts)
elementos$inicio <- start(datosts)
elementos$fin <- end(datosts)
return(elementos)
}
#' Tratamiento de fechas TSR
#'
#' Tratamiento para que una fecha tenga el formato Date. De no tener un formato convertible se retorna NA
#' valúes
#'
#' @param fecha_vector Vector de valores con las fechas en formato carácter o Date
#' Formatos soportador: "\%d/\%m/\%Y", "\%Y-\%m-\%d" y "\%d-\%m-\%Y"
#'
#' @return Vector en formato fecha estándar de R
#' @export
#' @encoding UTF-8
#'
#' @importFrom readr parse_date
#'
#' @examples
#' tiempo=seq(Sys.Date(),by="days",length=20)
#' tratamiento.fechas.TSR(tiempo)
#'
tratamiento.fechas.TSR <- function(fecha_vector){
#En caso de que el vector ya sea fecha, es necesario abortar
if(class(fecha_vector) == "Date"){
return(fecha_vector)
} else if(class(fecha_vector) != "Character"){
fecha_vector <- as.character(fecha_vector)
}
fecha_vector_tratamiento <- readr::parse_date(fecha_vector,"%d%.%m%.%Y") #Y debe ser mayuscula para 4 digitos de año
if(is.na(fecha_vector_tratamiento[1])){
fecha_vector_tratamiento <- readr::parse_date(fecha_vector,"%Y%.%m%.%d") #Y debe ser mayuscula para 4 digitos de año
}
if(is.na(fecha_vector_tratamiento[1])){
fecha_vector_tratamiento <- readr::parse_date(fecha_vector,"%m%.%Y%.%d") #Y debe ser mayuscula para 4 digitos de año
}
if(is.na(fecha_vector_tratamiento[1])){
message("\n Si la variable tiempo es fecha, use el formato year-month-day: \n
\"%Y-%m-%d\" ")
}
return(fecha_vector_tratamiento)
}
#' Mensaje de validación de los datos en rutina TSR
#'
#' Imprime un head de los datos y te pregunta si todo esta bien.
#'
#' Ideal para usarse en las rutinas de TSRutina que deben preguntar si se interpretaron
#' bien los datos.
#'
#' @param datos Data.Frame de 2 columnas, fecha y valores respectivamente. (se hace tratamiento de fechas con \code{\link{tratamiento.fechas.TSR}})
#' @param frecuencia Frecuencia de la serie de tiempo
#' @param inicio Inicio de la serie de tiempo
#'
#' @return Una lista \code{\link{list}} que contiene dos elementos, la base de datos tratada y su versión en un objeto TimeSeries
#' @export
#' @encoding UTF-8
#'
#' @importFrom crayon green red yellow
#'
#' @examples
#' base=data.frame(tiempo=seq(Sys.Date(),by="days",length=20),valores=(rexp(50)+1)*sin(1:50))
#' checar_datos(datos=base,frecuencia=4,inicio=2010)
#'
checar_datos <- function(datos,frecuencia,inicio,msg=TRUE) {
names(datos) <- c("x","y")
datos$x <- tratamiento.fechas.TSR(datos$x)
if(msg){
print(head(datos))
message("\n ¿Están bien los datos a usar? \n
Si hay un error [Esc] \n De lo contrario [Enter] para continuar")
continuar<-readline(": \t")
if(!continuar=="")
stop("Corrige el error")
}
#creando el objeto time series y si inicio no es el minimo de fecha
inicio <- min(datos$x)
datosts <- ts(data = datos$y,frequency = frecuencia,start=inicio)
return(list(datos=datos,datosts=datosts))
}
#' Pausa la consola al imprimir resultados
#'
#' Realiza una pausa tanto indeterminada como por una cantidad de segundos
#'
#' Cuando \bold{duracion} no es definido la pausa es indefinida y se espera la entrada de "stop" o Esc para terminar la pausa
#'
#' @param duracion Tiempo de pausa en segundos, default Inf (infinito)
#'
#' @return
#' @export
#' @encoding UTF-8
#'
#' @importFrom crayon green red yellow
#'
#' @examples
#'
#' #para 10 segundos
#' pausa(10)
#'
pausa <- function(duracion = Inf){
#parche para control de pausa
a <- try(get("pausa_off",envir = parent.frame()),1)
if(a==0){
return(invisible(NULL))
}
if (is.infinite(duracion)) {
arg <- "*"
while (arg != "") {
cat("\n")
arg <- readline("[Intro] to continue | [stop/Esc] to exit... ")
if (arg == "stop") {
stop("El programa finalizo", call. = FALSE)
}
}
} else {
cat("Pause of", duracion, "seconds\n")
Sys.sleep(duracion)
}
invisible(NULL)
}
#' @title Rutinas en una serie de tiempo
#'
#' @description Realiza una rutina para una base de datos en la que se
#' incluyen ajustes por:
#' \itemize{\item{Regresion lineal}{}}
#' \itemize{\item{Promedio movil simple}{}}
#' \itemize{\item{Promedio movil ponderado}{}}
#' \itemize{\item{Suavizamiento exponencial simple}{}}
#' \itemize{\item{Suav. expon. doble o de Holt}{}}
#' \itemize{\item{Suavizamiento Holt-Winter}{}}
#'
#' @param datos Dataframe de no mas de 2 columnas, en el orden primero tiempo y
#' luego valor el tiempo va en formato fecha, y tiene que ser en el orden
#' dia-mes-year. La versión 2.1 soporta objetos de la clase TS (time series)
#' @param frecuencia Este es el periodo de la serie, trimestral = 3, cuatrimestral = 4
#' ,mensual = 12, etc.
#' @param inicio Este es el year a iniciar la serie de tiempo
#'
#' @param pausa_off Logical ¿Desea quitar las pausas en la rutina?
#'
#' @param msg Logical ¿Desea ignorar las preguntas de la rutina y tomar los valores por defecto?
#'
#' @return La salida no es como tal un objeto, si no una serie de impresiones de varios
#' varios analisis.
#' \itemize{\item{\bold{Plost}}{ Arroja una lista de plots que ayudan a ver el comportamiento de la serie y como ciertos ajustes se aproximan mejor a ella}}
#' \itemize{\item{\bold{Resumenes}}{ Arroja ciertos resumenes de ciertos ajustes o pruebas que se hacen}}
#' \itemize{\item{\bold{Modelo}}{ Modelo con el menor MSE(Error cuadratico medio)}}
#'
#' @details Se usa una salida interactiva en la que el usuario debe agregar ciertos
#' datos o tomar ciertas decisiones durante la rutina.
#'
#' @author Cristhian Diaz
#' @export
#' @encoding UTF-8
#'
#' @import ggplot2
#' @importFrom forecast ggseasonplot
#' @importFrom pracma movavg
#' @importFrom forecast ses hw holt forecast
#' @importFrom greybox MSE
#' @importFrom crayon green red yellow
#' @importFrom magrittr `%$%`
#'
#' @examples
#' serie_tiempo_rutina(sunspot.year,5)
#'
#' base=data.frame(tiempo=seq(Sys.Date(),by="days",length=20),valores=(rexp(50)+1)*sin(1:50))
#' serie_tiempo_rutina(datos=base,frecuencia=4,inicio=2010)
#'
serie_tiempo_rutina <- function(datos,frecuencia=NULL,inicio=NULL,validar_=FALSE,pausa_off=1,msg = TRUE){
if(validar_){
paquetes.tsrutina()
pausa()
conditional.tsrutina(datos)
}
if(is.ts(datos)){
elementos <- tratamiento.ts_set(datos)
datos <- elementos$data
frecuencia <- ifelse(is.null(frecuencia),elementos$frecu,frecuencia)
inicio <- ifelse(is.null(inicio),elementos$inicio,inicio)
}
#verificar si los elementos se ven bien
elementos <- checar_datos(datos,frecuencia,inicio)
datos <- elementos$datos
datosts <- elementos$datosts
reporte_data <- list()
#Gráficos para ver si es estacional
p <- ggplot(datos, aes(x,y)) +
geom_point()+geom_line()+
ggtitle("Serie de tiempo Visualización")+
theme(plot.title = element_text(color = "Black",hjust = 0.5)) -> reporte_data$plots$General
print(p)
pausa()
if(frecuencia==1){
message("La frecuencia de la serie de tiempo es 1,\nusaremos frecuencia 12 para los siguientes 2 graficos")
frecuencia <- 12
datosts <- ts(datos$y,frequency = 12,start = start(elementos$datosts))
}
p <- try(
forecast::autoplot(stl(datosts, s.window = "periodic"),
ts.colour="blue",main="Ruido + Estacionalidad + Tendencia + SerieTemporal "),
silent = T
) -> reporte_data$plots$descomposicion
print(p)
pausa()
p <- forecast::ggseasonplot(datosts,col=rainbow(length(datos$y)/frecuencia),year.labels=TRUE,xlab="Tiempo",
ylab="Serie de tiempo",main = "Grafico Estacional de la Serie Temp.") -> reporte_data$plots$ggseason
print(p)
pausa()
### Abajo remplazo en formato ggplot
# boxplot(datosts~cycle(datosts),xlab = "Frecuencias",ylab = "Valores",main="Boxplot por cada valor de la frecuencia")
if(frecuencia == 12){
p <- ggplot(fortify(datosts) %>%
dplyr::mutate(mes = months(Index,F),
mes_num = lubridate::month(Index))) +
geom_boxplot(aes(x = reorder(mes,mes_num), y = Data, color = mes)) -> reporte_data$plots$boxplots
}else {
minimo <- which.min(c(end(datosts)[1] - start(datosts)[1],frecuencia))
if(minimo == 1L){
datos_by <- datos %>%
mutate(year = lubridate::year(x))
p <- ggplot(datos_by,aes(x = year, y = y, color = year,group = year))+
geom_boxplot(show.legend = F) +
labs(main = "Boxplot por año",x ="Año",y="") -> reporte_data$plots$boxplots
print(p)
}else{
datos_by <- datos %>% arrange(x) %>%
mutate(grupo = as.numeric(cycle({datosts})))
p <- ggplot(datos_by,aes(x = grupo, y = y, color = grupo,group = grupo))+
geom_boxplot(show.legend = F) +
labs(main = "Boxplot por frecuencia",x ="Frecuencia",y="") -> reporte_data$plots$boxplots
print(p)
}
}
print(p)
pausa()
separador()
#Modelo de regresion lineal
datos$periodos <- 1:length(datos$x)
datos_rl <- lm(y~x, data=datos)
reporte_data$lm <- summary(datos_rl)
print(reporte_data$lm)
#Se puede ver cuales variables son significativas en el modelo
pausa()
separador()
# Grafico del modelo de regresion lineal
# plot(datos$periodos,datos$y,type = "l",
# xlab="Periodos",ylab="Valor de la serie",
# main="Regresión Lineal")
# lines(datos_rl$fitted.values, col="blue")
p <- ggplot(datos,aes(x=periodos, y=y)) +
geom_line() + geom_smooth(method = "lm",formula = y ~ x) +
labs(x = "Periodos", y = "Valor de la serie", title = "Regresión Lineal") -> reporte_data$smooth$lm
print(p)
pausa()
separador()
# Promedio Movil simple
promo <- pracma::movavg(datosts, n=frecuencia, type="s")
# plot(datos$periodos,datos$y,type = "l",
# xlab="Periodos",ylab="Valor de la serie",
# main="Promedio movil simple")
# lines(promo,col="blue")
p <- ggplot(datos,aes(x=periodos, y=y)) +
geom_line()+
geom_line(aes(y=promo),color = "#7171f5")+
labs(x="Periodos", y= "Valor de la serie", title = "Promedio movil simple") -> reporte_data$smooth$promedio_movil
print(p)
pausa()
separador()
#Promedio Movil ponderado
promopo <- pracma::movavg(datosts, n=frecuencia, type="w")
plot(datos$periodos,datos$y,type = "l",
xlab="Periodos",ylab="Valor de la serie",
main="Promedio movil ponderado")
lines(promopo,col="blue")
p <- ggplot(datos,aes(x=periodos, y=y)) +
geom_line()+
geom_line(aes(y=promopo),color = "#7171f5")+
labs(x="Periodos", y= "Valor de la serie", title = "Promedio movil Ponderado") -> reporte_data$smooth$promedio_movil_ponderado
print(p)
pausa()
separador()
#Exponential Smoothing
pesoses <- forecast::ses(datos$y)
summary(pesoses)
# plot(datos$periodos,datos$y,type = "l",
# xlab="Periodos",ylab="Valor de la serie",
# main="Suavizamiento Exponencial")
# lines(datos$periodos,pesoses$fitted, col="blue")
p <- ggplot(datos,aes(x=periodos, y=y)) +
geom_line()+
geom_line(aes(y=pesoses$fitted),color = "#7171f5",type = "dashes")+
labs(x="Periodos", y= "Valor de la serie", title = "Suavizamiento exponencial") -> reporte_data$smooth$smooth_exponencial
print(p)
pausa()
separador()
#Holt's Exponential Smoothing
pesoholt <- forecast::holt(datos$y)
summary(pesoholt)
# plot(datos$periodos,datos$y,type = "l",
# xlab="Periodos",ylab="Valor de la serie",
# main="Suavizamiento Exponencial holt")
# lines(datos$periodos,pesoholt$fitted, col="blue")
p <- ggplot(datos,aes(x=periodos, y=y)) +
geom_line()+
geom_line(aes(y=pesoholt$fitted),color = "#7171f5")+
labs(x="Periodos", y= "Valor de la serie", title = "Suavizamiento Exponencial de Holt") -> reporte_data$smooth$holt
print(p)
pausa()
separador()
#Holt-Winters' Exponential Smoothing
datos$Ajustadohw <- tryCatch({
pesohw <- forecast::hw(datosts)
summary(pesohw)
#asignamos valores ajustados a una columna
datos$Ajustadohw <- as.numeric(pesohw$fitted)
},
error = function(e) {
cat(crayon::red("El ajuste por Holt-Winter fracasó por el método 1 ¿desea intentar usar ajuste por el método 2?"))
otro_metodo <- readline("\nRespuesta [TRUE/FALSE]: ")
otro_metodo <- ifelse(otro_metodo=="",TRUE,otro_metodo)
if(as.logical(otro_metodo)){
mensaje("\nIniciando proceso Holt-Winter por metodo 2: \n")
pesohw <- stats::HoltWinters(datosts)
print(pesohw)
datos$Ajustadohw <- pesohw$fitted %>% as.data.frame() %$% c(rep(NA,frecuencia),xhat)
return(datos$Ajustadohw)
}
})
#
# # plot(datos$periodos,datos$y,type = "l",
# # xlab="Periodos",ylab="Valor de la serie",
# # main="Suavizamiento Exponencial holt-winter")
# # lines(datos$periodos, datos$Ajustadohw, col="red", type="l")
# # lines(datos$periodos, pesohw$fitted, col="blue", type="l")
p <- ggplot(datos,aes(x=periodos, y=y)) +
geom_line()+
geom_line(aes(y=Ajustadohw),color = "#7171f5")+
labs(x="Periodos", y= "Valor de la serie", title = "Suavizamiento Exponencial Holt-Winter") -> reporte_data$smooth$holt_winter
print(p)
pausa()
#Ahora, para elegir el mejor modelo de "suavizamiento",
#usaremos el MSE (error cuadratico medio).
separador()
cat(crayon::green("\n MSE del modelo de Regresion Lineal \n"))
cat(MSE(datos$y, datos_rl$fitted.values))
cat(crayon::green('\n MSE del modelo de Promedios Moviles \n'))
cat(crayon::green('\n \t Simple \n\t'))
cat(MSE(datos$y, promo))
cat(crayon::green('\n \t Ponderado \n\t'))
cat(MSE(datos$y, promopo))
cat(crayon::green('\n MSE del modelo de Exponential Smoothing \n'))
cat(MSE(datos$y, pesoses$fitted))
cat(crayon::green('\n MSE del modelo de Holt \n'))
cat(MSE(datos$y, pesoholt$fitted))
cat(crayon::green('\n MSE del modelo Holt-Winters \n'))
cat(MSE(datos$y, datos$Ajustado))
cat(crayon::yellow('\n El error minimo se tiene con \n'))
a <- which.min(c(MSE(datos$y, datos_rl$fitted.values),
MSE(datos$y, promo),
MSE(datos$y, promopo),
MSE(datos$y, pesoses$fitted),
MSE(datos$y, pesoholt$fitted),
MSE(datos$y, datos$Ajustado)))
if(msg){#si no hay mensaje entonces predeterminado
h_pronostico <- readline('Inserte el número de periodos a pronosticar: \n(intro para 12 periodos) \n')
h_pronostico <- ifelse(h_pronostico != "",as.integer(h_pronostico),12)
}else{
h_pronostico <- 12L
}
switch(a,
'1' = {cat(crayon::green('Regresion lineal'))
pausa()
pronostico <- forecast::forecast(datos_rl$fitted.values,
h=h_pronostico,level=c(80,95))
},
'2' = {cat(crayon::green('Promedio movil simple'))
pausa()
pronostico <- forecast::forecast(promo,h=h_pronostico,level=c(80,95))
},
'3' = {cat(crayon::green('Promedio ponderado'))
pausa()
pronostico <- forecast::forecast(promopo,h=h_pronostico,level=c(80,95))
},
'4' = {cat(crayon::green('Exponencial simple'))
pausa()
pronostico <- forecast::ses(datos$y,h=h_pronostico,level=c(80,95))
},
'5' = {cat(crayon::green('Suavizamiento de Holt'))
pausa()
pronostico <- forecast::holt(datos$y,h=h_pronostico,level=c(80,95))
},
'6' = {cat(crayon::green('Suavizamiento de Holt-Winter'))
pausa()
pronostico <- tryCatch({
pesohw <- forecast::hw(datosts,h = h_pronostico,level=c(80,95))
pesohw
},
error = function(e) {
pesohw <- stats::HoltWinters(datosts)
pesohw_forecast <- stats:::predict.HoltWinters(pesohw,n.ahead = h_pronostico,prediction.interval = T,level = 0.95)
return(pesohw_forecast)
})
}
)
# plot(pronostico)
reporte_data$pronosticos$modelo <- pronostico
p <- autoplot(pronostico) +
lims(x = c(end(pronostico$x)[1]-1,NA))+
labs(title = paste0("Ajuste por ",pronostico$method)) -> reporte_data$pronosticos$plot
print(p)
separador()
pausa()
cat(crayon::green('Supuesto de Normalidad'))
prueba <- shapiro.test(pronostico$residuals)
print(prueba)
if(prueba$p.value>0.05){
message("No se puede rechazar H0:Hay Normalidad")
}else{
message("Se rechaza H0, se obta por H1: No hay normalidad")
}
separador()
return(invisible(reporte_data))
}
#' Depurador de dispositivos graficos
#'
#' Depura el espacio de dispositivos graficos atravez de dev.off() y dev.new(). Checando dev.list().
#'
#' Equivalente a un reinicio de la dev.list
#'
#' @return NULL return
#' @export
#' @encoding UTF-8
#'
#' @examples
#' dev.TRS()
#' #Se limpia los datos, consulte dev.list()
dev.TRS <- function(){
n <- length(dev.list())
if(n < 2){
dev.new()
dev.TRS()
}
if(n > 2){
for (i in 1:(n-2)){
dev.off(i+3)
}
}
}
#' Grafica de ajuste de TimeSeries en PNG (low computational cost)
#'
#' Obtén archivos .png de tus gráficos de serie de tiempo en tu directorio
#'
#' Se usa una salida interactiva para el proceso
#'
#' @param datos Dataframe de no mas de 2 columnas, en el orden primero tiempo y
#' luego valor el tiempo va en formato fecha, y tiene que ser en el orden
#' dia-mes-year. La versión 2.1 soporta objetos de la clase TS (time series)
#' @param frecuencia Este es el periodo de la serie, trimestral = 3, cuatrimestral = 4
#' ,mensual = 12, etc.
#' @param inicio Este es el year a iniciar la serie de tiempo
#' @param validar_ Boleano, True/False indica si se vericarán los datos
#' @param img_width Ancho de la imagen a guardar en pixeles
#' @param img_height Altura de la imagen a guardar en pixeles
#'
#' @return La salida no es como tal un objeto, si no una serie de impresiones de varios
#' analisis.
#' \itemize{\item{\bold{Plost}}{ Arroja una lista de plots que ayudan a ver el comportamiento de la serie y como ciertos ajustes se aproximan mejor a ella}}
#' \itemize{\item{\bold{PNG}}{ Imagenes grabadas en el directorio de trabajo}}
#'
#' @author Cristhian Diaz
#'
#' @export
#' @encoding UTF-8
#'
#' @import ggplot2
#' @importFrom forecast ggseasonplot seasonplot
#' @importFrom pracma movavg
#' @importFrom forecast ses hw holt forecast
#' @importFrom greybox MSE
#' @importFrom crayon green red yellow
#'
#'
#' @examples
#' serie_tiempo_plots(sunspot.year,15)
#'
#' base=data.frame(tiempo=seq(Sys.Date(),by="days",length=20),valores=(rexp(50)+1)*sin(1:50))
#' serie_tiempo_plots(datos=base,frecuencia=4,inicio=2010)
#'
serie_tiempo_plots <- function(datos,frecuencia=NULL,inicio=NULL,
validar_=FALSE,pausa_off=1,
img_width = 1024, img_height = 600){
if(validar_){
paquetes.tsrutina()
pausa()
conditional.tsrutina(datos)
}
if(is.ts(datos)){
elementos <- tratamiento.ts_set(datos)
datos <- elementos$data
frecuencia <- ifelse(is.null(frecuencia),elementos$frecu,frecuencia)
inicio <- ifelse(is.null(inicio),elementos$inicio,inicio)
}
#verificar si los elementos se ven bien
elementos <- checar_datos(datos,frecuencia,inicio)
datos <- elementos$datos
datosts <- elementos$datosts
#tratamiento para errores en mostrar graficos
dev.TRS()
#Graficos para ver si es estacional
png("serie_de_tiempo.png",width = img_width,height = img_height,units = "px") #se guarda una imagen, esta se plotea luego luego
print(
ggplot(datos, aes(x, y)) +
geom_point()+geom_line()+
ggtitle("Serie de tiempo Visualización")+
theme(plot.title = element_text(color = "Black",hjust = 0.5))
)
dev.off(4)
cat(crayon::green("Serie de tiempo serie_de_tiempo.png fue creada y guardada \n"))
message("¿Deseas ver en el grafico o seguir generando los siguientes graficos? \n")
imprime <- readline(" [Sí/Intro]:")
if(imprime %in% c("si","Sí","SI","yes","YES","Si","Yes","s","1")){
cat(crayon::yellow("Graficando..."))
dev.TRS() #Regulador de device graphics
print(
ggplot(datos, aes(x, y)) +
geom_point()+geom_line()+
ggtitle("Serie de tiempo Visualización")+
theme(plot.title = element_text(color = "Black",hjust = 0.5))
)
pausa()
}
png("st_descomposicion.png",width = img_width,height = img_height,units = "px") #se guarda una imagen, esta se plotea luego luego
print(
forecast::autoplot(stl(datosts, s.window = "periodic"), ts.colour="blue",
main = "Descomposición de Serie de tiempo")
)
cat(crayon::green("Descomposición de Serie de tiempo st_descomposicion.png fue creada y guardada \n"))
dev.off(4)
message("¿Deseas ver en el grafico o seguir generando los siguientes graficos? \n")
imprime <- readline(" [Sí/Intro]:")
if(imprime %in% c("si","Sí","SI","yes","YES","Si","Yes","s","1")){
cat(crayon::yellow("Graficando..."))
print(
forecast::autoplot(stl(datosts, s.window = "periodic"), ts.colour="blue",
main = "Descomposición de Serie de tiempo")
)
pausa()
}
png("st_seasonplot.png",width = img_width,height = img_height,units = "px") #se guarda una imagen, esta se plotea luego luego
print(forecast::seasonplot(datosts,col=rainbow(length(datos$y)/frecuencia),year.labels=TRUE,
xlab="Periodos",ylab="Serie de tiempo",
main = "Grafico Estacional de la serie de tiempo"))
cat(crayon::green("Temporalidad de Serie de tiempo st_seasonplot.png fue creada y guardada \n"))
dev.off()
message("¿Deseas ver en el grafico o seguir generando los siguientes graficos? \n")
imprime <- readline(" [Sí/Intro]:")
if(imprime %in% c("si","Sí","SI","yes","YES","Si","Yes","s","1")){
cat(crayon::yellow("Graficando..."))
print(
forecast::seasonplot(datosts,col=rainbow(length(datos$y)/frecuencia),year.labels=TRUE,
xlab="Periodos",ylab="Serie de tiempo",
main = "Grafico Estacional de la serie de tiempo")
)
pausa()
}
#Se aprecia que no hay estacionalidad
png("st_boxes_plot.png",width = img_width,height = img_height,units = "px") #se guarda una imagen, esta se plotea luego luego
print(boxplot(datosts~cycle(datosts),xlab = "Frecuencias",ylab = "Valores",main="Boxplot por cada valor de la frecuencia"))
cat(crayon::green("Temporalidad de Serie de tiempo st_boxes_plot.png fue creada y guardada \n"))
dev.off()
message("¿Deseas ver el grafico o seguir generando los siguientes graficos? \n")
imprime = readline(" [Sí/Intro]:")
if(imprime %in% c("si","Sí","SI","yes","YES","Si","Yes","s","1")){
cat(crayon::yellow("Graficando..."))
print(
boxplot(datosts~cycle(datosts),xlab = "Frecuencias",ylab = "Valores",main="Boxplot por cada valor de la frecuencia")
)
pausa()
}
#Modelo de regresion lineal
datos$periodos <- 1:length(datos$x)
datos_rl <- lm(y~periodos, data=datos)
#Grafico del modelo de regresion lineal
png("st_regresion_lineal.png",width = img_width,height = img_height,units = "px") #se guarda una imagen, esta se plotea luego luego
datos$lineal <- datos_rl$fitted.values
print(
ggplot(datos, aes(x,y)) +
geom_point()+geom_line()+
ggtitle("Serie de tiempo + Regresión lineal")+
theme(plot.title = element_text(color = "Black",hjust = 0.5))+
geom_line(aes(x,lineal), col="blue")
)
cat(crayon::green("Ajuste de regresión lineal st_regresion_lineal.png fue creada y guardada \n"))
dev.off()
message("¿Deseas ver en el grafico o seguir generando los siguientes graficos? \n")
imprime <- readline(" [Sí/Intro]:")
if(imprime %in% c("si","Sí","SI","yes","YES","Si","Yes","s","1")){
cat(crayon::yellow("Graficando..."))
print(
ggplot(datos, aes(x,y)) +
geom_point()+geom_line()+
ggtitle("Serie de tiempo + Regresión lineal")+
theme(plot.title = element_text(color = "Black",hjust = 0.5))+
geom_line(aes(x,lineal), col="blue")
)
pausa()
}
#abline(pesorl,col="blue")
#Promedio Movil simple
promo <- pracma::movavg(datosts, n=frecuencia, type="s")
png("st_prom_movil_simple.png",width = img_width,height = img_height,units = "px") #se guarda una imagen, esta se plotea luego luego
datos$promo <- promo
print(
ggplot(datos, aes(x,y)) +
geom_point()+geom_line()+
ggtitle("Serie de tiempo + Promedio movil simple")+
theme(plot.title = element_text(color = "Black",hjust = 0.5))+
geom_line(aes(x,promo), col="blue")
)
cat(crayon::green("Promedio Movil simple st_prom_movil_simple.png fue creada y guardada \n"))
dev.off()
message("¿Deseas ver en el grafico o seguir generando los siguientes graficos? \n")
imprime <- readline(" [Sí/Intro]:")
if(imprime %in% c("si","Sí","SI","yes","YES","Si","Yes","s","1")){
cat(crayon::yellow("Graficando..."))
print(
ggplot(datos, aes(x,y)) +
geom_point()+geom_line()+
ggtitle("Serie de tiempo + Promedio movil simple")+
theme(plot.title = element_text(color = "Black",hjust = 0.5))+
geom_line(aes(x,promo), col="blue")
)
pausa()
}
#Promedio Movil ponderado
promopo <- pracma::movavg(datosts, n=frecuencia, type="w")
png("st_prom_movil_ponderado.png",width = img_width,height = img_height,units = "px") #se guarda una imagen, esta se plotea luego luego
datos$promopo <- promopo
print(
ggplot(datos, aes(x,y)) +
geom_point()+geom_line()+
ggtitle("Serie de tiempo + Promedio ponderado")+
theme(plot.title = element_text(color = "Black",hjust = 0.5))+
geom_line(aes(x,promopo), col="blue")
)
cat(crayon::green("Promedio Movil ponderado st_prom_movil_ponderado.png fue creada y guardada \n"))
dev.off()
message("¿Deseas ver en el grafico o seguir generando los siguientes graficos? \n")
imprime <- readline(" [Sí/Intro]:")
if(imprime %in% c("si","Sí","SI","yes","YES","Si","Yes","s","1")){
cat(crayon::yellow("Graficando..."))
print(
ggplot(datos, aes(x,y)) +
geom_point()+geom_line()+
ggtitle("Serie de tiempo + Promedio ponderado")+
theme(plot.title = element_text(color = "Black",hjust = 0.5))+
geom_line(aes(x,promopo), col="blue")
)
pausa()
}
#Exponential Smoothing
pesoses <- forecast::ses(datos$y)
#summary(pesoses)
png("st_ajuste_exponencial.png",width = img_width,height = img_height,units = "px") #se guarda una imagen, esta se plotea luego luego
datos$sess <- pesoses$fitted
print(
ggplot(datos, aes(x,y)) +
geom_point()+geom_line()+
ggtitle("Serie de tiempo + Ajuste exponencial")+
theme(plot.title = element_text(color = "Black",hjust = 0.5))+
geom_line(aes(x,sess), col="blue")
)
cat(crayon::green("Suavizamiento Exponencial st_ajuste_exponencial.png fue creada y guardada \n"))
dev.off()
message("¿Deseas ver en el grafico o seguir generando los siguientes graficos? \n")
imprime <- readline(" [Sí/Intro]:")
if(imprime %in% c("si","Sí","SI","yes","YES","Si","Yes","s","1")){
cat(crayon::yellow("Graficando..."))
print(
ggplot(datos, aes(x,y)) +
geom_point()+geom_line()+
ggtitle("Serie de tiempo + Ajuste exponencial")+
theme(plot.title = element_text(color = "Black",hjust = 0.5))+
geom_line(aes(x,sess), col="blue")
)
pausa()
}
#Holt's Exponential Smoothing
pesoholt <- forecast::holt(datos$y)
#summary(pesoholt)
png("st_holt_exponencial.png",width = img_width,height = img_height,units = "px") #se guarda una imagen, esta se plotea luego luego
datos$holtt <- pesoholt$fitted
print(
ggplot(datos, aes(x,y)) +
geom_point()+geom_line()+
ggtitle("Serie de tiempo + Ajuste de Holt")+
theme(plot.title = element_text(color = "Black",hjust = 0.5))+
geom_line(aes(x,holtt), col="blue")
)
cat(crayon::green("Suavizamiento Holt st_holt_exponencial.png fue creada y guardada \n"))
dev.off()
message("¿Deseas ver en el grafico o seguir generando los siguientes graficos? \n")
imprime <- readline(" [Sí/Intro]:")
if(imprime %in% c("si","Sí","SI","yes","YES","Si","Yes","s","1")){
cat(crayon::yellow("Graficando..."))
print(
ggplot(datos, aes(x,y)) +
geom_point()+geom_line()+
ggtitle("Serie de tiempo + Ajuste de Holt")+
theme(plot.title = element_text(color = "Black",hjust = 0.5))+
geom_line(aes(x,holtt), col="blue")
)
pausa()
}
#Holt-Winters' Exponential Smoothing
pesohw <- forecast::hw(datosts,seasonal = "additive",h=10,level = 95)
#pesohw$model
#asignamos valores ajustados a una columna
png("st_holt-winter_exponencial.png",width = img_width,height = img_height,units = "px") #se guarda una imagen, esta se plotea luego luego
datos$Ajustado <- as.numeric(pesohw$fitted)
print(
ggplot(datos, aes(x,y)) +
geom_point()+geom_line()+
ggtitle("Serie de tiempo + Holt winter")+
theme(plot.title = element_text(color = "Black",hjust = 0.5))+
geom_line(aes(x,Ajustado), col="blue")
)
cat(crayon::green("Suavizamiento Holt-Winters st_holt-winter_exponencial.png fue creada y guardada \n"))
dev.off()
message("¿Deseas ver en el grafico o seguir generando los siguientes graficos? \n")
imprime <- readline(" [Sí/Intro]:")
if(imprime %in% c("si","Sí","SI","yes","YES","Si","Yes","s","1")){
cat(crayon::yellow("Graficando..."))
print(
ggplot(datos, aes(x,y)) +
geom_point()+geom_line()+
ggtitle("Serie de tiempo + Holt Winter")+
theme(plot.title = element_text(color = "Black",hjust = 0.5))+
geom_line(aes(x,Ajustado), col="blue")
)
pausa()
}
a <- which.min(c(MSE(datos$y, datos_rl$fitted.values),
MSE(datos$y, promo),
MSE(datos$y, promopo),
MSE(datos$y, pesoses$fitted),
MSE(datos$y, pesoholt$fitted),
MSE(datos$y, datos$Ajustado)))
cat(crayon::green('El ajuste con menor MSE es: '))
switch(a,
'1' = {cat(crayon::yellow('Regresion lineal'))
pausa()
pronosticado<-datos_rl$fitted.values},
'2' = {cat(crayon::yellow('Promedio movil simple'))
pausa()
pronosticado<-promo},
'3' = {cat(crayon::yellow('Promedio ponderado'))
pausa()
pronosticado<-promopo},
'4' = {cat(crayon::yellow('Exponencial simple'))
pausa()
pronosticado<-pesoses$model},
'5' = {cat(crayon::yellow('Ajuste de Holt'))
pausa()
pronosticado<-pesoholt$model},
'6' = {cat(crayon::yellow('Holt-Winter'))
pausa()
pronosticado<-pesohw$model}
)
pausa()
#Promedio Movil ponderado
pronostico <- forecast::forecast(pronosticado,h=20,level=c(80,95))
png("st_pronostico_ponderado_80-90.png",width = img_width,height = img_height,units = "px") #se guarda una imagen, esta se plotea luego luego
print(
forecast::autoplot(pronostico)
)
cat(crayon::green("Pronostico para el mejor ajuste st_pronostico_ponderado_80-90.png fue creada y guardada \n"))
dev.off()
message("¿Deseas ver en el grafico o seguir generando los siguientes graficos? \n")
imprime <- readline(" [Sí/Intro]:")
if(imprime %in% c("si","Sí","SI","yes","YES","Si","Yes","s","1")){
cat(crayon::yellow("Graficando..."))
print(
forecast::autoplot(pronostico)
)
}
}
#' Recomendaciones a un grafico ACF o PACF
#'
#' Obten un recomendación a tu ajuste de series de tiempo ACF o PACF
#'
#' El valor IC se calcula siguiendo la metodologia de la paqueteria Stats en a función acf(). Se retorna el valor positivo talque el intervalo se forma con (IC,-IC)
#'
#' @param objeto_acf objeto ACF o PACF
#' @param print_IC Indicador True/False para mostrar valor absoluto de intervalo de confianza ver details
#'
#' @return Valor del ultimo lag significativo de la funcion de autocorrelacion
#' @export
#' @encoding UTF-8
#'
#' @importFrom crayon green red yellow
#'
#' @details El intervalo de confianza (IC) de un gráfico ACF/PACF son las bandas a partir de las que una auto-correlación se considera significativas. Cuando las autocorrelaciones para cada lag resultan no significativas, se dice que no hay auto-correlación
#'
#' @examples
#' base=data.frame(x=seq(Sys.Date(),by="days",length=200),y=(rexp(50)+1)*sin(1:50))
#' recomendacion_autocorrelaciones(acf(base$y,plot = FALSE))
#'
recomendacion_autocorrelaciones <- function(objeto_acf,print_IC=FALSE) {
llamada <- match.call()
ruta <- match(c("objeto_acf"),names(llamada))
if(ruta!=2){#chequeo de que se agregaron bien los parametros
#stopifnot(ruta!=2)
message("Objeto_acf no encontrado o hay más de un parametro en la función")
stop()
}
#Zona de plot not True
a <- llamada[[2]]
a <- as.character(a)
if(length(a)!=3 | as.logical(a[3])){
message("EL objeto debe ser un ACF o PACF con parametro plot = FALSE")
message("Ver el ejemplo en la documentación")
?recomendacion_autocorrelaciones
}
#si la salida es un vector de 3 elementos entonces hay dos parametros
#serie=tryCatch(get(objeto_acf$series),error= function(e){message(e," \nSe busca otra entrada..."); return(NULL)})
# serie=get("base",envir = parent.frame())
#se corrigio el uso de envir
serie <- tryCatch(expr = get("base",envir = parent.frame()), error = function(e) {
message("\n Se intentara conectar con el objeto de llamado: ",
as.character(llamada[[2]][[2]]),"\n")#llamada[[2]][[2]] contempla que hay un objeto a llamar dentro de acf(*objeto)
get(as.character(llamada[[2]][[2]]),envir = parent.frame())
})
#en caso de serie NULL
# if(is.null(serie)){
# vec=strsplit(a[2],split = "$",fixed = TRUE)
# message("\nSe encontro la base de datos llamada: ",vec[[1]][1])
# message("\nDentro de ella se encontro el vector llamado: ",vec[[1]][2])
# serie=eval(str2lang(a[2]))
# cat("\nLos primeros 6 valores de este vector son:\n")
# print(head(serie))
# }
try({serie <- ts(serie$y)},silent = 1)
stopifnot(class(serie) == class(ts()))
order_ <- NULL
if(objeto_acf$type=="partial"){
matriz <- matriz_eacf(serie,ar.max = 1, ma.max = 15,
print_matrix = FALSE)
matriz <- matriz$symbol=="o"
for(i in 1:15){
if(matriz[1,i]==1){
order_ <- i
cat(crayon::yellow("\nProponemos MA(q) con q="),order_)
break
}
}
#parche para casos donde no hay valor a proponer
if(is.null(order_)){
cat(crayon::green("El valor de la q propuesta es mayor a 15..."))
order_ <- 16
}
}
if(objeto_acf$type=="correlation"){
matriz=matriz_eacf(serie,ar.max = 15,ma.max = 1,print_matrix = FALSE)
matriz=matriz$symbol=="o"
for(i in 1:15){
if(matriz[i,1]==1){
order_ <- i
cat(crayon::yellow("\nProponemos AR(p) con p="),order_)
break
}
}
if(is.null(order_)){
cat(crayon::green("El valor de la p propuesta es mayor a 15..."))
order_ <- 16
}
}
#obtener los intervalos de confianza dando el objeto
IC <- intervalo_confianza_acf(objeto_acf)
#mayores=abs(objeto_acf$acf)>IC
#cat("\nLos siguientes elementos son propuestas de r: ")
#posibles_lags=objeto_acf$lag[mayores]
#cat(posibles_lags)
#cat("\nProponemos que r sea:",posibles_lags[length(posibles_lags)])
if(print_IC){
cat("\nEl IC de modelo es: ",IC)
}
return(invisible(order_))
}
#' Recomendación de modelo ARMA
#'
#' Usando una matriz eacf de TSA paqueteria se propone un posible vector con los valores de p y q de un modelos ARMA(p,q)
#'
#' Se utiliza un metodo de busqueda de esquinas para proponer el valor de
#' p,q para ARMA(p,q).
#'
#' @param time_series Objeto Serie de tiempo
#' @param print_matrix Indicador, imprimir o no matriz de eacf
#' @param msg Mostrar mensaje en lugar de vector
#'
#' @return Vector de longitud 2, primera entrada valor de \bold{p}, segunda valor de \bold{q}
#' @export
#' @encoding UTF-8
#'
#' @importFrom crayon green red yellow
#'
#' @examples
#' recomendaciones_arma(AirPassengers)
recomendaciones_arma <- function(time_series,print_matrix=TRUE,msg=FALSE) {
x <- time_series
modelo_arma <- matriz_eacf(x,7,7,print_matrix)
matriz_true_false <- modelo_arma$symbol=="o"
matriz_true_false=matriz_true_false[-1,][,-1]
for (i in 1:7) {
if(sum(matriz_true_false[,i])>0){
for (j in 1:7) {
zz <- matriz_true_false[,i]
if(zz[j]==1){
#se analizan vecinos
izquierda <- matriz_true_false[j,i+1]
abajo <- matriz_true_false[j+1,i]
diagonal <- matriz_true_false[j+1,i+1]
#condicion algun vecino o diagonal no null
if(izquierda+abajo+diagonal >2 | diagonal>0){
vec <- c(i,j)
if(msg){
cat(crayon::yellow("Se recomienda ARMA("),vec[1],",",vec[2],crayon::yellow(")"))
return(invisible(NULL))
}
return(vec)
}
}
}
}
}
}
#' Ajuste de Modelo ARIMA a tu Serie de tiempo
#'
#' Ajusta un modelo arima usando un asistente que te mostrar si tu serie se ajusta a un modelo ARIMA en particular
#'
#' A diferencia de otras funciones esta utiliza el conocimiento del experto para ajustar el mejor modelo ARIMA a la serie de tiempo estudiada, mientras propone modelos que el asistente considera utiles.
#'
#' @param datos Dataframe de no mas de 2 columnas, en el orden primero tiempo y
#' luego valor el tiempo va en formato fecha, y tiene que ser en el orden
#' dia-mes-year. La versión 2.1 soporta objetos de la clase TS (time series)
#' @param frecuencia Este es el periodo de la serie, trimestral = 3, cuatrimestral = 4
#' ,mensual = 12, etc.
#' @param inicio Este es el year a iniciar la serie de tiempo
#' @param validar_ Boleano, True/False indica si se vericarán los datos
#'
#' @return La salida no es como tal un objeto, si no una serie de impresiones de varios
#' analisis. El mejor basando en criterio AIC.
#' @export
#' @encoding UTF-8
#'
#' @importFrom crayon green red yellow
#'
#' @examples
#' serie_tiempo_ARIMA(sunspot.year,5)
#'
#' base=data.frame(tiempo=seq(Sys.Date(),by="days",length=50),valores=(rexp(50)+1)*sin(1:50))
#' serie_tiempo_ARIMA(datos=base,frecuencia=4,inicio=2010)
#'
#' serie_tiempo_ARIMA(wineind)
#'
serie_tiempo_ARIMA<-function(datos,frecuencia=NULL,inicio=NULL,
validar_=FALSE,msg=TRUE,pausa_off=1){
if(validar_){
paquetes.tsrutina()
pausa()
conditional.tsrutina(datos)
}
if(is.ts(datos)){
elementos <- tratamiento.ts_set(datos)
datos <- elementos$data
frecuencia <- ifelse(is.null(frecuencia),elementos$frecu,frecuencia)
inicio <- ifelse(is.null(inicio),elementos$inicio,inicio)
}
if(!is.data.frame(datos)){
stop("El objeto debe ser un data frame con dos elementos")
}
#verificar si los elementos se ven bien
elementos <- checar_datos(datos,frecuencia,inicio)
datos <- elementos$datos
#datosts <- elementos$datosts #No se usa en la rutina
separador()
base <- datos
ban <- TRUE
numero_diferenciaciones <- 0
while(ban){
prueba<-tseries::adf.test(base$y)
print(prueba)
if(msg){
p_valor<-readline('\nInserte un p valor (intro para p=0.05): ')
}else{
p_valor<-0.05
}
if(p_valor==""){
p_valor<-0.05
}else{
cat("\n",crayon::red(sprintf("El valor de p= %s",p_valor)),"\n")
p_valor <- as.numeric(p_valor)
}
if(prueba$p.value>p_valor){
cat(crayon::yellow("\nNo se puede rechazar H0:Hay presencia de una raiz unitaria\n"))
cat(crayon::yellow("\nNo es estacionaria"))
pausa()
differenciado <- diff(base$y,lag = 1,differences = 1)
base <- base[-nrow(base),]
base$y <- differenciado
numero_diferenciaciones <- numero_diferenciaciones+1 #contador de diferenciaciones
cat(crayon::red('\nSe ha diferencio la base de datos para obtener estacionariedad\n'))
}else{
cat(crayon::yellow("\nSe rechaza H0, se obta por H1: La serie de tiempo es Estacionaria\n"))
ban <- FALSE
}
pausa()
}
#plotea el acf y analizas
acf(base$y,main="Autocorrelación, Analiza el valor de p en AR(p)")
#función de recomendación
rec <- recomendacion_autocorrelaciones(acf(base$y,plot = FALSE))
if(msg){
ar<-readline('Que AR(p) sospechas?, inserte el valor de p: ')
}else{
ar<-""
}
ar<-if(ar==""){
c(as.numeric(rec),numero_diferenciaciones,0)
}else{
c(as.numeric(ar),numero_diferenciaciones,0)
}
pausa()
#plotea el pacf
pacf(base$y,main="Autocorrelación Parcial,Analiza el valor de q en MA(q)")
rec <- recomendacion_autocorrelaciones(pacf(base$y,plot = FALSE))
if(msg){
ma <- readline('Que MA(q) sospechas?, inserte el valor de q: ')
}else{
ma <- ""
}
ma <- if(ma==""){
c(0,numero_diferenciaciones,as.numeric(rec))
}else{
c(0,numero_diferenciaciones,as.numeric(ma))
}
pausa()
#imprime arimas
cat(crayon::green("\nResumen estadistico de modelo AR(",ar[1],")"))
print(arima(base$y,order = ar))
pausa()
cat(crayon::green("\nResumen estadistico de modelo MA(",ma[3],")"))
print(arima(base$y,order = ma))
pausa()
#compara arimas
mama <- arima(base$y,order = ma)
rara <- arima(base$y,order = ar)
if(mama$aic<rara$aic){
cat(crayon::yellow(sprintf('\nEl modelo con menor AIC es el MA(%s)',ma[3])))
}else{
cat(crayon::yellow(sprintf('\nEl modelo con menor AIC es el AR(%s)',ar[1])))
}
pausa()
cat(crayon::green("\nPrueba de Box-Pierce and Ljung-Box Test"))
for (i in 1:2) {
modelo <- list(mama,rara)
if(modelo[[i]][["call"]][["order"]]=="ma"){
cat(crayon::green(sprintf("\nAnalisis de correlación en el modelo para Ma(%s)",ma[3])))
}else{
cat(crayon::green(sprintf("\nAnalisis de correlación en el modelo para AR(%s)",ar[1])))
}
box_test <- Box.test(modelo[[i]]$residuals, type ="Ljung-Box")
print(box_test)
cat(crayon::yellow("\nBox.test(), el p_valor > 0.05 entonces no hay correlacion ruido blanco"))
if(msg){
p_valor<-readline('Inserte un p valor, (intro para p=0.05): \n')
}else{
p_valor<-0.05
}
if(p_valor==""){
p_valor <- 0.05
}else{
cat("\n",crayon::red(sprintf("El valor de p= %s",p_valor)))
p_valor <- as.numeric(p_valor)
}
if(prueba$p.value>0.05){
cat(crayon::yellow("\nNo se puede rechazar H0 = No hay presencia de autocorrelacion "))
}else{
cat(crayon::yellow("\nSe rechaza H0, se obta por H1 = Hay presencia de autocorrelacion"))
}
sprintf("\n \n")
}
pausa()
#analisis ARMA si diferecias es mayor a cero, ARMA=ARIMA
datosts <- ts(base$y)
rec <- recomendaciones_arma(datosts) #checar como el objeto obtiene los symbol de la lista
cat(crayon::red("\nSe recomienda el modelo ARMA(p,q) con p=",rec[1]," q=",rec[2]))
if(msg){
cat("\nQue ARMA(p,q) sospechas?, inserte el valor de p,q separado por comas: ")
arma_pq<-readline('Ejemplo: 3,4 \t')
}else{
arma_pq<-""
}
arma_order<-if(arma_pq==""){
c(as.numeric(rec[1]),numero_diferenciaciones,as.numeric(rec[2]))
}else{
#dividir la entrada
w <- strsplit(arma_pq,",")
arma_p <- as.numeric(w[[1]][1])
arma_q <- as.numeric(w[[1]][2])
c(arma_p,numero_diferenciaciones,arma_q)
}
#se imprime el modelo arma
if(numero_diferenciaciones>0){
cat(crayon::red("\nSe detecto que la base a sido diferenciada al menos una vez\n
El modelo ARIMA final es de la forma (",arma_order,")"))
}
cat(crayon::green("\nA continuación su resumen estadistico"))
print(arima(base$y,order = arma_order))
pausa()
#ajuste ARIMA drif
return(NULL)
}
#' Rutina Principal de TSRutina
#'
#' Realiza todas las rutinas de la paqueteria TSRutina, ver detalles para más información
#'
#' Se hace un ajuste por metodo no ARIMA y se determina cual de estos modelos es el de menor
#' Error en los residuos. Ver \code{\link{serie_tiempo_rutina}}
#'
#' Se realiza las rutinas que verifican si la serie es o no estacionaria, si existe correlación
#' , luego se ajusta un modelo ARIMA y se guardan los graficos importantes en el directorio
#'
#' @param datos Datos para el analisis de serie de tiempo suport (data.frame and ts class)
#' @param frecuencia Frecuencia de la serie de tiempo, sirve para
#' reescribir la frecuencia cuando datos es un objeto ts
#' @param inicio Inicio de la serie de tiempo, igual que frecuencia
#' sobreescribe valores de objetos ts
#' @param validar_ (True or False) validar los datos
#' @param ... Not work
#'
#' @return La salida no es como tal un objeto, si no una serie de impresiones de varios
#' analisis. La siguiente lista detalla alguno de ellos:
#' \itemize{\item{\bold{Plost}}{ Arroja una lista de plots que ayudan a ver el comportamiento de la serie y como ciertos ajustes se aproximan mejor a ella}}
#' \itemize{\item{\bold{Resumenes}}{ Arroja ciertos resumenes de ciertos ajustes o pruebas que se hacen}}
#' \itemize{\item{\bold{Modelo}}{ Modelo con el menor MSE(Error cuadratico medio)}}
#'
#' @export
#' @encoding UTF-8
#'
#' @importFrom crayon green red yellow
#'
#'
#' @examples
#' base=data.frame(tiempo=seq(Sys.Date(),by="days",length=20),valores=(rexp(50)+1)*sin(1:50))
#' init(datos=base,frecuencia=4,inicio=2010)
init <- function(datos,frecuencia=NULL,inicio=NULL,validar_=FALSE,
msg=TRUE,pausa_off=1,...){
# Parte de validar integrada directamente en init() function
paquetes.tsrutina()
pausa()
conditional.tsrutina(datos)
if(is.ts(datos)){
elementos <- tratamiento.ts_set(datos)
datos <- elementos$data
frecuencia <- ifelse(is.null(frecuencia),elementos$frecu,frecuencia)
inicio <- ifelse(is.null(inicio),elementos$inicio,inicio)
}
cat(crayon::red("\n Inicio de rutina para tratamiento de una Serie de tiempo \n"))
serie_tiempo_rutina(datos = datos,frecuencia = frecuencia,
inicio = inicio,validar_ = validar_,
pausa_off = pausa_off,...)
cat(crayon::red("\n Inicio de pruebas para tratamiento de una Serie de tiempo \n"))
serie_tiempo_pruebas(datos = datos,frecuencia = frecuencia,
validar_ = validar_,
msg = msg,pausa_off = pausa_off)
cat(crayon::red("\n Ajuste de un modelo ARIMA para tratamiento de una Serie de tiempo \n"))
serie_tiempo_ARIMA(datos = datos,frecuencia = frecuencia,inicio = inicio,validar_ = validar_,msg = msg,pausa_off = pausa_off)
cat(crayon::red("\n Varios suavizamientos de una Serie de tiempo creación en workdir \n"))
serie_tiempo_plots(datos = datos,frecuencia = frecuencia,inicio = inicio,validar_ = validar_,pausa_off = pausa_off)
}
#' Realiza análisis de manera directa
#'
#' Esta función considera que las recomendaciones de la paquetería serán tomadas como los valores a usar
#'
#' Esta función comparte los mismos detalles que \code{\link{init}}
#'
#' @param datos Datos para el analisis de serie de tiempo suport (data.frame and ts class)
#' @param frecuencia Frecuencia de la serie de tiempo, sirve para
#' reescribir la frecuencia cuando datos es un objeto ts
#' @param inicio Inicio de la serie de tiempo, igual que frecuencia
#' sobreescribe valores de objetos ts
#' @param validar_ (True or False) validar los datos
#' @param ... Not work
#'
#' @return La salida no es como tal un objeto, si no una serie de impresiones de varios
#' analisis. La siguiente lista detalla alguno de ellos:
#' \itemize{\item{\bold{Plost}}{ Arroja una lista de plots que ayudan a ver el comportamiento de la serie y como ciertos ajustes se aproximan mejor a ella}}
#' \itemize{\item{\bold{Resumenes}}{ Arroja ciertos resumenes de ciertos ajustes o pruebas que se hacen}}
#' \itemize{\item{\bold{Modelo}}{ Modelo con el menor MSE(Error cuadratico medio)}}
#'
#' @export
#' @encoding UTF-8
#'
#' @examples
#' base=data.frame(tiempo=seq(Sys.Date(),by="days",length=20),valores=(rexp(50)+1)*sin(1:50))
#' Ajuste_rapido(datos=base,frecuencia=4,inicio=2010)
#'
Ajuste_rapido <- function(datos,frecuencia=NULL,inicio=NULL,
validar_=FALSE,msg=FALSE,pausa_off=1,...){
#Función para hacer el proceso de forma directa con las sugerencias como respuestas
#Probablemente se deba agregar en el futuro la opción de reporte
init(datos,frecuencia,inicio,validar_,msg,pausa_off)
}
#' Ajuste ARIMA con pruebas estadísticas
#'
#' @param datos Data.frame o objeto TS a analizar
#' @param frecuencia Frecuencia de los datos, en caso de TS sobrescribe los valores
#' @param inicio Inicio de la serie de tiempo, igual que frecuencia
#' sobrescribe valores de objetos ts
#' @param validar_ Indicador para verificar los datos [True/False]
#' @param ...
#'
#' @return La salida no es como tal un objeto, si no una serie de impresiones de varios
#' análisis. La siguiente lista detalla alguno de ellos:
#' \itemize{\item{\bold{Plost}}{ Arroja una lista de plots que ayudan a ver el comportamiento de la serie y como ciertos ajustes se aproximan mejor a ella}}
#' \itemize{\item{\bold{Resumenes}}{ Arroja ciertos resumenes de ciertos ajustes o pruebas que se hacen}}
#' \itemize{\item{\bold{Modelo}}{ Modelo con el menor MSE(Error cuadratico medio)}}
#'
#' @export
#' @encoding UTF-8
#'
#' @examples
#' base=data.frame(tiempo=seq(Sys.Date(),by="days",length=20),valores=(rexp(50)+1)*sin(1:50))
#' Ajuste_ARIMA_rapido(datos=base,frecuencia=4,inicio=2010)
#'
#' Ajuste_ARIMA_rapido(AirPassengers)
#'
Ajuste_ARIMA_rapido <- function(datos,frecuencia=NULL,inicio=NULL,
validar_=FALSE,msg=FALSE,pausa_off=1,...){
## validador integrado
paquetes.tsrutina()
pausa()
conditional.tsrutina(datos)
if(is.ts(datos)){
elementos <- tratamiento.ts_set(datos)
datos <- elementos$data
frecuencia <- ifelse(is.null(frecuencia),elementos$frecu,frecuencia)
inicio <- ifelse(is.null(inicio),elementos$inicio,inicio)
}
#solo se hace ajuste arima y prueba de estacionariedad
serie_tiempo_pruebas(datos,frecuencia,validar_,msg,pausa_off)
serie_tiempo_ARIMA(datos,frecuencia,inicio,validar_,msg,pausa_off)
}
#' De lista a reporte PDF
#'
#' Haz un reporte de tus objetos lista para poder analizar en formato de lectura tu serie de tiempo
#'
#' En desarrollo probablemente en uso hasta la versión 3.1< de la tabby verse
#'
#' @param datos Time Series a analizar
#'
#' @return PDF report
#' @export
#' @encoding UTF-8
#'
#' @import rmarkdown
#' @import knitr
#'
#' @examples
#'
#' a=TSRutina::Ajuste_ARIMA_rapido(AirPassengers,pausa_off = 0)
#' hacer_reporte(a)
#'
hacer_reporte <- function(datos){
#actual directorio del archivo
actual <- getwd()
file_direccion <- dirname(rstudioapi::getSourceEditorContext()$path)
setwd(file_direccion)
readr::write_csv(datos,"datos_temporales.csv")
#jejeje no tengo perra idea de como hacer esto :c
# rmarkdown::render("R/reporte.Rmd", "pdf_document")
# rmarkdown::render("R/reporte.Rmd", "pdf_document")
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.