Finding Numerical Results to Large Scale Economic Models Using Path-Following Algorithms: A Vintage Capital Example

Board of Governors of the Federal Reserve System International Finance Discussion Papers Number 728 June 2002 FINDING NUMERICAL RESULTS TO LARGE SCALE ECONOMIC MODELS USING PATH-FOLLOWING ALGORITHMS: A VINTAGE CAPITAL EXAMPLE Brett D. Berger NOTE: International Finance Discussion Papers are preliminary materials circulated to stimulate discussion and critical comment. References to International Finance Discussion Papers (other than an acknowledgment that the writer has had access to unpublished material) should be cleared with the author or authors. Recent IFDP’s are available on the Web at www.federalreserve.gov/pubs/ifdp/

FINDING NUMERICAL RESULTS TO LARGE SCALE ECONOMIC MODELS USING PATH-FOLLOWING ALGORITHMS: A VINTAGE CAPITAL EXAMPLE Brett D. Berger* Abstract: This paper describes the numerical optimization methods used in Berger (2001) to find the complete time paths of key economic variables in neoclassical vintage capital models. An interior and a non-interior point method are discussed. Both of the methods are part of the general class of “path-following” algorithms. These algorithms can be efficiently applied to convex programming problems; and due to the standard shape of production and utility functions, many economic problems can be written as convex programming problems. Vintage capital models add scale and complexity to standard growth models because one must now handle the dynamics of multiple capital stocks. This increase in complexity will often prevent the discovery (or existence) of closed form solutions, making numerical solutions of the type found in Berger (2001) necessary. Keywords: optimization, productivity, technology *The author is an economist in the Division of International Finance, Board of Governors of the Federal Reserve System. He can be reached at Mail Stop 42B, Federal Reserve Board, Washington, DC 20551; brett.d.berger@frb.gov. The majority of this research was done while the author was a graduate student at the University of Washington. The author is grateful for the comments from his supervisory committee, Stephen Turnovsky, Charles Engel, Richard Hartman, and Terry Rockafellar, and for comments from colleagues at the Board. The views of this paper are solely the responsibility of the author and should not be interpreted as reflecting the views of the Board of Governors of the Federal Reserve System.

1. Introduction To more accurately mirror the world, economic models have become increasingly complex. In many cases, it is necessary to find numerical solutions to parameterized versions of the models if one is to obtain information on transitional dynamics. Optimization versions of vintage capital models are an example of such models. There are two general methods for modeling technology, “embodied” and “disembodied.” Vintage capital or “embodied technology” models are models in which technology is inherent in the capital stock of the economy (i.e. the physical equipment, buildings, etc.). The embodiment of technology implies a heterogeneous capital stock, and that it is necessary to invest in new capital if the economy is to reap the benefits of new technology. It is obvious that much of the technological growth in the world is of the embodied variety, however, most economic models assume technology is disembodied owing to the difficulty of tracking multiple capital stocks. This paper discusses numerical optimization methods that are capable of handling the large-scale nature of vintage capital models, and their application to the three basic neoclassical vintage capital models: putty-putty, clay-clay, and putty-clay.1 The optimization methods used in this paper have the potential for application to a wide array of economic models. Diminishing returns in both productive factors and utility often implies that models can be posed as convex programming problems, which mathematical theory states will have global solutions (assuming a solution exists). The algorithms 1 A detailed description of the results can be found in Berger (2001). 1

discussed in this paper are guaranteed to converge when applied to convex programming problems. 2. Optimization Method 1: Interior Point Method The algorithms used to obtain the numerical solutions of Berger (2001) are part of the general class of “path-following” algorithms. An interior point method was utilized to obtain approximations to the solutions of the putty-putty and clay-clay models. Interior point methods are so named because the choice variables lie in the strictly ( ) feasible set of the optimization problem, P , at each iteration. This is achieved by the addition of penalty functions, f , to the objective function, f . The penalty functions are 1 0 functions of the constraint values. The penalty functions increase in value as the boundary of the feasible set is approached, and take on infinite values at the boundary, thereby creating a “barrier” to exiting the feasible set. The penalty functions are multiplied by a positive scalar, µ, which is monotonically decreasing in the iterations, k. k The solutions to this barrier problem, ( P(µ) ) , as µ →0, comprise the “central path” k k ( ) which is followed to the solution of P . Below is the formal statement of the standard ( ) problem P and the resulting optimality conditions. 2

