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][1] 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")


  [1]: https://stats.stackexchange.com/questions/229624/fitting-truncated-distributions-using-fitdistrplus-with-a-lower-bound-of-zero