A Computer program to minimize a function with many variables using computer evaluated exact higher-order derivatives

A Computer Program to Minimize a Function with Many Variables Using Computer Evaluated Exact Higher-Order Derivatives R. Kalaba and A. Tishler* Department of Economics University of Southern California Los Angeles, Califmia 90089-0035

ABSTRACT

This paper describes a computer program for an optimization algorithm using firstup to the rth-order derivatives to find the optimum of rcontinuouslydifferentiable functions of many variables. The algorithm, developed by Kalaba and Tishler [6, 71, obtains the exact values of the derivatives required for the optimization from the “table algorithm” presented in [4] and [5]. The optimization algorithm described here collapses to the well-known Newton-Raphson algorithm when only first- and secondorder derivatives are used.

1.

INTRODUCTION

This paper describes a computer program of an optimization algorithm using first- up to the &-order derivatives to find the optimum of r-continuously-differentiable functions of many variables. The algorithm, developed by Kalaba and Tishler [6, 71, obtains the exact values of the derivatives required for the optimization from the “table algorithm” presented in [4] and [5]. The optimization algorithm described here collapses to the well-known NewtonRaphson algorithm when only first- and second-rder derivatives are used. Thus, we call the algorithm a “generalized Newton algorithm of order r,” GN( r ), when the highest derivatives to be used in the optimization process are of order r + 1. The integration of the GN(?') and the table algorithms into a combined algorithm, denoted AUTOGN(T), is derived in [7]. *On leave from the Faculty of Management,Tel Aviv University, Israel APPLIED MATHEMATICS AND COMPUTATION 13:143-172 0 Elsevier Science Publishing Co., Inc., 1983 52 Vanderbilt Ave., New York, NY 10017

(1983) 143 00963003/83/$03.00

R. KALABA AND A. TISHLER

144

The computer program presented here can calculate the exact values of derivatives up to order four; hence, the AUTOGN(T) algorithm uses derivatives up to order two, three, or four, depending on the researcher’s choice. Generalization of AUTOGN(T) to use derivatives of order five or higher is possible along the same lines discussed in this paper. This article proceeds as follows: Section 2 presents the GN(~) optimization algorithm, and Section 3 describes its programming. Sections 4 and 5 detail the programming issues of the table algorithm using the valuation formula of Black and Scholes [2] as an example. Section 6 is devoted to further discussion and summary. The complete computer program is listed in the appendix. 2.

THE GN(T) ALGORITHM:

A REVIEW

The optimization problem is given by min F(x), XEP

(2.1)

where F(x) is r-continuously-differentiable. respect to xi be denoted by yi; that is,

Let the derivative of F(x)

-aF(x) ax = (Yl,Y2,...'YK)"

with

(2.2)

Let the “super-Hessian” matrix A be given by A= [(l);(y,)>...,

(YK);(Y12),(YlY2)...~(y~);...;

(y;>,(Y;-‘Y,),...,(Y~)l~

(2.31

where the expressions in the parentheses are the following column vectors:

(z)=

Iz;g )...)

Yg$.&Jg

)...)

1

1 I

a22

ax;

-;...,

.- a*2 a*;>

TZ

iTi2 .&_‘.&,

‘..‘)

ax;,

.

(2.4)

A GN(~) algorithm for finding the optimal r in problem (2.1) is defined as

A Computer Program to Minimize a Function

145

follows: a.0 a.i

Choose x0 and form A0 (A at x = x0) according to Equation (2.3). i=1,2,3 ,... a.i.1 Calculate xi = xi-’ + a’-‘, where S’-’ includes the second up to the K + 1 elements in the first row of A-’ evaluated at xi-‘. Co to a.i.2. a.i.2 If [(xi - xi-l]] < E for some tolerance E > 0, go to a.i.3.; otherwise, compute A’ using Equation (2.3) at x = xi and go to a.i.1. a.i.3 xi is a local minimum for F(X) if 8F(x’)/& = 0 and the Hessian of F is positive definite.

3.

DESCRIPTION

OF THE GN(r) ALGORITHM

3.1

General Remarks As is evident from the algorithmic scheme in Section 2, the computation of the matrix A [see (2.3)] is the major step in each iteration of the GN(T) algorithm. Once A is given, the next approximation for the vector x is given by step a.i.1. The construction of A for any r > 1 is given in [7]. In this paper we restrict ourselves to r = 2,3,4 according to the program available in the appendix. Note that for any GN(T), A is a square matrix with number of columns equal to one plus the number of distinct first- up to r + lst-order derivatives of F with respect to X. We first describe the elements in the first row of A. For GN(~), the first element of the first row of A is equal to one. Positions 2 up to K + 1 are occupied by the K first-order derivatives of F [in (2.1)] with respect to xi, x s,“.,xK; that is, y,, yz,. . . , yK. If c~(2) is in use, the first K + 1 elements in the first row of A are identical to those used for GN(1). In addition, positions K + 2 up to 1 + K + K( K + 1)/2 include all distinct products of pairs of yi, yj, i, j= 1,. . . , K, in increasing order of i and j where i < j; that is, y;, y, y,, . . . , ylyK; y$, ysy3,. . . , y,y,;...;yi_,,yK_,yK;yi. For~~(3)positions 1 up to l+K+K(K+1)/2 in the first row of A are identical to those used in GN(2). Positions 2 + K + K(K + 1)/2 up to the last one (which is equal to 1 + K + K(K + 1)/2 + Cf= 1 ZyciIF_ j) include all distinct products of triples yi yjyl,in increasing orderofi,jl=l,..., K,wherei
146

R. KALABA AND A. TISHLER

element in the first position (row) of that column. The order of the higher-order derivatives in each column in A is identical to that in Equation (2.4). Table 1 presents the matrix A for GN(~). Its (K + l)X(K + 1) upper left block gives the matrix A for GN(~), and its [l + K + K(K + 1)/2] X [l + K + K(K + 1)/2] upper left block gives the matrix A for GN(~). To facilitate further discussion we use the following definition.

DEFINITION 3.1. The vector z is the vector form of the function Z(X), x = (x,, x2,. . . , x,)‘, if it includes in its first position the value of z(x), and then the first- up to the rth-order (distinct) derivatives of Z(X) with respect to Xl> x B’...>XK, in increasing order of indexing as is given in (2.4). At this stage we note that applying the table algorithm (see Section 4) to any function of x, say g(x), results in a vector containing the function g(r), the first-order derivatives Jg(x)/&,,.. .,ag(x)/&,, and the second- up to rth-order derivatives, all in precisely the same order which is described in forming the columns of A. Thus, we apply the table algorithm to each of the components of the first row in A to compute the full matrix A. Clearly, each column of A is the vector form of the element in its first position (row). 3.2.

Description

of the Subroutines

Which

Perform the Optimization

a. Subroutine GNR. In iteration i, subroutine GNR (see the program listed in the appendix) gets as inputs the initial conditions xi- ’ (in the vector x), the function F and all its derivatives (w.r.t. x) up to order r + 1 (in the vector F) and an index stating what is the value of r + 1 (denoted I R in our program). Subroutine GNR prints xi- ‘, the value of F( xi _ ’ ) and its gradient (the values of yr(x’-‘),...,yK(x’+‘)). Th en it calls subroutine OGEN. Subroutine OGEN gets the vector F and the value of r as inputs, and returns the vectors (yK) [see (2.3) and (2.4)] in the columns of the matrix denoted (Y,),(Y,),..*, Y. The matrix A is then constructed according to the value of IR (which could be 2, 3, or 4 in our program) and inverted using the IMSL subroutine L INV2 F. The vector xi is calculated as described in step a.i.1 in the algorithmic scheme in Section 2. Finally, subroutine GNR prints ri and the vector ri - x1-l. The output of GNR is the vector ri. b. Subroutine OGEN. The inputs for of r and the vector F containing the value rth order dist #catderivatives with respect in increasing order of indexing; that is,

subroutine OGEN include the value of the function F and its first up to to r. The derivatives are organized if we minimize a function of two

axf

~a2Y1

ax;

‘For [l +

ax;

ax;

ax,

a3(d)

ax;

.'.

ax;

ax;

a3(YB) a3(dY2)

ax;

ax;

ax;

.,. a3&)

ax;

x

aYyi?) axfax

aYy3 ax;

a$

axlax

az(d >

K + K(K + 1)/2]

ax;

a$l") a3(Y,2Y2)

...

a"(yl") a3(y12y2). ax; ax;