Minimize f (x) 0 over x∈Rn subject to: G(x)≤0 ( ) P where: f :Rn (cid:54) R 0 G:Rn (cid:54) Rm G(x)=( g (x),g (x),…,g (x) ) ' with g :Rn (cid:54) R 1 2 m i ( ) P is a general problem formulation because equality constraints can be included by ( ) splitting the constraint into two inequality constraints. The lagrangian for P is: L(x;λ)= f (x)+λ'G(x) (1) 0 ( ) This leads to the optimality conditions of P : ∇f +∇G(x)'λ=0 (2) 0 G(x)≤0 (3) λ≥0 (4) ΛG(x)=0 (5) where Λ is a diagonal matrix with λ on the diagonal. (5) is the complementarity condition stating that either λ, g , or both must equal zero for all i. Let the penalty i i function be the standard log barrier function, such that for z∈Rm:   ∑m ln(z ) if z >0 ∀i ln(z)=  i=1 i i (6)  ∞ otherwise Thus, ( P(µ) ) ,the log barrier version of the standard problem, is: k 3

Minimize f 0 (x)−µ k ln (−G(x) ) ( P(µ) ) k over x∈Rn where µ is a positive scalar. This leads to the optimality conditions: k ∇f +∇G(x)'µ−G(x) −1 e=0 (7) 0 k G(x)<0 (8) where z =diag(z)∈Rm×m and e=(1,…,1)'∈Rm. (8) must hold at a solution because otherwise the value of the objective is infinity by the definition of the penalty function. Convex programming problems are minimization problems in which the objective function and the inequality constraints (written as g ≤0) are convex functions, and the i equality constraints (written as g =0) are affine functions. Theory states that all i solutions to convex programming problems are global solutions. Further, when a solution exists, convex programming problems with strictly convex objective functions have a unique, global solution. For convex programming problems, Proposition 4.1.1 of Bertsekas (1999) states { } that “every limit point of a sequence x generated by a barrier method is a global k { } minimum of the original constrained problem.” In other words, the sequence x , the k solutions to the problems ( P(µ) ) as µ →0, goes to a global minimum, x , of ( P ) . k k Both ( P ) and ( P(µ) ) can be rewritten in equivalent equilibrium formats where k slack variables, restricted to be non-negative, are added to the inequality constraints to make them equality constraints. This is advantageous because of the resulting structure of the Jacobian of the optimality conditions, and readily obtained dual variables. By two 4

problems being “equivalent” it is meant that a stationary point of one problem is also a stationary point of the other. ( ) The equilibrium formatted problem, E and its corresponding lagrangian follow: Minimize f (x) 0 over x∈Rn ( ) subject to: E s≥0 G(x)+s =0 L(x,s;λ)= f (x)+λ' ( G(x)+s )+γ'(−s) (9) 0 ( ) (9) leads to the optimality conditions of E : ∇f +∇G(x)'λ=0 (10) 0 λ−γ=0 ⇒ λ=γ (11) G(x)+s =0 (12) −s≤0 (13) γ≥0 (14) Γs =0 where Γ=diag(λ) (15) ( ) ( ) It is simple to show that the optimality conditions of P are equivalent to those of E . ( E(µ) ) , is formed by using the log barrier function to force s to lie in the positive k orthant. Minimize f (x)−µln(s) 0 k over x∈Rn s∈Rm ( E(µ) ) k subject to: G(x)+s =0 5

It is not necessary to have an explicit non-negativity constraint for s since that will necessarily be the case at a solution. The lagrangian for ( E(µ) ) is given by (16) and k followed by the corresponding optimality conditions: L(x,s;λ)= f (x)−µln(s)+λ' ( G(x)+s ) (16) 0 k ∇f +∇G(x)'λ=0 (17) 0 −µS−1e+λ=0 ⇒ Sλ=µe (18) k k G(x)+s =0 (19) s≥0 (20) ( ) where S =diag(s). (18) and (19) imply that s =−G(x) and λ=µ −G(x) −1 e. Then k (7) is equivalent to (17); and (8) is equivalent to (19) combined with (20). Therefore ( P(µ) ) and ( E(µ) ) are equivalent problems. Hence, one can work with the k k equilibrium form of the problem and be assured that the theoretical results corresponding to convex programming and barrier problems will apply. Most applications of the interior-point algorithm will not strictly follow the central path to the solution but instead reduce µ by a variable amount at each iteration. k This was true in obtaining the results of Berger (2001). For that paper, the variables also ( ) were not restricted to lie in the strictly feasible region of P for all iterations. Instead, using formulation ( E(µ) ) , only the restrictions that x, s, and λ are positive are enforced k at each iteration. ( G(x)+s )→0 with s >0, but the solution at termination must be checked to make sure that G(x)≤0 since s can converge to zero along with G(x) and 6

