Calculate the Area under a Curve

Question

I would like to calculate the area under a curve to do integration without defining a function such as in integrate().

My data looks as this:

Date          Strike     Volatility
2003-01-01    20         0.2
2003-01-01    30         0.3
2003-01-01    40         0.4
etc.

I plotted plot(strike, volatility) to look at the volatility smile. Is there a way to integrate this plotted "curve"?

Answer 1

The AUC is approximated pretty easily by looking at a lot of trapezium figures, each time bound between x_i, x_{i+1}, y{i+1} and y_i. Using the rollmean of the zoo package, you can do:

library(zoo)

x <- 1:10
y <- 3*x+25
id <- order(x)

AUC <- sum(diff(x[id])*rollmean(y[id],2))

Make sure you order the x values, or your outcome won't make sense. If you have negative values somewhere along the y axis, you'd have to figure out how exactly you want to define the area under the curve, and adjust accordingly (e.g. using abs() )

Regarding your follow-up : if you don't have a formal function, how would you plot it? So if you only have values, the only thing you can approximate is a definite integral. Even if you have the function in R, you can only calculate definite integrals using integrate(). Plotting the formal function is only possible if you can also define it.

Answer 2

Just add the following to your program and you will get the area under the curve:

require(pracma)
AUC = trapz(strike,volatility)

From ?trapz:

This approach matches exactly the approximation for integrating the function using the trapezoidal rule with basepoints x.

Answer 3

Three more options, including one using a spline method and one using Simpson's rule...

# get data
n <- 100
mean <- 50
sd <- 50

x <- seq(20, 80, length=n)
y <- dnorm(x, mean, sd) *100

# using sintegral in Bolstad2
require(Bolstad2)
sintegral(x,y)$int

# using auc in MESS
require(MESS)
auc(x,y, type = 'spline')

# using integrate.xy in sfsmisc
require(sfsmisc)
integrate.xy(x,y)

The trapezoidal method is less accurate than the spline method, so MESS::auc (uses spline method) or Bolstad2::sintegral (uses Simpson's rule) should probably be preferred. DIY versions of these (and an additional approach using the quadrature rule) are here: http://www.r-bloggers.com/one-dimensional-integrals/

Answer 4

OK so I arrive a bit late at the party but going over the answers a plain R solution to the problem is missing. Here goes, simple and clean:

sum(diff(x) * (head(y,-1)+tail(y,-1)))/2

The solution for OP then reads as:

sum(diff(strike) * (head(volatility,-1)+tail(volatility,-1)))/2

This effectively calculates the area using the trapezoidal method by taking the average of the "left" and "right" y-values.

NB: as @Joris already pointed out you could use abs(y) if that would make more sense.

Answer 5

In the pharmacokinetics (PK) world, calculating different types of AUC is a common and fundamental task. The are lots of different AUC calculations for pharmacokietics, such as

• AUC0-t = AUC from zero to time t
• AUC0-last = AUC from zero to the last time point (may be same as above)
• AUC0-inf = AUC from zero to time infinity
• AUCint = AUC over a time interval
• AUCall = AUC over the whole time period for which data exists

One of the best packages which does these calculations is the relatively new package PKNCA from the folks at Pfizer. Check it out.

Answer 6

Joris Meys's answer was great but I struggled to remove NAs from my samples. Here is the little function I wrote to deal with them :

library(zoo) #for the rollmean function

######
#' Calculate the Area Under Curve of y~x
#'
#'@param y Your y values (measures ?)
#'@param x Your x values (time ?)
#'@param start : The first x value 
#'@param stop : The last x value
#'@param na.stop : returns NA if one value is NA
#'@param ex.na.stop : returns NA if the first or the last value is NA
#'
#'@examples 
#'myX = 1:5
#'myY = c(17, 25, NA, 35, 56)
#'auc(myY, myX)
#'auc(myY, myX, na.stop=TRUE)
#'myY = c(17, 25, 28, 35, NA)
#'auc(myY, myX, ex.na.stop=FALSE)
auc = function(y, x, start=first(x), stop=last(x), na.stop=FALSE, ex.na.stop=TRUE){
  if(all(is.na(y))) return(NA)
  bounds = which(x==start):which(x==stop)
  x=x[bounds]
  y=y[bounds]
  r = which(is.na(y))
  if(length(r)>0){
    if(na.stop==TRUE) return(NA)
    if(ex.na.stop==TRUE & (is.na(first(y)) | is.na(last(y)))) return(NA)
    if(is.na(last(y))) warning("Last value is NA, so this AUC is bad and you should feel bad", call. = FALSE) 
    if(is.na(first(y))) warning("First value is NA, so this AUC is bad and you should feel bad", call. = FALSE) 
    x = x[-r]
    y = y[-r]
  }
  sum(diff(x[order(x)])*rollmean(y[order(x)],2))
}

I then use it with an apply onto my dataframe : myDF$auc = apply(myDF, MARGIN=1, FUN=auc, x=c(0,5,10,15,20))

Hope it can help noobs like me :-)

EDIT : added bounds

Answer 7

You can use ROCR package, where the following lines will give you the AUC:

pred <- prediction(classifier.labels, actual.labs)
attributes(performance(pred, 'auc'))$y.values[[1]]