Paqueterias.

library(forecast)
library(ggplot2)
library(zoo)
library(gridExtra)
library(imputeTS)
library(zoo)
library(TSstudio)
library(gt)
library(tidyverse)
library(glue)
library(dplyr)
library(lubridate)
library(highcharter)
library(tseries) ; library(TSA)
library(timeSeries)
library(rugarch)
library(fGarch)
library(MASS)
library(astsa)
library(nortest)

Preprocesamiento de la serie

Base de datos Energy Market Price recuperada de: https://www.kaggle.com/pranavkasela/energymarket-price-time-series/metadata

El conjunto de datos comprende 3287 observaciones de precios diarios del mercado de energía que comprenden desde el 01-Enero-2010 al 31-Diciembre-2018. La base está compuesta por dos columnas:

  1. Fecha [yyyy-mm-dd]
  2. Precio [\(\$\)] - (Value)

Los datos corresponden a los precios diarios del mercado de energía, en el periodo comprendido de 01 de enero del 2010 al 31 de diciembre del 2018, los cuales se representan en la siguiente gráfica

dataset <- read.csv("/cloud/project/time_series_dataset.csv", sep=";")

value = dataset[,2]
inds <- seq(as.Date("2010-01-01"), as.Date("2018-12-31 "), by = "day")
serie = zoo(value,inds)
##
   ## Alterntiva
tsData = ts(value, start = c(2010,1),frequency = 56)
##
hchart(as.ts(value))

De la serie original, se puede suponer que se tiene una varianza constante, sin embargo, lo ideal sería buscar una transformación que permita garantizar la homocedasticidad de los errores. Además,se observa que la tendencia no puede ser considerada como creciente, constante o decreciente ya que se visualizan periodos donde la tendencia oscila entre los precios. Por ende,puede tratarse de una serie estacionaria y con posibles ciclos estacionales.

Analisis Descriptivo

## Toma de parametros una serie. 

descrip = function(s){
  ini = start(s)
  end = end(s)
  f = frequency(s)
  nas = sum(is.na(s))
  Data_f = data.frame(ini,end,f,nas) 
  colnames(Data_f) = c("Inicio","Fin","Frecuencia","Datos Faltantes")
  # Imprimimos los datos
  print(Data_f)
}

descrip(tsData)
##   Inicio  Fin Frecuencia Datos Faltantes
## 1   2010 2068         56               0
## 2      1   39         56               0

En este caso al tratar de ajustar una curva para la componente estacional se puede notar que ninguna es capaz de ajustarse adeacuadamente a la serie por lo que es posible necesitar hacer modificaciones a la serie para obtener un mejor tratamiento de la serie de tiempo.

## Tendencia

f_tendencia = function(s){

library(tidyverse)
library(ggplot2) ; library(xts);library(zoo)
  
data = data.frame(1:length(coredata(s)),coredata(s))
colnames(data) = c("x","y")


exp.model <- lm(log(data$y)~data$x,data = data) 
a = min(length(data$x),length(fitted(exp.model)))
exp.model.df <- data.frame(x=data$x[1:a],
                           y=exp(fitted(exp.model)[1:a]))

ggplot(data, aes(x = x, y = y)) + geom_line() +
  stat_smooth(method = 'lm', aes(colour = 'lineal'), se = FALSE) +
  stat_smooth(method = 'lm', formula = y ~ poly(x,2), aes(colour = 'cuadratico'), se= FALSE) +
  stat_smooth(method = 'lm', formula = y ~ poly(x,3), aes(colour = 'cubico'), se = FALSE)+
  stat_smooth(data=exp.model.df, method = 'loess',aes(x,y,colour = 'exponencial'), se = FALSE) +
  ggtitle("Ajuste de Curvas",subtitle = "Eleccion de curva por coeficiente de determinacion" ) + theme_bw()
}

f_tendencia(tsData)

Una tendencia polinomial no se ajusta adecuadamente a la serie de tiempo debido a que el coeficiente de correlacion no es suficientemente grande.

library(purrr)
library(tidyverse)
library(ggplot2) ; library(xts);library(zoo)