therefore the solution at termination of the algorithm could lie just outside the boundary of feasible set (by the value of the tolerance for the optimality conditions). For many applications, the level of this discrepancy is of little practical value (such as in the vintage capital case). The iterations use the Newton method and it remains to show that the method is defined for each iteration. We will do this by examining the Jacobian of the optimality conditions of ( E(µ) ) . In certain cases it is preferable, for numerical stability reasons, to k multiply (17) by a positive definite diagonal matrix D. Without loss of generality (i.e. D=I), (17) can be rewritten as: D (∇f +∇G(x)'λ)=0 (21) 0 The Jacobian of the optimality functions can be written as: M M 0 1 2   J = M 0 I (22)  3   0 S Λ   where ( ) M = D ∇2f (x)+∇2G(x)'λ, M ∈Rn×n (23) 1 0 1 M = D (∇G(x)' ) , M ∈Rn×m (24) 2 2 M =∇G(x), M ∈Rm×n (25) 3 3 with ∇2G(x)'λ=λ∇2g (x)+λ∇2g (x)+…+λ∇2g (x) 1 1 2 2 m m The Newton direction is d =−J−1g where g , the vector of optimality k k k k functions, and J are both evaluated at the kth iteration values of x, s, and λ. From the k 7

definition of d , there is a unique Newton direction for each iteration if J is nonsingular. k J is nonsingular if only if there is a unique solution to the following system: M M 0u a 1 2      M 0 I v = b (26)  3      0 S Λ  w   c       From (26), v and w can be solved in terms of u by the following method: w=b−M u (27) 3 v=S−1(c−Λw) (28) Substituting (27) and (28) into (26): ( ) M S−1ΛM +M u =a−M S−1(c−Λb) (29) 2 3 1 2 From (29), the Newton step has a unique solution if and only if the matrix ( ) M S−1ΛM +M is nonsingular. It will be shown that under certain conditions 2 3 1 ( ) M S−1ΛM +M is positive definite and therefore nonsingular. 2 3 1 Proposition 1. Let x∈Rn (non-negative orthant) and ( P ) be a convex programming + problem. Then the Newton direction has a unique solution. Proof. As shown above, the Newton direction has a unique solution if and only if ( M S−1ΛM +M ) is non-singular. Suppose x∈Rn and ( P ) is a convex programming 2 3 1 + problem. First will be shown that M S−1ΛM is positive definite. Then it will be shown 2 3 that M is positive semi-definite. The sum of a positive definite and positive semi- 1 8

definite matrix is positive definite and non-singular. Hence, this is sufficient for the proof. From (24) and (25), M S−1ΛM = DM 'S−1ΛM = DH'H , where D is positive 2 3 3 3 ( )1 ( )1 definite and H = S−1Λ 2 M = S−1Λ 2∇G(x). 3 With no abstract constraints in ( P ) , x∈Rn implies that ∇G(x)∈Rm×n has rank + n, with m≥n. Then nul ( G(x) )=0. Since at each iteration S−1 and Λ are positive definite diagonal matrices, this implies that nul ( H )=0. It is obvious that ( Hw ) 'Hw≥0, and nul ( H )=0 implies that ( Hw ) 'Hw=0 if and only if w=0. Let Q= H'H . Then w'H'Hw=w'Qw. Since w'Qw>0 ∀w≠0, Q is positive definite. Since D is a diagonal positive definite matrix, DQ= DH'H = DM 'S−1ΛM =M S−1ΛM is positive 3 3 2 3 definite. ( ) M = D ∇2f (x)+∇2G(x)'λ . So, M is positive semi-definite if and only if 1 0 1 ∇2f (x)+∇2G(x)'λ is positive semi-definite. If ( P ) is a convex programming problem 0 then ∇2f (x), the Hessian of the objective function, is positive semi-definite. Similarly, 0 ∇2G(x)'λ=λ∇2g (x)+λ∇2g (x)+…+λ∇2g (x), and if ( P ) is a convex 1 1 2 2 m m programming problem, then each g is convex and ∇2g is positive semi-definite. i i Therefore, by the fact that a sum of positively weighted positive semi-definite matrices is positive semi-definite, ∇2G(x)'λ is positive semi-definite. By the same fact, 9

∇2f (x)+∇2G(x)'λ is positive semi-definite. Therefore, since D is a diagonal positive 0 definite matrix, M is positive semi-definite. (cid:44) 1 The method of reducing µ in Berger (2001) was adapted from Wright (1997). k At each iteration, µ is defined to be the average value of the complementarity k conditions (i.e. µ =(λ's)/m). The initial guess for λ and s give µ. A parameter k 0 σ∈(0,1) determines by what fraction the next iteration should attempt to decrease µ k (µ =σµ) . There will be a speed tradeoff in the choice of σ. A lower value of σ may k+1 k allow for fewer iterations if at each iteration the step stays close to the central path. However, too small a value of σ may cause the step to diverge far enough from the central path that the Newton method does not converge well on subsequent iterations. 3. Non-interior Path Following Algorithm Solutions to the putty-clay model in Berger (2001) were obtained using a noninterior path-following algorithm based on Burke and Xu (2000). Burke and Xu, however, present a non-interior path following algorithm for a linear complementarity problem, whereas the putty-clay model is a general non-linear problem. As such a problem, multiple local solutions are possible. The following is the format of a general non-linear optimization problem (except for the abstract constraint on x(cid:4)): 10

