Stochastic Newton-like methods for computing equilibria in general equilibrium models

Stochastic Newton-like methods for computing equilibria in general equilibrium models

Nataša Krejić^{I, *}; Zorana Lužanin^{I, *}; Zoran Ovcin^II

^IDepartment of Mathematics and Informatics, Faculty of Science University of Novi Sad, Trg Dositeja Obradovica 4, 21000 Novi Sad, Serbia, [e-mail] / [e-mail]

^IIFaculty of Technical Sciences, University of Novi Sad, Trg Dositeja Obradovica 6 21000 Novi Sad, Serbia [e-mail]

ABSTRACT

Calculating an equilibrium point in general equilibrium models in many cases reduces to solving a nonlinear system of equations. Taking model parameter values as random variables with a known distribution increases the level of information provided by the model but makes computation of equilibrium points even more challenging. We propose a computationally efficient procedure based on application of the fixed Newton method for a sequence of equilibrium problems generated by simulation of parameters values. The convergence conditions of the method are derived. The numerical results presented are obtained using the neoclassic exchange model and the spatial price equilibrium model. The results show a clear difference in the quality of information obtained by solving a sequence of problems if compared with the single equilibrium problem. At the same time the proposed numerical procedure is affordable.

Mathematical subject classification: 65H10, 91B50.

Key words: nonlinear system of equations, Newton-like method, general equilibrium model.

1 Introduction

The problem we consider in this paper is

where x ∈ Rⁿ is an unknown vector, W is a m-dimensional vector of parameters and F : Rⁿ ®Rⁿ.

The most simple case of this type of problems arises when W ∈ W" is constant. Many economic models can be written in the form of (1) with W being a vector of parameters that are estimated. Estimation of the parameter values based on real data is a nontrivial task given that W is most likely a random vector. As F is often highly nonlinear and n can be large, solving (1) becomes a challenging issue so W is often replaced by a constant real valued vector, say by the expected value. Such modeling simplification makes solving (1) easier but also introduces additional uncertainty into the model. An important part of the information available from the real data might be lost with the constant value parameters approach and hence the model properties might deteriorate. Therefore the question of finding a computationally affordable procedure for solving (1) with W being more general than just a constant vector naturally arises.

Let us assume that W in (1) is estimated in such a way that its probability distribution is known. So W is a random vector on some probability space (Ω, F, P)• This kind of parameter estimation provides significantly larger amount of information than constant value and thus implies better properties of the model. In that case (1) becomes a stochastic problem and its solution can be a random variable. Analytical solution of (1) is not available in general and thus an affordable numerical procedure is of great interest.

Solving stochastic nonlinear systems of equations numerically is a difficult problem. A number of numerical procedures is suggested and analyzed in literature if the parameter vector W is a constant, [10, 6, 7]. The proposed numerical procedures are in fact designed for solving the deterministic nonlinear problems of the form (1) and they are mainly based on the interior-point or homotopy methods. In this paper we propose a different approach aiming to preserve as much information about the vector W as possible but keeping the computational cost at a reasonable level. If W is a random vector with known distribution then for a given sample {w¹, ..., w^N} one could solve the sequence of problems

getting a sequence of solutions (x*¹, ..., x*^N) ∈ R^nxN. If N is large enough then (x *¹, ..., x *^N) might provide a good idea of the true properties of the random vector x* which is a solution of (1). However solving the sequence of nonlinear systems (2) for large values of N might be computationally unaf-fordable. Therefore we are proposing a numerical procedure based on the fixed Newton idea that generates a sequence of approximate solutions of (2) at an affordable cost and thus provides significantly more information that a solution of (1) for a constant value (i.e. the expected value) of the vector W.

Section 2 of this paper contains the algorithm we are proposing and its convergence analysis. Two equilibrium models of the form (1) are described in Section 3 and the numerical results are presented in Section 4.

2 The Algorithm

Let us assume that W is a random vector with a known distribution and the sample (w¹, ..., w^N) is drawn. We are interested in solving the sequence of problems

and thus obtaining the sequence of optimal solutions

The final goal is to obtain an approximation of x *, which is a random variable that solves (1). So for N large enough we can use the sequence (x*¹, ..., x*^N) to estimate the distribution of each component of x^*.

