Using copulas to fit hourly observations to daily data

Question

I've been trying to automatically find low periods in some data that I have. The data is structured as hourly observations across a period of two years. Thus far I've experimented with a number of approaches (in particular, Twitter's AnomalyDetection package), without much success. My current approach is to first fit a distribution to each hour, and then to use a copula to combine those distributions into a joint, daily, distribution. However, no matter what distributions I pick, I wind up with the same problem - the daily data all appear to be very unlikely to occur (see the final graph in the code below).

I have a number of questions (how to select the right copula, whether this approach is even valid, whether another approach might be better suited to this - in particular, Renato suggested here that I try a Generalized Additive Model), but mostly would like to know what (if anything) I'm doing wrong in the code below.

Thank you for your help! Any thoughts or suggestions are much appreciated.

# first load the data 
library(repmis)
dropFile <- 
    "https://dl.dropboxusercontent.com/u/428109976/testData.csv"
testData <- repmis::source_data(dropFile, sep = ",", header = TRUE, 
    stringsAsFactors = FALSE)
testData$V1 <- NULL
library(lubridate)
testData$time <- lubridate::mdy_hm(testData$time)

# now let's try fitting a normal to demonstrate
library(fitdistrplus)
    
for (i in unique(hour(testData$time))){
  thishour <- which(hour(testData$time)==i)
  fw <- fitdist(testData[thishour, ]$obs, "norm", start = c(1,1))

  par(mfrow = c(2, 2), oma=c(0,0,2,0))
  plot.legend <- "Normal"
  denscomp(list(fw), legendtext = plot.legend)
  qqcomp(list(fw), legendtext = plot.legend)
  cdfcomp(list(fw), legendtext = plot.legend)
  ppcomp(list(fw), legendtext = plot.legend)
  title(paste("For Hour", i), outer = TRUE)
}

# clearly a normal isn't the way to go - better would probably be a 
# weibull or GPD, but this serves as a good example
# maybe a truncated distribution
# from here we'll contruct a joint distribution using copulas

library(copula)
    
# now capture the distributions fitted to each hour
groupDist <- function(x, time.colname="time", data.colname="obs", 
    distribution, timeFUN = hour, startlist = NULL) {
  timeFUN <- match.fun(timeFUN)
  distList <- list()
  for (i in unique(timeFUN(x[,time.colname]))){
    timeiter <- which(timeFUN(x[,time.colname])==i)
    fg <- fitdist(x[timeiter, data.colname], distribution, 
                  start = startlist)
    distList[[i + 1]] <- fg
  }
  distList
}
    
fitList <- groupDist(testData, distribution = "norm", 
                     startlist = c(1,1))
    # 
    param_list <- list()
    for (i in seq_along(fitList)){
      param_list[[i]] <- 
        list(mean = as.numeric(fitList[[i]]$estimate[1]),
                          sd = as.numeric(fitList[[i]]$estimate[2]))
    }
    
# now let's reorder the data so that it fits well into copula functions
testData$hour <- hour(testData$time)
testData$date <- date(testData$time)
testData$time <- NULL

library(tidyr)

wide_data <- spread(testData, hour, obs)
# for now we'll remove all rows with missing data
wide_data <- wide_data[complete.cases(wide_data[, 
                            names(wide_data) != "time"]),]
# this data doesn't fit in [0,1] as needed for `copula`, so we use pobs()
pData <- pobs(wide_data[, names(wide_data) != "time"])

# now we want to fit a copula to the data
# I'm not sure how to select which copula to use, 
# but I suspect I need an Archimedean
# as dim > 2
# let's arbitrarily choose a Frank copula for now
library(copula)
frank_model <- archmCopula("frank", dim = ncol(pData))
fit.f <- fitCopula(frank_model, pData, method = 'ml')
coef(fit.f)
library(copula)
frank_model <- archmCopula("frank", dim = ncol(pData))
fit.f <- fitCopula(frank_model, pData, method = 'ml')
coef(fit.f)

# creating our joint distribution given the marginal distributions 
# and our copula
myMvd.frank <- mvdc(copula = frankCopula(
    param = as.numeric(coef(fit.f)), dim = ncol(pData)), 
    margins = rep("norm", ncol(pData)), paramMargins = param_list)
# finally, we have the CDF of the joint distribution 

cdf_frank <- pMvdc(pData, myMvd.frank)

# let's plot to see what things look like  

library(ggplot2)  

# why are so many of the days so unlikely? this seems wrong
# I've seen this no matter my selectin of distribution and 
# (archimedean) copula
# so far I've tried fitting truncated normal, weibull, truncated 
# lognormal, truncated logis,
# gamma and exponential distributions but nothing seems to work 

ggplot(as.data.frame(cdf_frank), aes(cdf_frank, ..density..)) +
  stat_bin(bins = 100) +
  ggtitle("CDF of Frank Copula over all daily observations")