Minimize f (x(cid:4)) 0 over x(cid:4)∈Rn (non-negative orthant) + subject to: G(x(cid:4))≤0 ( ) P where: gnl f :Rn (cid:54) R 0 G:Rn (cid:54) Rm ( ) G(x(cid:4))= g (x(cid:4)),g (x(cid:4)),…,g (x(cid:4)) ' with g :Rn (cid:54) R 1 2 m i ( ) The following is the lagrangian corresponding to P : gnl L(x(cid:4);λ)= f (x(cid:4))+λ'G(x(cid:4)) (30) 0 The first order optimality conditions are: ∂L =∇f (x(cid:4))+∇G(x(cid:4))'λ≥0 (31) ∂x(cid:4) 0 ∂L =G(x(cid:4))≤0 (32) ∂λ  ∂L ∂L Let x=(x(cid:4),λ)'. Let M(x)= ,− . Then the Karush-Kuhn-Tucker (KKT)    ∂x(cid:4) ∂λ  conditions are equivalent to and are satisfied for an (x,v) combination of variables if: M(x)−v=0 (33) x≥0 (34) v≥0 (35) Vx=0 (36) where V =diag(v). 11

For non-interior path-following algorithms, penalty functions are not used. Instead, the scalar µ ≥0 is a parameter in a function φ(x,v,µ):R(2n+2m+1) (cid:54) R(n+m), k k [ ] where φ(x,v,µ)= φ(x ,v ,µ),…,φ(x ,v ,µ) , v = M(x) and: k 1 1 k n+m n+m k φ(x ,v ,µ)= x +v − (x −v )2 +4µ2 (37) i i k i i i i k ( ) The optimality conditions of P will hold if and only if φ=0 and µ =0. The gnl k central path is defined as the (x,v) combinations such that φ=0 for a sequence of nonnegative µ. The Newton method is once again used to find a stationary point of the k optimization problem. Reductions in µ are done in a predictor-corrector step k framework. A predictor step, which attempts a reduction of µ to zero, is first used in k the iteration. A corrector step is then utilized in which only a partial reduction in µ is k attempted. It is called a “corrector” step because it is used to guide values closer to the central path. Note 1. φ(x ,v ,µ)=0 if and only if x ≥0, v ≥0 and xv = µ2 i i k i i i i k Proof. Let φ(x ,v ,µ)=0. i i k Then x +v = (x −v )2 +4µ2 . i i i i k Squaring both sides: x2 +v2 +2xv =(x −v )2 +4µ2 ⇒ 4xv =4µ2 ⇒ i i i i i i k i i k xv = µ2 (38) i i k By (38), µ ≠0 implies x and v are non-zero and the same sign. k i i But x ,v <0 implies that φ(x ,v ,µ)<0. Therefore: i i i i k 12

µ ≠0 ⇒ x ,v >0 (39) k i i By (38), µ =0 implies x and/or v equals zero. By (37), if φ=0 and either x or v is k i i i i zero, then the other is non-negative. Therefore x≥0, v≥0 and xv = µ2. It is trivial to k see that if xv =µ2, then φ(x ,v ,µ)=0. (cid:44) i i k i i k Let:  M(x)−v   φ(x ,v ,µ) 1 1 k     F(x,v,µ)= φ(x,v,µ) with φ(x,v,µ)= (cid:35) (40) k  k  k    µ  φ(x ,v ,µ)     k n n k Note 1 shows that (33)-(36) hold if and only if F =0. Newton's method leads to solving (41) for the predictor step:  ∆x    F(x,v,µ)+∇F(x,v,µ) ∆v =0 (41) k k   ∆µ   k where:  ∇M(x) −I 0    ∇F(x,v,µ)= ∇ φ(x,v,µ) ∇ φ(x,v,µ) ∇ φ(x,v,µ) k  x k v k µ k  k  0 0 1    Then solve the following for the corrector step:  ∆x   0      F(x,v,µ)+∇F(x,v,µ) ∆v = 0 (42) k k     ∆µ (1−σ)µ     k k Proposition 2. The Newton direction has a unique solution if and only if ∇ φ(x,v,µ)+∇ φ(x,v,µ)∇M(x) is nonsingular. This will occur if ∇M(x) is positive x k v k 13