ax;

A is the ( K + 1) X (K + 1) upper left block. For GN(~), A is the [ 1 + GN(~), A is the fd matrix.

ax;

ax;

~a(Yif) ax,

3

a2(Y;) a$?) a2(y?y2)... az(Y;) 37 Yi >

... ~a3YK a"(Yl">aG,Y,)

.

...

a(yi$2) .

2

ax: ax; ax; a2(y,2y2) aYy;) az(yl") ... axlax ax,ax, axlax

ax;

a%hY,)

ax,

a:;)

3

axfax axfax, axfax '.. axtax axfax, axfax

a% a%3

T

. . a%& a2(d) a2(Y12Y2) . . . a%3

ax,

a%hY2)

ax;

K + K(K + 1)/2] upper left block. For

GN( l),

ax;

axfax

- a%,

. . . 0 -@Y,

. . .

0

ax;

a2(YlY2)

_a% .. a% aY;l”)aYYlY2) ...

. . .

ax;

2

GN( r)”

... a(&

.. .

a2yK a"(~:)aYy,y,) . . . ax,ax2 ax,ax, ax,ax2

a"(d)

a~;

ax,

a2YK

a(Y;Y,)

ax,

8x1

qY,Yz)

YlY2

. .. ~4Yti‘i) ax,

TABLE 1 A FOR

MATRIX

&l")

ax,

~(Y!h

o _a3y, ... a%, a%3 ax; ax; ax;

o

. . .

ax,ax,

..

.,.

..’

YK 8Y1 aYK .. ax, 8x1 . . * aYK ax, ... ax,

YI

o a%

0

0%

0 . . .

1 z

THE

148

R. KALABA AND A. TISHLER

variables, say, the vector F contains (for r = 3)

where F. denotes the value of the function and other subscripts denote the relevant derivatives. The output of subroutine OGEN consists of the variables yi,. . . ,yK in vector form as defined by (2.3)-(2.4). Given the vector F, subroutine OGEN constructs yr, . . . ,yK by a search over the indices of the elements of F. That is, to construct yr, start at the first component of F (which is F,), and proceed to the second, the third, and up to the last, where every time an index (an element of the subscript of a component of the vector F) equals 1 (in any position in the subscript), this component is added to the vector y,, in the order of the search described above. For ya, start at F,, repeat the search, and every time a subscript containing the index 2 arises, add the component of F with that subscript-in that order-to the vector ya. The same procedure is applied to y,, . . . , yK. c. Subroutine MULTGN. When IR > 2, GNR calls subroutine MULTGN to generate the products yiyj (and also yiyjyl when I R = 4) in vector fom for the matrix A (see Section 2). Subroutine MULTGN is similar to MULT (see description in Section 5), but is restricted to evaluate yiyj (and yiyjyl when I R = 4) and all its first-, second-, and third-order (if I R = 4) derivatives, as is required for the matrix A. The precise description of the computation of the vector forms of yiyj and yiyjyl is given in Section 4. d. The MAIN. The MAIN of the program (see appendix) sets up the number of variables in the function F (denoted NVAR), the value of r + 1 (denoted I R), and the order of derivatives to be calculated (denoted by I D, where I D> I R). Next, it calculates the number of derivatives of second, third, and fourth (if required) order to be calculated (denoted NSEC, NTHR, and NFOR respectively). Initial conditions for x = x0 are given (the vector x in the program). Starting at x = x0, the program successively updates the vector X, by first calling subroutine FUN CT (to be described in Section 5) which produces the vector F given x (using the table algorithm) and then calling subroutine GNR which produces the next approximation to the vector x. This loop is repeated until a maximal number of iterations is reached (denoted by NIT). Clearly, we can introduce, after subroutine GNR, a subroutine to check whether convergence is achieved. Since many convergence criteria are possible, we don’t present any of them in this paper. However, a simple possible approach is to check whether I/x- xi- ‘11~ E, where i is the iteration number and E > 0 is some tolerance. Note that in the above process we call the subroutine SEDER before calling FUNCT in each

149

A ComputerProgram to Minimize a Function

iteration. This is required in order to transfer the values of the vector x, at iteration i, to the table algorithm in the computation of the function F (see Section 5).

4.

THE “TABLE

ALGORITHM”-BRIEF

DESCRIPTION

The precise description and properties of the table algorithm are given in [5] and [4]. In this paper we shall describe it only briefly. In principle, a table algorithm can be constructed to calculate the value and partial derivatives through order r of any admissible (r, n) function defined as follows:

A function F: L) + R, D c R”, n > 1, will be referred to DEFINITION 4.1. as an admissible (r, n) fin&on if, given any x = (xi,. . . ,xn) in D, the value F(x) can be sequentially calculated in a finite number of steps by means of the n initial conditions sl=xl,...,

8,

