R/series_de_tiempo.R

Defines functions hacer_reporte Ajuste_ARIMA_rapido Ajuste_rapido init serie_tiempo_ARIMA recomendaciones_arma recomendacion_autocorrelaciones serie_tiempo_plots dev.TRS serie_tiempo_rutina pausa checar_datos tratamiento.fechas.TSR tratamiento.ts_set conditional.tsrutina serie_tiempo_pruebas separador paquetes.tsrutina

Documented in Ajuste_ARIMA_rapido Ajuste_rapido checar_datos conditional.tsrutina dev.TRS hacer_reporte init paquetes.tsrutina pausa recomendacion_autocorrelaciones recomendaciones_arma separador serie_tiempo_ARIMA serie_tiempo_plots serie_tiempo_pruebas serie_tiempo_rutina tratamiento.fechas.TSR tratamiento.ts_set

# 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")

}
crissthiandi/TSRutina documentation built on Dec. 3, 2024, 8:57 p.m.