# 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")

# 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")

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.

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.