Chapter 2

Numerical Integration

INTRODUCTION

Integration plays a central role in the mathematical models that have been developed to price options. Integration usually involves finding the antiderivative of a function and applying the Fundamental Theorem of Calculus, according to which the antiderivative is evaluated at the endpoints of the integral only. This well-known analytical procedure yields the area under the curve defined by the function, in between the two endpoints.

Unfortunately, while many option pricing models require evaluating the integral of a function to find the area under the curve, it is often impossible to find an antiderivative for the function. One example is the option pricing model of Heston and Nandi (2000) reviewed in Chapter 6, for which they find a semiclosed form solution for an option price, on an underlying asset that has volatility defined by a generalized autoregressive conditional heteroskedasticity (GARCH) process. The solution is in semiclosed form because finding the option price involves a complex integral that must be evaluated by numerical integration. We will need also to apply numerical integration to calculate model-free volatility, skewness, and kurtosis, in later chapters.

This chapter introduces numerical integration, presents algorithms that are used to perform numerical integration, and illustrates the VBA code to implement these algorithms. There are two main families of algorithms that have been developed for numerical integration: Newton-Cotes formulas, also called quadratures, and Gaussian quadratures. The main difference between these algorithms is the way in which the function being integrated is divided into subintervals along the horizontal axis. Newton-Cotes formulas divide the axis into subintervals at each integration point, or abscissa. Lagrange polynomials of varying complexity are used to approximate the areas under the function between successive pairs of abscissas. These areas are then aggregated to produce an approximation to the original integral. The subintervals are usually chosen to be of equal width.

Gaussian quadratures use a much more complex algorithm to approximate the integral. The chief difference is that the abscissas are determined not by the analyst, but by the algorithm itself, as the roots of approximating polynomials. Gaussian quadratures involve the sum of the product of the function evaluated at the optimal abscissas and a weighting function. The fundamental theorem of Gaussian quadratures defines the abscissas as the roots of the orthogonal polynomial for the same interval and weighting function. Although more complicated than Newton-Coates formulas, Gaussian quadratures are more accurate because they fit polynomials up to degree 2n − 1 exactly, where n is the number of abscissas. With today’s high computer speeds, however, a very large number of points can be chosen, and Newton-Cotes approximations can be computed very quickly. In light of this, for many financial analysts the increased accuracy of the Gaussian quadratures does not justify the extra programming effort required to implement them. It is therefore not surprising that the more straightforward Newton-Cotes formulas are the most widely used of numerical integration algorithms. Consequently, in this chapter we focus mostly on Newton-Cotes formulas. In later chapters, however, the superiority of Gaussian quadratures over Newton-Cotes formulas will become evident.

NEWTON-COATES FORMULAS

In this section, seven popular Newton-Coates numerical integration formulas are introduced: the left-, right-, and midpoint rules, the trapezoidal rule, Simpson’s rule and Simpson’s three-eighths rule, and Boole’s rule. We present the VBA code for implementing each rule, and compare the accuracy of each using a simple example. Finally, we assess how the accuracy of each approximation increases as the number of integration points increases.

The Left-, Right- and Midpoint Rules

Finding the integral of a function is equivalent to finding the area under the curve traced by that function. One way to approximate that area is by dividing the area into rectangles and adding the rectangles together. Suppose we wish to approximate the area under the curve traced by the function f(x) = x3 + 10. This is, of course, a simple function for which an antiderivate exists. It can be integrated exactly, but it is used throughout this chapter to illustrate the methodology of each rule. To approximate the area under the curve within the interval [a, b] on the horizontal axis (x-axis), n + 1 equally spaced integration points, or abscissas, are used to divide the interval into a set of n subintervals of equal width:

where x0 ≡ a and xn ≡ b are the endpoints. A rectangle is constructed in each subinterval and its area is obtained. Aggregating the areas of the rectangles produces the approximation. The left-, right-, and midpoint rules differ only by how the height of the each rectangle is chosen.

According to the right-point rule, the height of each rectangle is defined as the value of the function f(x) at the rightmost endpoint of each subinterval. Figure 2.1 illustrates the right-point rule using our function as an example, with a = −2 and b = 2. This rule can be expressed mathematically as the approximation

FIGURE 2.1 Right-Point Rule Numerical Integration

(2.1)

where

Δx = (b − a)/n = the width of each subinterval

xi = the rightmost point of each subinterval

a, b = endpoints of the interval over which the integral is evaluated

f (xi) = the value of the function evaluated at xi

n = the number of subintervals.

Note that each abscissas xi can be expressed in terms of the leftmost endpoint a and the subinterval width Δx as xi = a + iΔx, for i = 0, 1, 2, …, n. This convenient representation will be exploited later in this chapter, when we implement numerical integration with VBA.

The left-point rule sets the height of each rectangle to the leftmost value of f(x) in each subinterval, so that the approximation is

(2.2)

Figure 2.2 illustrates our example using the left-point rule approximation rule. Figures 2.1 and 2.2 demonstrate the shortcomings of the left- and right-point rules. Although these rules each produce a numerical approximation to the integral, a substantial amount of information is lost between the n + 1 points defined by the subintervals. It easy to see that when the function is monotone nondecreasing (as in our example), the right-point rule produces an approximation that is too large, while the left-point rule produces one that is too small. The reverse is true for a nonincreasing monotone function.

FIGURE 2.2 Left-Point Rule Numerical Integration

There are two ways to incorporate information ignored by the left- and right-point rules. The first method requires new function outputs, namely the y-values defined at the midpoint of each subinterval. This is aptly named the midpoint rule. The midpoint rule is identical to the left- and right-point rules, except that it uses the value of f(x) obtained at the midpoint of each subinterval for the height of the rectangles.

Mathematically, the midpoint rule can be expressed as

(2.3)

where

As illustrated by Figure 2.3, the increase in accuracy brought on by choosing the midpoint rule over either the right-point or left-point rule is substantial. Nevertheless, these three methods are relatively crude compared to the methods that follow.

FIGURE 2.3 Midpoint Rule Numerical Integration

The Trapezoidal Rule

The trapezoidal rule does not require additional information other than the value of f(x) at each abscissa. Implementing this rule is no more complicated than implementing the left-, right-, or midpoint rules. The interval [a, b] is again divided into a set of n subintervals of equal widthΔx. The goal of this approach, however, is to approximate the rate at which the function f(x) is increasing or decreasing within each subinterval. This is achieved by forming a line segment to join f(x) at the endpoints of each subinterval, which produces a trapezoid instead of a rectangle. The area of the trapezoid defined by the subinterval (xi−1, xi) is

Aggregating this quantity for all subintervals produces the trapezoidal rule, which can be expressed mathematically as

(2.4)

Because trapezoids are used instead of rectangles, the accuracy of the trapezoidal approximation is much higher than that achieved with the midpoint rule. Figure 2.4 illustrates our example with the trapezoidal rule...