Tarea 4
Sebastián Lillo
2023-10-15
library(data.table)
library(lattice)
library(ggplot2)
library(caret)
library(jtools)
library(scales)
library(lattice)
library(MASS)
library(corrplot)## corrplot 0.92 loadedlibrary(SimCorMultRes)
library(fitdistrplus)## Loading required package: survival##
## Attaching package: 'survival'## The following object is masked from 'package:caret':
##
## clusterlibrary(dplyr)##
## Attaching package: 'dplyr'## The following object is masked from 'package:MASS':
##
## select## The following objects are masked from 'package:data.table':
##
## between, first, last## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(nortest)
library(MASS)
library(boot)##
## Attaching package: 'boot'## The following object is masked from 'package:survival':
##
## aml## The following object is masked from 'package:lattice':
##
## melanomalibrary(goftest)##
## Attaching package: 'goftest'## The following objects are masked from 'package:nortest':
##
## ad.test, cvm.testrm(list = ls())Pregunta 1
Pregunta a)
Bondad de ajuste plots
# Your data
prices <- c(100.28, 99.89, 100.13, 99.38, 100.11, 102.88, 101.67, 100.29, 102.03, 105.08,
103.34, 104.36, 104.35, 103.58, 102.16, 100.48, 99.38, 96.46, 96.95, 97.68,
96.92, 98.30, 99.55, 98.17, 97.86, 97.80)
# Calculate the ratio
ratio <- c(prices[-1] / prices[-length(prices)], 1)
# Create a data frame with the ratios
df <- data.table(OriginalPrices = prices, ShiftedPrices = c(prices[-1], 1), Ratio = ratio)
df[,log_test:=log(df$Ratio)+1]
df[,log_ratio:=log(df$Ratio)]# Eliminando columnas usadas para el test inicial
df <- df[, c('OriginalPrices', 'Ratio', 'log_ratio', 'log_test')]descdist(df$log_ratio)
## summary statistics
## ------
## min: -0.02982247 max: 0.02945508
## median: -0.001888046
## mean: -0.0009631422
## estimated sd: 0.01407979
## estimated skewness: 0.3732178
## estimated kurtosis: 3.020156hist(df$log_ratio)
#Sin considerara la traslación
b_norm <- fitdist(df$log_ratio, "norm")
b_logis <- fitdist(df$log_ratio, "logis")
b_unif <- fitdist(df$log_ratio, "unif")
# Considerando la tralación
a_norm <- fitdist(df$log_test, "norm")
a_logis <- fitdist(df$log_test, "logis")
a_lnorm<-fitdist(df$log_test,"lnorm")
a_gamma<-fitdist(df$log_test,"gamma")
# Plot trasladado
plot.legend <- c("Normal","Logistic","Log Normal", "Gamma")
denscomp(list(a_norm,a_logis,a_lnorm,a_gamma), fitlty=1, legendtext=plot.legend)
# Plot general
plot.legend <- c("Normal","Logistic","Uniforme")
denscomp(list(b_norm,b_logis,b_unif), fitlty=1, legendtext=plot.legend)
cdfcomp(list(a_norm,a_logis,a_lnorm,a_gamma), fitlty=1, legendtext=plot.legend, verticals=T,
do.points=F)
ppcomp(list(a_norm,a_logis,a_lnorm,a_gamma),legendtext=plot.legend)
qqcomp(list(a_norm,a_logis,a_lnorm,a_gamma),legendtext=plot.legend)
cdfcomp(list(a_norm,a_logis,a_lnorm,a_gamma), fitlty=1, legendtext=plot.legend, verticals=T,
do.points=F)
Para la resolución de esta sección se realiza un cálculo manual de las ratios solicitadas, guardando la variable dentro de la tabla para poder utilizarla. En este caso se puede apreciar que la distribución de los datos es difícil de determinar de forma visual. Para el ploteo de más distribuciones se implementa una traslación de la distribución de los satos para segura su no negatividad. Esto permite incorporar una distribución log normal y una distribución gamma al análisis. Pese a que la distribución log normal es redundante para el análisis general, se presenta que tanto esta distribución como la gamma presentan una distribución similar a la normal. Adicionalmente, por grafico de Cullen and Frey podemos determinar una cercanía de la distribución a nuestras distribuciones tentativas. En este caso se realiza una comparación entre la log normal, uniforme, normal, logística y gamma. Junto a esto, se realizan los test visuales respectivos llegando a determinar que la distribución normal tiene un mejor comportamiento en las colas y en el centro de su distribución.
Pregunta b)
Bondad de ajuste test estadisticos
gofval <- gofstat(list(a_norm,a_logis,a_lnorm,a_gamma),fitnames=c("norm","logis","lnorm","gamma"))
gofval## Goodness-of-fit statistics
## norm logis lnorm gamma
## Kolmogorov-Smirnov statistic 0.10271436 0.09797345 0.10115442 0.10154523
## Cramer-von Mises statistic 0.03759276 0.03364960 0.03563911 0.03617506
## Anderson-Darling statistic 0.27633499 0.25711658 0.26387125 0.26761488
##
## Goodness-of-fit criteria
## norm logis lnorm gamma
## Akaike's Information Criterion -144.9117 -144.3040 -145.0317 -144.9932
## Bayesian Information Criterion -142.3955 -141.7878 -142.5156 -142.4770gofval_1 <- gofstat(list(b_norm,b_logis,b_unif),fitnames=c("norm","logis","unif"))
gofval_1## Goodness-of-fit statistics
## norm logis unif
## Kolmogorov-Smirnov statistic 0.10271436 0.09791911 0.1829550
## Cramer-von Mises statistic 0.03759276 0.03363013 0.1894371
## Anderson-Darling statistic 0.27633499 0.25706076 Inf
##
## Goodness-of-fit criteria
## norm logis unif
## Akaike's Information Criterion -144.9117 -144.3040 -142.9273
## Bayesian Information Criterion -142.3955 -141.7878 -140.4111Dada la implementación del test estadístico de AD-test más el resto de análisis estadísticos, podemos notar que el mejor ajuste a nivel estadístico esta dado por la distribución logística, no obstante dado que la representación gráfica de la distribución logística muestra un comportamiento muy alejado al de los datos en el caso del análisis de colas por lo cual podemos determinar que la distribución normal es la mejor distribución para implementar dado que presenta un comportamiento razonable a lo largo de todos los test implementados.
pv.chi.norm <- gofval$chisqpvalue["norm"]
pv.chi.logis <- gofval$chisqpvalue["logis"]
pv.chi.gamma <- gofval$chisqpvalue["gamma"]
pv.chi.lnorm <- gofval$chisqpvalue["lnorm"]
pv.cvm <- nortest::cvm.test(prices)$p.value #Solamente para Normal
pv.ad <- nortest::ad.test(prices)$p.value #Solamente para Normal
pv.ks <- nortest::lillie.test(prices)$p.value #Solamente para Normal
pv.pea <- nortest::pearson.test(prices)$p.value #Solamente para Normal
pvalue_norm <- data.frame(ChiC=pv.chi.norm,CvM=pv.cvm,AD=pv.ad,KS=pv.ks,PEA=pv.pea)
pvalue_gamma <- data.frame(ChiC=pv.chi.gamma,CvM=1,AD=1,KS=1,PEA=1)
pvalue_lnorm <- data.frame(ChiC=pv.chi.lnorm,CvM=1,AD=1,KS=1,PEA=1)
if (!is.na(pv.chi.logis))
{ pvalue_norm <- rbind(pvalue_norm,pvalue_gamma,pvalue_lnorm,data.frame(ChiC=pv.chi.logis,CvM=NaN,AD=NaN,KS=NaN,PEA=NaN))
} else {pvalue_norm <- rbind(pvalue_norm,pvalue_gamma,pvalue_lnorm,data.frame(ChiC=NaN,CvM=NaN,AD=NaN,KS=NaN,PEA=NaN))}
row.names(pvalue_norm) <- c("norm","logis",'gamma','lnorm')
pvalue_normEn esta sección se incorpora el valor p y el valor del test AD para comparar los valores. Pero el único valor consistente para nuestro análisis es el valor de la normal.
Confirmación con valor previo desplazamiento de media
pv.chi.norm <- gofval_1$chisqpvalue["norm"]
pv.cvm <- nortest::cvm.test(prices)$p.value #Solamente para Normal
pv.ad <- nortest::ad.test(prices)$p.value #Solamente para Normal
pv.ks <- nortest::lillie.test(prices)$p.value #Solamente para Normal
pv.pea <- nortest::pearson.test(prices)$p.value #Solamente para Normal
pvalue_norm <- data.frame(ChiC=pv.chi.norm,CvM=pv.cvm,AD=pv.ad,KS=pv.ks,PEA=pv.pea)
row.names(pvalue_norm) <- c("norm")
pvalue_norm# mean.data == mu_i
# sd.data == sigma_i
# Media y sigma trasladados
mu_a <- a_norm$estimate[1]
sigma_a <- a_norm$estimate[2]
# Media y sigma estandar
mu_b <- b_norm$estimate[1]
sigma_b <- b_norm$estimate[2]
est_a <- data.frame(mu=mu_a,sigma=sigma_a)
est_b <- data.frame(mu=mu_b,sigma=sigma_b)
est <- rbind(est_a,est_b)
row.names(est) <- c("Media trasladada","Media 0")
estPregunta c)
Union de pasos (i),(ii),(iii),(vi).
library(fitdistrplus)
library(ggplot2)
set.seed(1234)
mu_b <- mu_b # Replace with your desired mean
sigma_b <- sigma_b # Replace with your desired standard deviation
Z.norm <- rnorm(length(df$log_ratio), mu_b, sigma_b)
# Create empirical cumulative distribution functions (ECDF)
X <- df$log_ratio
eX <- ecdf(X)
eZN <- ecdf(Z.norm)
# Plot empirical cumulative distribution functions
plot(eX, ylab="Fn(x)", main="Empirical Distribution", col='purple',
verticals = TRUE, do.points = FALSE, xlim = c(min(X), max(X)))
legend(0.75*max(X), 0.1, leg='Datos', col = 'black', lty = 1)
par(new = TRUE)
plot(eZN, ylab="Fn(x)", col = 'darkgreen', verticals = TRUE,
do.points = FALSE, xlim = c(min(X), max(X)))
legend(0.75*max(X), 0.2, leg='simulado', col = 'black', lty = 1)
# Procedimiento de simulación
# Step (i)
sim <- Z.norm
# Media y sigma estandar
c_norm <- fitdist(sim, "norm")
mu_c <- c_norm$estimate[1]
sigma_c <- c_norm$estimate[2]
# Paso (iii)
M <- 1000 # Número de repeticiones
# Crear un vector para almacenar las estadísticas de prueba
estadisticas_T <- numeric(M)
# Generar M veces las muestras y calcular la estadística de prueba T
for (i in 1:M) {
muestras <- rnorm(1000, mean = mu_c, sd = sigma_c)
estadisticas_T[i] <- mean(muestras)
}
# Calcular p-value y el intervalo de confianza del 95%
p_values <- 1 - pnorm(estadisticas_T)
p_mean <- mean(p_values)
p_sd <- sd(p_values)
lower_ci <- p_mean - 1.96 * (p_sd / sqrt(M))
upper_ci <- p_mean + 1.96 * (p_sd / sqrt(M))
# Print p-value and confidence interval
cat("Estimación de p-value:", p_mean, "\n")## Estimación de p-value: 0.5019693cat("Intervalo de confianza al 95% para p-value: [", lower_ci, ",", upper_ci, "]\n")## Intervalo de confianza al 95% para p-value: [ 0.5019593 , 0.5019793 ]# Plot the last simulated samples and test statistics
plot(muestras, main = "Simulated Samples", ylab = "Values")
plot(estadisticas_T, main = "Test Statistics", ylab = "Values")
En esta sección se generan los pasos en conjunto de la simulación con las condiciones solicitadas. En este caso de resguarda el test T para poder realizar la comparación con estadístico de AD. Principalmente el test nos otorga una mirada a la normalidad, que al parecer es normal significativamente considerando los valores obtenidos y la visualización de la distribución empírica.
Pregunta d)
Comparación del valor P obtenido.
cat("Simulation p-value:", p_mean, "\n")## Simulation p-value: 0.5019693cat("p-value goftest:", as.numeric(pvalue_norm[1]), "\n")## p-value goftest: 0.1854036Mediante la comparación de estos dos valores, podemos ver que claramente con el goftest se tiene una mayor probabilidad de fallar en el ajuste de una distribución a los datos. Por otro lado, dado el intervalo de confianza y la generación aleatoria de muestras podemos determinar que el mayor avance que se muestra en estas secciones viene de la mano con la comparación de los valores p antes mencionados.
Pregunta e)
set.seed(1234)
Z.norm <- rnorm(length(df$log_ratio),mu_b,sigma_b)
X <- df$log_ratio
leg <- c("Datos","Normal")
color <- c(1,4)
eX <- ecdf(X)
eZN <- ecdf(Z.norm)
plot(eX, ylab="Fn(x)", main="Distribución empirica",col=color[1], verticals = TRUE,do.points=FALSE,xlim=c(min(X),max(X)))
legend(0.75*max(X),0.1, leg[1], col=color[1], lty=1)
par(new=T)
plot(eZN, ylab="Fn(x)", main="Distribución empirica",col=color[2],verticals = TRUE,do.points=FALSE,xlim=c(min(X),max(X)))
legend(0.75*max(X),0.2, leg[2], col=color[2], lty=1)
La distribución empírica presentada muestra un comportamiento cercano al
de los datos, pese al volumen reducido de información que manejamos
podemos decir que la generación aleatoria se presenta correctamente.
# Procedimiento de simulación
# Paso (i)
sim <- Z.norm
# Media y sigma estandar
c_norm <- fitdist(sim, "norm")
mu_c <- c_norm$estimate[1]
sigma_c <- c_norm$estimate[2]
mu_c## mean
## -0.004941907sigma_c## sd
## 0.01266194En esta sección realizamos la simulación principal y se guardan los valores necesarios del promedio y la desviación estándar. En particular, se generan un total de datos igual al numero de filas en la base de datos, para rescatar de esta forma los valores estadísticos muestrales de la simulación.
# Paso (ii)
ad.test(Z.norm, "pnorm")##
## Anderson-Darling test of goodness-of-fit
## Null hypothesis: Normal distribution
## Parameters assumed to be fixed
##
## data: Z.norm
## An = 9.7672, p-value = 3.353e-05Dado el test de AD nos podemos percatar que tiene un valor mucho mayor al previamente calculado. Así como también se presenta un valor p mucho más reducido, dando a entender que la distribución no se ajusta necesariamente bien a los datos. Por lo tanto, en este caso particular vemos que la distribución normal no cumple con los test de forma satisfactoria en este caso particular.
#Paso (iii)
M <- 1000 # Número de repeticiones
# Crear un vector para almacenar las estadísticas de prueba
estadisticas_AD <- numeric(M)
# Generar M veces las muestras y calcular la estadística de prueba
for (i in 1:M) {
muestras <- rnorm(1000, mean = mu_c, sd = sigma_c)
estadisticas_AD[i] <- ad.test(muestras[i], "pnorm")$statistic[1]
}
plot(muestras)
plot(estadisticas_AD)
Esta repetición permite generar un mayor numero de datos con los que trabajar, así como también guardar las estadísticas del test AD, donde se puede apreciar que los datos se concentran en una bondad de ajuste alta dado el bajo valor del test.
# Calcular el intervalo de confianza para p para cada estadística de prueba
intervalos_confianza_p <- numeric(M)
for (i in 1:M) {
intervalos_confianza_p[i] <- 1 - pchisq(estadisticas_AD[i], df = 1)
}
# Calcular el intervalo de confianza del 95% para p
intervalo_confianza_p_95 <- quantile(intervalos_confianza_p, c(0.025, 0.975))
# Imprimir el intervalo de confianza del 95% para p
print("Intervalo de confianza del 95% para p:")## [1] "Intervalo de confianza del 95% para p:"print(intervalo_confianza_p_95)## 2.5% 97.5%
## 0.5339573 0.5342537# Calcular la probabilidad p correspondiente a cada estadística de prueba
probabilidades_p <- 1 - pchisq(estadisticas_AD, df = 1)
# Calcular la media de las probabilidades p
valor_p <- mean(probabilidades_p)
# Imprimir el valor p
print("Valor estimado de p:")## [1] "Valor estimado de p:"print(valor_p)## [1] 0.5341918Por lo tanto, utilizando los M valores rescatados para la estadística del test AD se puede concluir que la distribución de las simulaciones reiteradas tiene un valor p elevado para asegurar un buen ajuste de una distribución normal en los datos, así como también un intervalo de confianza acotado a valores cercanos. Por lo tanto, la probabilidad de que la distribución de los datos sea normal es elevada considerando la diferencia que muestra con los test anteriormente trabajados.
Pregunta 2
library(data.table)
library(lattice)
library(ggplot2)
library(caret)
library(jtools)
library(scales)
library(lattice)
library(MASS)
library(corrplot)
library(SimCorMultRes)
library(fitdistrplus)
library(dplyr)
library(nortest)
library(MASS)
library(boot)
library(goftest)
library(purrr)##
## Attaching package: 'purrr'## The following object is masked from 'package:scales':
##
## discard## The following object is masked from 'package:caret':
##
## lift## The following object is masked from 'package:data.table':
##
## transposelibrary(tidyr)rm(list = ls())#Ajustes para gráficos via ggplot
mytheme <- theme(plot.title = element_text(size = 22),
axis.title.y=element_text(size = 20),
axis.text.y=element_text(size = 15),
axis.title.x=element_text(size = 20),
axis.text.x=element_text(size = 15),
legend.position="bottom",
legend.text=element_text(size = 15),
legend.title=element_text(size = 15))df <- fread("energy_dataset.csv")
df1 <- fread("weather_features.csv")
df1 <- df1 %>%
rename(time = dt_iso)# Merge
df_time <- full_join(df, df1, by = 'time')
# Separación de tiempo
df_time$year <- year(df_time$time)
df_time$month <- month(df_time$time)
df_time$day <- mday(df_time$time)
df_time$hour <- hour(df_time$time)
df_time$hour[df_time$hour == 0] <- 24
# Selección de Barcelona
df_barcelona <- df_time %>%
filter(city_name == 'Barcelona')
df_barcelona <- df_barcelona[complete.cases(df_barcelona$`generation solar`), ]En esta sección se realiza la limpieza de los datos y se filtra la base considerando la ciudad de Barcelona. La elección de ciudad se dio de forma arbitraria ya que Barcelona es la ciudad más importante a nivel arquitectónico considerando los aportes de Gaudi.
Sección a)
(i) Pronostico por regresión.
# Selección de la data
data7<-df_barcelona[year(time)<=2018 &
month(time)==7 & mday(time)<31,
.(anio=year(time),
dia= mday(time),
hora= hour(time),
solar=`forecast solar day ahead`,
generation=`generation solar`,
eol=`forecast wind onshore day ahead`,
dem=`total load forecast`,
precio=`price actual`,
time= as.POSIXct(time))] #Guarda el campo "time" para uso posterior
test7<-df_barcelona[year(time)==2018 &
month(time)==7 & mday(time)==31,
.(anio=year(time),
dia= mday(time),
hora= hour(time),
solar=`forecast solar day ahead`,
generation=`generation solar`,
eol=`forecast wind onshore day ahead`,
dem=`total load forecast`,
precio=`price actual`,
time= as.POSIXct(time))]
generation7_test <- df_barcelona[year(time)==2018 &
month(time)==7 & mday(time)==31,.(generation=`generation solar`)]Sea realiza una división de la data considerando el proceso de test y generación de un modelo bajo los criterios vistos en clases.
mod7<-with(data7,lm(formula=generation ~ hora+eol+dem))
summary(mod7)##
## Call:
## lm(formula = generation ~ hora + eol + dem)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4337.0 -1050.6 -227.1 944.3 4866.6
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.335e+03 1.975e+02 -21.95 < 2e-16 ***
## hora -2.623e+01 4.423e+00 -5.93 3.4e-09 ***
## eol -1.333e-01 1.162e-02 -11.47 < 2e-16 ***
## dem 2.414e-01 6.548e-03 36.86 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1495 on 2890 degrees of freedom
## Multiple R-squared: 0.3692, Adjusted R-squared: 0.3686
## F-statistic: 563.9 on 3 and 2890 DF, p-value: < 2.2e-16Este modelo de regresión lineal muestra que los valores asociados a los parámetros de cada una de las variables incorporadas son significativo para la interpretación del modelo. En este caso se incorpora la información atmosférica pertinente seleccionada desde la base de datos considerando la variable objetivo que se busca predecir y las características de las demás variables, por lo que se excluyen las variables de generación eléctrica y solo se incorpora la información relacionada con el clima y su forcast.
pron_horario <- predict(mod7,newdata=test7,level=0.95,int="p") #Añade intervalo de predicción
#head(fore_hour)
#Arma el dataframe con los datos pronosticados
types <- c("generation"="real", "fit"="pronóstico","lwr"="inf","upr"="sup")
horas <- data.frame(hora=0:23)
pron.df <- cbind(horas,pron_horario)
pron <- pron.df %>% tidyr::gather("type", "generation", 2:4)
plot.pronostico <- ggplot(pron,aes(hora,generation,group=type,color=types[type]))+
geom_line()+
labs(x = 'Hora', y = 'generation',
title='Pronóstico 95% para el día 31/7/2018',
color="generation")+
mytheme
real.df <- cbind(horas,generation7_test)
real <- real.df %>% tidyr::gather("type", "generation", 2)
plot.pronostico+geom_line(data=real,aes(hora,generation,color=types[type]))
res7<-resid(mod7)
sim.df <- cbind(horas,sim=predict(mod7,newdata=test7) + rnorm(24,0,sd(res7)))
sim <- sim.df %>% tidyr::gather("type", "generation", 2)
plot.pronostico+geom_line(data=sim,aes(hora,generation,color=types[type]))+
labs(title='Simulación para el día 31/7/2018')
Este grafico muestra la distribución y los márgenes de predicción con un 95% de confianza. Los resultados se fundamentan en el modelo regresivo lineal, por lo que podemos notar que existe cierta diferencia entre la predicción y los datos reales considerando que el modelo utilizado tiene un R^2 bastante bajo en un 30% y un test estadístico F relativamente alto, por lo que es posible que se estén excluyendo algunas variables significativas para la predicción. No obstante, es claro que este modelo al ser tan básico en su predicción no es capas de representar de forma fehaciente la generación solar.
(ii) Regreción categorica
# generación de variables categoricas
binary_vars <- data.frame(matrix(0, nrow = nrow(df_barcelona), ncol = 24))
for (i in 1:24) {
binary_vars[, i] <- ifelse(df_barcelona$hour == i, 1, 0)
}
df_barcelona <- cbind(df_barcelona, binary_vars)Para realizar la regresión categórica se considero el proceso descrito en el enunciado. Se genero un matriz vacía para contener las variables y se implementa un cambio en la numeración de las horas colocando la hora 0 como la hora 24 para facilitar el cálculo y generación de las variables.
# Selección de la data
data7 <- df_barcelona[year(time) <= 2018 & month(time) == 7 & mday(time) < 31,
.(anio = year(time),
dia = mday(time),
hora = hour(time),
solar = `forecast solar day ahead`,
generation = `generation solar`,
eol = `forecast wind onshore day ahead`,
dem = `total load forecast`,
precio = `price actual`,
time = as.POSIXct(time), X1,X2,X3,X4,X5,X6,X7,X8,X9,X10,X11,X12,X13,X14,X15,X16,X17,X18,X19,X20,X21,X22,X23,X24)] #Guarda el campo "time" para uso posterior
test7<-df_barcelona[year(time)==2018 &
month(time)==7 & mday(time)==31,
.(anio = year(time),
dia = mday(time),
hora = hour(time),
solar = `forecast solar day ahead`,
generation = `generation solar`,
eol = `forecast wind onshore day ahead`,
dem = `total load forecast`,
precio = `price actual`,
time = as.POSIXct(time), X1,X2,X3,X4,X5,X6,X7,X8,X9,X10,X11,X12,X13,X14,X15,X16,X17,X18,X19,X20,X21,X22,X23,X24)]
generation7_test <- df_barcelona[year(time)==2018 &
month(time)==7 & mday(time)==31,.(generation=`generation solar`)]mod7<-with(data7,lm(formula=generation ~ hora+eol+dem+ X1+X2+X3+X4+X5+X6+X7+X8+X9+X10+X11+X12+X13+X14+X15+X16+X17+X18+X19+X20+X21+X22+X23))
summary(mod7)##
## Call:
## lm(formula = generation ~ hora + eol + dem + X1 + X2 + X3 + X4 +
## X5 + X6 + X7 + X8 + X9 + X10 + X11 + X12 + X13 + X14 + X15 +
## X16 + X17 + X18 + X19 + X20 + X21 + X22 + X23)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3553.3 -243.2 127.4 506.0 1414.7
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.375e+02 1.494e+02 2.928 0.003439 **
## hora 1.104e+00 4.538e+00 0.243 0.807876
## eol -9.086e-02 6.397e-03 -14.204 < 2e-16 ***
## dem 1.867e-02 5.009e-03 3.727 0.000198 ***
## X1 -3.080e+01 1.018e+02 -0.303 0.762212
## X2 -6.570e+01 9.987e+01 -0.658 0.510721
## X3 -1.188e+02 9.802e+01 -1.212 0.225790
## X4 -2.337e+02 9.602e+01 -2.434 0.015004 *
## X5 -3.111e+02 9.551e+01 -3.257 0.001139 **
## X6 7.741e+01 9.575e+01 0.808 0.418914
## X7 1.291e+03 9.755e+01 13.233 < 2e-16 ***
## X8 2.585e+03 9.857e+01 26.228 < 2e-16 ***
## X9 3.288e+03 9.908e+01 33.190 < 2e-16 ***
## X10 3.647e+03 9.949e+01 36.651 < 2e-16 ***
## X11 3.841e+03 9.996e+01 38.422 < 2e-16 ***
## X12 3.874e+03 9.880e+01 39.214 < 2e-16 ***
## X13 3.774e+03 9.769e+01 38.633 < 2e-16 ***
## X14 3.498e+03 9.754e+01 35.862 < 2e-16 ***
## X15 3.029e+03 9.792e+01 30.934 < 2e-16 ***
## X16 2.456e+03 9.845e+01 24.948 < 2e-16 ***
## X17 1.762e+03 9.936e+01 17.730 < 2e-16 ***
## X18 8.719e+02 9.988e+01 8.729 < 2e-16 ***
## X19 1.494e+02 1.015e+02 1.473 0.140897
## X20 -1.699e+01 1.021e+02 -0.166 0.867847
## X21 -2.197e+01 1.015e+02 -0.217 0.828577
## X22 -3.082e-02 1.023e+02 0.000 0.999760
## X23 NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 806.6 on 2868 degrees of freedom
## Multiple R-squared: 0.8179, Adjusted R-squared: 0.8163
## F-statistic: 515.3 on 25 and 2868 DF, p-value: < 2.2e-16En este modelo se puede apreciar mejor que contamos con un R^2 mucho mayor que en el modelo anterior, pero también se sufre de tener incluidas variables poco significativas o con colinealidad. Lo más interesante es notar que las horas de la tarde poseen una mayor significancia que las variables que se asocian a horas de la madrugada o de la noche. Por lo tanto, este modelo presenta un mejor ajuste para la predicción particular.
pron_horario <- predict(mod7,newdata=test7,level=0.95,int="p") #Añade intervalo de predicción## Warning in predict.lm(mod7, newdata = test7, level = 0.95, int = "p"):
## prediction from a rank-deficient fit may be misleading#head(fore_hour)
#Arma el dataframe con los datos pronosticados
types <- c("generation"="real", "fit"="pronóstico","lwr"="inf","upr"="sup")
horas <- data.frame(hora=0:23)
pron.df <- cbind(horas,pron_horario)
pron <- pron.df %>% tidyr::gather("type", "generation", 2:4)
plot.pronostico <- ggplot(pron,aes(hora,generation,group=type,color=types[type]))+
geom_line()+
labs(x = 'Hora', y = 'generation',
title='Pronóstico 95% para el día 31/7/2018',
color="generation")+
mytheme
real.df <- cbind(horas,generation7_test)
real <- real.df %>% tidyr::gather("type", "generation", 2)
plot.pronostico+geom_line(data=real,aes(hora,generation,color=types[type]))
res7<-resid(mod7)
sim.df <- cbind(horas,sim=predict(mod7,newdata=test7) + rnorm(24,0,sd(res7)))## Warning in predict.lm(mod7, newdata = test7): prediction from a rank-deficient
## fit may be misleadingsim <- sim.df %>% tidyr::gather("type", "generation", 2)
plot.pronostico+geom_line(data=sim,aes(hora,generation,color=types[type]))+
labs(title='Simulación para el día 31/7/2018')
Considerando el modelo generado y la representación grafica que tenemos se puede concluir que este modelo genera una mejor simulación de la generación solar considerando las variables de la hora de forma particular en ves de agrupadas en una sola variable. Como se presenta en los gráficos se reduce el margen de error asociado a las predicción y se realiza un mejor ajuste de la línea predictiva.
(iii) 5 simulaciones
set.seed(1234) # Set seed for reproducibility
num_simulations <- 5
simulated_generations <- matrix(NA, ncol = 24, nrow = num_simulations)
for (i in 1:num_simulations) {
simulated_generations[i,] <- predict(mod7, newdata = data.frame(hora = 0:23,
eol = sample(data7$eol, 24, replace = TRUE),
dem = sample(data7$dem, 24, replace = TRUE),
X1 = sample(data7$X1, 24, replace = TRUE),
X2 = sample(data7$X2, 24, replace = TRUE),
X3 = sample(data7$X3, 24, replace = TRUE),
X4 = sample(data7$X4, 24, replace = TRUE),
X5 = sample(data7$X5, 24, replace = TRUE),
X6 = sample(data7$X6, 24, replace = TRUE),
X7 = sample(data7$X7, 24, replace = TRUE),
X8 = sample(data7$X8, 24, replace = TRUE),
X9 = sample(data7$X9, 24, replace = TRUE),
X10 = sample(data7$X10, 24, replace = TRUE),
X11 = sample(data7$X11, 24, replace = TRUE),
X12 = sample(data7$X12, 24, replace = TRUE),
X13 = sample(data7$X13, 24, replace = TRUE),
X14 = sample(data7$X14, 24, replace = TRUE),
X15 = sample(data7$X15, 24, replace = TRUE),
X16 = sample(data7$X16, 24, replace = TRUE),
X17 = sample(data7$X17, 24, replace = TRUE),
X18 = sample(data7$X18, 24, replace = TRUE),
X19 = sample(data7$X19, 24, replace = TRUE),
X20 = sample(data7$X20, 24, replace = TRUE),
X21 = sample(data7$X21, 24, replace = TRUE),
X22 = sample(data7$X22, 24, replace = TRUE),
X23 = sample(data7$X23, 24, replace = TRUE)))
}par(mfrow=c(3, 2)) # Set up the layout for subplots
for (i in 1:num_simulations) {
plot(simulated_generations[i,], type = 'l', xlab = 'Hora', ylab = 'Generacion',
main = paste("Simulacion ", i), col = i)+ mytheme
}
# Set up the layout for subplots
par(mfrow=c(1, 1))
# Plot all simulations on one plot
plot(simulated_generations[1,], type = 'l', xlab = 'Hora', ylab = 'Generacion',
main = "Simulaciones", col = 1, ylim = c(min(simulated_generations), max(simulated_generations))+ mytheme)
for (i in 2:num_simulations) {
lines(simulated_generations[i,], col = i)
}
# Add a legend
legend("topright", legend = paste("Simulacion", 1:num_simulations), col = 1:num_simulations, lty = 1)
Considerando las 5 simulaciones generadas podemos notar que, si bien el modelo realiza un buen ajuste a la predicción de los datos, sufre de una gran variabilidad considerando que todas las simulaciones generadas poseen valores distintos y muy alejadas unas de otras. De esta forma, el modelo presenta una gran variabilidad en la forma en la que genera sus predicciones y esto puede afectar en las conclusiones extraídas del mismo.
Sección b)
(i) Modelo Auto-regresivo.
# Selección de la data
data7<-df_barcelona[year(time)<=2018 &
month(time)==7 & mday(time)<31,
.(anio=year(time),
dia= mday(time),
hora= hour(time),
solar=`forecast solar day ahead`,
generation=`generation solar`,
eol=`forecast wind onshore day ahead`,
dem=`total load forecast`,
precio=`price actual`,
time= as.POSIXct(time))] #Guarda el campo "time" para uso posterior
test7<-df_barcelona[year(time)==2018 &
month(time)==7 & mday(time)==31,
.(anio=year(time),
dia= mday(time),
hora= hour(time),
solar=`forecast solar day ahead`,
generation=`generation solar`,
eol=`forecast wind onshore day ahead`,
dem=`total load forecast`,
precio=`price actual`,
time= as.POSIXct(time))]
generation7_test <- df_barcelona[year(time)==2018 &
month(time)==7 & mday(time)==31,.(generation=`generation solar`)]Se vuelve a incorporar la data solicitada para mantener constancia con la pregunta y la comparación pertinente.
library(fitdistrplus)
library(data.table)
library(nortest)
library(ggplot2)
library(tidyr)generation7_train <- data7[anio==2018,generation]
armod7<-ar(generation7_train)
pron_ar <- predict(armod7,n.ahead=24)
cat(sprintf("Valores pronosticados\n\n"))## Valores pronosticadospron_ar$pred## Time Series:
## Start = 722
## End = 745
## Frequency = 1
## [1] 756.5094 805.8843 808.0238 757.5847 583.2555 513.1999 1140.7329
## [8] 2507.3076 3839.5688 4582.0724 4887.9199 5035.2651 5112.0617 5070.3631
## [15] 4826.2521 4411.7658 3880.3068 3092.7160 2008.1416 1064.4981 649.2585
## [22] 576.5472 546.1445 535.8009cat(sprintf("\n\n"))cat(sprintf("Errores estándar\n\n"))## Errores estándarpron_ar$se## Time Series:
## Start = 722
## End = 745
## Frequency = 1
## [1] 240.2539 446.1470 599.0785 704.5765 778.5476 829.3511 861.6895 880.0018
## [9] 889.1991 893.0485 895.0069 895.6937 895.7150 895.7652 895.8991 896.2686
## [17] 896.8708 897.3456 897.7267 898.2602 898.8160 899.4132 900.2848 900.3695Este modelo presenta una consideración de la generación solar en el año 2018, por lo que incluye disantos valores que se presentan en el código anterior.
Gráfico Comparado
library(forecast)## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zootypes <- c("pron.Point.Forecast"="pronóstico","pron.Lo.95"="inf","pron.Hi.95"="sup", "generation"="real", "sim"="sim")
horas <- data.frame(hora=0:23)
pron_ar_conf <- forecast(armod7,24,level=95) #Level es el nivel de confianza (95%)
pron.ar.df <- data.frame(hora=horas,pron=pron_ar_conf)
#head(pron.ar.df)
pron.ar2 <- pron.ar.df %>% tidyr::gather("type", "generation", 2:4)
#head(pron.ar2)
plot.pronostico <- ggplot(pron.ar2,aes(hora,generation,group=type,color=types[type]))+
geom_line()+
labs(x = 'Hora', y = 'Generation',
title='Pronóstico 95% para el día 31/7/2018',
color="Precios:")+
mytheme
real.df <- cbind(horas,generation7_test)
#head(real.df)
real <- real.df %>% tidyr::gather("type", "generation", 2)
#head(real)
plot.pronostico <- plot.pronostico+geom_line(data=real,aes(hora,generation,color=types[type]))
plot.pronostico
El grafico comparado presentado por el modelo auto regresivo presenta datos muy similares a la data real, por lo que hasta el momento parece ser el mejor modelo para poder predecir la generación solar de forma fidedigna.
m.ar <- armod7$x.mean
s.ar <- sqrt(armod7$var.pred)
p.ar <- armod7$order
cat("mean=", m.ar, " sd=",s.ar, " order=",p.ar)## mean= 1996.083 sd= 240.2539 order= 27cat("Coeficientes AR:")## Coeficientes AR:coef <- as.numeric(armod7$ar)
coef## [1] 1.564729491 -0.784264800 0.166867899 0.007332011 -0.046152510
## [6] -0.004426049 -0.004346847 0.007969413 -0.008969440 0.060564680
## [11] -0.099061411 0.029879977 0.050703381 -0.043929044 -0.021809380
## [16] 0.045086212 0.004018969 -0.042666776 0.000799198 0.037447872
## [21] -0.040126162 -0.001791610 0.258179204 0.090237300 -0.290114000
## [26] -0.083636136 0.114646500sim.ar <- vector(mode="numeric",length=24)
ndata <- length(generation7_train)
prev <- generation7_train[ndata:(ndata-p.ar+1)]-m.ar
for (i in 1:24){
sim.ar[i] <- m.ar + coef%*%prev + rnorm(1,0,s.ar)
prev <- c(sim.ar[i]-m.ar,prev[-p.ar])
}sim.df <- cbind(horas,sim=sim.ar)
#head(sim.df)
sim.plot <- sim.df %>% tidyr::gather("type", "generation", 2)
#head(sim)
plot.pronostico+geom_line(data=sim.plot,aes(hora,generation,color=types[type]))+
labs(title='Simulación via AR para el día 31/7/2018')
Adicionalmente se considera una simulación para poder comparar los resultados obtenidos. De esta forma, es rango generado es el mejor de entre los tres modelos generados, así como también en la generación de una simulación para los datos en este rango particular. Si se considera una simulación adicional se puede apreciar que existe una variabilidad mayor pero aun cercana al modelo ajustado.
(ii) 5 Simulaciones.
simulated_generations <- matrix(NA, ncol = 24, nrow = 5)
for (i in 1:5) {
simulated_generations[i,] <- predict(armod7, n.ahead = 24)$pred
}
# Print the simulated_generations matrix
print(simulated_generations)## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,] 756.5094 805.8843 808.0238 757.5847 583.2555 513.1999 1140.733 2507.308
## [2,] 756.5094 805.8843 808.0238 757.5847 583.2555 513.1999 1140.733 2507.308
## [3,] 756.5094 805.8843 808.0238 757.5847 583.2555 513.1999 1140.733 2507.308
## [4,] 756.5094 805.8843 808.0238 757.5847 583.2555 513.1999 1140.733 2507.308
## [5,] 756.5094 805.8843 808.0238 757.5847 583.2555 513.1999 1140.733 2507.308
## [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16]
## [1,] 3839.569 4582.072 4887.92 5035.265 5112.062 5070.363 4826.252 4411.766
## [2,] 3839.569 4582.072 4887.92 5035.265 5112.062 5070.363 4826.252 4411.766
## [3,] 3839.569 4582.072 4887.92 5035.265 5112.062 5070.363 4826.252 4411.766
## [4,] 3839.569 4582.072 4887.92 5035.265 5112.062 5070.363 4826.252 4411.766
## [5,] 3839.569 4582.072 4887.92 5035.265 5112.062 5070.363 4826.252 4411.766
## [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]
## [1,] 3880.307 3092.716 2008.142 1064.498 649.2585 576.5472 546.1445 535.8009
## [2,] 3880.307 3092.716 2008.142 1064.498 649.2585 576.5472 546.1445 535.8009
## [3,] 3880.307 3092.716 2008.142 1064.498 649.2585 576.5472 546.1445 535.8009
## [4,] 3880.307 3092.716 2008.142 1064.498 649.2585 576.5472 546.1445 535.8009
## [5,] 3880.307 3092.716 2008.142 1064.498 649.2585 576.5472 546.1445 535.8009par(mfrow=c(3, 2)) # Set up the layout for subplots
for (i in 1:5) {
plot(simulated_generations[i,], xlab = 'Hora', ylab = 'Generacion',
main = paste("Simulacion ", i), col = i)+ mytheme
}
Como se ha presentado a lo largo de este trabajo todos los modelos
tienen sus características propias y aportan a la realización de un
modelo en conjunto. Considerando el caso del modelo AR se puede apreciar
un caso muy particular en el que las distintas simulaciones presentan
una distribución similar. Este caso es producido dada la característica
de la data y del modelo en sí, el cual considera los datos del día
anterior para ir avanzando en la generación del modelo. Por esta razón,
la construcción de los modelos y su análisis son completamente distintos
y sus resultados se pueden ver claramente diferenciados por la variación
en simulaciones múltiples y en la bondad de ajuste que presentan en su
modelo particular. Con respecto a los pronósticos obtenidos en a(iii),
las 5 simulaciones obtenidas en a son claramente diferenciables de las
simulaciones en b, dado que se presenta una mayor variabilidad en a y
diferencias entre las simulaciones. Por el otro lado en b se presenta
que las 5 simulaciones son iguales por lo que no incluye ninguna
variabilidad lo cual puede afectar al modelo dada su rigidez. En
conclusión, es necesario genera run modelo medio que considera cierta
variabilidad para poder generar un buen pronóstico y no sobre ajustar
respecto a los valores que presentamos en los modelos.