f_tendencia_r = function(s){
  
  data = data.frame(1:length(coredata(s)),coredata(s))
  colnames(data) = c("x","y")
  if(sum(is.na(data$y))==0){
  ### Todos los que se me ocurrieron.
  model_exponential <- lm(data = data,log(data$y)~data$x)
  ###
  z=sqrt(sum(residuals(model_exponential)^2) / df.residual(model_exponential))
  model_linear <- lm(data = data,data$y~data$x)
  adj_r_squared_linear <- summary(model_linear) %>% 
  .$adj.r.squared %>% 
  round(4)
  
  y_predicted <- model_exponential$fitted.values %>% 
  map_dbl(~exp(.+z^2/2))
  #########
  r_squared_exponential <- (cor(y_predicted,data$y))^2

  adj_r_squared_exponential <- 1-(1-r_squared_exponential)*
  ((nrow(data)-1)/(nrow(data)-1-1))
  ########
  model_quadratic <- lm(data = data,data$y~poly(data$x,2))
  model_cubic <- lm(data = data,data$y~poly(data$x,3))
  adj_r_squared_quadratic <- summary(model_quadratic) %>% 
  .$adj.r.squared

  adj_r_squared_cubic <- summary(model_cubic) %>% 
  .$adj.r.squared
  adj_r_squared_quadratic <- summary(model_quadratic) %>% 
  .$adj.r.squared

  adj_r_squared_cubic <- summary(model_cubic) %>% 
  .$adj.r.squared

  adj_r_squared_all <- c(linear=round(adj_r_squared_linear,4),
                       exponential=round(adj_r_squared_exponential,4),
                       quadratic=round(adj_r_squared_quadratic,4),
                       cubic=round(adj_r_squared_cubic,4))
  print(adj_r_squared_all)
  }else{
    print("En carguese de los NA´s primero")
  }

}
f_tendencia_r(tsData)
##      linear exponential   quadratic       cubic 
##      0.1728      0.1840      0.2585      0.4068
## Ordenar serie por temporalidad
  # s: serie
  # t: Temporalidad
df <- data.frame(date = inds,
                 value = value)

Agrupacion por ciclos

Como lo observamos en la serie, algunas parte de la serie tiene rangos mas altos que otros aunque la variabilidad nos muestra que el precio se mantienen constante.

El boxplot anterior nos muestra que los precios no siguen un comportamiento al alza o a la baja, en realidad siempre se mantienen en constante cambio , en términos financieros, estamos frente a un mercado con comportamiento latente.

Descomposicion

Al continuar con el análisis, podemos observar en el decompose que la tendencia no sigue un comportamiento lineal en realidad es muy cambiante, también podemos observar ciclos estacionales aproximadamente cada inicio de año, además de que la parte aleatoria no sigue el comportamiento de un ruido blanco.

ts_decompose(tsData,type = "both")

La varianza y la media se muestran constantes con respecto al tiempo por lo que hay indicios de estacionariedad. Aunque al ser precios podriamos considerar la serie de los rendimientos.

### ACF sin diferenciar

Aqui vemos que la serie necesita ser diferenciada pues tiene valores cercanos a 1 a pesar de tener varios retrasos en el acf posiblemente estamos hablando de una serie no estacionaria, esto era evidente ya que en la descomposicion vimos una patron estacional.

En el diagrama del ACF se puede observar que los datos inicialmente tenían una persistencia de valores altos, lo cual representaba una tendencia positiva a largo plazo, como se ha podido observar en la gráfica, sin embargo, al final tenemos el caso opuesto, el inicio de una tendencia bajista a corto plazo, que posteriormente toma un repunte a la alza. Por el lado del PACF se observa muy poca correlación entre los datos

Diferenciando la serie y viendo el ACF de los rendimientos al cuadrado vemos que la parte estacional parece que no tiene efectos sobre la serie de tiempo, quiza se deba a que el numero de observaciones perteneciente a los ciclos era de 7.

Espectro de la serie

Se requirió hacer un ajuste en la serie original para agrupar los datos por promedio de cada 7 observaciones, es decir, se realizó un promedio semanal y la nueva serie nos arroja 469 datos.

Se construye una función que devuelve el período de la frecuencia dominante de una serie de tiempo basándose en su análisis espectral, el análisis espectral es una técnica que nos permite descubrir periodicidades subyacentes ya que en este caso era muy difícil interpretar el acf y pacf.

