Working on seasonality; very much a WIP

R/clusterviz.R

@@ -1,6 +1,9 @@
 library(reticulate)
 library(ggplot2)
 library(dplyr)
+library(tidyr)
+library(TSA)
+library(forecast)
 theme_set(theme_bw())
 use_virtualenv("../venv/")
 
@@ -11,6 +14,7 @@ aggdf <- p$read_pickle("../data/9-clusters.agg.pkl")
 aggdf <- as.data.frame(aggdf)
 aggdf$cluster <- factor(aggdf$cluster)
 str(aggdf)
+clusters = levels(aggdf$cluster)
 
 ggplot(aggdf, aes(y = kwh_tot_mean, x = cluster)) + geom_boxplot()
 
@@ -91,3 +95,54 @@ ggsave("jul-9-fre.png", julfre, path = "../img/", dpi = "retina", width = 40, he
 ggsave("oct-9-fix.png", octcon, path = "../img/", dpi = "retina", width = 40, height = 25, units = "cm")
 ggsave("oct-9-fre.png", octfre, path = "../img/", dpi = "retina", width = 40, height = 25, units = "cm")
 
+
+# ----
+
+
+cacf <- list()
+perd <- list()
+
+for (c in clusters) {
+    cagg <- filter(aggdf, cluster == c)
+    cacf[[c]] <- acf(cagg$kwh_tot_mean, lag.max = 48*365, plot = FALSE)$acf
+    per <- periodogram(cagg$kwh_tot_mean, plot = FALSE)
+    perd[[c]] <- data.frame(freq = per$freq, spec = per$spec) %>% mutate(cluster = c, period = (1/freq)) 
+    #%>% arrange(desc(spec)) %>% head(5)
+}
+
+acfm <- sapply(cacf, as.numeric) %>% as.data.frame() %>% mutate(hour = ((1:length(`1`)) - 1)/2) %>% 
+    gather(key = "cluster", value = "acorr", clusters) %>% mutate(day = hour / 24, week = hour / (24 * 7))
+
+
+fcorr <- ggplot(acfm, aes(x = week, y = acorr, color = cluster)) + geom_line(size = 1.5) +
+    scale_color_manual(values = cbp) + facet_grid(cluster ~ .) + coord_cartesian(expand = FALSE) +
+    theme(legend.position = "none") + labs(title = "Autocorrelation plot (full year)", 
+                                           y = "Autocorrelation", x = "lag (weeks)")
+
+wcorr <- ggplot(acfm, aes(x = day, y = acorr, color = cluster)) + geom_line(size = 1.5) +
+    scale_color_manual(values = cbp) + facet_grid(cluster ~ .) +
+    scale_x_continuous(breaks = unique(floor(acfm$day / 7)) * 7) +
+    theme(legend.position = "none") + coord_cartesian(xlim = c(0, 15), expand = FALSE) +
+    labs(title = "Autocorrelation plot (two weeks)", y = "Autocorrelation", x = "lag (days)")
+
+ggsave("full-autocorr.png", fcorr, path = "../img/", dpi = "retina", width = 40, height = 25, units = "cm")
+ggsave("week-autocorr.png", wcorr, path = "../img/", dpi = "retina", width = 40, height = 25, units = "cm")
+
+perd <- bind_rows(perd)
+
+ggplot(perd, aes(x = freq, y = spec)) + 
+    geom_line() + facet_grid(cluster ~ ., scales = "free") +
+    scale_x_continuous(breaks = 1 / (c(48, 48*7, 48*7*4, 48*365)))
+
+c1ts <- filter(aggdf, cluster == "1")$kwh_tot_mean
+cts <- ts(c1ts, frequency = 48, start = c(1, 1))
+# carima <- auto.arima(cts, trace = TRUE)
+# plot(forecast(carima, h = 480))
+ctsnp <- msts(c1ts, c(48, 48*7))
+ctbats <- tbats(ctsnp)
+plot(forecast(ctbats, h = 48 * 7 * 4))
+
+c9ts <- filter(aggdf, cluster == "9")$kwh_tot_mean
+ctsnp <- msts(c9ts, c(48, 48*7))
+ctbats <- tbats(ctsnp)
+plot(forecast(ctbats, h = 48 * 7 * 4))