help mata _bcsf_bracketing-------------------------------------------------------------------------------

Title

_bcsf_bracketing() -- Bracketing and bracket convergence functions usedby censlope

Syntax

real scalar_bcsf_bracketing(yfun,y,yfarleft,yfarright,mbracket,scalefactor,bracketmat,leftindex,rightindex)

real scalar_bcsf_bisect(objfun,x0,x1,y0,y1,itcount,iterate,tolerance[,log])

real scalar_bcsf_regula(objfun,x0,x1,y0,y1,itcount,iterate,tolerance[,log])

real scalar_bcsf_ridders(objfun,x0,x1,y0,y1,itcount,iterate,tolerance[,log])where

yfun:pointer (real scalar function) scalary:real scalaryfarleft:real scalaryfarright:real scalarmbracket:real scalarscalefactor:real scalarbracketmat:real matrixleftindex:real scalarrightindex:real scalarobjfun:pointer (real scalar function) scalarx0:real scalarx1:real scalary0:real scalary1:real scalaritcount:real scalariterate:real scalartolerance:real scalarlog:real scalar

DescriptionThe prefix

_bcsfstands for "bracket convergence for step functions". The programcenslopeuses these Mata functions for the numerical solution of inequalities involving step functions. If a scalar object functionG()is a step function (instead of being a continuous function), then there may not be a unique solution to the equationG(x)=0. Instead, there may be either no solution or a nonempty interval of solutions. However, ifG()is monotonically nonincreasing, then there may be a supremum for the set of values ofxsuch thatG(x)is positive (or nonnegative), which is also the infimum for values ofxsuch thatG(x)is nonpositive (or negative). Similarly, ifG()is monotonically nondecreasing, then there may be a supremum for the set of values ofxsuch thatG(x)is negative (or nonpositive), which is also the infimum for values ofxsuch thatG(x)is nonnegative (or positive). Examples of such step functionsG()include functions derived as continuous monotonic transformation of cumulative distribution functions (CDFs). The suprema and infima will then be percentiles, percentile differences, percentile ratios, or percentile slopes. The function_bcsf_bracketing()is used to find pairs of X values that bracket the suprema and infima. When such a pair of X values is found, the functions_bcsf_bisect(),_bcsf_regula(), and_bcsf_ridders()can be used to redefine these X brackets progressively, until the relative difference between the upper and lowerx-brackets has been reduced to a level at or below a user-defined tolerance level. The brackets are then said to have converged. Bracket convergence is achieved using variants of standard numerical methods used to find solutions for scalar equations involving continuous functions, such as the methods described in chapter 9 of Press et al. (1992). All the_bcsffunctions return an integer-valued scalar return code, indicating the outcome of their attempts to carry out their bracketing or convergence tasks.

_bcsf_bracketing(yfun,y,yfarleft,yfarright,mbracket,scalefactor,bracketmat,leftindex,rightindex)inputs a Y value and finds a pair of X values whose corresponding Y values bracket the input Y value. The argumentsyfun,y,yfarleft,yfarright,mbracket, andscalefactorare input, the argumentsleftindexandrightindexare output, and the argumentbracketmatis a two-column bracket matrix, which must have at least two rows on input and may be extended with additional rows on output. A bracket matrix is a two-column matrix whose first column is assumed to contain a set of X values in ascending order and whose second column is assumed to contain the corresponding Y values. The argumentyfuncontains a pointer to a scalar function, which is used to calculate a Y value from the corresponding X value, and is assumed to be bounded and monotonically nonincreasing or nondecreasing. The argumentycontains the Y value to be bracketed. The argumentsyfarleftandyfarrightare assumed to contain the limits of (*yfun)(x) asxtends to minus infinity and to plus infinity, respectively. The argumentmbracketcontains the maximum number of rows that the bracket matrix is allowed to have at the end of execution. The argumentscalefactorcontains a scale factor, which must be nonmissing and strictly greater than 1, and is used to calculate X values for additional rows of the bracket matrix. It is assumed that the lowest X value in the bracket matrix is negative and that the highest X value in the bracket matrix is positive, so that an additional X value on the left (or right) can be calculated by multiplying the lowest (or highest) existing X value byscalefactor. The output argumentsleftindexandrightindexare set during execution to contain the highest index of the bracket matrix that bracketsyon the left and the lowest index of the bracket matrix that bracketsyon the right. This bracketing must be strict, except ify==yfarleft(in which case the left bracket may be partial) or ify==yfarleft(in which case the right bracket may be partial). The assumptions mentioned are assumed and are not necessarily tested. The terminology is detailed inDefinitions and methodsbelow.

