[spcopula-commits] r118 - in pkg: . R man

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Thu Dec 5 22:39:46 CET 2013


Author: ben_graeler
Date: 2013-12-05 22:39:45 +0100 (Thu, 05 Dec 2013)
New Revision: 118

Added:
   pkg/R/tailDependenceFunctions.R
   pkg/R/tawn3pCopula.R
   pkg/man/tailDepFun.Rd
   pkg/man/tawn3pCopula-class.Rd
   pkg/man/tawn3pCopula.Rd
Modified:
   pkg/DESCRIPTION
   pkg/NAMESPACE
   pkg/R/linkingVineCopula.R
   pkg/R/spVineCopula.R
   pkg/R/vineCopulas.R
Log:
- added 3 parameter Tawn Copula
- added parametric and empirical tail dependence functions

Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION	2013-11-27 15:21:31 UTC (rev 117)
+++ pkg/DESCRIPTION	2013-12-05 21:39:45 UTC (rev 118)
@@ -39,4 +39,6 @@
   returnPeriods.R
   spatialGaussianCopula.R
   tawnCopula.R
+  tawn3pCopula.R
+  tailDependenceFunctions.R
   zzz.R

Modified: pkg/NAMESPACE
===================================================================
--- pkg/NAMESPACE	2013-11-27 15:21:31 UTC (rev 117)
+++ pkg/NAMESPACE	2013-12-05 21:39:45 UTC (rev 118)
@@ -13,6 +13,7 @@
 export(surGumbelCopula, r90GumbelCopula, r270GumbelCopula)
 export(tawnT1Copula, surTawnT1Copula, r90TawnT1Copula, r270TawnT1Copula)
 export(tawnT2Copula, surTawnT2Copula, r90TawnT2Copula, r270TawnT2Copula)
+export(tawn3pCopula)
 export(spCopula, stCopula)
 export(vineCopula, spVineCopula, stVineCopula)
 export(neighbourhood, stNeighbourhood)
@@ -38,6 +39,8 @@
 
 # fitting
 export(fitCorFun, loglikByCopulasLags, fitSpCopula, composeSpCopula)
+export(tailDepFun, lowerTailDepFun, upperTailDepFun)
+export(empTailDepFun, lowerEmpTailDepFun, upperEmpTailDepFun)
 
 # MRP functions
 export(genEmpKenFun, genInvKenFun)

Modified: pkg/R/linkingVineCopula.R
===================================================================
--- pkg/R/linkingVineCopula.R	2013-11-27 15:21:31 UTC (rev 117)
+++ pkg/R/linkingVineCopula.R	2013-12-05 21:39:45 UTC (rev 118)
@@ -172,12 +172,7 @@
   return(matrix(tmp, ncol = 2))
 }
 
-## fit a Copula through fitCopula with BiCopEst
-## -> BB8Copula (collate issues)
 
-
-
-
 # ## transform a fit from VineCopula to a list of copula objects
 # castVineCopula <- function(cdvEst) {
 #   copulas <- NULL

Modified: pkg/R/spVineCopula.R
===================================================================
--- pkg/R/spVineCopula.R	2013-11-27 15:21:31 UTC (rev 117)
+++ pkg/R/spVineCopula.R	2013-12-05 21:39:45 UTC (rev 118)
@@ -132,7 +132,7 @@
   }
   
   return(new("fitCopula", estimate = spVineCop at parameters, var.est = matrix(NA), 
-             method = method, 
+             method = sapply(method,paste,collapse=", "), 
              loglik = sum(l0)+loglik,
              fitting.stats=list(convergence = as.integer(NA)),
              nsample = nrow(data at data), copula=spVineCop))

