Calculate the Area under a Curve

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.

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.

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/

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.

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.