_bcsf_bisect(objfun,x0,x1,y0,y1,itcount,iterate,tolerance,log)carries out bracket convergence by using a revised version (for step functions) of the bisection method, described in Press et al. (1992). The argumentobjfunis input and contains a pointer to a scalar object function, assumed to be a monotonically nonincreasing or nondecreasing step function. The argumentsx0andx1contain two X values on entry, which are assumed to be boundaries for an interval containing the supremum or infimum X value that we need to estimate, and will be revised inward during each iteration to contain boundaries for another smaller interval containing the same supremum or infimum. The argumentsy0andy1are assumed to contain, on entry, the values&objfun(x0)and&objfun(x1), respectively, and, if so, then they will still be equal to&objfun(x0)and&objfun(x1)on output. The argumentitcountcontains a running total of iterations and will be increased on exit by the number of iterations carried out bybcsf_bisect(). The argumentiteratecontains the maximum number of iterations thatbcsf_bisect()may carry out. The argumenttoleranceis input and specifies that no more iterations are to be carried out after the value ofreldif(x0,x1)is equal to or less than the value oftolerance. The argumentlogis input, and, if it is not equal to zero, thenbcsf_bisect()will output an iteration log, giving the values ofx0,x1,y0, andy1at the end of each iteration. The argumentsy0andy1must partially bracket a Y value of zero on entry if iterations are to proceed, withy0as the partial bracket andy1as the strict bracket. If this is the case on entry, then it will also be the case for the revised values ofy0andy1on exit. The terminology is detailed inDefinitions and methodsbelow.

_bcsf_regula(objfun,x0,x1,y0,y1,itcount,iterate,tolerance,log)carries out bracket convergence by using a revised version (for step functions) of the regula falsi (or false position) method, described in Press et al. (1992). The arguments have the same names and functions as those forbcsf_bisect().

_bcsf_ridders(objfun,x0,x1,y0,y1,itcount,iterate,tolerance,log)carries out bracket convergence by using a revised version (for step functions) of the Ridders method, described in Press et al. (1992). The arguments have the same names and functions as those forbcsf_bisect().

Definitions and methodsA real number value

yis said to be strictly bracketed by two other valuesy0andy1if and only if(

y0-y) * (y1-y) < 0or (in other words) if (

y0-y) and (y1-y) have opposite nonzero signs. If the Y values are defined using a step function, then we may use a more relaxed definition of bracketing. If we definesign0 = sign(

y0-y) sign1 = sign(y1-y)then we say that

yis partially bracketed byy0andy1if and only ifsign1!=0 & sign1!=sign0

This definition implies that

y0may be equal toybuty1may not. We then say thaty1is the strict bracket andy0is the partial bracket.Whether the bracketing is strict or partial, there may be a function

F()such thatF(x0)=y0andF(x1)=y1for somex0andx1. We then say that they-bracket associated with the lower X value is the left bracket and that they-bracket associated with the higher X value is the right bracket.The function

_bcsf_bracketing()finds a pair ofy-brackets for a Y value using a bracket matrix, which is a matrixbracketmatwith two columns and at least two rows, sorted in ascending order of the first column, with the property that, for each indexi,

bracketmat[i,2] = (*yfun)(bracketmat[i,1])In other words, the first column contains an ascending series of X values, and the second column contains the corresponding Y values.

_bcsf_bracketing()aims to identify a pair of indicesleftindexandrightindex, having the property thatbracketmat[leftindex,2] andbracketmat[rightindex,2] bracket the input valuey. This bracketing must be strict, except ify==yfarleft(in which case the left bracket may be partial) or ify==yfarright(in which case the right bracket may be partial). If the first row ofbracketmatdoes not contain a potential left bracket, thenbracketmatis extended by adding additional rows before the first row, until the first row of the extended matrix contains a potential left bracket. If the last row ofbracketmatdoes not contain a potential right bracket, thenbracketmatis extended by adding additional rows after the last row, until the last row of the extended matrix contains a potential right bracket. In both cases, an additional row is defined by multiplying the X value in the existing end row by the factorscalefactorto derive a new X valuexnew, and by defining the corresponding Y value as

