Numerical Integration

Highlights

  • Numerical integration
    • Numerical univariate integration
    • Implemented as a Mata class
  • Additive quadrature model
    • Adaptive Gauss—Kronrod method
    • Implemented as a Mata class
  • Robustness to singular points

Stata provides statistical solutions developed by StataCorp, and it provides programming tools for those who want to develop their own solutions. There are two Stata programming languages: ado, which is easy to use, and Mata, which performs numerical heavy lifting. And Stata is integrated with Python.

Mata's new Quadrature() class provides adaptive Gaussian quadrature for numerically integrating univariate functions. It approximates the integral from a to b of f(x), where a can be minus infinity or finite and b can be finite or positive infinity.

Quadrature() uses the adaptive Gauss—Kronrod method. It also provides the Simpson method for use in teaching.

Let's see it work

To use Quadrature(), you

  1. define the function to be integrated,
  2. set the lower and upper limits of integration, and
  3. call the Quadrature() to perform the integration.

For instance, imagine we want to integrate the normal density function even though Mata has a built-in cumulative normal. Define the function to calculate f(x) by typing

real scalar myfunc(real scalar x)
{
      return(normalden(x))
}

or, if you prefer,

real scalar myfunc(real scalar x)
{
      return( (1/sqrt(2*pi()))*exp(-(x^2/2)) )
}

To obtain the integral, you might write

real scalar Integrate_myfunc(real scalar a, real scalar b)
{
    class Quadrature scalar q

    q.setEvaluator(&myfunc())
    q.setLimits((a, b))

    return(q.integrate())
}

To obtain the integral from 0 to infinity, code Integrate_myfunc(0, .). Returned will be 0.49999999999936, a result accurate to 12 digits.

If you needed a result accurate to 13 digits, you could set an absolute or relative tolerance. You could call q.setAbstol(1e-13) before calling q.integrate(). Or you could call q.setReltol(1e-13). Because the answer is 0.5 in this case, absolute and relative tolerances are about the same.

Post your comment

Timberlake Consultants