Added: pkg/R/tailDependenceFunctions.R
===================================================================
--- pkg/R/tailDependenceFunctions.R	                        (rev 0)
+++ pkg/R/tailDependenceFunctions.R	2013-12-05 21:39:45 UTC (rev 118)
@@ -0,0 +1,39 @@
+# adopted from http://www.r-bloggers.com/copulas-and-tail-dependence-part-1/, 04.11.2013
+
+lowerEmpTailDepFun <- function(u) {
+  empFun <- function(x) sum((u[,1]<=x)&(u[,2]<=x))/sum(u[,1]<=x)
+  function(x) sapply(x,empFun)
+}
+
+upperEmpTailDepFun <- function(u) {
+  empFun <- function(x) sum((u[,1]>=x)&(u[,2]>=x))/sum(u[,1]>=x)
+  function(x) sapply(x,empFun)
+}
+
+empTailDepFun <- function(u) {
+  function(z) {
+    res <- z
+    res[z>0.5] <- upperEmpTailDepFun(u)(z[z>0.5])
+    res[z<=0.5] <- lowerEmpTailDepFun(u)(z[z<=0.5])
+    return(res)
+  }
+}
+
+##
+
+lowerTailDepFun <- function(copula) {
+  function(z) pCopula(cbind(z,z),copula)/z
+}
+
+upperTailDepFun <- function(copula) {
+  function(z) (1-2*z+pCopula(cbind(z,z),copula))/(1-z)
+}
+
+tailDepFun <- function(copula) {
+  function(z) {
+    res <- z
+    res[z>0.5] <- upperTailDepFun(copula)(z[z>0.5])
+    res[z<=0.5] <- lowerTailDepFun(copula)(z[z<=0.5])
+    return(res)
+  }
+}
\ No newline at end of file

Added: pkg/R/tawn3pCopula.R
===================================================================
--- pkg/R/tawn3pCopula.R	                        (rev 0)
+++ pkg/R/tawn3pCopula.R	2013-12-05 21:39:45 UTC (rev 118)
@@ -0,0 +1,128 @@
+#######################################
+## tawn copula with all 3 parameters ##
+#######################################
+
+setClass("tawn3pCopula", representation(exprdist = "expression"),
+         contains = "evCopula")
+
+Atawn3p <- function(t, param = c(0.9302082, 1, 8.355008)) {
+  alpha <- param[1]
+  beta  <- param[2]
+  theta <- param[3]
+  (1-beta)*(t) + (1-alpha)*(1-t) + ((alpha*(1-t))^theta+(beta*t)^theta)^(1/theta)
+
+}
+
+ATawn <- function(copula, w) {
+  Atawn3p(w,copula at parameters)
+}
+
+setMethod("A",signature("tawn3pCopula"),ATawn)
+
+dAduTawn <- function(copula, w) {
+  alpha <- copula at parameters[1]
+  beta  <- copula at parameters[2]
+  theta <- copula at parameters[3]
+
+  # 1st derivative
+  p1 <- (alpha*(alpha*(-(w-1)))^(theta-1)-beta*(beta*w)^(theta-1)) 
+  p2 <- ((alpha*(-(w-1)))^theta+(beta*w)^theta)^(1/theta-1)
+  
+  # 2nd derivative
+  p3 <- (alpha*(-(w-1)))^(theta-2)
+  p4 <- (beta*w)^(theta-2)
+  p5 <- ((alpha*(-(w-1)))^theta+(beta*w)^theta)^(1/theta-2)
+  
+  data.frame(der1=alpha-beta-p1*p2,
+             der2=alpha^2*beta^2*(theta-1)*p3*p4*p5)
+}
+
+setMethod("dAdu",signature("tawn3pCopula"),dAduTawn)
+
+tawn3pCopula <- function (param = c(0.5, 0.5, 2)) {
+  # A(t) = (1-beta)*t + (1-alpha)*(1-t) + ((alpha*(1-t))^theta+(beta*t)^theta)^(1/theta)
+  # C(u1,u2) = exp(log(u1*u2) * A(log(u2)/log(u1*u2)))
+  #          = u1*u2 + exp(A(log(u2)/log(u1*u2)))
+
+  cdf <- expression(exp(log(u1*u2)*((1-beta)*(log(u2)/log(u1*u2)) +
+                                    (1-alpha)*(1-log(u2)/log(u1*u2)) +
+                                    ((alpha*(1-log(u2)/log(u1*u2)))^theta+(beta*log(u2)/log(u1*u2))^theta)^(1/theta))))
+  dCdU1 <- D(cdf, "u1")
+  pdf <- D(dCdU1, "u2")
+  
+  new("tawn3pCopula", dimension = 2L, exprdist = c(cdf = cdf, pdf = pdf),
+      parameters = param, param.names = c("alpha", "beta", "theta"), 
+      param.lowbnd = c(0,0,1), param.upbnd = c(1,1,Inf), 
+      fullname = "Tawn copula family with three parameters; Extreme value copula")
+}
+
+dtawn3pCopula <- function(u, copula, log=FALSE, ...) {
+  dim <- copula at dimension
+  for (i in 1:dim) {
+    assign(paste("u", i, sep=""), u[,i])
+  }
+  alpha <- copula at parameters[1]
+  beta  <- copula at parameters[2]
+  theta <- copula at parameters[3]
+  
+  val <- c(eval(copula at exprdist$pdf))
+  ## improve log-case
+  if(log) 
+    log(val) 
+  else 
+    val
+}
+
+setMethod("dCopula", signature("matrix", "tawn3pCopula"), dtawn3pCopula)
+setMethod("dCopula", signature("numeric", "tawn3pCopula"),dtawn3pCopula)
+
+ptawn3pCopula <- function(u, copula, ...) {
+  dim <- copula at dimension
+  for (i in 1:dim) {
+    assign(paste("u", i, sep=""), u[,i])
+  }
+  alpha <- copula at parameters[1]
+  beta <-  copula at parameters[2]
+  theta <-copula at parameters[3]
+  
+  val <- c(eval(copula at exprdist$cdf))
+}
+
+setMethod("pCopula", signature("matrix", "tawn3pCopula"),  ptawn3pCopula)
+setMethod("pCopula", signature("numeric", "tawn3pCopula"), ptawn3pCopula)
+
+setMethod("rCopula", signature("numeric", "tawn3pCopula"),
+          function(n, copula, ...) {
+            copula:::revCopula(n, copula, ...)[,2:1]
+          })
+
+fitTawn3pCop <- function(copula, data, method = c("mpl", "ml", "itau", "irho"), 
+                         start = NULL, lower = NULL, upper = NULL, optim.method = "BFGS", 
+                         optim.control = list(maxit = 1000), estimate.variance = FALSE, 
+                         hideWarnings = TRUE) {
+  if(is.null(start))
+    start <- copula at parameters
+  
+  # copied from copula::fitCopula
+  if (!is.matrix(data)) {
+    warning("coercing 'data' to a matrix.")
+    data <- as.matrix(data)
+    stopifnot(is.matrix(data))
+  }
+  
+  cat(match.arg(method))
+  
+  res <- switch(match.arg(method), 
+               ml = copula:::fitCopula.ml(copula, data, start, lower, upper, 
+                                          optim.method, optim.control, 
+                                          estimate.variance, hideWarnings),
+               mpl = copula:::fitCopula.mpl(copula, data, start, lower, upper, 
+                                            optim.method, optim.control, 
+                                            estimate.variance, hideWarnings),
+               itau = copula:::fitCopula.itau(copula, data, estimate.variance),
+               irho = copula:::fitCopula.irho(copula, data, estimate.variance))
+  cat(match.arg(method))
+  return(res)
+}
+            
+setMethod("fitCopula", signature("tawn3pCopula"), fitTawn3pCop)
\ No newline at end of file