The problem (1) is considered under a set of standard assumptions. Let D ⊂ Rⁿ such that the following assumptions are satisfied.

(A1) For all w ∈ Ω there exists x^w* ∈ D such that F(x^w*, w) = 0.

(A2) For all w ∈ Ω the Jacobian matrix F'(x^w*, w) is nonsingular.

(A3) For all x, y ∈ D and w ∈ Ω there exists y > 0 such that

(A4) For all w_t ,wj ∈ Ω and x ∈ D there exists y_W > 0 such that

One obvious possibility is to use the Newton method to solve each of the systems (2). The Newton method is a reliable choice given that it yields quadratic convergence for a good initial approximation. So given x^i,0 the sequence

will converge to the solution xⁱ* if F(x, wⁱ) satisfies the standard assumptions A1-A3. However such procedure would be quite expensive for large N if the dimension n of the problems (2) is even a modest number. On the other hand the Jacobian matrices F'(x^i,k, wⁱ) could be expected to have very similar structure for different sample realizations w^l. Therefore a numerical procedure should take advantage of such similarity. A natural way to employ the favorable properties of the sequence of problems we are solving would be consider the procedure based on the fixed Newton approach. The fixed Newton method calculates the Jacobian only once and proceeds with iterations using the initial Jacobian during the whole process. The method is convergent if the initial approximation is close enough to the solution of the considered nonlinear system. So we will consider the following algorithm for solving (2).

Algorithm FNM

StepO: Let x⁰ ∈ Rⁿ and (w¹,...,w^N) be given. Calculate the mean value

and solve F (x ,w) = 0 by the Newton method. Denote the solution by x ^w and define A = F'(x ^w, w).

Step 1: For i = 1, N

Step 1i: Set x^i,0 = x^w

Step 2i: Repeat until convergence

Step 3i: Set xⁱ* = x^i,k

Remark. The suggested matrix A = F'(x^w, w) is one possible choice for the constant matrix used in Step 2i. The main computational burden in (3) comes from the calculation of the Jacobian at each iteration and for each sample realization wⁱ, followed by solving the corresponding linear systems that define the iteration increments s^i,k. Therefore the main advantage of the suggested procedure is that all steps are obtained by solving a system of linear equation with the same matrix A. Naturally A should be a reasonably good approximation of the true Jacobians at all iterations. Furthermore A can be factored only once and used in each Step 2i thus decreasing the linear algebra costs significantly. So calculating A with the sample mean w appears natural given that such A should be reasonably close to F'(x, wⁱ) for all wⁱ and x close to xⁱ*. But there are several other possibilities, for example A = F'(x⁰, w^j), with j being an arbitrary number 1 < j < N. The numerical testing we performed did not indicate that the algorithm exhibits large differences depending on the choice of w or any w^j when defining A. More details are reported in Section 4.

The algorithm above will generate the sequence (x¹*, ..., x^N*) if Step 2i is well defined i.e. if every single iterative procedure for i = 1, ..., N ends up with some x^i,k that satisfies the convergence conditions and could be taken as x^i* in Step 3i. One of the advantages of the above algorithm is that it could be easily used in parallel environment [2].

We will prove that the standard assumptions (A1)-(A4) imply that Step 2i is well defined. The convergence theorem presented below is based on the Banach contraction principle. The convergence condition might be seen as a generalization of two neighbourhod theorem as it consists of an inequality connecting the distance between the initial point and the solution and the variance of the model parameters W. Such condition might seem rather strong as common in local convergence statements. But the numerical results given in Section 4 will demonstrate the effectiveness of the algorithm despite the strong convergence condition.

Let us denote B = B(x, δ) = {y ∈ Rⁿ : ||x — y || < δ}. Assuming that Assumptions A1-A3 are satisfied one can easily see that the Newton method for solving F(x, w) = 0 is locally convergent i.e. there exists δ₁> 0 such that for x⁰∈ B(x^w*,δ₁) Step 0 is well defined and generates a convergent sequence. Let x ^w be the value obtained after s iterations of the Newton method in Step 0 and let s be large enough such that x^w is close enough to x^w*. Then the matrix F'(x^w, w) is nonsingular and there exists M > 0 such that ||F'(x^w,w)^-1||= M. The convergence of FNM Algorithm is stated in the theorem below.

