- Email: [email protected]
A numerical search technique for optimal capital accumulation programs
Socio-Econ. Plan. Sci. Vol.9, pp. 4%47.PergamonPress 1975.Printedin Great Britain.
NOTE A NUMERICAL
SEARCH TECHNIQUE FOR OPTIMAL ACCUMULATION PROGRAMS
CAPITAL
JAMESA. YUNKER Department of Economics, Western Illinois University, Macomb, Illinois 61455 (Receit:ed 8 J u n e 1974)
Finding an explicit analytical solution to optimal capital accumulation problems usually requires fairly unrealistic simplifying assumptions with respect to the utility, production, and other functions employed. At the same time. if these functions are left algebraically unspecified, the obtainable results are limited and usually depend on the steady-state growth assumption. This note proposes a method for finding a numerical solution to the optimal capital accumulation problem. The idea is to set an initial value for capital accumulation, extrapolate a numerical solution through the Euler-Lagrange equation, and on the basis of the error between the actual and desired terminal capital stock, set a new initial capital accumulation value and extrapolate again. This is continued until the actual terminal capital is sufficiently close to the desired. An example of application is offered which suggests that the method could be successfully applied, without undue costs in computer or programming time, to a wide range of problems that would not be susceptible to solution using other methods.
where K = capital stock, L = labor (in general exogenously determined by time), F is a concave production function (in K and L), K = capital accumulation, c = per capita consumption, U is a concave utility function, d(t) is the discount function, Tis the time horizon, Ko and K T are respectively initial and terminal capital stock, (3) says that capital stock cannot decline faster than normal, proportional depreciation 6 K , and (4) says that per capita consumption cannot drop below the subsistence minimum c,,. The optimal time path of K , K*(t), must satisfy the Euler Lagrange equation of the classical calculus of variations, which may be written: FK
subject to the boundary conditions: K(T)= K T
(2)
and subject to the constraints: f ; >_ - 6 K
(3)
c > c,,
(4)
a
u"
L
d
U'
b.
(5)
A solution of this second-order differential equation which satisfies the boundary conditions and the constraints will necessarily (owing to the concavity of the production and utility functions) result in a maximum value of (1) subject to the conditions and constraints.:[: In order to obtain an explicit analytical solution to this problem, it is necessary to specify algebraic forms for U(c), F ( K , L , t ) , L(t) and d(t). The complexities of differential equations being what they are, in order to be able to obtain an analytical solution, the forms of these functions have to be rather carefully selected, and they have to be fairly elementary.§ For example, even using a Cobb-Douglas production function will complicate the problem to the level of insoluability.II Some analysis can be undertaken without approaching the question of an explicit analytical solution, but the key simplifying assumption usually required here. that the system be in steady-state equilibrium, may be unrealistic near the boundary of the planning horizon if steady-state growth does not bring K near to the desired K r at time T. The purpose of this note is to suggest a method of ascertaining a numerical solution to the problem stated in (1V(4), given algebraically and numerically specified U, F, Land d functions. There is no requirement that these functions be of any particular type--they can be as simple or as complicated as is necessary to achieve the desired level of realism. While it cannot be guaranteed that there are no functions sufficiently complicated to make the method impracticable, whether because of a failure to converge or excessive requirements in computer lime, successful experiments with a sample model strongly suggest that the method could probably be used in a wide variety of realistic problems of practical interest. The basic differential equation of the capital accumulation problem is (5). The numerical extrapolation (also
I. THE METHOD The standard finite horizon optimal capital accumulation problem may be stated as follows.t Maximize the functional:
K(0) = K o
£
t For example: See Ref. [2], p. 363; See Ref. [3], p. 267. :~ See Ref. [4], Part IV, for the relationship between the calculus of variations and the maximum principle approaches to the dynamic control problem. § For example, Burmeister and Dobell (1970) utilize F = b K in a worked example (p. 374), as do Hadley and Kemp [3], p. 267. II For example, even if we simplify to the extent of letting U = V - e -'~ (so that U " / U ' = - v), d = 1, L = 1, if we make F = A K ~ (0 < ~ < 1), then equation (5) becomes: /< + (1 - Ifl )~ a K z - t = 0 which is impervious to analytical solution. 45
46
Note
called simulation or "solution") of any paramatrically specified second order differential equation, no matter how difficult or impossible it might be in terms of an explicit analytical solution, is easily accomplished given two initial conditions Ko and K,. But the logic of the calculus of variations approach requires the solution to satisfy two boundary conditions K(0) = K~,, K ( T ) = Kr. The basic idea of the proposed method is to set /4. (Ko is given), numerically extrapolate the solution using the discrete-period analogue to (5) for Tperiods (applying (3) and (4) if necessary), and compare K ( T ) with the desired Kr. On the basis of the error, a revised K . is set on which another extrapolation is based. By means of adjustments in /~o, the K ( T ) value is adjusted to the desired K-r level. In other words, we are considering K ( T ) as a function of /~,. The fnnction is given by the discrete-period analogue of (5), taking into account also (3) and (4). K0 is varied so as to set K ( T ) to the particular desired K r level. When this has been accomplished, we have a numerical time path in K that satisfies (2~ (5). and which therefore gives a m a x i m u m of the functional (1) within the constraints. The idea is closely related to numerical search as used in, for example, singlevariate optimization problems. Say y =.f{x) is a very complicated but numerically specified function for which we would like to know the value x* that maximizes it. The complexity of the function precludes setting./'(x) - 0 and solving explicitly for x. The numerical approach is to vary x in some systematic way and calculate the corresponding y values, v is moved in whatever direction raises y, and when y is finally observed to decline, a converging procedure is commenced to bring x to the x* level (to the desired degree of precision)at which y is maximized. Any search technique used to bring x to the x* level at which the function of x, y, is at a peak, may also be used in the proposed method to bring Ko to the /£* level at which the function of /£,, K(T), is at K,rA" The only difference is that the converging procedure starts when the given K r is passed by K ( T ) rather than a m a x i m u m value of K ( T ~ w e are not. of course, looking for the m a x i m u m value of K ( T ) (it is doubtful whether a finite m a x i m u m even exists). A convenient method of organizing the search is ~ follows. The lower and upper boundaries respectively of a feasible initial value for capital accumulation (K~} are set by - 6K0 (the m i n i m u m capital decline) and F(Ko, Lo, O) - cmL0 (the m a x i m u m capital increase given the subsistence minim u m in consumption). Divide the interval between - 6 K , and Fo - c,~L~ into n equal sub-intervals and extrapolate for Tperiods according to (5), (4) and (3) (n + l) times using the K,>'s equal to the sub-interval boundaries. This is likely to suggest that the function relating/£o and K ( T ) is monotonic. At any rate take those intervals (if the function is monotonic there will be just one such) for which the corresponding K ( T ) values bracket the desired K r, and use a convergence procedure to a s c e r t a i n / ~ , to the desired level of accuracy. For example, say K ( T ) appears that it would be a monotonic increasing function of/{(~ in the interval. The/~,> on iteration i is related to the/(,) on iteration (i - 1) by the following rule: /~'],i -g,l.i
--
/~0.i I + ri where ri = ri- 1 if K ( T ) i_ t < K r /'~f).i 2 or- "i where r i = ½ri_ i if K(T)i-1 > K r
+ Numerical search methods of this sort are covered in detail by Beveridge and Schechter [1] Chapter 6. {; Numerical extrapolation of differential equations is covered by Kaplan [5] Chapter 10.
where r,, the first iteration adjustment factor, is cqual t(~ thc length oftheinterval, and Ko.~> the initial /~0, is equal to the lefl-hand boundary of the interval in which K~ is included. This is continued until iK(T) - K l i < ~, where represents the desired level of accuracy. 2. ~.N EXAMPLE OF APPLICATION The method was tried out on the [bllowing model: U(c)-=
I
c
"
(6)
F(K.L} - ,4K'le
(7)
l , [ t l - L,e"'
18)
d(t) = e ;':
19)
This utility function is especially convenient when utilized in equation (5), since in this case U ' / U ' = - v. The labor growth and discounting functions are standard. The main "complicating factor" here is the use of the C o b b - D o u g l a s production function. This would preclude an analytical solution of(5). Equation [5) can be written, using (6)-(9), as: /~ = (n + FK)I~ + 07 + p - FD(L/r) + n [ F c L - F).
(I0)
Given initial values K(0t and /~(0), /~(0) is calculated from this equation. Then /£(I) = /~(0) + K(0), K(t) = K(O) + K(0L and the values of K{l), /~(01 are used in (10i to obtain /~(2).~ This is continued for the desired number of periods. thus generating a numerical approximation to an analytical solution. On the basis of the difference between K ( T ) and the desired Kr, a new /<(0) is set according to the method described above. The problem was programmed in Fortran and run on an IBM S/360 m50computer. Experience with this problem indicated virtually insignificant m a n - h o u r costs in writing and debugging the program, and in compilation and execution computer time. Double precision was required to bring about convergence. In all. twelve different sets of parameter values were used, and in each case convergence was achieved within 50 iterations to within five {terminal capital stock being measured in thousands), each iteration requiring extrapolation over twenty-five periods. In each case. K ( T ) was a monotonic increasing function of/(,, and there was thus a unique value of K o that produced KI. The non-varied parameters were A = 500, :< = (t.4, fl = 0.6, K 0 = 7000, L0 = 1, T = 25. Systematic variation was tried in the utility parameter (v = 0.1, 0.5, 1.0), the discount rate (p = 0.01, 0.05, 0.1), the labor growth rate (n = 0, 0.01, 0.02, 0.04), a n d the terminal capital stock ( K r = 7000, 800(t. 9000, 10000). In general, the results indicated that of the c, p, and t¢ parameters, only the labor growth rate has any significant effect on the optimal solution. The tendency is for the capital stock to grow at the same proportionate rate as labor {given the linear homogeneity of the production function this implies a complete steady-state system) for as much of the planning period as possible, with the last few periods being devoted to as quick an approach of K to K r as is possible given the depreciation or subsistence consumption constraints (the depreciatkm conslraint is operative if the n-rate steady-state path of K(t) passes considerably above K,r and the subsistence constraint if it passes considerably below it). Figure I is a sketch of different solution paths in K for which all parameters were constant except for n. Table I shows K(t) values at intervals for three different runs where only the discount rate p was varied. This shows the insignificant effect of p.
Note
47
K(T}
/
1400C
1
n =O.O4
1300C
Kr
t
1
]¢~
I
1
f~
n
1
I
t
IU- V, ,
1200(
IIOOC
IO00C
(A)
(B)
(C)
(D)
K0 900C
Fig. 2:
K~8ooc
7000
I0
20
25
t
Fig. 1.
3. E V A L U A T I O N
In view of the success of the experiment, it would seem that the proposed method of numerical solution could well have a great deal to offer in the area of practical calculation of optimal capital accumulation programs. This is because it can be applied to many problems which are not amenable to analytical solution because of the complexity of the functions involved.
Table 1. Optimal capital accumulation for various discount rates
t
K,(p = 0.01)
K,(p = 0.05)
Kt(p = 0.1)
4
7007
7007
7007
9
7017
7017
7016
14
7027
7026
7025
19 20 21 22 23 24 25
7036 7038 7038 7037 7033 7024 7003
7035 7036 7036 7035 7032 7022 7001
7033 7034 7034 7033 7029 7019 6998
Parameters: A, et, fl, Ko, Lo, Tas stated in text, n = 0, v = 0"5, KT = 7000.
The most obvious potential problem with the method is that there might be some cases where the implied K(T) = ~(/
1. G. S. G. Beveridge and R. S. Schechter, Optimization: Theory and Practice. McGraw-Hill, New York (1970). 2. E. Burmeister and A. R. Dobell, Mathematical Theories of Economic Growth. Macmillan, New York (1970). 3, G. Hadley and M. C. Kemp, Variational Methods in Economics. Elsevier, New York (197 l). 4, M. D. Intriligator, Mathematical Optimization atld Economic Theory. Prentice-Hall, Englewood Cliffs, N.J. (1971). 5. W. Kaplan, Ordinary Differential Equations. AddisonWesley, Reading, Mass. (1958).
NOTE A NUMERICAL
SEARCH TECHNIQUE FOR OPTIMAL ACCUMULATION PROGRAMS
CAPITAL
JAMESA. YUNKER Department of Economics, Western Illinois University, Macomb, Illinois 61455 (Receit:ed 8 J u n e 1974)
Finding an explicit analytical solution to optimal capital accumulation problems usually requires fairly unrealistic simplifying assumptions with respect to the utility, production, and other functions employed. At the same time. if these functions are left algebraically unspecified, the obtainable results are limited and usually depend on the steady-state growth assumption. This note proposes a method for finding a numerical solution to the optimal capital accumulation problem. The idea is to set an initial value for capital accumulation, extrapolate a numerical solution through the Euler-Lagrange equation, and on the basis of the error between the actual and desired terminal capital stock, set a new initial capital accumulation value and extrapolate again. This is continued until the actual terminal capital is sufficiently close to the desired. An example of application is offered which suggests that the method could be successfully applied, without undue costs in computer or programming time, to a wide range of problems that would not be susceptible to solution using other methods.
where K = capital stock, L = labor (in general exogenously determined by time), F is a concave production function (in K and L), K = capital accumulation, c = per capita consumption, U is a concave utility function, d(t) is the discount function, Tis the time horizon, Ko and K T are respectively initial and terminal capital stock, (3) says that capital stock cannot decline faster than normal, proportional depreciation 6 K , and (4) says that per capita consumption cannot drop below the subsistence minimum c,,. The optimal time path of K , K*(t), must satisfy the Euler Lagrange equation of the classical calculus of variations, which may be written: FK
subject to the boundary conditions: K(T)= K T
(2)
and subject to the constraints: f ; >_ - 6 K
(3)
c > c,,
(4)
a
u"
L
d
U'
b.
(5)
A solution of this second-order differential equation which satisfies the boundary conditions and the constraints will necessarily (owing to the concavity of the production and utility functions) result in a maximum value of (1) subject to the conditions and constraints.:[: In order to obtain an explicit analytical solution to this problem, it is necessary to specify algebraic forms for U(c), F ( K , L , t ) , L(t) and d(t). The complexities of differential equations being what they are, in order to be able to obtain an analytical solution, the forms of these functions have to be rather carefully selected, and they have to be fairly elementary.§ For example, even using a Cobb-Douglas production function will complicate the problem to the level of insoluability.II Some analysis can be undertaken without approaching the question of an explicit analytical solution, but the key simplifying assumption usually required here. that the system be in steady-state equilibrium, may be unrealistic near the boundary of the planning horizon if steady-state growth does not bring K near to the desired K r at time T. The purpose of this note is to suggest a method of ascertaining a numerical solution to the problem stated in (1V(4), given algebraically and numerically specified U, F, Land d functions. There is no requirement that these functions be of any particular type--they can be as simple or as complicated as is necessary to achieve the desired level of realism. While it cannot be guaranteed that there are no functions sufficiently complicated to make the method impracticable, whether because of a failure to converge or excessive requirements in computer lime, successful experiments with a sample model strongly suggest that the method could probably be used in a wide variety of realistic problems of practical interest. The basic differential equation of the capital accumulation problem is (5). The numerical extrapolation (also
I. THE METHOD The standard finite horizon optimal capital accumulation problem may be stated as follows.t Maximize the functional:
K(0) = K o
£
t For example: See Ref. [2], p. 363; See Ref. [3], p. 267. :~ See Ref. [4], Part IV, for the relationship between the calculus of variations and the maximum principle approaches to the dynamic control problem. § For example, Burmeister and Dobell (1970) utilize F = b K in a worked example (p. 374), as do Hadley and Kemp [3], p. 267. II For example, even if we simplify to the extent of letting U = V - e -'~ (so that U " / U ' = - v), d = 1, L = 1, if we make F = A K ~ (0 < ~ < 1), then equation (5) becomes: /< + (1 - Ifl )~ a K z - t = 0 which is impervious to analytical solution. 45
46
Note
called simulation or "solution") of any paramatrically specified second order differential equation, no matter how difficult or impossible it might be in terms of an explicit analytical solution, is easily accomplished given two initial conditions Ko and K,. But the logic of the calculus of variations approach requires the solution to satisfy two boundary conditions K(0) = K~,, K ( T ) = Kr. The basic idea of the proposed method is to set /4. (Ko is given), numerically extrapolate the solution using the discrete-period analogue to (5) for Tperiods (applying (3) and (4) if necessary), and compare K ( T ) with the desired Kr. On the basis of the error, a revised K . is set on which another extrapolation is based. By means of adjustments in /~o, the K ( T ) value is adjusted to the desired K-r level. In other words, we are considering K ( T ) as a function of /~,. The fnnction is given by the discrete-period analogue of (5), taking into account also (3) and (4). K0 is varied so as to set K ( T ) to the particular desired K r level. When this has been accomplished, we have a numerical time path in K that satisfies (2~ (5). and which therefore gives a m a x i m u m of the functional (1) within the constraints. The idea is closely related to numerical search as used in, for example, singlevariate optimization problems. Say y =.f{x) is a very complicated but numerically specified function for which we would like to know the value x* that maximizes it. The complexity of the function precludes setting./'(x) - 0 and solving explicitly for x. The numerical approach is to vary x in some systematic way and calculate the corresponding y values, v is moved in whatever direction raises y, and when y is finally observed to decline, a converging procedure is commenced to bring x to the x* level (to the desired degree of precision)at which y is maximized. Any search technique used to bring x to the x* level at which the function of x, y, is at a peak, may also be used in the proposed method to bring Ko to the /£* level at which the function of /£,, K(T), is at K,rA" The only difference is that the converging procedure starts when the given K r is passed by K ( T ) rather than a m a x i m u m value of K ( T ~ w e are not. of course, looking for the m a x i m u m value of K ( T ) (it is doubtful whether a finite m a x i m u m even exists). A convenient method of organizing the search is ~ follows. The lower and upper boundaries respectively of a feasible initial value for capital accumulation (K~} are set by - 6K0 (the m i n i m u m capital decline) and F(Ko, Lo, O) - cmL0 (the m a x i m u m capital increase given the subsistence minim u m in consumption). Divide the interval between - 6 K , and Fo - c,~L~ into n equal sub-intervals and extrapolate for Tperiods according to (5), (4) and (3) (n + l) times using the K,>'s equal to the sub-interval boundaries. This is likely to suggest that the function relating/£o and K ( T ) is monotonic. At any rate take those intervals (if the function is monotonic there will be just one such) for which the corresponding K ( T ) values bracket the desired K r, and use a convergence procedure to a s c e r t a i n / ~ , to the desired level of accuracy. For example, say K ( T ) appears that it would be a monotonic increasing function of/{(~ in the interval. The/~,> on iteration i is related to the/(,) on iteration (i - 1) by the following rule: /~'],i -g,l.i
--
/~0.i I + ri where ri = ri- 1 if K ( T ) i_ t < K r /'~f).i 2 or- "i where r i = ½ri_ i if K(T)i-1 > K r
+ Numerical search methods of this sort are covered in detail by Beveridge and Schechter [1] Chapter 6. {; Numerical extrapolation of differential equations is covered by Kaplan [5] Chapter 10.
where r,, the first iteration adjustment factor, is cqual t(~ thc length oftheinterval, and Ko.~> the initial /~0, is equal to the lefl-hand boundary of the interval in which K~ is included. This is continued until iK(T) - K l i < ~, where represents the desired level of accuracy. 2. ~.N EXAMPLE OF APPLICATION The method was tried out on the [bllowing model: U(c)-=
I
c
"
(6)
F(K.L} - ,4K'le
(7)
l , [ t l - L,e"'
18)
d(t) = e ;':
19)
This utility function is especially convenient when utilized in equation (5), since in this case U ' / U ' = - v. The labor growth and discounting functions are standard. The main "complicating factor" here is the use of the C o b b - D o u g l a s production function. This would preclude an analytical solution of(5). Equation [5) can be written, using (6)-(9), as: /~ = (n + FK)I~ + 07 + p - FD(L/r) + n [ F c L - F).
(I0)
Given initial values K(0t and /~(0), /~(0) is calculated from this equation. Then /£(I) = /~(0) + K(0), K(t) = K(O) + K(0L and the values of K{l), /~(01 are used in (10i to obtain /~(2).~ This is continued for the desired number of periods. thus generating a numerical approximation to an analytical solution. On the basis of the difference between K ( T ) and the desired Kr, a new /<(0) is set according to the method described above. The problem was programmed in Fortran and run on an IBM S/360 m50computer. Experience with this problem indicated virtually insignificant m a n - h o u r costs in writing and debugging the program, and in compilation and execution computer time. Double precision was required to bring about convergence. In all. twelve different sets of parameter values were used, and in each case convergence was achieved within 50 iterations to within five {terminal capital stock being measured in thousands), each iteration requiring extrapolation over twenty-five periods. In each case. K ( T ) was a monotonic increasing function of/(,, and there was thus a unique value of K o that produced KI. The non-varied parameters were A = 500, :< = (t.4, fl = 0.6, K 0 = 7000, L0 = 1, T = 25. Systematic variation was tried in the utility parameter (v = 0.1, 0.5, 1.0), the discount rate (p = 0.01, 0.05, 0.1), the labor growth rate (n = 0, 0.01, 0.02, 0.04), a n d the terminal capital stock ( K r = 7000, 800(t. 9000, 10000). In general, the results indicated that of the c, p, and t¢ parameters, only the labor growth rate has any significant effect on the optimal solution. The tendency is for the capital stock to grow at the same proportionate rate as labor {given the linear homogeneity of the production function this implies a complete steady-state system) for as much of the planning period as possible, with the last few periods being devoted to as quick an approach of K to K r as is possible given the depreciation or subsistence consumption constraints (the depreciatkm conslraint is operative if the n-rate steady-state path of K(t) passes considerably above K,r and the subsistence constraint if it passes considerably below it). Figure I is a sketch of different solution paths in K for which all parameters were constant except for n. Table I shows K(t) values at intervals for three different runs where only the discount rate p was varied. This shows the insignificant effect of p.
Note
47
K(T}
/
1400C
1
n =O.O4
1300C
Kr
t
1
]¢~
I
1
f~
n
1
I
t
IU- V, ,
1200(
IIOOC
IO00C
(A)
(B)
(C)
(D)
K0 900C
Fig. 2:
K~8ooc
7000
I0
20
25
t
Fig. 1.
3. E V A L U A T I O N
In view of the success of the experiment, it would seem that the proposed method of numerical solution could well have a great deal to offer in the area of practical calculation of optimal capital accumulation programs. This is because it can be applied to many problems which are not amenable to analytical solution because of the complexity of the functions involved.
Table 1. Optimal capital accumulation for various discount rates
t
K,(p = 0.01)
K,(p = 0.05)
Kt(p = 0.1)
4
7007
7007
7007
9
7017
7017
7016
14
7027
7026
7025
19 20 21 22 23 24 25
7036 7038 7038 7037 7033 7024 7003
7035 7036 7036 7035 7032 7022 7001
7033 7034 7034 7033 7029 7019 6998
Parameters: A, et, fl, Ko, Lo, Tas stated in text, n = 0, v = 0"5, KT = 7000.
The most obvious potential problem with the method is that there might be some cases where the implied K(T) = ~(/
1. G. S. G. Beveridge and R. S. Schechter, Optimization: Theory and Practice. McGraw-Hill, New York (1970). 2. E. Burmeister and A. R. Dobell, Mathematical Theories of Economic Growth. Macmillan, New York (1970). 3, G. Hadley and M. C. Kemp, Variational Methods in Economics. Elsevier, New York (197 l). 4, M. D. Intriligator, Mathematical Optimization atld Economic Theory. Prentice-Hall, Englewood Cliffs, N.J. (1971). 5. W. Kaplan, Ordinary Differential Equations. AddisonWesley, Reading, Mass. (1958).