SECTION 12.

INTEGRATION METHODS

The present section describes a variety of methods for numerical integration in 1, 2 and 3 dimensions. All input and output variables are REAL*8.

SECTION 12.1 INTEGRATION OF THE NORMAL DENSITY

12.1.1 FUNCTION BIV2 AND BIV2U

BIV2 evaluates the lower tail integral of the bivariate normal
density which has been normalized to have zero means and unit
variances, i.e., from (-infinity,-infinity) to the point (AX,AY).
BIV2U evaluates the upper tail integral, i.e., from (AX,AY) to
(infinity,infinity). For description of the method see Owen,
D.B., "Tables for Computing Bivariate Normal
Probabilities," *Annals of Mathematical Statistics*,
27 (1956), 1075-1090. BIV2 and BIV2U are called as

ANSLOW=BIV2(AX,AY,R,IER) ANSHIGH=BIV2U(AX,AY,R,IER)

where AX and AY are the upper limits of integration along the X- and Y- axes and R is the correlation between X and Y.

IER = 0 is normal return.

IER =-3 indicates input error (i.e., R is greater than or equal
to 1.0 in absolute value).

12.1.2 FUNCTIONS DUTT2 AND DUTT3

These compute upper tail integrals for the bivariate and
trivariate normal densities normalized to have zero means and
unit variances. The method is described in Dutt, J.E.,
"Numerical Aspects of Multivariate Normal Probabilities in
Econometric Methods," *Annals of Economic and Social
Measurement*, 17 (1976). 547-561. The calls are

ANSWER=DUTT2(AX,AY,R12) ANSWER=DUTT3(AX,AY,AZ,R12,R13,R23)

respeectively, where AX, AY, AZ are the lower limits of integration and R12, R13, R23 are the correlations X and Y, X and Z and Y and Z respectively. The correlation matrix is not checked for positive definiteness in DUTT3.

12.1.3 FUNCTIONS CLARK2 AND CLARK3

These compute approximations to the upper tail integrals for
the bivariate and trivariate normal densities normalized to have
zero means and unit variances. The computations are based on
Clark, C.E., "The Greatest of a Finite Set of Random
Variables," *Operation Research*, 9 (1961), 145-162.
See also Quandt, R.E., *The Econometrics of Disequilibrium*,
Blackwell. The calls are

ANSWER=CLARK2(AX,AY,R12) ANSWER=CLARK3(AX,AY,AZ,R12,R13,R23)

(see 12.1.2). The correlation matrix is not checked for positive in CLARK3.

SECTION 12.2 NEWTON-COTES INTEGRATION OF GENERAL FUNCTIONS

Newton-Cotes integration in 1, 2, and 3 dimensions over rectangular regions is accomplished by

CALL NEWCT1(AX,BX,NX,FUNC,ANS,IER) CALL NEWCT2(AX,BX,AY,BY,NX,NY,FUNC,ANS,IER) CALL NEWCT3(AX,BX,AY,BY,AZ,BZ,NX,NY,NZ,FUNC,ANS,IER)

where

AX, BX, AY, BY, AZ, BZ = lower and upper limits of integration
along the X, Y, Z axes

NX, NY, NZ = number of points of integration (2.LE.NX,NY,NZ.LE.9)
along the three axes

FUNC = the name of the FUNCTION subprogram that evaluates the
integrand, defined as FUNCTION FUNC(X), or FUNCTION FUNC(X,Y), or
FUNCTION FUNC(X,Y,Z). FUNC must be declared in an EXTERNAL
statement in the calling routine.

ANS = the answer will be returned here.

IER = 0 if normal return,

IER =-3 if input error (i.e., if the limits of integration do not
define intervals of positive width).

12.2.1 WEDDLE'S RULE (WEDDLE)

Weddle's Rule is a Newton-Cotes type integral formula (see F.
B. Hildebrand, *Introduction to Numerical Analysis*,
McGraw-Hill, 1956). WEDDLE is a univariate integration formula
and a 7-point Weddle Rule is used over each of n subintervals.
This is accomplished by

CALL WEDDLE(AX,BX,N,FUNC,ANS,IER)

where AX and BX are the lower and upper limits of integration respectively, N is the number of subintervals to be used from AX to BX, FUNC is the name of the function that evaluates the integrand, and ANS will contain the answer. IER = 0 normal return = -3 AX is greater than BX

SECTION 12.3 GAUSS-LEGENDRE INTEGRATION

12.3.1 N-POINT INTEGRATION

N-point integration over rectangular regions in 1, 2, 3 dimensions is accomplished by

