Harmonic stuff

R/clusterviz.R

@@ -193,3 +193,85 @@ ggplot(nft, aes(x, res)) + geom_line() + scale_x_datetime(date_breaks = "1 week"
 
 ggplot(nft, aes(y = res)) + geom_boxplot()
 ggplot(nft, aes(x = res)) + geom_histogram()
+
+
+c9ts <- filter(aggdf, cluster == "9")$kwh_tot_mean
+c9ft <- as.numeric(abs(fft(c9ts)))
+c9ftdf <- data.frame(x = 1:length(c9ft), ft = c9ft)
+# %>% arrange(desc(ft)) %>% mutate(xt = 1/x)
+
+c9df <- data.frame(x = 1:length(c9ts), y = c9ts, ly = log(c9ts + min(c9ts) + 1))
+c9harm. <- harmonic(ts(c9ts, frequency = 17544), 2)
+sinmod <- lm(c9df$y ~ c9harm.)
+c9df$fit <- fitted(sinmod)
+c9df$res <- resid(sinmod)
+
+ggplot(c9df, aes(x = x, y = ly)) + geom_point() +
+    geom_line(aes(y = fit), size = 2, colour = "blue")
+
+plot(c(c9df$res))
+hist(log(c9df$y), breaks = 30)
+
+n <- 1.5
+ytra <- c9df$y / c9df$fit^n - 1/c9df$fit^(n-1)
+plot(ytra)
+
+
+yts <- msts(ytra, seasonal.periods = c(48, 48*7))
+autoplot(mstl(yts))
+
+ytsd <- ts(ytra, frequency = 48)
+ydharm. <- harmonic(ytsd, 10)
+sdmod <- lm(ytra ~ ydharm.)
+sddf <- data.frame(x = 1:length(ytra), y = ytra, f = fitted(sdmod), r = resid(sdmod))
+ggplot(sddf[1:(48*7*2), ], aes(x = x, y = y)) + geom_point() +
+    geom_line(aes(y = f), colour = "blue", size = 2)
+
+plot(resid(sdmod)[1:(48*7*2)])
+
+
+
+
+
+
+c9tsy <- ts(c9ts, frequency = floor(48 * 365.25))
+c9tsw <- ts(c9ts, frequency = floor(48 * 7))
+c9tsd <- ts(c9ts, frequency = floor(48))
+
+harm.y <- harmonic(c9tsy, 2)
+harm.w <- harmonic(c9tsw, 2)
+harm.d <- harmonic(c9tsd, 6)
+
+cmod <- lm(c9ts ~ harm.y + harm.w + harm.d)
+summary(cmod)
+plot(fitted(cmod))
+plot(resid(cmod))
+
+
+cmodi <- lm(c9ts ~ harm.y * harm.w * harm.d)
+
+summary(cmodi)
+
+plot(fitted(cmodi))
+plot(c9ts)
+plot(resid(cmodi))
+hist(resid(cmodi), breaks = 100)
+# Heavy tails
+qqnorm(resid(cmodi))
+
+
+clust = "6"
+
+cts <- filter(aggdf, cluster == clust)$kwh_tot_mean
+ctsy <- ts(cts, frequency = floor(48 * 365.25))
+ctsw <- ts(cts, frequency = floor(48 * 7))
+ctsd <- ts(cts, frequency = floor(48))
+harm.y <- harmonic(ctsy, 5)
+harm.w <- harmonic(ctsw, 3)
+harm.d <- harmonic(ctsd, 3)
+cmodi <- lm(cts ~ harm.y * harm.w * harm.d)
+summary(cmodi)
+
+plot(fitted(cmodi))
+plot(cts)
+plot(resid(cmodi))