( ) semi-definite. ∇M(x) is positive semi-definite if P is a convex programming gnl problem. Proof. The Newton direction has a unique solution if and only if there is a unique solution to the following system of equations:  ∇M(x) −I 0 u a      ∇ φ(x,v,µ) ∇ φ(x,v,µ) ∇ φ(x,v,µ) v = b (43)  x k v k µ k     k      0 0 1  w c      By (43): w=c (44) v=∇Mu−a (45) By (43), (44) and (45): (∇ φ+∇ φ∇M ) u =b+∇ φa−∇ φc (46) x v v µ k Hence, there is a unique Newton direction if (∇ φ+∇ φ∇M ) is non-singular. x v x −v ∇ φ(x ,v ,µ)=1− i i (47) x i i i k (x −v )2 +4µ2 i i k x −v ∇ φ(x ,v ,µ)=1+ i i (48) v i i i k (x −v )2 +4µ2 i i k ( ) ∇ φ(x,v,µ)=diag ∇ φ(x ,v ,µ),…,∇ φ(x ,v ,µ) (49) x k x 1 1 k x n n k 1 n ( ) ∇ φ(x,v,µ)=diag ∇ φ(x ,v ,µ),…,∇ φ(x ,v ,µ) (50) v k v 1 1 k v n n k 1 n It can be seen from (47)-(50) that µ >0 implies: k ∇ φ(x ,v ,µ)>0 and ∇ φ(x,v,µ) is positive definite (51) x i i k x k i 14

∇ φ(x ,v ,µ)>0 and ∇ φ(x,v,µ) is positive definite (52) v i i k v k i Therefore, (∇ φ+∇ φ∇M ) is non-singular if ∇M is positive semi-definite. x v  M M  ∇M = 1 2 (53)   −M ' 0   2 where: M =∇2f (x)+∇2G(x)'λ 1 0 M =∇G(x)' 2 However: [ ] u u v ∇M =u'M u (54)   v 1 By (54), ∇M is positive semi-definite if and only if M is positive semi-definite. 1 ( ) M is positive semi-definite if P is a convex programming problem (see proof of 1 gnl Proposition 1). (cid:44) 4. General Model Descriptions Berger (2001) found solutions to discrete time, finite horizon, parameterized, optimization versions of the three main types of vintage capital models. These models are referred to as putty-putty, clay-clay, and putty-clay models and vary according to the substitutability of factors. Putty refers to the ability to pair labor and capital in any production ratio. Clay refers to a fixed capital-labor ratio. The first word in each pair refers to substitutability at the time of installation of the capital; and the second word refers to the substitutability for all time after the installation. For example, the term 15

putty-clay refers to a model in which the capital-labor ratio for capital of a particular vintage may be chosen at the time of installation, but for all time thereafter labor must be used with that vintage of capital according to the chosen proportion. Each of the three models is based on the basic consumption-savings model. A central planner of the economy chooses a time path of consumption with the objective of maximizing the sum of discounted utility of the population, which is solely a function of consumption. Utility is represented by a finite, constant-elasticity-of-intertemporalsubstitution utility function: 1 ( ) U(C )= C1−γ−1 γ≠1 t 1−γ t (55) =ln(C ) γ=1 t 1 The elasticity of intertemporal substitution equals − , and for fixed C this is a γ t continuous function of γ. Output not consumed each period is invested and becomes capital that can be used to increase production in future periods. Each period in the putty-putty and clayclay models a static optimization takes place in which the central planner maximizes total output given the capital stock, labor stock, and given parameters. For non-vintage models this allocation is trivial. However, for vintage capital models this is a non-trivial problem in which scarce labor must be assigned to the different vintages of capital. The three models are based on a Cobb-Douglas production function: f (N ,K )=d A N1−αKα (56) v v v t v v v 16

where v represents the vintage, d is the disembodied technology parameter, A is the t v embodied technology parameter, and 0<α<1. Total output for period t is: V+t−1 f =d ∑ A N1−αKα (57) t t v v v v=1 where V is the number of vintages available in the initial period.. One advantage of the putty-putty model with this production function is that the static optimization problem of maximizing output giving existing resources can be written as a convex programming problem with a strictly convex objective. Therefore, since the problem has a nonempty feasible set, there exists a unique global solution to maximize output each period. This solution, the indirect production function, can be written as a function of time t parameters: f (N ,Q )= N1−αQ1−α (58) t t t t t where Q is an aggregated capital stock which is a weighted sum of the stocks of each t vintage (Appendix 1). The clay-clay and putty-clay models maximize the output of the same Cobb-Douglas function, but subject to a fixed capital-labor ratio, r . Therefore, v output corresponding to each vintage each period is: f (N ,K )=min(a N ,d K ) (59) tv tv tv tv tv tv tv where a =d Arα and d =d Arα−1 measure the productivity of labor and capital, tv t v v tv t v v respectively.2 In the clay-clay model, a and d are parameters, whereas in the puttytv tv 2 d and d represent distinct concepts. d is disembodied technology corresponding to time period t t tv t ( ) whereas d is the productivity of capital and includes both embodied A and disembodied technology. tv v 17