ynew=(*yfun)(xnew)and defining the new row as (

xnew,ynew). The argumentscalefactoris usually set to 2, implying successive doubling of the outermost X value inbracketmat. This method requires that the first X value and the last X value inbracketmatare initialized to have opposite signs. Whenbracketmatis known to contain at least one potential left bracket and one potential right bracket,_bcsf_bracketing()returns the highest possible index for a left bracket inleftindexand returns the lowest possible index for a right bracket inrightindex.The functions

_bcsf_bisect(),_bcsf_regula(), and_bcsf_ridders()input a pointerobjfun, assumed to point to a nonmonotonically increasing or decreasing function, and aim to estimate the supremum (or infimum) of the set of X values with positive (or negative) values of (*objfun)(x). The functions start with a pair of X valuesx0andx1, with corresponding Y valuesy0=(*objfun)(x0) andy1=(*objfun)(x1), which partially bracket a Y value of zero, withy1as the strict bracket andy0as the partial bracket. The functions then carry out a sequence of iterations whose maximum number is given by the argumentiterate, in which one or the other of the (x,y) pairs is replaced with a new (x,y) pair, reducing the difference abs(x0-x1) but preserving the bracketing properties. The aim of the iterations is to reach a point where the relative difference reldif(x0,x1) is no more than the value of thetoleranceargument, and therefore either of the X values can be used as the estimate of the supremum or infimum that we aim to estimate. The methods used by the functions are versions of those described in Press et al. (1992), modified for use with step functions. The function_bcsf_bisect()uses a bisection method. The function_bcsf_regula()uses the bisection method if the value ofy0is zero at the start of the iteration and otherwise uses the regula falsi method. The function_bcsf_ridders()uses the bisection method if the value ofy0is zero at the start of the iteration and otherwise uses the method of Ridders (1979).

Return codes

_bcsf_bracketing()returns one of the following return codes:0: Bracketing successful > 1: Invalid row number for

bracketmaton input > 2: Missing values in input arguments or inbracketmaton input > 3: First and last X values inbracketmatdo not have opposite signs > 4:scalefactor<= 1 on input > 5: Left-extension ofbracketmatunsuccessful > 6: Right-extension ofbracketmatunsuccessful > 7: Left bracket index could not be located > 8: Right bracket index could not be located >

_bcsf_bisect(),_bcsf_regula(), and_bcsf_ridders()return one of the following return codes:0: Convergence successful > 1: Maximum iterations completed without convergence > 2: Missing X value calculated > 3: Missing Y value calculated > 4: Missing value for

iterateon input > 5:y0andy1do not partially bracket zero on input >

Conformability

_bcsf_bracketing(yfun,y,yfarleft,yfarright,mbracket,scalefactor,bracketmat,leftindex,rightindex):yfun: 1x1y: 1x1yfarleft: 1x1yfarright: 1x1mbracket: 1x1scalefactor: 1x1bracketmat:M x2 whereM>= 2leftindex: 1x1rightindex: 1x1

_bcsf_bisect(objfun,x0,x1,y0,y1,itcount,iterate,tolerance,log),_bcsf_regula(objfun,x0,x1,y0,y1,itcount,iterate,tolerance,log),_bcsf_ridders(objfun,x0,x1,y0,y1,itcount,iterate,tolerance,log):objfun: 1x1x0: 1x1x1: 1x1y0: 1x1y1: 1x1itcount: 1x1iterate: 1x1tolerance: 1x1log: 1x1

Diagnostics

_bcsf_bracketing(yfun,y,yfarleft,yfarright,mbracket,scalefactor,bracketmat,leftindex,rightindex)can fail only if the function indicated byyfunfails.

_bcsf_bisect(objfun,x0,x1,y0,y1,itcount,iterate,tolerance,log),_bcsf_regula(objfun,x0,x1,y0,y1,itcount,iterate,tolerance,log)and_bcsf_ridders(objfun,x0,x1,y0,y1,itcount,iterate,tolerance,log)can fail only if the function indicated byobjfunfails.

Source code_bcsf_bracketing.mata, _bcsf_bisect.mata, _bcsf_regula.mata, _bcsf_ridders.mata

AuthorRoger Newson, Imperial College London, UK. Email: r.newson@imperial.ac.uk

ReferencesPress, W. H., S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery. 1992.

Numerical Recipes in C: The Art of Scientific Computing. 2nd ed. Cambridge, UK: Cambridge University Press.Ridders, C. J. F. 1979. A new algorithm for computing a single root of a real continuous function.

IEEE Transactions on Circuits and Systems, vol. CAS-26(11): 979-980.

Also seeManual:

[M-0] introOnline:

mata,somersd(if installed)