Modified: pkg/R/vineCopulas.R
===================================================================
--- pkg/R/vineCopulas.R	2013-11-27 15:21:31 UTC (rev 117)
+++ pkg/R/vineCopulas.R	2013-12-05 21:39:45 UTC (rev 118)
@@ -43,9 +43,9 @@
 
 ## density ##
 
-dRVine <- function(u, copula, log=F) {
+dRVine <- function(u, copula, log=FALSE) {
   RVM <- copula at RVM
-  vineLoglik <- RVineLogLik(u, RVM, separate=T)$loglik
+  vineLoglik <- RVineLogLik(u, RVM, separate=TRUE)$loglik
   if(log)
     return(vineLoglik)
   else
@@ -89,14 +89,18 @@
 # fitting using RVine
 fitVineCop <- function(copula, data, method) {
   stopifnot(copula at dimension==ncol(data))
+  if(!is.null(method[["familyset"]]))
+    familyset <- method[["familyset"]]
+  else
+    familyset <- NA
   if("StructureSelect" %in% method)
-    vineCop <- vineCopula(RVineStructureSelect(data, indeptest="indeptest" %in% method))
+    vineCop <- vineCopula(RVineStructureSelect(data, familyset, indeptest="indeptest" %in% method))
   else
-    vineCop <- vineCopula(RVineCopSelect(data, Matrix=copula at RVM$Matrix, 
+    vineCop <- vineCopula(RVineCopSelect(data, familyset, copula at RVM$Matrix, 
                                          indeptest="indeptest" %in% method))
   
   return(new("fitCopula", estimate = vineCop at parameters, var.est = matrix(NA), 
-             method = method, 
+             method = sapply(method,paste,collapse=", "), 
              loglik = RVineLogLik(data, vineCop at RVM)$loglik,
              fitting.stats=list(convergence = as.integer(NA)),
              nsample = nrow(data), copula=vineCop))

Added: pkg/man/tailDepFun.Rd
===================================================================
--- pkg/man/tailDepFun.Rd	                        (rev 0)
+++ pkg/man/tailDepFun.Rd	2013-12-05 21:39:45 UTC (rev 118)
@@ -0,0 +1,70 @@
+\name{tailDepFun}
+\alias{tailDepFun}
+\alias{lowerTailDepFun}
+\alias{upperTailDepFun}
+
+\alias{empTailDepFun}
+\alias{lowerEmpTailDepFun}
+\alias{upperEmpTailDepFun}
+
+\title{
+Tail dependence functions
+}
+\description{
+Functions returning a (empirical) tail dependence function for a copula (sample). The tail dependence functions can be upper, lower or joint, where values below 0.5 are calculated from the lower tail dependence function and values larger 0.5 for the upper tail dependence function.
+}
+\usage{
+tailDepFun(copula)
+lowerTailDepFun(copula)
+upperTailDepFun(copula)
+
+empTailDepFun(u)
+lowerEmpTailDepFun(u)
+upperEmpTailDepFun(u)
+}
+
+\arguments{
+  \item{copula}{
+an object of class \code{\linkS4class{copula}}
+}
+  \item{u}{
+a bivariate sample on (0,1)
+}
+}
+
+\value{
+A function taking arguments from the unit interval (0,1) and returning the corresponding tail index.
+}
+\references{
+Inspired by: \url{http://www.r-bloggers.com/copulas-and-tail-dependence-part-1/}
+}
+\author{
+Benedikt Graeler
+}
+
+\examples{
+data(simulatedTriples)
+X <- rankTransform(triples[,c(1,3)])
+  
+tdfEmp <- empTailDepFun(X)
+plot(tdfEmp,ylim=c(0,1),
+     ylab="tail dependence index")  
+abline(v=0.5, col="grey")
+
+smplTau <- cor(X,method="kendall")[1,2]
+
+# Gauss
+tdfGauss <- tailDepFun(normalCopula(sin(smplTau*pi/2)))
+curve(tdfGauss,add=TRUE,col="blue")
+
+# survival Gumbel
+tdfGumbel <- tailDepFun(surGumbelCopula(1/(1-smplTau)))
+curve(tdfGumbel,add=TRUE,col="darkgreen")
+
+# survival BB6 copula
+tdfBB6 <- tailDepFun(surBB6Copula(c(4.65,2.28)))
+curve(tdfBB6,add=TRUE,col="red")
+
+legend("bottomleft",c("empircal","Gauss","surv. Gumbel","surv. BB6"),
+       col=c("black","blue","darkgreen","red"),lty=1)
+}
\ No newline at end of file

Added: pkg/man/tawn3pCopula-class.Rd
===================================================================
--- pkg/man/tawn3pCopula-class.Rd	                        (rev 0)
+++ pkg/man/tawn3pCopula-class.Rd	2013-12-05 21:39:45 UTC (rev 118)
@@ -0,0 +1,60 @@
+\name{tawn3pCopula-class}
+\Rdversion{1.1}
+\docType{class}
+\alias{tawn3pCopula-class}
+\alias{A,tawn3pCopula-method}
+\alias{dAdu,tawn3pCopula-method}
+\alias{dCopula,matrix,tawn3pCopula-method}
+\alias{dCopula,numeric,tawn3pCopula-method}
+\alias{fitCopula,tawn3pCopula-method}
+\alias{pCopula,matrix,tawn3pCopula-method}
+\alias{pCopula,numeric,tawn3pCopula-method}
+\alias{rCopula,numeric,tawn3pCopula-method}
+
+\title{Class \code{"tawn3pCopula"}}
+\description{
+S4-Class representing the three parameter Tawn Copula family. 
+}
+\section{Objects from the Class}{
+Objects can be created by calls of the form \code{new("tawn3pCopula", ...)} or more likely with the help of the constructor \code{\link{tawn3pCopula}}.
+}
+\section{Slots}{
+  \describe{
+    \item{\code{exprdist}:}{Object of class \code{"expression"} defining the cdf and pdf of the copula.}
+    \item{\code{dimension}:}{Object of class \code{"integer"} defining the dimension (currently fixed at 2.) }
+    \item{\code{parameters}:}{Object of class \code{"numeric"} providing the three parameters. }
+    \item{\code{param.names}:}{Object of class \code{"character"} giving the names of the parameters. }
+    \item{\code{param.lowbnd}:}{Object of class \code{"numeric"} defining the lower bounds of the three parameters. }
+    \item{\code{param.upbnd}:}{Object of class \code{"numeric"} defining the upper bounds of the three parameters. }
+    \item{\code{fullname}:}{Object of class \code{"character"} providing a descriptive name for the copula family.}
+  }
+}
+\section{Extends}{
+Class \code{"\linkS4class{evCopula}"}, directly.
+Class \code{"\linkS4class{copula}"}, by class "evCopula", distance 2.
+Class \code{"\linkS4class{Copula}"}, by class "evCopula", distance 3.
+}
+\section{Methods}{
+  \describe{
+    \item{A}{\code{signature(copula = "tawn3pCopula")}: ... }
+    \item{dAdu}{\code{signature(copula = "tawn3pCopula")}: ... }
+    \item{dCopula}{\code{signature(u = "matrix", copula = "tawn3pCopula")}: ... }
+    \item{dCopula}{\code{signature(u = "numeric", copula = "tawn3pCopula")}: ... }
+    \item{fitCopula}{\code{signature(copula = "tawn3pCopula")}: ... }
+    \item{pCopula}{\code{signature(u = "matrix", copula = "tawn3pCopula")}: ... }
+    \item{pCopula}{\code{signature(u = "numeric", copula = "tawn3pCopula")}: ... }
+    \item{rCopula}{\code{signature(n = "numeric", copula = "tawn3pCopula")}: ... }
+	 }
+}
+
+\author{
+Benedikt Graeler
+}
+
+\seealso{
+\code{\linkS4class{tawnCopula}} for a symmetric one parameter, \code{\linkS4class{tawnT1Copula}} and \code{\linkS4class{tawnT2Copula}} for asymmetric two-parameter sub-families of the Tawn copula family. 
+}
+\examples{
+showClass("tawn3pCopula")
+}
+\keyword{classes}

Added: pkg/man/tawn3pCopula.Rd
===================================================================
--- pkg/man/tawn3pCopula.Rd	                        (rev 0)
+++ pkg/man/tawn3pCopula.Rd	2013-12-05 21:39:45 UTC (rev 118)
@@ -0,0 +1,42 @@
+\name{tawn3pCopula}
+\alias{tawn3pCopula}
+
+\title{
+Tawn Copula Family constructor using all three parameters
+}
+\description{
+Constructor of the Tawn Copula family with three parameters allowing for asymmetric copula families.
+}
+\usage{
+tawn3pCopula(param = c(0.5, 0.5, 2))
+}
+
+\arguments{
+  \item{param}{
+  The vector defining the three parameters of the Tawn Copula.
+}
+}
+\details{
+Pickand's dependece function is given through: 
+A(t) = (1-beta)*t + (1-alpha)*(1-t) + ((alpha*(1-t))^theta+(beta*t)^theta)^(1/theta)
+}
+\value{
+An instance of the \code{\linkS4class{tawn3pCopula}} class.
+}
+
+\author{
+Benedikt Graeler
+}
+
+\seealso{
+\code{\link{tawnCopula}} for a symmetric one parameter, \code{\link{tawnT1Copula}} and \code{\link{tawnT2Copula}} for asymmetric two-parameter sub-families of the Tawn copula family. 
+}
+
+\examples{
+persp(tawn3pCopula(c(0.4,0.9,4)), dCopula)
+persp(tawn3pCopula(c(0.4,0.9,4)), pCopula)
+
+}
+
+\keyword{ copula }
+\keyword{ distribution }



More information about the spcopula-commits mailing list