CALL GAULG1(AX,BX,NX,FUNC,ANS,IER) CALL GAULG2(AX,BX,AY,BY,NX,NY,FUNC,ANS,IER) CALL GAULG3(AX,BX,AY,BY,AZ,BZ,NX,NY,NZ,FUNC,ANS,IER)

where the calling parameters have the same meaning as in SECTION 12.2. Note that here NX, NY, NZ must satisfy 2.LE.NX,NY,NZ.LE.11.

12.3.2 ITERATED GAUSS-LEGENDRE INTEGRATION

The calls are

CALL GAUSA(AX,BX,N,FUNC,ANS,IER) CALL BIINT(AX,BX,AY,BY,N,FUNC,ANS,IER) CALL TRINT(AX,BX,AY,BY,AZ,BZ,N,FUNC,ANS,IER) CALL QUINT(AX,BX,AY,BY,AZ,BZ,AW,BW,N,FUNC,ANS,IER)

perform in 1, 2, 3, 4 dimensions respectively, N times iterated 2-point Gauss-Legendre integration over the rectangular region given by AX,...,BZ. All calling parameters have the same meaning as in SECTION 12.2, except that no separate N may be specified for the X, Y, and Z dimensions. Warning: N must be even!! Note: In the Lahey compiler compatible version QUINT is called QUINTA

SECTION 12.4 BIVARIATE INTEGRATION WITH VARIABLE LIMITS

1. BIVAR2. The subroutine BIVAR2 performs 10-point Gaussian integration of the function FUNCA(X,Y) between limits AX and BX along the X-axis and limits AY(X) and BY(X) along the Y-axis. Hence it computes

The call to BIVAR2 is

CALL BIVAR2(AX,BX,FUNCA,ANS,IER)

where

AX = lower limit along X-axis

BX = upper limit along X-axis

FUNCA = name of function that evaluates the integrand

ANS = contains the answer

IER = 0 if AX .LT. BX and AY(X) .LT. BY(X) and

IER = -3 otherwise. In this latter case no integration is
attempted and the value of ANS is meaningless. The function FUNCA
may be named by the user, but the name must be declared in an
EXTERNAL statement in the MAIN program. It should be defined as

FUNCTION FUNCA(X,Y) IMPLICIT REAL*8 (A-H,O-Z) . . . . FUNCA= RETURN END

If the user needs to pass parameters to this function, labelled COMMON should be utilized for the purpose. In addition, the user needs to define two additional functions AY(X) and BY(X) which give the lower and upper limits of integration corresponding to each X-value. These functions should be

FUNCTION AY(X) IMPLICIT REAL*8 (A-H,O-Z) AY= RETURN END FUNCTION BY(X) IMPLICIT REAL*8 (A-H,O-Z) BY= RETURN END

As before, if these functions require parameters, the user can
pass them via labelled COMMON. As an illustration, consider the
FUNCTION H(X,Y) which is to be integrated from AX to BX along the
X-axis and from 0 to exp(-ax) along
the Y-axis, where H is given by exp(bx)y^{2}:

Since the function and the upper Y-limit depend on parameters a and b, they will be passed through labelled COMMON/USER1/A,B. The MAIN program allows the user the input various values of A and B as well as various values of AX and BX. The analytic answer is also computed and is printed along the numerical integral. The listing follows.

IMPLICIT REAL*8 (A-H,O-Z) COMMON/USER1/A,B EXTERNAL H OPEN (UNIT=9,FILE='BIVAR.DAT',STATUS='OLD') READ (9,1000) A,B,AX,BX 1000 FORMAT(4F6.0) WRITE (6,1001) A,B,AX,BX 1001 FORMAT(' A,B,AX,BX ',4F8.3) TEMP=B-3*A C=AX D=BX TRUTH=1./(3.*TEMP)*(DEXP(TEMP*D)-DEXP(TEMP*C)) CALL BIVAR2(AX,BX,H,ANS,IER) WRITE (*,1002) TRUTH,ANS,IER 1002 FORMAT(' TRUTH,ANS,IER ',2F12.6/1X,I3) CLOSE(9) STOP END FUNCTION AY(X) IMPLICIT REAL*8 (A-H,O-Z) AY=0. RETURN END FUNCTION BY(X) IMPLICIT REAL*8 (A-H,O-Z) COMMON/USER1/A,B BY=DEXP(-A*X) RETURN END FUNCTION H(X,Y) IMPLICIT REAL*8 (A-H,O-Z) COMMON/USER1/A,B H=DEXP(B*X)*Y**2 RETURN END

where the file BIVAR.DAT might contain the line 1.0 1.0 0.0 3.0

