help for solve

Find a solution to an equation of the form F(x) = y.

Syntax

solve retname
=exp
,
program(progname)
argname(argname)
arg_low(#)
arg_high(#)
tolerance(#)
[
maxiter(#)
hole1(#)
hole2(#)
hole3(#)
hole4(#)
nudgeval(#)
other_prog_options(options_to_prog)
verbose
]

Options

program(progname) Specify the name of a program that implements the function of interest.
This is required.

argname(argname) Specify the name of the targeted argument for progname;
it is also used as the name of the returned result from solve. This is required.

arg_low(#) One bound of the search interval. This is required.

arg_high(#) The other bound of the search interval. This is required.

tolerance(#) The maximal acceptable difference between the produced value and exp.
This is required.

maxiter(#) The maximal number of iterations allowed; default=60.

hole1(#) etc. Specify an argument for which the function has a "hole" (explained below).

nudgeval(#) The amount by which an argument should be incremented to avoid a hole; default=0.00001.

other_prog_options(options_to_prog) Other options to progname.

verbose Display the sequence of intermediary results.

Preliminary Remarks

Suppose F is a real-valued function of one real variable. Let y be a fixed value.
The goal is to find an x such that F(x) = y .

As a practical matter, the goal is to find x such that F(x) is within tolerance of y.
Also, it would be desirable for the solution to be unique, but let's put that aside for a while.

In the syntax,

exp corresponds to y;

argname corresponds to x.

progname corresponds to F;

retname corresponds to F(x);

It is necessary to create, or you may already have, a program that implements the function.
The name of the program is specified in program(progname).

progname must...

be declared as rclass;

include a syntax command, with real-typed option named argname
to serve as the argument of interest
(see optional_real_value in syntax);

return its resulting value in retname, preferably as a scalar.

Thus, the program progname must be established in advance of calling solve.
Its "input" is argname, and its "output" is retname. You inform
solve of these names via the argname option and retname feature.

It is possible for progname and retname to be the same. You may find this
convenient - or confusing.

The returned value is best produced as a scalar. Thus in progname, there should be
the command

return scalar retname = some_expression

(return macro is a possibility, but you may lose precision.)

You also need to specify a targeted interval
in which to search for the solution; this interval is bounded by
arg_low(#) and arg_high(#). Usually, it is necessary that F(arg_low) and
F(arg_high) occur on opposite sides of exp. That is, they straddle exp.

If progname requires that you specify other options in addition to
argname, place them in the other_prog_options option. Thus, you will
specify options inside an option specification. See the examples.

Examples

. prog def square, rclass

. syntax, x(real)

. return scalar sq = `x' * `x'

. end

. solve sq= 2, prog(square) argname(x) arg_low(0) arg_high(3) tol(.0001)

[output omitted]

solve result: x = 1.4142151

. prog def raise, rclass

. syntax, x(real) p(real)

. return scalar r = `x' ^`p'

. end

. solve r = 2, prog(raise) argname(x) arg_low(0) arg_high(3) tol(.0001) other_prog_options(p(2))

[output omitted]

solve result: x = 1.4142151

. solve r = 2, prog(raise) argname(x) arg_low(0) arg_high(3) tol(.0001) other_prog_options(p(3))

[output omitted]

solve result: x = 1.259903

. ret list x

scalar r(x) = 1.259902954101563

Of course, there are more effective ways to obtain a square or cube root, but this
demonstrates the use of solve. Note that the prior result is not particularly
precise. We can redo it to greater precision:

. solve r = 2, prog(raise) argname(x) arg_low(0) arg_high(3) tol(.00000001) other_prog_options(p(3))

[output omitted]

solve result: x = 1.259921

. disp 2^(1/3)

1.259921

Following is the original usage that motivated the creation of solve.
The program calc_vsl is not shown, as it is lengthy. but suffice it to say that it has
many options in its syntax, including sigma.
It does some of its work in data, and calls other programs.
The vsl option refers to a variable name; the returned value is also named
vsl (it is the sum of the vsl variable).

. #delimit ;

. solve vsl = 2757032.2, prog(calc_vsl) argname(sigma)

arg_low(.5) arg_high(6.5)

tolerance(1e-8) maxiter(150)

hole1(1)

other_prog_options(

age(age) z0_z_ratio(z0_z_ratio) y_f_hat(y_f_hat) c_f_hat(c_f_hat)

discount_factor(`discount_factor')

ref_age_index(`ref_age_index') max_age_index(`max_age_index')

survprob(`survprob_modified') phi(`phi') integrandA(`integrandA')

vsl(`vsl')

)

;

[output omitted]

solve result: sigma = 1.3616268

Note the presence of the hole1(1) option. calc_vsl has a hole at
sigma=1; it yields a missing value. (In an earlier version, it yielded 0, which was
out-of-line with its neighboring values.) See below for an explanation of holes.

Additional Remarks

In a more general case, F may have a multitude of arguments, that is, it has the form

F(x, other_arguments)

but your interest is focused on only x; you are concerned only with how F behaves as x varies,
while other_arguments are held at some particular fixed value (or tuple of values).
Thus, F can be viewed as a function of x alone, while
other_arguments are considered as parameters that adjust the behavior of F.

This would correspond to progname having additional options (options_to_prog), though there are
other avenues for influencing the result, as will be discussed shortly.

Similarly, progname may return other values in addition to retname, but
for the present purpose, you would be interested in only retname.

As noted, there are many potential factors that could influence the resulting value of
retname. These include...

argname,

options_to_prog (if progname has such options),

data,

scalars,

global macros,

system parameters (see creturn),

characteristics and notes,

any similar influencing factors to programs that are called by progname.

solve will issue multiple invocations of progname while varying argname
as it seeks a solution.
It should be the case (it is the user's responsibility) that all other influencing factors
(other than argname) remain constant as solve makes these invocations of progname.
This way, for the duration of running solve, progname can be viewed as a
function of argname only; it should yield consistent results in repeat call with the same
value of argname (and options_to_prog if applicable).

(options_to_prog may contain names of variables or other named entities.
So, while its textual value is fixed for the duration of running solve, the
values in the referent entities could possibly vary; it is important to not have such
variation.)

progname may alter data or other entities in the process of calculating its result,
but the data elements that it relies on should not be altered. Thus, it might
generate or replace some data variables in the process of calculating a result,
but such variables should not be referenced as a source of information in the calculation
of retname. For example, say that a, b, and c are variables in the data. You then
may generate d as some calculation that depends on a, b, c, and argname, followed by
summing d, and returning the value of that sum in retname. Note that, while you
have made changes to the data, you have not changed a, b, or c; the changes are
only in d, which has no effect on a subsequent call. (You might drop d, or
it may have been a tempvar.)

Usually, progname should be such that it does not alter the set of observations;
it should not drop or create new observations, or collapse the data.
The exception would be if it reloads the (same) data in each invocation.

exp is evaluated once, prior to seeking a solution. Thus, it is not a moving target.
Typically it would be a constant, but generally may be an expression. As such, it
would be evaluated as a scalar quantity, meaning that any variable references that
are not specifically indexed default to index 1.

maxiter(#)} specifies the maximum number of iterations that will be performed.
It is not the same as the maxiter system parameter c(maxiter), which is used in
maximum likelihood estimation procedures. (See creturn and maximize.)

Note that the result is not necessarily exact;
generally it is "close enough", where tolerance sets the standard of what
that means. You may specify a tolerance of 0, signifying that you desire an exact
solution, but you are likely to not achieve such a solution in general.

The result is returned in r(argname). Thus, argname
serves two purposes: the name of the argument to progname, and the returned
result from solve.

It is important that the function be continuous within the specified range (except for possible holes, see below).
It is also somewhat important that it be monotonic. If it is not continuous,
then solve may fail to yield a result, even if a solution does exist.
If is is not monotonic, then solve may find one of possibly many solutions.

solve will attempt to detect and issue warnings about non-monotonicity.
However, it cannot detect all instances of non-monotonicity.

The function may have "holes". A hole is an argument that yields a value that is
out-of-line with its neighbors, or yields a missing value. Thus, a hole is a point
of benign discontinuity. The function is "almost continuous", in that its
right and left limits are equal, but not equal to its value at that argument.
A classic example is F(x) = x^2/x, which has a hole at x=0 (in Stata terms, it yields missing),
but otherwise, is the same as F(x) = x, which is continuous.

Another possible type of hole is one where the function
has a non-missing value that jumps out of line at a particular point. This may seem hypothetical,
but it happened in the original application that motivated the creation of solve, and
thus motivated the inclusion of the hole1() (etc.) options.

If one of the candidate arguments should land on a designated hole, it will be
altered a slight bit to avoid the hole; it is moved out of the way by
incrementing it by nudgeval.

Note that, if the value at a hole is non-missing and equals
exp, then the given argument is a solution, but will not be detected.

It is the user's responsibility to determine the location of holes, if they exist.

The syntax allows up to four holes to be specified. Contact the author if more are needed.

Of course, a jump discontinuity may prevent a solution from being found, even if
one exists.

arg_low(#) and arg_high(#) are the bounds of the search interval.
They may be regarded as initial estimates for
the solution, where the solution is believed to lie between these two arguments.
It is usually necessary that a solution exists between these arguments; the exeption
would be if one of them qualifies as a solution.

arg_low(#) and arg_high(#) may occur in either order; it is not necessary that
arg_low(#) < arg_high(#).

Note that the resulting soultion may be sensitive to the choice of these arguments, as
well as to tolerance.

solve does not allow for progname to have varlist, namelist, anything, if, in, using,
=exp, or weights. Contact the author if these limitations present a problem.

Methodology

solve uses the bisection method, which is just one of several established methods
for numerically solving such equations.
It starts by evaluating progname at two arguments, call them A and B,
which may be regarded as preliminary guesses to the solution.
Initially, A and B are arg_low and arg_high.

Let...

valA be the result of evaluating progname at argname=A;

valB be the result of evaluating progname at argname=B.

If valA is within tolerance of exp, then A is the solution;

If valB is within tolerance of exp, then B is the solution.

If both valA and valB are within tolerance of exp, then choose the one
that is closer to exp, giving preference to valA, and a solution of A, if they are equally close.

Otherwise, bisect the interval, evaluate at the midpoint, and turn your
attention to the half-interval that straddles exp. Repeat the process on this
new, smaller interval; repeat until a solution is found.

Also See

Bisection is regarded as the least sophisticated among several methods for solving
such equations. Other methods include

Newton-Raphson method;

Brents's method.

See the Wikipedia articles on those methods.

Also, see mm_finvert in the Moremata module, available on ssc.
This is intended for functions defined in Mata. Type "ssc desc moremata".

Thus, mm_finvert is useful for functions defined in Mata, while
solve is useful for functions defined in Stata. But it is presumably possible to
use mm_finvert for Stata functions, or solve for Mata functions, though this may
involve some convoluted programming techniques.

Author

David Kantor; email kantor.d@att.net