(4.0

=x,,

the twovariable algebraic special functions f: R2 + R given by

w=u+u,

w=u-u,

and arbitrary, one-variable, &-order tions of the form h: M+R,

w = uu,

w =

u/v,

w =

uv,

(4.2)

continuously differentiable special func-

McR.

(4.3)

Given any admissible (r, n) function F with domain point x, the general table algorithm for F(x) is constructed as follows: First, form the admissible Zist for F(x), i.e., a-list of initial conditions (4.1) and special functions (4.2) and (4.3) whose sequential evaluation yields F(x). Second, replace (using a calculus subroutine) the initial conditions si = xl in the admissible list for F(x) by a one-dimensional array (vector) containing the value and partial derivatives through order r of the function g : R” + R defined by g(x, ,..., x,)= xl. Speciry analogous arrays (vectors) for the functionsg(x, ,..., rn)Eri,i=2 ,..., n, to replace the remaining initial conditions si = xi, i = 2,. . . , n, in the admissible list for F(x).

150

R. KALABA AND A. TISHLXR

Third, replace each special function (4.2) and (4.3) in the admissible list for F(x) by a calculus subroutine for outputting the value of the special function, together with all of its partial derivatives with respect to x through order r. For any two-variable special function w = f(u, u), the input to the calculus subroutine will be the previously calculated one-dimensional arrays U = (u, u,), u,~, . . .) and V = (0, ux,, ur2,. . . ), and the output will be a one-dimensional array W = (w, wz,, wX,, . . . ). For any one-variable special function z = h( y ), the input to the calculus subroutine will be the previously calculated one-dimensional array Y = (y, yx,, yx2,. . .), and the output will be a onedimensional array Z = (z, z,,, z,~, . . . ). The list of calculus subroutines resulting from these three steps is called the table algorithm for F at x. To simplify the exposition we use the following example:

F(x)

= eax +

Px2 + y,

(4.4)

where o = - 1.0, p = 0.5, and y = 10.0 are known constants. F(x) is computed sequentially in Table 2. All right-hand terms in the expressions in (4.5)-(4.13) involve operations on X, on a known constant, or on previously [in the sequence (4.5)-(4.13)] calculated variables for any given values of x. The sequential evaluation of steps (4.5) to (4.13) yields the corresponding value F for the function defined by (4.4) for initial values (conditions) of x = 1~‘. Moreover, the output at each step in the sequence (4.5)-(4.13) includes the value of the function computed in the second column of Table 2, and all its derivatives up to order 4 in columns 3-6. Thus, each operation on a variable z in Table 2 is done in a vector form (see Definition 3.1). The important feature in the computation of the derivatives in Table 2 lies in the simplicity of the computation in each step and in the sequential calculation of its entries. For example, to calculate Q, Q,, Q,,, Q,,,, and Q xxxxin step (4.11) one needs to have D, D,, OX,, DXX,, D,._, P, I’,, Pz,, PXz,, and LX7 all of which were calculated in the previous steps (4.8) and (4.10). Note that in each of the steps (4.5) to (4.13) we actually compute the five component vectors Z = (2, zX, zXX, zxxx, .z,,,, ) of the relevant left-hand-side variable, which is defined by a very simple operation on the right-hand side of the equality sign. Since the number of these operations is very small (as we shall show below) and most functions F use identical types of operations, the software given in the appendix can be applied (possibly with some modifications) to any function F.

P = /?H Q=P+D S = yB

(4.10) (4.11) (4.12)

(4.13)

Q+ S

H=A.A

(4.9)

F=

A = x0 B = 1.0 C= aA D=&

(4.5)

(4.6) (4.7) (4.8)

Function

Step

a

= YB,,

L

F,, = Qrr + L

Sx = YB,

Fx= Qx + S,

Pm = BH,, Q_=Px,+Dx,

Px.xX= BH,,, tzx = Pm + D,,, xxx = YB,,, Fz+x= Q,,x + L

Cm = aA,,, Dx,, = D&z +2DxCm + DGx* H,,, = A,,,A +3A,,A, + 3ArAm + AA,,,

C,, = aA,, Dm = D&x + DC,, H,, = A,,A +2A,A, + A&

A xx+= 0 B xxx = 0

a3 ax3

A,,=0 B,, = 0

ax2

a2

P, = BH, Qx=Pz+Dx

+AAz

H,=A,A

A,=1 B, = 0 C, = aA, D,=DC,

z

SEQUENTIAL

TABLE 2 COMPUTATION OF F, F,, F,,, F,,,, ANDF,,,,

= $:: = xxx* = Fzxxr=

Hz,,, =

A IXIX= Bx**x = CXXXI= D zxxx =

Q-

+ Lx,

aA.... D&x + 3D&, + 3D,c,,, + DL, A,,,,A +4A,,,A, + 3A,,A,, + 3A,,A,, +4A,Ax,, + AA,,,, BH,,,, Pm,, + DYB-

0 0

a4 ax4

s 3 s

k

t.

ii 2.

3

;

il.

R. KALABA AND A. TISHLEZR

152

The operations in Table 2 (or for that matter, for any other function F) can be classified to one of the following three types: (a) defining a variable (“vectorizing”), e.g., A = x0, (b) operation on one variable, e.g., D = ec, (c) operation on two variables, e.g., Q = P + D. It is obvious that all relevant functions can be sequentially evaluated using only these three types of operations. One could, however, introduce special functions of three or more variables (see [S]). Clearly, the order of the sequence need not be unique. For example, (4.8) and (4.9) can change positions in the order of evaluation, without affecting the outcome of F in (4.13). Thus, one can choose, arbitrarily, any possible ordering for the sequence of calculations, provided it yields the correct value of F in the last entry of Table 2. To implement the computation one needs to supply routines which calculate the five-component vectors 2 = (z, zx, zrx, z_, z,_), where 2 stands for each of the left-hand-side variables in the second column of Table 2. It is easy to see that this amounts to a fairly small number of routines. More specifically, the required routines for each type of operations are as follows. (a): A “vectorizing” routine is needed for each variable. In our example a vectorizing routine for step (4.5) produces the vector

In general, we order the variables ri, x2,. , . ,xK [see (2.1)] in some arbitrary, but fixed, order. “ Vectorizing” a variable is executed by specifying CALL VECTtJ,

XJ,

YJ)

where J is the index of the Jth variable (of rj), XJ is the value of xj (initial condition for r j), and YJ is the output vector from VECT given by YJ=(XJ,O

,..., O,l,O ,..., 0)‘. Jy1

We recommend, but it is not necessary, that the first K operations in the computation of F(x), where x is a K X 1 vector be to “vectorize” xi, i=l ,.**, K. (b): Each operation of one variable requires a special routine. Three special routines are needed for Table 2. That is, in steps (4.7), (4.10), and (4.12) we

A Computer Program to Minimize a Function

153

need a routine which multiplies a variable by a constant. This routine is denoted MULCON in the appendix and calculates the first-, second-, third-, and fourthorder (if required) derivatives of arA, /3H, and yB respectively. Actual execution is done by specifying [for step (4.10) for example] CALL

MULCON(j3,

H, P).

where the inputs are B, a known constant, and H, a variable (in vector form). The output is given in the vector I’, as is described in Table 2. The step (4.6) is executed by CALL

UNIT

(B).

The output is B = (l,O,O,O,O)‘. Step (4.8) requires a routine which produces the value and all relevant derivatives of ec, where C is a function of the vector x. This is done by specifying CALL

EXPP( C, D),

where the input vector is C [which was computed in (4.7)], and the output vector is D. The construction of D in subroutine EXPP is similar to that in Table 2, step (4.8). Other one-variable operations in our program (which are not used in the above example) are executed as follows:

(i)

CALL

NULL(A),

which produces the null vector in A.

6)

CALL

RECP(Z,Y),

where Z is an input variable (in vector form). The output is l/Z form.