Theorem 1.Let F satisfy assumptions A1-A4. If for s_W = diam(Q) and M = ||F'(x^w,w)^-1|| there exists δ > 0 such that α = M(yδ + y_Wε_W) < 1 and || F '(x^w, w )^-1 F (x^w, wⁱ )|| < δ(1 — α) then for every wⁱ ∈ Ω sequence {x^wi,k }^∞_k=0, i = 1, 2, ..., N defined by Step 1 of Algorithm FNM converges linearlyto solution F(x, wⁱ) = 0.

Proof. Let G(x) = x — F'(x^w, w)^—1F(x, wⁱ). Then, for x ∈ B

From the mean-value theorem we obtain

Therefore, G is contraction on B. Further more

and G maps B into itself. Therefore the contraction mapping theorem applies and there exists a unique fixed point in B and that point is the solution of F(x,wⁱ) = 0 in B.

3 Equilibrium models

Neoclassical exchange economy

The model of neoclassical exchange economy is one of the oldest and most classical economic models. It was developed by Walras, later extended and formalized by Arrow and Debreu, [7, 11] and since then extensively analyzed in the literature, [11, 3, 9]. One survey of the general equilibrium models is presented in [6]. The most common approaches for solving equilibrium problems include simplex methods based on the constructive proof of the Brouwer fixed point theorem, tâtonnement approach, Newton method, Smale's method, interior-point method (see [6] and references therein). Homotopy method is analyzed in [4].

Let Rⁿ₊ be the set of nonnegative real numbers. The model we will consider assumes that there are n commodities in X c Rⁿ₊ and m economic agents. Each agent has a utility function u^j, and an initial endowment w^j for j = 1, ..., m. Therefore each agent has the budget π • w^j where π ∈ Rⁿ₊ is the price vector of commodities. The choice of utility function is defining the model. We will consider the following three possibilities.

• Cobb-Douglas:

• fixed proportions demand function:

• CES demand function:

In any of these utility functions the set of parameters a₁, ..., a_n, b has to be estimated. The models considered in [6] assume that these values are constant, while we consider a more general case in which these parameters are random variables with known distribution.

All agents are profit maximizers, the goal is to maximize the chosen utility function for each agent i.e. to determine

For the considered utility functions we have the following maximizers

• Cobb-Douglas demand function:

• fixed proportions demand function:

• CES demand function:

If the excess demand function is defined as

then any strictly positive π * such that

is an equilibrium price. Spatial price equilibrium model

The class of market equilibrium problems where supply, demand and transportation cost functions are nonlinear and separable is considered in [10] and the references cited therein. The distinguishing characteristic of spatial price equilibrium models lies in their recognition of the importance of space and transportation costs associated with shipping a commodity from a supply market to a demand market. These models are perfectly competitive partial equilibrium models. The model we consider assumes that there are many producers and consumers involved in the production and consumption, respectively, of one or more commodities.

We will use the following notation.

• Yi - total supply (output) at origin (market) i ∈ I = {1, ..., m}

• Zi - total demand (input) at destination (market) j ∈ J = {1, ..., n}

• x_ij - quantity shipped (traded) from origin i ∈ I to destination j ∈ J

• Si (Yi) - inverse supply function at origin i

• Ij (Zj) - inverse demand function at destination j

• c_ij (x_ij) - transportation (transaction) cost between origin i and destination j

The equilibrium can be described as a triple {(Y_i)_i ∈_I, (Zj) _j ∈_J, (x_ij)_i ∈_Ij ∈ _j} that solves the system of equalities and inequalities

The model is developed assuming that the following assumptions hold.

M1 The supply, demand and transportation cost functions are nonnegative C ¹(R₊) functions. The supply and transportation functions are nondecreasing while the demand price is nonincreasing.

M2

If M1-M2 hold then the equilibrium conditions (5) are the KKT conditions for the minimization problem

Supply, demand and transportation costs are model parameters that need to be estimated while x_ij are the unknowns we want to compute at the equilibrium.

If a shipping quantity x_ij- is zero, the corresponding inequality in (5) can be strict. As in ([10]) we can assume the strict complementarity hypothesis:

in the solution (x*, Y*, Z*). Otherwise we have degeneracy. Most common situation for degeneracy is when the transportation cost functions c_ij (x_ij) are constant. Under the strict complementarity assumption, we can apply some strategy for determining set of indices P with nonzero shipping quantities. That gives us system of equations from (5) for (i, j) ∈ P. However, in our computations it was not needed, as the problems we considered had all x_ij positive (P = I x J).

Denoting by x the components xij, i ∈ I, j ∈ J and by W all model parameters we can write the conditions (5) as F(x, W) = 0.

4 Numerical results

In a Neoclassical exchange economy we find the equilibrium price by solving the system (4) for a price vector π*. As a consequence of the homogeneity of degree zero of the demand and supply functions in the prices, without loss of generality, we can define prices and look for the solution on the positive simplex

In numerical calculations we shall substitute the last equation in the system (4) with the equation

Statistics

In order to see the effects of our Algorithm, we perform several calculations with different parameter values. In all examples the values of parameters are generated as pseudo-random numbers having normal or uniform distribution with given deviation centered about the value that occurs in the standard testing problems. The first phase of the algorithm yields x^w where w is the sample mean. Given that all mappings are nonlinear that point is just an approximation to the solution mean. Although this approximation was relatively close to the sample mean in all tested examples Step 1 of FNM algorithm further improves the approximate solution x*ⁱ for each sample wⁱ and thus generates a sample of solution which allow us more detailed statistical analysis.

For each solution some descriptive statistics are presented. We also perform a statistical test to verify or reject the hypothesis that the solution sample is normally distributed as the parameters samples. The test for normality we use is the Anderson-Darling test ([1]) with significance level α = 5%. The Anderson-Darling statistics is

where N is the sample size, and z_i is the i-th value in the sorted sample. We use the modified Anderson-Darling statistics

whose critical value for significance level α = 0.05 is 0.752, and critical value for α = 0.01 is 1.035.

The values for kurtosis and skewness of a component of the solution sample π *ⁱ, i = 1, ..., N are presented in the Figures. These values describe the shape of the probability density function of the sample component distribution. The formulas for skewness and kurtosis are given by

respectively, where

Kurtosis is 0 and skewness is 3 for normal distribution.

The values of β₁², β₂, histogram, statistical moments m₁ and u_k, the Anderson-Darling statistics and other descriptive statistics data are examples of valuable information on the behavior of the equilibrium prices under Gaussian noise i.e. under the assumption of normal distribution of parameter values.

Example 1

We consider a Neoclassical exchange economy with 2 goods and 2 agents with initial endowments w¹ = (3, 1) and w² = (1, 2). The first agents' utility function is the Cobb-Douglas function The second agents utility function is the fixed proportions utility function

The excess demand functions in this case are the following

The initial value for the Newton method (NM) is taken as π₀ = (0.1, 0.9) for all sample values and the exit criterion for both the Newton method and the fixed Newton method (FNM) ||ξ|| < ∈ = 10^—6.

Taking parameters a₁¹ = 0.4, a₁² = 2, a₂² = 3, we obtain the solution π^* = (0.2087, 0.7913) by the Newton method in 6 iterations.

Next, we make samples of N = 500 values for parameters from normal distribution: a¹ : N(0.4, 0.05), a₁² : N(2.0, 0.05), a²₂ : N(3.0, 0.05). Let W denotes the triple of parameters W = (a¹, a², a²₂)and wⁱ, i = 1, ..., N be the set of the sample values.