clay model the central planner can choose r thereby affecting the levels of a and d . v tv tv Unlike the putty-putty model, there are no closed-form solutions for the indirect production functions in the clay-clay and putty-clay models. The putty-putty and clay-clay models can be written as convex programming problems with objective functions that are strictly convex in consumption. Therefore there is a unique global solution for consumption. Consumption determines the time path of output in these two models and therefore there is also a unique global solution to the time path of output. The putty-clay model, however, cannot be written as a convex programming problem. While the objective functions in all three models are the same, the introduction of the capital-labor ratio choice variables in the putty-clay model leads to a non-convex feasible set. Therefore the same claim as to globalness and uniqueness of solutions cannot be made for the putty-clay model. 5. Initialization of Algorithms Initialization of the models in Berger (2001) is simple. For the putty-clay model, first choose arbitrary capital-labor ratios for each vintage. Consequently for all three models, given a level of capital, labor, and parameter values in any given period, one can determine the maximum level of output. For the clay-clay and putty-clay models, maximizing output involves an algorithm which first allocates labor to the most productive capital until the capital-labor ratio is reached, and then allocating labor to the second most productive capital, continuing down the capital spectrum until either labor or capital run out. For the initial period, output can be set to a fraction of maximum output, 18

and consumption to a fraction of the chosen output. Output and consumption determine the following period's capital and therefore the procedure can be repeated for all of the periods. The slack variables should be set such that the optimality conditions are nearly satisfied given the initialization of the other primal variables. There is often no good intuition for the initial guess of the dual variables except for which will be zero (set these positive, but close to zero) and the relative scale of the others. 6. Putty-Putty Model The putty-putty model examined in Berger (2001) is a discrete time version of the Phelps (1962) model with a finite time horizon. The Phelps model uses a constant rate of savings, and is divorced from dynamic optimization unlike this model. Technological growth is exogenous. The following model already incorporates the static production maximization, so that the production function given is the indirect production function (Appendix). The formal optimization problem is: T 1 ( ) Minimize −∑βt−1 C1−γ−1 1−γ t t=1 over C ,Y,Q for t =1…T t t t subject to: −C ≤0 ( ) t P pp Y −d N1−αQα≤0 t t t t C −Y ≤0 t t Q −Q≤0 1 Q −Q + A 1 α(Y −C )≤0 for t =1…T −1 t+1 t t t t C is the chosen level of consumption in period t. Y is the chosen level of output in t t period t. Q is the level of the aggregate capital stock at time t. N is the given labor t t 19

stock each period. Q is the given level of aggregate capital in the initial period. d is t the given level of disembodied technology at time t. A is the given level of embodied t technology corresponding to the vintage produced at time t. T is the given total number of time periods. ( ) ( ) P is a convex programming problem. It can be shown that P is equivalent pp pp to a problem written in standard economic form, a maximization with the last four constraints written as equalities. The inequality constraints can be written as equalities because there are no transaction costs associated with the use of capital for production, or in carrying over capital and investment from previous periods into subsequent periods. Therefore, at a solution, the production and dynamic constraints will hold with equality. Results for the putty-putty model were obtained for T =200 so that the barrier form of the problem, as actually solved, had 1,800 primal variables (including 1,200 slack variables) and 1,200 dual variables (the number of constraints, including 600 nonnegativity constraints). The value of the constant structure of the Jacobian of the optimality conditions to the equilibrium barrier formulation (22) is apparent in that the Newton direction can be found by inverting a 600 by 600 matrix, rather than the entire 3,000 by 3,0000 Jacobian. 7. Clay-Clay Model The clay-clay model, also known as the Leontief or fixed factor production model, is a discrete time variation of Solow, Tobin, Weizsäcker, Yaari (1966) with a finite time horizon. Unlike in the putty-putty model, there is no aggregate capital stock 20

because the embedded production maximization each period does not have a closed form solution. Thus, each vintage must be tracked separately. The capital stock in period t is a vector of V+t-1 elements, where V is the number of vintages available in the initial period ( t =1 ) . The savings of period t becomes the quantity of capital of vintage V+t. The production function, Y =min(d K ,a N ), where N is labor in period t tv tv tv tv tv tv assigned to capital of vintage v, can be equivalently written as two inequality constraints. Y ≤a N (60) tv tv tv Y ≤d K (61) tv tv tv The equivalency of the “min” function to the two inequality constraints arises because production will be maximized each period, for each vintage, given levels of labor and capital assigned to each vintage, and therefore (60) and/or (61) must hold with equality at a solution. What does it mean if (60) does not hold strictly for time period t and vintage v, at a solution? It means that there is unproductive labor applied to vintage v. This can only be the case at a solution if there are no other vintages in that period to which the labor could be applied and be productive. If this is the case, then N , the labor available in t period t, could be increased arbitrarily without changing the optimal time path of consumption and hence the optimal value of the optimization problem. Since the labor constraint each period, ∑V+t−1 N ≤ N , is only relevant when (60) is binding, we can v=1 tv t Y substitute tv for N in the labor constraint, eliminating N as a choice variable. This a tv tv tv yields the clay-clay model (a convex programming problem) solved in Berger (2001). 21