(iii)

CALL

in vector

LOGG( U, RI,

where U is the input variable in vector form, and the output (stored in R) is log(U) in vector form.

(iv)

CALL

CONTTU(b,U,

Z),

R. KALABA AND A. TISHLEX

154

where the inputs are b, a known constant, and U, a variable (in vector form). The output is Z = bU in vector form. CALL

(v> where /“,

Z

is

the

input

( 1/fi)e-it2dt

NOMC(Z,Y),

variable

(in

vector

form).

The

output

is

in vector form.

Other one-variable operations like sin(x), tanh(x), etc., are easily available using the same structure of subroutines like RECP, EXPP, LOGG, etc. For almost all functions of one variable, we use the chain rule in the computation of the first, second, third, and fourth derivatives. Assume that the function is given by (4.14)

U=h(w)

where w is a function of two variables rr and x2. Then (for simplicity in exposition let r = 3) (4.15)

U= h(w), U, = h’(w)w,,

(4.16)

U, = h’( w)w,,

(4.17)

Urr= h”(w)wf+

(4.18)

h’(w)w,,,

U,, = h”(w)wlwz + h’(w)w12,

(4.19)

U,=

(4.20)

h”(w)w;+

h’(w)w,,,

(4.21)

Un1=h”‘(w)w;+3h”(w)wrwrr+h’(w)wlll, u,,,=

h”‘(w)wzwf+

h”(w)(2wlzwl+

wZql)+h’(w)wm, (4.22)

L&,=h”‘(w)wrw;+

h”(w)(2wlzwz+

~~w~)+h’(w)w~~~, (4-B)

u,,

= h “‘(w)w;

+3h”(w)w,w,

+ h’(w)w,,,.

(4.24)

A Computer Progam to Minimize a Function

155

Thus, it is sufficient, for any h, to compute once h(w), h’(w), h”(w), and h”‘(w) and then form (4.15)-(4.24) in the same way for all functions h. Specifically, we form a subroutine (DER in the appendix) that gets as input h(w), h’(w), h”(w), h “‘(w), and the variable w (in a vector form), which was computed in a previous step, to form the vector (U, Vi, V,, Uii,. . . , U,,) in (4.15)-(4.24). Extending the above procedure to evaluate fourth-order (or higher) derivatives of U [in (4.14)] is straightforward. (c): There are only five possible functions of two variables in the table algorithm: addition, subtraction, multiplication, division, and raising a variable to the power of another variable. All these operations are done in vector form. In Table 2 we use multiplication in (4.9) and addition in (4.11) and (4.13). Operations of functions of two variables are executed as follows (letting A and B be the input vectors and C be the output vector): CALL

6)

MULTtA,

B, Cl,

where C= A.B; (ii)

CALL

ADD(A,

B,C),

CALL

SUB(A,

B, Cl,

CALL

DIV(A, B,C),

where C = A + B; (iii) whereC=A-B;

(4 where C = A /B;

CALL

(4

UTTVtA,

B, Cl,

where C = AB. Finally, for a given function F, the sequence of CALL statements which sequentially calculates F (in vector form), as is exemplified in Table 2, is organized in subroutine FUNCT. 5.

AN APPLICATION

Examples of the use of GN(T) algorithms when r = 1 and 2 are found in [5, 61. Here we apply the AIJTOGN(T) algorithm for T = 1, 2, 3 to the

R. KALABA AND A. TISHLER

156

Black-&holes

[2] formula describing the value of an option; that is, let (5.1)

where

d

h(Yi/ci)+(R+ia2)7i

,=

II

=

d 2r

(5.2)

uri

h(yi/ci)+(R-!fiU2)7i

ari

>

(5.3)

and

Equation (5.1) gives the value of the option at time i, where yi and ci are the stock price and the execution price of the option at that time. We let +ribe the number of periods remaining (at time i) until maturity (execution) of the option. The short-term interest rate, denoted R, is assumed to be constant through time, as is u2 (and hence a), the variance rate (“volatility”) of the return on the stock. Thus, given observations on yi, ci, ri and the actual observed (in the market) value of the option (denoted si), for i = 1,. . . , N, we shah estimate R, u, and (Ywhich minimize

where (~--a constant-stands for the transaction cost of buying an option (see [2, p. 6531). The subroutine FUNCT1 in Section A.3 of the appendix was used in the AUTOCN(T) to optimize (5.5). We use synthetic data (15 sample points) as follows: yi = ci = 1, ri = 0.1+ (i - 1)X0.01, and Gi are generated according to (5.1). The “true” parameters are R = 0.1, u = 0.5, and (Y= 0.3. Estimation results for r = 2,3,4 are presented in Table 3. The estimation of (5.5) is quite fast for r = 1,2,3. GN(~) and GN(~) are almost identical in terms of number of iterations, while GN(~) is a close third.

0 1 2 3 4 5

0 1 2 3 4

2

3

0.07 0.06688 0.08706 0.09932 0.10000

0.07 0.06823 0.07725 0.09963 0.09997 0.10000

0.07 0.06900 0.07066 0.10316 0.10222 0.10002 O.lOoOO

R

0.54772 0.53902 0.51307 0.50071 0sOOOu

0.54772 0.54209 0.52258 0.50047 0.50003 0.50000

0.54772 0.54455 0.52888 0.49755 0.49776 0.49998 0sOOoO

u

- 0.47 0.76x10-* o.25x1o-4 0.43x10-5 0.14xlo~s

-0.94 o.80x10~4 0.48x10-4 0.82~10~~ 0.28x10-s

- 5.84 -0.93xlo~s 0.29x10p3 o.49x10p4 0.16~10~’

-5.84 0.40~10-~ o.11x10~3 0.17x10-3 0.72~10~’ o.22x10~‘2

- 5.84 -0.19x10-” 0.70~10-~ o.13x10~2 0.70x10-6 0.67~10~~ o.55x1o-g

F,

0.57 0.29x10-6 0.52x10-s 0.55 x 10 - lo 0.15x 10 - l6

0.57 0.76~10-~ o.11x1o-7 053x10~9 0.21~10~‘~ 0.88 x 10 - 27

0.57 0.13x10-5 0.18~10~’ 0.31x10-’ 0.11x10-~ 0.84X10~‘2 0.22x10-”

F

at the University of Southern California. Iteration 0 gives the initial

0.1 0.29641 0.29897 0.29994 03OoOO

- 0.94 0.14x10-3 o.19x10~4 0.28x10-4 0.12x10-7 o.36x10~‘3

F,

- 0.47 0.1 0.13x10-3 0.29581 0.29823 O.89x1O-5 0.29996 0.15x10-4 0.3oOOO 0.54x10-s 0.30000 o.19x10-‘3

FR -0.94 0.18~10-~ o.llx1o-5 o.22x1o-3 0.83x10 -’ o.llxlo-5 0.89 x 10 - lo

a

- 0.47 0.1 0.17x10-3 0.29536 0.29774 -o.58x10~6 0.12x10-3 0.30018 0.13xlo~s 0.30018 0.59x10-6 0.30000 0.59 x 10 - lo 0.3OoOO

“Calculations done on IBM 370/H% conditions.

0 1 2 3 4 5 6

Iteration No.

1

r

TABLE 3 ESTIMATION BESUL-~~%OR BLACKANDSCHOLES’S FORMULA USINGAUTOCN(r), r = 1,2,3

R. KALABA AND A. TISHLER

158

The CPU time required for GN(~) is quite large-24 seconds, as compared to about 2 seconds for GN(~) and about 0.5 seconds for GN(~). Thus, it seems that the use of GN(~) in its current version can be useful for limited cases in problems with l-2 variables only. The GN(~) is quite efficient and proved to be superior to GN(~) in various examples (see also [5]). The superiority of GN(~) over algorithms which use first-order derivatives only is discussed in detail by Belsley [l] and hence omitted here. 6.

DISCUSSION

AND FURTHER

CLARIFICATIONS

At this stage of software development, one needs to prepare a sequence of statements (see subroutine F UNCT in the appendix) in order to execute the algorithm as described in this paper. It is possible however to combine our algorithm with a program that “translates” the sequence as given in the second column of Table 2 into the appropriate sequence of CALL statements. This becomes useful when the sequence to evaluate F is very long. We have actually experimented with a simple “interpreter” and feel that it can be very beneficial for a user who does not want to “know” the program but intends to use it. Further improvement will be to combine our program with a formal algebraic package which, when given the function in the form of (4.4) produces the sequence of the CALL statement corresponding to the second column in Table 2 directly. This, however, makes the user dependent on large-scale programs that he usually cannot modify by himself. CALL

APPENDIX A. 1. Zntroduction The complete FORTRANprogram used in this paper is listed below. The important parameters in the program are as follows: NVAR

NS EC NTH R NFOR

IR

ID

the number the number the number the number a parameter

of elements in X. of second order derivatives. of third order derivatives. of fourth order derivatives. which determines the order of the GN: IR=2

==+. GN(I)

IR=3

3

IR=4

==+. GN(3)

GN(2)

the order of derivatives to be calculated. I D 2 I R.

A Computer Program to Minimize a Function

ND IN

X F

159

the number of elements in each vector: 1 +NVAR+NSEC

for

l+NVAR+NSEC+NTHR

for CN(2),

1 +NVAR+NSEC+NTHR+NFOR

forGN(3).

CN( 1))

the vector of variables. the vector which includes F in (2.1) in vector form.

Subroutine SEDER was not discussed in the paper. Its role is to transfer the values of the X(I), I=1 ,...,NVAR, from the vector X to the scalars X1,X2,X3, . . . . This is done to allow the user to vectorize all the NVAR variables using only one VECT subroutine. It was not possible to transfer X( I) in the CALL VECT statements, which creates the need for SEDER. The appropriate values for X 1, X2, . . . are transferred to subroutine FUNCT through a COMMON statement. The dimensions in the program allow the use of GN(~) or GN(2) for functions with 1 up to 10 variables and GN(~) with 1 up to 5 variables. Note that IR = 3 and NVAR = 10 imply NDIM= 286, which appears in the dimension of all relevant vectors. Thus, if one wants to use GN(~) with NVAR > 10, wherever 286 appears, one should replace it with the appropriate number computed by 1 +NVAR+NSEC+NTHR. To clarify, note that the calculation of NSEC, NTHR, and NFOR is performed in the MAIN. In addition, the dimension of X should be NVAR, and the dimensions of A, Hl, R, W, U, V, and Y in GNR should be at least ACm, m), HI Cm, m), R(m), W(m),UCm), V(m), and Y(m,NVAR), where m=l+NVAR for GN(~), m=l+NVAR+NSEC for GN(~), and m = 1 +NVAR+NSEC+NTHR for GN(~). Finally, we use the IMSL routine LINV2F to invert A in GNR. This requires specifying appropriately the dimension of W5 to be at least m2 + Sm. The statement COMMON/All/ in SEDER, GNR, FUNCT, and MAIN should include NVAR names (xi, x2, . . . ), and subroutine SEDER should be adjusted appropriately (see SEDER in Section A.2). A.2 C

The Program MAIN PROGRAM IMPLICIT REAL*8(A-H,O-2) COMMON/A1O/NVAR,NDIM,NSEC,NTHR,NFOR COHMON/A11/X1,X2,X3,X4,X5,X6,X7,X8,X9,X1O COMMON/A13/ID,IR DIMENSION X(lO),F(286) NVAR=3

R. KALABA AND A. TISHLJZR

160 IR=3 ID=IR C C IR--GN ORDER C--ID--DERIVATIVE ORDER C C NIT=20 DO 5 I=l,NVAR 5 X(I)=1 .DO IFCIR.GT.ID) ID=IR NSEC=CNVAR+l)*NVAR/Z NDIW=l+NVAR+NSEC IFCID.EQ.2) GOT0 1001 KK=0 DO 100 I=l,NVAR DO 100 J=I,NVAR DO 100 K=J,NVAR 100 KK=KK+l NTHR=KK NDIM=NDIM+NTHR IFCID.EQ.3) GOT0 1001 KK=0 DO 200 I=l,NVAR DO 200 J=I,NVAR DO 200 K=J,NVAR DO 200 L=K,NVAR C

200

1001

1 11

KK=KK+l NFOR=KK NDIM=NDIM+NFOR CONTINUE XC1 )=O.O7DO XCZ)=DSQRT(O.3DO) XC3)=0.1DO DO 1 IJK=l,NIT IJKl=IJK-1 YRITEC6,ll) IJKl CALL SEDER(X) CALL FUNCTC F) CALL GNRCIR,X,F) CONTINUE FORMATCZX,//,2X,‘ITERATION STOP END