De esta manera es como se puede afirmar que la serie presentaba dos tipos de estacionalidad, aquella que se veía reflejada al inicio de cada año y la que se veía reflejada al analizar los acf y pacf, ya que se mostraba la dependencia cada 7 observaciones atrás.

fantasma <- function(x)
{
  n <- length(x)
  spec <- spec.ar(c(x),plot=FALSE)
  if(max(spec$spec)>10)
  {
    period <- round(1/spec$freq[which.max(spec$spec)])
    if(period==Inf) 
    {
      j <- which(diff(spec$spec)>0)
      if(length(j)>0)
      {
        nextmax <- j[1] + which.max(spec$spec[j[1]:500])
        period <- round(1/spec$freq[nextmax])
      }
      else
        period <- 1
    }
  }
  else
    period <- 1
  return(period)
}

fantasma(value)
## [1] 7

Agrupacion por promedios.

### Agrupamos por dias de la Semana con domingo como dia de inicio. 
df <- data.frame(date = inds,
                 value = value)
## Quite 2 observaciones que pertenecian a viernes y sabado para iniciar con 
## Semanas enteras. 
agrupacion = df %>%
mutate(month = format(date, "%w"), year = format(date, "%Y")) %>%
group_by(month, year)
### 
yupi = length(value)

## Ya empieza en Domingo
agrupacion = agrupacion[3:(yupi-2),]

Una vez ajustando las observaciones para promediarlas con sus respectivos dias de la semana se procede a hacer promedios semanales.

value_we = c()
Indi_we = c()

### Aqui empezamos a calcular promedios

for (i in 1:(length(agrupacion$value)/7)){
  k = 7*i
  if(i ==1 ){
  value_we = append(value_we,mean(agrupacion$value[i:k]))
  }else{
    value_we = append(value_we,mean(agrupacion$value[k-6:k]))
  }
}

for(i in 1:(length(agrupacion$value)/7)){
  k = i*7
  Indi_we = append(Indi_we,agrupacion$date[k])
}
### Serie de promedios semanales
Ser = zoo(value_we,Indi_we)
Ser = Ser[4:length(Ser)]
ts.plot(as.ts(value_we), main = "Promedios Semanales",
        ylab="Precio")
grid()

Ahora podemos notar que las curvas ajustadas a la componente estacional ya tiene mejor \(R^{2}\)

f_tendencia(Ser)

f_tendencia_r(Ser)
##      linear exponential   quadratic       cubic 
##      0.7384      0.7226      0.8433      0.9691

Holt Winters

## Base de Datos. 

### 
value = dataset$value
### Codificamos la serie 
winters_data = ts(value, start = c(2010,1), end = c(2018,12), frequency = 120)

Ajuste

newts = Ser
lambda <- BoxCox.lambda(newts) 

cox = BoxCox(newts, lambda)
plot(cox,main= "Con Transformacion",col="deeppink4");grid()

Quitamos la tendencia con el metodo de diff

tsData_diff <- diff(cox)
tsdisplay(tsData_diff, main = "No tendencia")

tsData_diff2 <- diff(tsData_diff, lag = 7, differences = 1)

tsdisplay(tsData_diff2, main = "No tendencia, No ciclos estacionales")

adf.test(cox, k = 7) 
## 
##  Augmented Dickey-Fuller Test
## 
## data:  cox
## Dickey-Fuller = -1.9472, Lag order = 7, p-value = 0.6003
## alternative hypothesis: stationary
arima <- auto.arima(cox) 
res <- residuals(arima)  # get residuals from arima
plot(res, main = "Residuales auto.arima")

Box.test(arima$residuals)
## 
##  Box-Pierce test
## 
## data:  arima$residuals
## X-squared = 10.675, df = 1, p-value = 0.001086
qqnorm(residuals(arima), col = "purple3", lwd = 2) 
qqline(residuals(arima), col = "red", lwd = 2)

ad.test(residuals(arima))
## 
##  Anderson-Darling normality test
## 
## data:  residuals(arima)
## A = 112.13, p-value < 2.2e-16
shapiro.test(residuals(arima))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(arima)
## W = 0.14518, p-value < 2.2e-16