T 1 ( ) Minimize −∑βt−1 C1−γ−1 1−γ t t=1 over C ,Y ,K for t =1…T for v=1…V +t−1 t tv tv subject to: C ≥0 t Y ≥0 tv Y −d K ≤0 tv tv tv ∑V+t−1 Y tv −N ≤0 ( ) v=1 a t P tv cc C −∑V+t−1 Y ≤0 t v=1 tv K = K 1 1   0 0     (cid:37) (cid:35) (cid:35) C  K =  K +  t for t =1…T −1 t+1  1 t  0 0  Y    t     0 (cid:34) 0  −1 1 In Berger (2001), results were obtained for T =45 and with no depreciation of capital. Depreciation can be easily incorporated into the model by changing the diagonal ones in the capital accumulation equation to fractional values. Consumption and output converge quickly to their equilibrium growth path in the clay-clay model, so 45 periods is ample time to observe the behavior of the system. Tracking each vintage each period greatly increases the number of variables in the model. The number of variables can be decreased somewhat by writing the problem in reduced form, substituting for K with a function of C, Y, and K. The reduced form, equilibrium barrier formulation, with T=45 and V =2, has 3,420 primal variables (including 2,295 slack variables), and 2,295 dual variables. 22

8. Putty-Clay The putty-clay model is similar to the clay-clay model. The primary difference is that r , the capital-labor ratio for each vintage v, is a choice variable. v T 1 ( ) Maximize ∑βt−1 C1−γ−1 1−γ t t=1 over C ,Y ,K ,r for t =1…T for v=1…V +t−1 t tv tv v subject to: C ≥0 t Y ≥0 tv A Y − tv K ≤0 tv r1−α tv v ∑V+t−1 Y tv −N ≤0 ( P ) v=1 A rα t pc tv v C −∑V+t−1 Y ≤0 t v=1 tv K = K 1 1   0 0     (cid:37) (cid:35) (cid:35) C  K =  K +  t for t =1…T −1 t+1  1 t  0 0  Y    t     0 (cid:34) 0  −1 1 Unlike the putty-putty and clay-clay models, the putty-clay model can not be written as a convex programming problem (the feasible set of is not a convex set). Therefore, there can be multiple solutions to the model, and they need not be global. Berger (2001) does in fact find multiple solutions to the putty-clay problem. Finding the results for the putty-clay model is more difficult than for the previous models. Because the problem is not a convex programming problem, there is no guarantee that the non-interior path-following algorithm will converge. The algorithm is highly sensitive to the initial guess, and often gets mired away from a solution. In order 23

to find solutions, it was necessary to use a bootstrapping method in which a solution was found for a low value of T, and this solution was used as the basis for an initial guess to a slightly longer problem. This procedure was repeated until a solution to the desired length problem was found. Also, because r is to a fractional power, it must be restricted to positive values in the algorithm. However, limiting the step size to keep the capitallabor ratio positive caused the program to stop progressing as the step size rapidly approached zero. So the primal variable r was replaced with s where r =re −s v , and r v v v is a user defined algorithm parameter. Two types of solutions to the model were found. One type has r constant from v an early time period through the terminal period. The results suggest that there may be a local minimum for every constant level of r up to some threshold level. To obtain v solutions of this type, it is possible to obtain a solution to a shorter horizon model, e.g. T=10, find the constant r , and then treat the problem as a clay-clay model using the v obtained capital-labor ratios. When a solution to the 45 period clay-clay model is found, it can be used as an initial guess for the putty-clay model (re-running the non-interior point method is necessary in order to get the correct early period values for the 45 period putty-clay model). The non-interior method then only takes a few iterations to arrive at a solution. When r is set sufficiently high, another solution type is obtained. This solution has r increasing over time. This solution type also required bootstrapping to attain v solutions to longer time period problems. Unlike for the constant capital-labor ratio solution however, there is no way to use the more stable clay-clay model to help obtain 24