NO.

‘,14,2X,//)

SUBROUTINE FUNCTCF) IMPLICIT REAL*8 CA-H,O-2) COMMON/A11/X1,X2,X3,X4,X5,X6,X7,X8,X9,X10 DIMENSION AC286),BC286),C(286),DC286),P(286),PC286)

A Computer Program to Minimize a Function S

100 4 3 2

200 8 7 6

161

4(286),SC286),FC286) ALFi=-1 . DO BETA=0.5D0 GAMMA=lO.DO CALL VECTCl,Xl,A) CALL UNIT(I)) CALL MULCONCALFA,A,C) CALL EXPPCC,D) CALL MULTCA,A,H) CALL MULCONCBETA,H,P) CALL ADDCP,D,Q) CALL MULCONCGAHMA,B,S) CALL ADDCQ,S,F) RETURN END SUBROUTINE OGENCIR,Y,F) IMPLICIT REAL*8 CA-H,O-2) CONMON/A10/NVAR,NDIM,NSEC,NTHR,NFOR DIRENSION YC71,1O),FC286) DO 1 I=l,NVAR YCl,I)=FCI+l) DO 2 IJK=l,NVAR KK=NVAR+l KJ=l DO 3 I=l,NVAR DO 4 J=I,NVAR KK=KK+l 1FCIJK.EP.I.OR.IJK.EQ.J) GOT0 GOT0 4 KJ =KJ +l YCKJ,IJK)=F(KK) CONTINUE CONTINUE CONTINUE IFCIR.EQ.2) GOT0 300 00 5 IJK=l,NVAR KK=l +NVAR+NSEC KJ=l +NVAR DO 6 I=l,NVAR DO 7 J=I,NVAR DO 8 K=J,NVAR KK=KK+l 1FCIJK.EQ.I.OR.IJK.EQ.J.OR.IJK.EQ.K) GOT0 8 KJ=KJ+l YCKJ,IJK)=FCKK) CONTINUE CONTINUE CONTINUE