2. BIVAR3. BIVAR3 is similar to BIVAR2 with two significant exceptions: (a) Instead of 10-point Gaussian integration, it uses iterated Gauss-Legendre Integration (Section 12.3), where the user can set the number of points of integration arbitrarily; (b) Either or both of the limits AX and BX can be infinite (AX -infinity and/or BX +infinity). The call to BIVAR3 is

CALL BIVAR3(AX,BX,NPT,ICASE,INORM,DELX,DELY,FUNCT,NTR,ANS,IER)

Before we describe the meaning of the various parameters, we need to say a word about how BIVAR3 works when either AX or BX is non-finite. If AX is -infinity, BIVAR3 first computes the bivariate integral over the stripe from BX-DELX to BX; then from BX-2*DELX to B-DELX, etc., until the marginal stripe's integral divided by the sum of the integrals over the previous stripes is less than EPPSS (set at 1.D-8). If AX is finite and BX is + infinity, integration begins over the stripe AX to AX+DELX, then from AX+DEL to AX+2*DELX, and proceeds analogously. If AX is - infinity and BX is +infinity, integration is over pairs of stripes; from 0 to DELX and from -DELX to 0; then from DELX to 2*DELX and from -2*DELX to -DELX, and proceeds analogously to the previous cases. For this to work in any of the three cases, it is important that the integrand be "single peaked;" say we integrate from 0 to +infinity, and DELX = 1; Assume that the integrand is positive from 0 to 1, 0 from 1 to 2, and again positive from 2 to 3. Then BIVAR3 will quit when the upper limit is 2, which is an error. The arguments in the call to BIVAR3 have the following meaning:

AX = lower limit (unless -infinity, in which case it need not
bet set by the user)

BX = upper limit (unless + infinity, in which case it need not be
set by the user)

NPT = Number of point of integration; MUST BE AN EVEN NUMBER

ICASE = 0 both AX and BX are finite and must be specified;

ICASE = 1 AX is -infinity and need not be specified

ICASE = 2 BX is + infinity and need not be specified

ICASE = 3 AX = -infinity, BX = + infinity, and neither needs to
be specified

INORM = 0 if the user will provide a function FUNCT (which must
be declared in an EXTERNAL statement;

INORM = 1 if the integrand is to ve the standard bivariate normal
density function

DELX = width of horizontal stripe (needs to be specified only if
AX or BX is non-finite)

DELY = irrelevant and need not be specified (but a place holder
must be provided for it in the CALL BIVAR3 statement); this is
present for future modifications or if the user wants to provide
for infinite limits in AY(X) or BY(X)

FUNCT = the name of the FUNCTION subprogram that provides the
value of the integrand (must be be defined as FUNCTION FUNCT(X,Y)
and parameters must be passed to it inm labelled COMMON)

NTR = Number of stripe width increases permitted in the initial
stripe until a nonzero integral is obtained (If in the initial
stripe the integral is .LE. 0, the stripe width is increased and
BIVAR3 tries again. A value of 10 might be reasonable.

ANS = the answer

IER = error return;

= 0 for normal return,

= -3 if AX .GE. BX, NPT is not even, or DELX .LE. 0. Note that
appropriate values of DELC and NTR (and also NPT) are problem
dependent, and experimenation with alternative values may be
useful.

SECTION 12.5 UNIVARIATE INTEGRATION OF N(0,1) FROM -INFINITY TO X (DERF)

To compute the integral of the standard normal density from minus infinity to a value x, use

RESULT=0.5+0.5*DERF(X/DSQRT(2.D0)

where both X and RESULT are REAL*8.

SECTION 12.6 ROMBERG INTEGRATION

To compute a univariate integral by Romberg integration, use

CALL ROMB(A,B,EPS,RES,FUNC,IER)

where

A = lower limit

B = upper limit

EPS = tolerance (e.g.,1.D-5)

RES = result

FUNC= name of function evaluating the integrand (must be declared
in an EXTERNAL statement)

IER = error code, = 0 for notmal return

IER =-3 input error (EPS .LE. 0 or A. GE. B)

All real variables must be declared REAL*8. See, for example,
Herbert S. Wilf, "Advances in Numerical Quadrature," in
*Mathematical Methods for Digital Computers*, Vol. II,
1966.

SECTION 12.7 GAUSS-LAGUERRE INTEGRATION

To integrate the function f(x) from A to infinity, write a Fortran FUNCTION F(X) which evaluates exp(x)f(x), declare it in an EXTERNAL statement in the calling routine, and

CALL GAUSLA(A,N,F,ANS,IER)

where A, N, F are inputs, ANS, IER are outputs, and

A = lower limit of integration (REAL*8)

N = integer number of points of evaluation (2.LE.N.LE.11)

F = REAL*8 function defined above

ANS = value of the integral (REAL*8)

IER = error code (= 0 normal return, = -3 input error)