the full 45 period solution. The putty-clay model has more variables than the clay-clay model because of the addition of 46 ( V +T −1 ) capital-labor ratio variables. 9. Conclusion In Berger (2001) numerical results were obtained for the three basic utilitymaximizing growth models with vintage capital. The fact that complete time paths of the key economic variables were obtained for these models demonstrates the potential of the algorithms discussed in this paper. Given the specified parameter values, the results for the putty-putty and clay-clay model are definitive since the models can be written as convex programming problems with strictly convex objectives, and therefore have unique global solutions. The putty-clay model, however, cannot be written as a convex programming problem and multiple local solutions were found. Numerical results had not been previously obtained for optimization, discrete time versions of these models. It was also shown that for convex programming problems the Newton method will be well defined for both the interior-point and non-interior point path-following algorithms. The results for the putty-clay model demonstrate both the difficulty, but also the possibility of attaining results for general non-linear models. The non-interior path following algorithm was unstable when applied to this model, but using bootstrapping techniques, solutions were found for the model. The putty-clay model is a large scale general non-linear model with several complications including fractional exponents. For smaller scale or convex programming problems, the non-interior method may be equally or more effective than the interior point method. 25

Models in economics are becoming increasingly complex. It will not be possible to find closed form theoretical solutions to many of these models. Analyzing models which more accurately reflect the world in which we live will require the use of algorithms such as those outlined in this paper. 26

Appendix 1: Discrete-Time, Vintage Capital Cobb-Douglas Production Maximization Each instant, a static optimization takes place in which given the stocks of technology, labor, and capital, labor must be allocated to the different vintages of capital in order to maximize production. The formal problem is as follows: V Maximize Y =∑A N1−αKα 0<α<1 v v v v=1 over N v subject to: (62) N ≥0 v V ∑N ≤ N v v=1 V is the newest vintage available. A is the level of technology associated with vintage v. v N , the sole choice variable, is the quantity of labor assigned to capital of vintage v. N is v the total amount of labor available. K is the amount of capital of vintage v available. v It is obvious that a solution will have N >0 because the marginal product of v labor of each vintage goes to infinity as N goes to zero. Also, the total labor constraint v will hold with equality since the objective function rewards the use of labor with each vintage and there is no cost of using labor. The lagrangian for this problem is: V  V  L(N ;λ)=∑A N1−αKα+λ N −∑N (63)   v v v v v v=1  v=1  A representative first order condition is: ∂L =(1−α)A N−αKα−λ=0 (64) ∂N v v v v 27

An equality can be used instead of an inequality, because N >0 at a solution. (64) v states that the marginal product of labor of each vintage must be equal to the shadow value of labor. Labor, being homogeneous, has a unique shadow value regardless of the vintage to which it is applied. Setting two first order conditions for different vintages equal to each other we obtain: ANα−1Kα= A Nα−1Kα (65) i i i j j j Solving for N i 1  A α K  N =  i   i N (66) i  A   K  j  j   j  This equation must hold true for every i given a fixed A , N , and K . Therefore, j j j summing over the N : i V N V 1 N =∑N = j ∑AαK (67) i 1 i i i=1 AαK i=1 j j Solving for N in terms of the given parameters we obtain a closed form solution for j each choice variable: 1 NAαK N* = j j (68) j V 1 ∑AαK i i i=1 Let: V 1 Q=∑AαK (69) i i i=1 28

Q can be viewed as an aggregated capital stock with the weights of each vintage a function of the technology inherent in each vintage. Substituting N* into the direct production function, we obtain the indirect j production function, which is a function only of given parameters. 1−α  1  Y*(A,N,K,σ ,σ ,σ )=∑ V A Kα NA v αK v  A N K v v  Q  v=1     V 1 = N1−αQα−1∑AαK (70) v v v=1 = N1−αQα 29

References Berger, Brett D. (2001). “Convergence in Neoclassical Vintage Capital Growth Models.” International Finance Discussion Papers, Board of Governors of the Federal Reserve System, 713, Novemeber 2001. Bertsekas, Dimitri. P. (1999). Nonlinear Programming, Second Edition. Athena Scientific. Burke, J., Xu, S. (2000). “A Non-Interior Predictor-Corrector Path Following Algorithm for the Monotone Linear Complementarity Problem.” Mathematical Programming 87, 113-130. Greenwood, J., Z. Hercowitz, and P. Krussel. (1997). “Long-Run Implications of Investment-Specific Technological Change.” The American Economic Review 87, 342-362. Phelps, Edmund S. (1962). “The New View of Investment: A Neoclassical Analysis.” Quarterly Journal of Economics 76, 548-567. Rockafellar, R. T. (1970). Convex Analysis. Princeton University Press. Solow, Robert M. (1956). “A Contribution to the Theory of Economic Growth.” Quarterly Journal of Economics 70, 65-94. Solow, Robert M. (1959). “Investment and Technical Progress.” Stanford Symposium on Mathematical Methods in the Social Sciences, 89-104. Solow, R., J. Tobin, C. von Weizsäcker, and M. Yaari. (1966). “Neoclassical Growth with Fixed Factor Proportions.” Review of Economic Studies 33, 79-155. Wright, Stephen J. (1997). Primal-Dual Interior-Point Methods. Society for Industrial and Applied Mathematics. 30