100

GOT0

200

162 5

210 20 19

300

R. KALABA AND A. TISHLER CONTINUE IF(IR.EQ.3) GOT0 300 DO 19 IJK=l,NVAR KKK=l+NVAR+NSEC+NTHR KJ=l+NVAR+NSEC DO 20 I=l,NVAR DO 20 J=I,NVAR DO 20 K=J,NVAR DO 20 L=K,NVAR KKK=KKK+l IF~IJK.EQ.I.OR.IJK.EQ.J.OR.IJK.EQ.K.OR.IJK.EQ.L~GOTO GOT0 20 KJ=KJ+l YCKJ,IJK)=F(KKK) CONTINUE CONTINUE RETURN END

210

SUBROUTINE SEDER(X) IMPLICIT REAL*8 CA-H,O-2) COMMON/All/Xl,X2,X3,X4,X5,X6,X7,X8,X9,X10 DIMENSION X(10) x1=x(1 1 x2=x(2) x3=x(3) x4=x(4) x5=x(5) x6=x(6) x7=x(7) x8=x(8) x9=x(9) x10=xC10) RETURN END

61 62

SUBROUTINE GNR(IR,X,F) IMPLICIT REAL*8 CA-H,O-2) COt4MON/A1O/NVAR,NDIM,NSEC,NTHR,NFOR DIMENSION X(10),FC286),YC71,10),R(71 ,W5C5260),AC71,71 ),H1(71,71 S NGN=l +NVAR IFCIR.EQ.3) NGN=NGN+NSEC IF(IR.EQ.4) NGN=NGN+NSEC+NTHR WRITEC6,61) WRITE(6,62)CXCI),I=l,NVAR) WRITEC6,61) FORMATCZX,////) FORMATC2X,5CD20.10,2X),/,2X,5CD20.10,2X)) WRITEC6,62) F(1)

),WC71 ),VC71)

),UC71)

A

Computer Program to Minimize a Function WRITEC6,61) WRITEC6,62)CFCI+l),I=T,NVAR) WRITEC6,61) CALL DO

OGENCIR,Y,F)

1 I=l,NGN

ACI,l)=O.DO DO 2 1

2 J=l,NVAR

ACI,J+l)=YCI,J) CONTINUE ACl,l)=l.DO IFCIR.EQ.2)

GOT0

200

KK=NVAR+l

13

DO

10

II=l,NVAR

DO

13

I=l,NGN

RCI)=YCI,II) DO

11

JJ=II,NVAR

KK=KK+l DO 20

20

I=l,NGN

WCI)=YCI,JJ) CALL DO

MULTGNCIR,R,W,U)

21

I=l,NGN

21

ACI,KK)=UCI)

11

CONTINUE

10

CONTINUE IFCIR.EQ.3)

GOT0

200

LL=l+NVAR+NSEC

130

131

DO

100

II=l,NVAR

DO

130

I=1 ,NGN

RCI)=YCI,II) DO

110

JJ=II,NVAR

DO

131

I=l,NGN

VCI)=YCI,JJ) DO

111

KK=JJ,NVAR

LL=LL+l DO 132

132

I=l,NGN

WCI)=YCI,KK) CALL

MULTGNCIR,W,V,U)

CALL

MULTGNCIR,U,R,W)

DO

210

I=l,NGN

210

ACI,LL)=WCI)

111

CONTINUE

110

CONTINUE

100

CONTINUE

200

CONTINUE CALL DO

30

30

LINV2FCA,NGN,i'1,H1,O,W5,0) I=l,NVAR

XCI)=XCI)+HlCl,I+l) WRITEC6,62)CXCI),I=l,NVAR) WRITEC6,62)CHlC1,I+l),I=l,NVAR) RETURN END

103

R. KALABA AND A. TISHLJIR

164

3 2

11 10

6 5 4 200

1

SUBROUTINE MULTGNCIR,X,Y,Z) IMPLICIT REAL*8 CA-H,O-2) COMMON/A1O/NVAR,NDIM,NSEC,NTHR,NFOR DIMENSION XC71 ),YC71),2(71) zC1)=xC1)*YC1) DO 1 I=l,NVAR zc1+1 )=x(1 )*Y(I+l )+Y(l )*x(1+1 1 KK=l +NVAR DO 2 I=l,NVAR DO 3 J=I,NVAR KK=KK+l ZCKK)=XCKK)*YCl )+X(1+1 )*Y(J+l )+XCJ+l )*Y(I+l CONTINUE IFCIR.EQ.3) GOT0 200 KKK=l+NVAR+NSEC KK=l +NVAR DO 4 I=l,NVAR DO 5 J=I,NVAR KK=KK+l DO 6 K=J,NVAR J J=NVAR+l DO IO II=l,J DO 11 IJ=II,NVAR JJ=JJ+l CONTINUE J J=J J-NVAR+K KKK=KKK+l J J J =KK+K-J ZCKKK)=XCKKK)*YCl )+X(KK)*YCK+l )+XCJJJ)*YCJ+l S +X(JJ)*Y(I+1)+X(J+1)*Y(JJJ)+X(K+1)*Y(KK)~X~l~*Y~KKK~ CONTINUE CONTINUE CONTINUE RETURN END SUBROUTINE VECTCK,ALFA,A) IMPLICIT REAL*8CA-H,O-Z) COMMON/AlO/NVAR,NDIM,NSEC,NTHR,NFOR DIMENSION AC2861 DO 1 I=l,NDIM ACI)=O.DO AC1 )=ALFA ACK+l )=I .DO RETURN END SUBROUTINE UNIT(R) IMPLICIT REAL*8CA-H,O-Z) COMMON/A1O/NVAR,NDIM,NSEC,NTHR,NFOR

)+X(1

)*Y(KK)

)+XCI+l)*YCJJ)

A Computer Program to Minimize a Function

1

1