We perform both methods (Newton's and FNM) for the parameter sample values obtaining the sample of solutions π*^l, i = 1, ..., N.

In Table 1 we compare the total number of iterations, number of function evaluations and Jacobian evaluations for both NM and FNM with the Jaco-bian A = ξ'(x^w, w) taken at the solution x^w of NM in the sample mean w.

The average number of iterations in each of 500 problems is 5.91 for NM and 9.11 for FNM but the CPU time clearly demonstrate the advantage of FNM approach as it is 0.156 for NM and 0.098 for FNM.

The histogram and descriptive statistics of the solution sample for both components π₁* and π₂* are shown in Figure 1. We get the same values of π₁* and π₂* as they are not independent (π₁ + π₂ = 1).

Problem (4) is a continuous mapping from the parameter set into the solution set. Thus the obtained solution sequence is a sample from a random variable. If it is a sample from a random variable, from the histogram we could say that the probability distribution of the underlying distribution is bell shaped, slightly slanted (β₁² = 0.1090, should be zero for symmetric random variable like normal distribution). The kurtosis parameter of the sample is β₂ = 2.9419, near 3, the kurtosis of normal distribution.

However, test for normality of the sample gives the value for adjusted AD statistics A²_m = 1.3263 that is larger than the critical value. So we reject the null hypothesis of the normality of the solution sample (P < 0.0019).

This tells us that if one wants to estimate confidence intervals for sample values, knowing the distribution of parameters, the limiting values of the confidence intervals for normal distribution cannot be used.

Example 2

We consider a Neoclassical exchange economy with n = 3 goods andm = 4 agents with the CES utility function. The Excess demand functions are

In this Example, parameters and initial endowments are given in the following matrices

and i-th row of a matrix represents i-th agents' data. In our notation, the elements in the i-th row and j-th column of matrices α and ω are α_ij and ωⁱ_j. We get ξ by substituting the last equation in (4) with

Given the initial value π₀= (0.4, 0.2, 0.4), the Newton method with exit criterion ||ξ|| < ∈ = 10^—6 in 7 iterations gives the solution π * = (0.2441, 0.5566, 0.1993).

The sample used in this Example is obtained by perturbing components of parameter b: bⁱ, i = 1, 2, 3, 4 with N = 500 normally distributed values, bⁱ : N(0.5, 0.1), i = 1, 2, 3, 4.

Like in Example 1, the total number ofiterations, number offunction evaluations and Jacobian evaluations for both NM and FNM are given in Table 2. The average numbers of iterations are 6.60 and 9.16 for NM and FNM respectively but the corresponding CPU times are 1.267 and 0.820.

Descriptive statistics of the solution sample are given in Figure 2.

We can easily see from the histogram that the sample is asymmetric (β₁²> 0.6) and the kurtosis is far from 3, the kurtosis of normal distribution. Even more, AD test tells us the underlying distribution certainly (P < 10^—5) is not normal. This should be taken in account if confidence intervals for equilibrium prices are sought.

Example 3

We consider single commodity, nonlinear, separable spatial price equilibrium model (5) from [10] with two sources and two destinations.

Numerical tests are performed assuming that the supply and demand functions are

and the cost functions is linear,

Parameters u_t, vj, θj and y_ij belong to the intervals [3, 5], [15, 20], [0.1, 0.5], [5, 10] respectively.

The initial value for the Newton method for all sample values is and the exit criterion for both methods is ||F|| < ∈ = 10^—6.

For the following values of the parameters

the solution by the Newton method is , and it is obtained in 4 iterations.

Using a sample of N = 500 random values for independent components of v: v₁ and v₂ , both normally distributed v_1,2 : N(17.5, 1) and other parameters as above, we solved the equilibrium with the Newton method. After that, as in Example 1 we used the fixed Newton method, with the Jacobian and starting point obtained from the Newton method computed with the mean of parameter sample. The number of calculations is given in Table 3.

The average numbers of iterations are 10.93 and 6.45 for NM and FNM respectively with the corresponding CPU times 1.272 and 0.289. The descriptive statistics of the solution sample is given in Figure 3.

In two out of four cases for the solution sample the null hypothesis is not rejected. In these cases, an even larger sample is needed to reach the real properties of the solution and in these cases FNM algorithm will be even more efficient.

CPU time and larger dimensions

The results shown in Table 1 and 2 clearly demonstrate the advantage of FNM method. The reported CPU time is the total time for all sample values in seconds. If the dimensions of problems are larger the advantages of the proposed method are even more evident. To support this claim we tested Example 2 and Example 3 for different values of m and n. Thus for Example 2 we generate the matrix α₀ using random numbers scaled to row sums equal to 1, with w₀ = 2 + U(0, 1), b = [0.9; 0.9; ...; 0.9]. Here U(0, 1) denotes independent random numbers with uniform distribution on [0, 1 ]. The tested sample size is 500 again and the results (the total number of iterations and the total CPU ¹ time) are reported in Table 4 for NM and FNM method.

For Example 3 parameters were u_i = 4, vj = 17.5 + 0.5 • N(0, 1), θj = 0.3 and y_ij= 7.5 for all i and j. N(0, 1) denotes independent normally distributed values with mean 0 and standard deviation 1. The initial value was

the exit criterion was ||F|| < ∈ = 10^-6. The results are given in Table 5 - the total number of iteration and the total CPU time in seconds.

Conclusion

Equilibrium problems are often formulated as systems of nonlinear equations that depend on parameters. Significantly more information on equilibrium point can be obtained by simulation and solving a sequence of nonlinear systems rather than by solving a single nonlinear system with constant parameter values. We introduce a perturbation of model parameters taking their values from a given distribution sample. For each sample value we solve the corresponding system of nonlinear equations using the same starting point and applying Fixed Newton method. Thus we obtain the sample of solutions for each problem using a relatively cheap procedure. The sample of solutions yields better description of the observed system. For example the presented descriptive statistics show that the assumption of Gaussian distribution of the solution due to the Gaussian distribution of parameters is not guaranteed.

The question of confidence interval formulas for equilibrium prices is still open and it should be investigated from problem to problem. The stochastic simulations (Monte-Carlo) could give more information. If tests find the underlying distribution was normal, one could use well known formulas for confidence intervals for normal distribution. However, the algorithm proposed in this paper can give easier access to larger solution sample.

Solving many nonlinear systems requires a lot of CPU power which is usually limited. In our examples we showed that using similarities within sample data, an affordable procedure can be implemented. In Example 1 and 2 the total number of function calculations has risen following the rise of the number of iterations, but the number of expensive Jacobian calculations was almost zeroed. In Example 3 even the number of iterations lowered. In order to demonstrate the effectiveness of the proposed algorithm we included some tests with larger dimension problems. The number of iterations in these case is in line with small dimensional cases but the CPU time definitely favors the use of FNM method instead of NM method.

Acknowledgement. We are grateful to the two referees whose insightful comments improved the quality of this paper.

Received: 15/08/10.

Accepted: 05/01/11.

#CAM-319/10.

• [1]  T.W. Anderson and D.A. Darling, Asymptotic Theory of Certain "Goodness of Fit" Criteria Based on Stochastic Processes. The Annals of Mathematical Statistics, 23(2) (1952), 193212.
• [2] J. Arnal, H. Migallón, V. Migallón and J. Penadés, Parallel Newton Iterative Methods Based on Incomplete LU Factorizations for Solving Nonlinear Systems. Lecture Notes in Computer Science, 3402 (2005), 716729.
• [3] A. De la Fuente, Mathematical Methods and Models for Economists. Cambridge University Press (2000).
• [4] B.C. Eaves and . Schmedders, General equilibrium models and homotopy methods. Journal of Economics & Control, 23 (1999), 12491279.
• [5] M.B.M. Elgindi, On the application of Newtonâs and chord methods to bifurcation problems. Internat. J. Math. & Math. Sci., 17(1) (1994), 147154.
• [6] M. Esteban-Bravo, Computing equilibria in general equilibrium models via interior-point methods. Computational Economics, 23 (2004), 147171.
• [7] M. Esteban-Bravo, An interior-point algorithm for computing equilibria in economies with incomplete asset markets. Journal of Economic Dynamics and Control, 32(3) (2008), 677694.
• [8] J. Focke, A Simplified Newton Method for Computing the Factor Loadings in Maximum Likelihood Factor Analysis. Biometrical Journal, 28 (1986), 441453.
• [9] K.L. Judd, Numerical Methods in Economics. MIT Press (1998).
• [10] P. Marcotte, G. Marquis and L. Zubieta, A Newton-SOR Method for Spatial Price Equilibrium. Transportation Science, 26(1) (1992), 3647.
• [11]  H.E. Scarf, The Computation of Equilibrium Prices: An Exposition. In Handbook of Mathematical Economics (K.J. Arrow and M.D. Intriligator (eds.) Vol. 2, North-Holland, Amsterdam, (1981), 10071061.
• [12] X. Chen and T. Yamamoto, A convergence ball for multistep simplified Newtonlike methods. Numerical Functional Analysis and Optimization, 14(12) (1993), 1524.

*

Research supported by Ministry of Science of Serbia, grant 144006.

Publication Dates

• Publication in this collection
  22 Mar 2011
• Date of issue
  2011

History

• Accepted
  05 Jan 2011
• Received
  15 Aug 2010