DIMENSION R(286) DO 1 I=l,NDIM R(I)=O.DO R(1 )=I .DO RETURN END SUBROUTINE NULL(A) IMPLICIT REAL*8(A-H,O-2) COMMON/A1O/NVAR,NDIM,NSEC,NTHR,NFOR DIMENSION A(2861 DO 1 I=l,NDIM A(I)=O.DO RETURN END SUBROUTINE RECP(X,Y) IMPLICIT REAL*8(A-H,O-2) DIMENSION X(286),Y(286),F(5) F(l)=l.DO/X(l) F(2)=-1 .DO/(X(l )*X(1 1) F(3)=2.DO/(X(l)*X(1)*X(l)) F(4)=-6.DO/(X(1)*X(1)*X(l~*X(l~) F(5)=-4.DO*F(4)/X(l) CALL DER(F,X,Y) RETURN END SUBROUTINE EXPP(V,C) IMPLICIT REAL*8(A-H,O-2) DIMENSION ‘.‘(286),C(286),F(5) A=V(l) F(1 )=DEXP(A) F(Z)=F(l) F(3)=F(l) F(4)=F(l) F(5)=F(l) CALL DER(F,V,C) RETURN END SUBROUTINE CONTTU(CON,U,Z) IMPLICIT REAL*8(A-H,O-Z) DIMENSION U(286),Z(286),F(5) A=U(l) B=DLOG(CON) Al =A*B F(1 )=DEXP(Al) F(2)=F(l )*B

165

R. KALABA AND A. TISHLER

lf36 FC3)=FCZ)*B FC4)=FC3)*B FCS)=FC4)*B CALL DERC F,U,Z) RETURN END SUBROUTINE LOGGCU,R) IMPLICIT REAL*8CA-H,O-Z) DIMENSION UC286),RC286),FCS) A=UCl) FC1 )=DLOGCA) FC2)=1 .DO/UCl) FC3)=-FC2)*FC2) FC4)=2.DO/CA*A*A) FC5)=-3.DO*FC4)/A CALL DERCF,U,R) RETURN END SUBROUTINE NOMCCZZ,Y) IMPLICIT REAL*8CA-H,O-Z) DIMENSION FCS),ZZC286),YC286) z=zzCl) D2=DSQRTC6.28318530717959DO) IFCZ.LT.-lS.DO)Z=-15.DO IFCZ.GT.15.DO)Z=15.D0 D’DEXPC-0.5DO*Z*Z)/D2 CALL MDNORDCZ,DD) FCI )=DD FC2)=D F(3)=-ZZCI )*D FC4)=CZZCl)*ZZCl)-l.DO)*D FC5)=ZZCl)*C3.DO-ZZCl)*ZZ(1))*D CALL DERC F,ZZ,Y) RETURN END SUBROUTINE MULCONCC,A,B) IMPLICIT REAL*8CA-H,O-Z) COMMON/A1O/NVAR,NbIM,NSEC,NTtlR,NFOR DIMENSION A(286),8(286)

DO 1 I=l,NDIM B(I)=C*A(I) 1

CONTINUE RETURN END SUBROUTINE DIVCX,Y,Z) IMPLICIT REAL*8CA-H.O-Z)

A Computer Program to Minimize a Function

167

DIMENSION xC286),YC286),Z(286),UC286) CALL RECPCY,U) CALL MULTCX,U,Z) RETURN END SUBROUTINE SUB(A,B,C) IMPLICIT REAL*8CA-H,O-Z) COMMON/A1O/NVAR,NDIM,NSEC,NTHR,NFOR DIMENSION A(286),8(286),C(286) DO 1 I=l,NDIM C(I)=A(I)-B(I)

CONTINUE RETURN END SUBROUTINE ADDCX,Y,Z) IMPLICIT REAL*8CA-H,O-Z) COMMON/A1O/NVAR,NDIM,NSEC,NTHR,NFOR DIMENSION X(286),YC286),2(286) DO 1 I=l,NDIM z(I)=x(I)+Y(I)

CONTINUE RETURN END SUBROUTINE UTTVCU,V,Z) IMPLICIT REAL*8CA-H,O-Z) DIMENSION UC286),VC286),2(286),AC286),C(286) CALL LOGGCU,A) CALL MULTCA,V,C) CALL EXPPC C,Z) RETURN END

SUBROUTINE MULTCX,Y,Z) IMPLICIT REAL*8CA-H,O-Z) COMMON/A13/ID,IR COMMON/A1O/NVAR,NDIM,NSEC,NTHR,NFOR DIMENSION XC286),YC286),2(286) zCl)=xCl)*YCl) DO 1 I=l,NVAR 1

3 2

zt1+1

)=X(1

)*Y(I+l

KK=NVAR+l DO 2 I=l,NVAR Do 3 J=I,NVAR KK=KK+l ZCKK)=XCKK)*YCl CONTINUE

)+Y(l)*x(I+l)

)+X(1+1

)*Y(J+l

)+XCJ+l)*YCI+l

)+X(1

)*YCKK)

R. KALABA AND A. TISHLER

168

11 10

6 5 4

30

31

41

IF(ID.EQ.2) GOT0 100 KKK=l+NVAR+NSEC KK=l +NVAR DO 4 I=l,NVAR DO 5 J=I,NVAR KK=KK+l DO 6 K=J,NVAR J J=NVAR+l DO 10 II=l,J DO 11 IJ=II,NVAR JJ=JJ+l CONTINUE JJ=JJ-NVARtK KKK=KKKtl JJJ=KKtK-J Z(KKK)=X(KKK)*Y(l)tX(KK)*Y(Ktl~tX(JJJ)*Y~Jt~~tX~Itl~*Y~JJ~ tX~JJ)*Y~It1~tX~Jt1~*Y~JJJ~tX~Kt~~*Y~KK~tX~~~*Y~KKK~ 0 CONTINUE CONTINUE IF(ID.EQ.3) GOT0 IJ=ltNVAR IJK=IJtNSEC IJKL=IJKtNTHR DO 20 I=l,NVAR DO 21 J=I,NVAR IJ=IJ+l DO 22 K=J,NVAR IK=IJtK-J Kl =I tNVAR DO 30 Il=l,J DO 30 IZ=Il,NVAR Kl =Kl tl JK=Kl -NVARtK IJK=IJKtl DO 23 L=K,NVAR IL=IKtL-K JL=JKtL-K Kl =I tNVAR DO 31 Il=l,K DO 31 IZ=Il,NVAR Kl=Kltl KL=Kl -NVARtL IJL=IJKtL-K Kl =I tNVARtNSEC DO 41 Il=l,I DO 41 12=Il,K DO 41 13=12,NVAR Kl =Kl tl IKL=Kl-NVARtL Kl =I tNVARtNSEC

100

A Computer Program to Minimize a Function

42

23

22 21 20 100

1

3 2

11 10

6 5

DO 42 Il=l,J DO 42 IZ=Il,K DO 42 13=12,NVAR Kl =Kl +l J KL=Kl -NVAR+L IJKL=IJKL+l Z(IJKL)=XCIJKL)*YCl)+X(IJK)*Y(L+1 )+XCIJL)*YCK+l) +X(IJ)*Y(KL)+X(IKL)*Y(J+1)+X(IK)*Y(JL)+X(IL)*Y(JK) f +XCI+l)*YCJKL)+XCJKL)*Y(I+1)+X(JK)*Y(IL) S +XCJL)*YCIK)+XCJ+l)*Y(IKL)+X(KL)*YCIJ) S +XCK+l )*YCIJL)+XCL+l )*YCIJK)+XCl )*YCIJKL) S CONTINUE CONTINUE CONTINUE RETURN END

SUBROUTINE DERCF,U,Z) IMPLICIT REAL*8CA-H,O-2) COMMON/Al S/ID, IR COMMON/A1O/NVAR,NDIM,NSEC,NTHR,NFOR DIMENSION FC5),2(286),UC286) ZCl )=FCl) DO 1 I=l,NVAR ZCI+l )=FCZ)*UCI+l) KK=NVAR*l DO 2 I=l,NVAR DO 3 J=I,NVAR KK=KK+l ZCKK)=FC3)*UCI+l )*UCJ+l )+FCZ)*UCKK) CONTINUE IFCID.EQ.2) GOT0 100 KK=l +NVAR KKK=l+NVAR+NSEC DO 4 I=l,NVAR DO 5 J=I,NVAR KK=KK+l DO 6 K=J,NVAR JJ=NVAR+l DO 10 II=l,J DO 11 IJ=II,NVAR JJ=JJ+l CONTINUE JJ=JJ-NVAR+K KKK=KKK+l J J J=KK+K-J ZCKKK)=FC4)*UCI+l )*U(J+l )*lJCK+l )+FCZ)*U(KKK) S +FC3)*CUCJJJ)*UCJ+1)+UCJJ)*U(I+1)+U(KK)*U~K+l~~ CONTINUE

169

R. KALABA AND A. TISHLEX

170 4

CONTINUE IF(ID.EQ.3)

GOT0

100

IJ=l+NVAR IJK=IJ+NSEC IJKL=IJK+NTHR DO

20

I=l,NVAR

DO

21

J=I,NVAR

IJ=IJ+l DO

22

K=J,NVAR

IK=IJ+K-J

30

31

41

42

23

22 21 20 100

Kl =l +NVAR DO 30 Il=l,J DO 30 IZ=Il,NVAR Kl =Kl +I JK=Kl -NVAR+K IJK=IJK+l DO 23 L=K,NVAR IL=IK+L-K JL=JK+L-K Kl =I +NVAR DO 31 Il=l,K DO 31 IZ=Il,NVAR Kl =Kl +I KL=Kl -NVAR+L IJL=IJK+L-K Kl =I +NVAR+NSEC DO 41 Il=l,I DO

41

12=Il,K

DO

41

13=12,NVAR

Kl =Kl +I IKL=Kl -NVAR+L Kl =I +NVAR+NSEC DO 42 Il=l,J DO 42 IZ=Il,K DO 42 13=12,NVAR Kl =Kl +I JKL=Kl -NVAR+L IJKL=IJKL+l Z(IJKL)=F(S)*U(L+l )*U(K+l )*U(J+l )*U(I+l) $ +F(4)*(U(KL)*U(J+l )*lJ(I+l )+U(K+l )*U(JL)*U(I+l) f +U(K+l )*lJ(J+l )*U(IL)+U(L+l )*U(JK)*U(I+l) 3 +U(L+l)*U(J+1)*U(IK)+U(L+1)*U(K+1)*U(IJ)~ t +F(3)*(U(JKL)*U(I+l)+U(JK)*U(IL)+U(JL)*U~IK~ f +U(J+1)*U(IKL)+U(KL)*U(IJ~+U(K+1)*U(IJL) 16 +U(L+1)*U(IJK))+F(2)*U(IJKL) CONTINUE CONTINUE CONTINUE RETURN END

A Computer Program to Minimize a Function

A.3

The Function for the Black-Scholes

171

Valuation Formula

SUBROUTINE FUNCTlCN,W,XX,C,T,SS) IMPLICIT REAL*8CA-H,O-Z) COMHON/A11/X1,X2,X3,X4,X5,X6,X7,X8,X9,X1O DIMENSION W~90~,xx~90~,C~90~,T~90~,sS~286~,Y1(286~,Y2~286~, $ R(286),A(286),8(286),P(286),D(286),E(286),F(286),G(286), t H(286),P(286),S(286),Dl(286),D2(286),V(286),U(286),Y3(286) CALL NULLCSS) CALL UNIT(R) ONE=1 .DO DO 1 I=l,N Al=XX(I)/C(I) Al =DLOGC Al 1 A2=0.5DO A3=TC I) A4=DSQRTCA3) AS=XXCI) A6=CC I) A7=-A3 A8=WC I) A200=DSQRTCU.O50DO) CALL VECT(l,Xl ,Yl 1 CALL VECTC2,X2,Y2) CALL VECTC3,X3,Y3) CALL RULCONC Al ,R,A) CALL MULTCY2,Y2,B) CALL MULCONCA2,B,P) CALL ADDCYl,P,D) CALL SUBCYl,P,E) CALL MULCONCA3,D,F) CALL MULCONCA3,E,G) CALL ADDCA,F,H) CALL ADDCA,G,Q) CALL MULCONCA4,YL,S) CALL DIVCH,S,Dl) CALL DIVCQ,S,DZ) CALL MULCONCA7,Yl,V) CALL EXPPCV,U) CALL MULCONCA6,U,A) CALL NOMCCDl,B) CALL NOMCCD2,D) CALL MULCONCA5,B,E) CALL MULTCA,D,F) CALL SUBCE,F,G) CALL MULCONCAI,R,E) CALL SUBCE,G,H) CALL SUBCH,Y3,D) CALL MULTCD,D,F) CALL ADDCSS,F,Q)

R. KALABA AND A. TISHLER

172

1

CALL MULCON(ONE,Q,SS) CONTINUE RETURN END

NOTES. The vectors U, XX, C, and T in FUNCTI include the data (generated in a subroutine which is not given here) corresponding to the values of ii?,, y,, ci. and rl, i = 1,. , N, in the discussion in Section 6. The authors warmly thank their colleagues H. Kagiwada, N.Rasakhoo, K. Spingarn, and L. Tesfatsion for generously providing ideas and inspiration for this study.

REFERENCES D. A. Belsley, On the efficient computation of the nonlinear full-information maximum-likelihood estimator, J. EConometrics 14203-255 (1980). F. Black and M. Scholes, The Pricing of Options and Corporate Liabilities, J. Pal. Ecott. 81637-659 (1973). S. M. Goldfeld and R. E. Quandt, Nonlinear Methods in Ecunumetrics, NorthHolland, Amsterdam, 1972. R. Kalaba, N. Rasakhoo, and A. Tishler, Nonlinear least squares via automatic derivative evaluation, Appl. Math. Cwmput., to appear. R. Kalaba and A. Tishler, Optimization in non-linear models via automatic derivative evaluation, Dept. of Economics Working Paper ##8212, Univ. of Southem California, 1982. R. Kalaba and A. Tishler, A generalized Newton algorithm using higher order derivatives, 1. Optim. Theory Appl., to appear. R. Kalaba and A. Tishler, A generalized Newton algorithm to minimize a function with many variables using computer evaluated exact higher order derivatives, J. Optim. Theory Appl., to appear. D. A. Pierre, Optimization Theory with Applications, Wiley, 1969.