RIA

laboria

Institut de Recherche
d'Informatique
et d'Automatique

Domaine de Voluceau
Rocquencourt
B. P. 105 78150 - Le Chesnay
France
Tél.: 954 90 20

laboratoire de recherche
en informatique
et automatique

AN AUTOMATIC METHOD
TO GET AN
ECONOMETRIC MODEL
IN A QUASI-
TRIANGULAR FORM

Pierre NEPOMIASTCHY
Alain RAVELLI
François RECHENMANN

Rapport de Recherche N° 313

Juin 1978
<figure><img src="image_1.png" /></figure>

AN AUTOMATIC METHOD TO GET AN ECONOMETRIC MODEL
IN A QUASI-TRIANGULAR FORM

Pierre Nepomiastchy*, Alain Ravelli**, François Rechenmann***

Résumé :

La résolution d'un modèle macroéconomique non linéaire de dimension n par une méthode de type Newton est présentée en utilisant la structure spécifique de la matrice de Jacobi du système d'équations qui décrit le modèle, de sorte que, si s des variables (s beaucoup plus petit que n) sont fixées, le système résultant est récursif et peut être résolu directement. Ces variables sont appelées les variables de bouclage du modèle. Un algorithme de réordonnancement des variables et des équations est proposé de manière à déterminer un nombre minimal de variables de bouclage par analyse directe de l'influence de telles permutations sur la structure de la matrice d'incidence du système global. Dans une approche entièrement automatique du problème, un préprocesseur peut calculer la matrice d'incidence par analyse syntaxique, équation par équation, et obtenir en sortie le code FORTRAN des équations réordonnées. L'efficacité de l'algorithme est démontrée sur quelques exemples.

Abstract :

A Newton-like method for solving non linear quasi-triangular econometric models with n equations is presented using a specific structure of the Jacobian in such a way that, as soon as a number s (much smaller than n) of variables, called loop variables, is given, the resulting model is recursive and can be solved directly. A renumbering algorithm, which analyzes the direct influence of varia- bles (and equations) permutations on the structure of the incidence matrix of the global system, is proposed. This algorithm computes a minimal set of those loop variables. An automatic approach is presented, where the incidence matrix is evaluated by a syntaxic analysis using a pre- processor, and FORTRAN code of the renumbered equations and variables is generated. Effecti- veness of the method is illustrated with some concrete examples.

* IRIA

** IRIA

*** Université de Grenoble
# 1. INTRODUCTION.

An econometric model is usually given as a set of dynamic equations, the unknowns of which are supposed to describe the behavior of an economy.

In mathematical terms, the simulation of a model, for a given time period, means the solution of a system with n equations and n unknowns (endogeneous variables), the coefficients of which are evaluated through real economic data (exogeneous variables), or desired estimated goals (control variables).

In recent years, a large increase of the dimension of models can be observed. For example, the DMS model [1] of INSEE (France), includes 1262 equations.

By the way, to use and even to construct such models (which means several simulation tests), asks for new specific technics which are to be different from the standard ones that have been performed for smaller models, or which can be found in general operation research packages.

For the present, seemly, the best way to take into account peculiarities of econometric systems is to adapt known solving methods using the sparsity of the Jacobi matrix of the system of equations, that describes the model.

There are few such iterative methods in the literature, which are closely related to other branches of mathematics : linear programming, graph theory, electrical network analysis (see §3.).

We recall in §2. a solving method that we have presented in [2]. This method uses the so-called quasi-triangular form of the Jacobi matrix J of the system (figure 1) in such a way that, as soon as a number of variables are given, the resulting system with n-s equations and n-s unknowns can be solved recursively (after a possible renumbering of variables and equations) :
<figure><img src="image_2.png" /></figure>

The non zero pattern of matrix J is hachured.

<caption>- Figure 1 -</caption>

Clearly, as it was noticed in [2], the use of such a structure is justified only when the number s of those "main" variables is much smaller than the total number of variables, and accordingly when the Jacobi matrix J is almost triangular. An optimal structure (i.e. a structure for which the number s takes a minimal value), when it has not been given by the constructors of the model, has to be determined by the users before any simulation.

Without looking systematically for this optimum, we bring up, in §4., an algorithm for renumbering equations and variables, in order to obtain a model as recursive as possible. This will be defined later as a minimal structure. In this paper, we are only studying symmetric renumerations of equations and variables ; in other words, we are looking for an n × n permutation matrix P, such that matrix P J P^T is the one in figure 1 with s << n (P^T denotes the transposed matrix of P).

Our algorithm gave good results for models (*) of different nationalities and with number of equations from 4 to 400. Through the variety of those models, this proves the great generality of methods which take into account the quasi-triangular structure as it seems to appear in most of econometric models.

In §5., we expose an entirely automatic approach in view to integrate the programs in a modular library.

<footnote>(*) ANDOMINI [5], PIMPON [6], FAIR [7], STAR [8], and model MPS from the Federal Reserve Bank.</footnote>
Finally, as it could be noticed, our study is uniquely based upon a qualitative analysis of interdependence of variables. Another method is a quantitative one which examines whether one variable is more important than others or not. This last problem is related to the construction of reduced models. A systematic study of the question is done in [3] by Dufour and Gilles. The authors measure the degree of interdependence within each couple of variables, using methods from information theory, in order to get rid of all connections beneath a "reasonable" level.
## 2. OUTLINE OF THE MODEL.

In this chapter, we briefly recall the terminology we used in [2]. For more details see [4].

### 2.1. *Structure of a quasi-triangular model.*

A discrete time dynamic model is considered with period t taking values from 1 to T.

Let $x_t = (x_t^1, \dots, x_t^n)$ be the vector in $\mathbb{R}^n$ of endogeneous variables at period t. The model simulation for this period t means the solution of a system with n equations, which can be written :

$$ (1) \qquad G_t(x_t) = 0 $$

where $G_t = (g_t^1, \dots, g_t^n)$ is a map $G_t : \mathbb{R}^n \to \mathbb{R}^n$.

Now, in most econometric models, a great number of equations (the behavioral equations) are estimated as relations between one variable (depicted variable) and other variables (depicting variables). It seems then natural to write such an equation in the form $x_t^i = f_t^i(x_t)$, and so we can entirely describe the system (1) writing it in the form :

$$ (2) \qquad K x_t = F_t(x_t) $$

where $F_t = (f_t^1, \dots, f_t^n)$ is a map $F_t : \mathbb{R}^n \to \mathbb{R}^n$, and K is an $n \times n$ diagonal matrix :

$$ K = \left( \begin{array}{c@{\quad}c} \begin{matrix} k_1 & \\ & k_n \end{matrix} & 0 \end{array} \right), \text{ such that :} $$

$$ \left\{ \begin{array}{ll} \bullet & k_i = 1 \text{ if and only if equation i is analytically solvable with respect to } x_t^i \text{ for every t (this equation is then written in (2) in its solved form),} \\ & k_i = 0 \text{ otherwise.} \end{array} \right. $$

The matrix K is therefore independent ($^*$) of the period t.

The solved form is very useful while simulating the model equation per equation, for we can obtain variable $x_t^i$ as an explicit function of other variables. If equation i is not solvable with respect to $x_t^i$ (i.e. $k_i = 0$), then we have to write it : $x_t^i = x_t^i + \alpha_i f_t^i(x_t)$ (with some $\alpha_i \in \mathbb{R}$ : weight ), and we need to know an initial value of $x_t^i$ in order to

<footnote>(*) In case we should do a dynamic simulation.</footnote>
solve equation i, using such an iterative method as :

$$x_t^{i,n} = x_t^{i,n-1} + \alpha_t f_t^i(x_t^{n-1}) \quad (\text{Gauss-Seidel}).$$

In figure 2, we represent the number of depicted variables and the total number of variables for each model.

<table><caption>- Figure 2 -</caption><thead><tr><th></th><th>ANDOMINI<br>(5)</th><th>PIMPON<br>(6)</th><th>FAIR<br>(7)</th><th>STAR<br>(8)</th><th>MPS</th></tr></thead><tbody><tr><th>Var.</th><td>4</td><td>14</td><td>83</td><td>139</td><td>400</td></tr><tr><th>Depicted Var.</th><td>4</td><td>14</td><td>82</td><td>136</td><td>400</td></tr></tbody></table>

**Definition 1** : A model is said to be put in a quasi-triangular form, if there are some integers $l$, $m$ and $s$ such that the associated system of equations (2) may be split into 4 successive blocks :

$$(3a) \qquad x_t^i = f_t^i(x_t^1, \dots, x_t^{i-1}), \quad i=1, \dots, l$$

$$(3b) \qquad x_t^i = f_t^i(x_t^1, \dots, x_t^{i-1}, x_t^{m-s+1}, \dots, x_t^m), \quad i=l+1, \dots, m-s$$

$$(3c) \qquad k_i x_t^i = f_t^i(x_t^1, \dots, x_t^m), \quad i=m-s+1, \dots, m$$

$$(3d) \qquad x_t^i = f_t^i(x_t^1, \dots, x_t^{i-1}), \quad i=m+1, \dots, n$$

Those blocks are respectively called : prologue, heart, loop equations, and epilogue.

Such a structure is clearly distinguished by the form of the Jacobi matrix $J$ of the function $F_t = (f_t^1, \dots, f_t^n)$ in system (3) (which will be called the Jacobian of system (3) in brief) :
<figure><img src="image_3.png" /><figcaption>- Figure 3 -</figcaption></figure>

where submatrices $J_1$ and $J_3$ are low triangular and the central block $J_2$ has a quasi-triangular structure like the one in figure 1.

It could be noticed that the general system (2) has a "rough" quasi-triangular structure setting $l=0$ and $m=s=n$. Nevertheless this "rough" structure does not take into account sparsity of matrix J as we shall do later.

## 2.2. Notion of loop variable.

Let us consider the general system with n equations (2).

**Définition 2** : The incidence matrix of system (2) at period *t* is the n × n matrix $E_t = (e_{ij}^t)$ such that :

$$ egin{cases} 
\begin{array}{l} 
\bullet \ i \neq j \quad e_{ij}^t = 0 \quad , \ i=1, \dots, n \quad , \ j=1, \dots, n, \text{ if and only if } \frac{\partial}{\partial x_t^j} (f_t^i(x_t)) \\ 
\text{is identically equal to } 0, \text{ and } e_{ij} = 1 \text{ otherwise} 
\end{array} \\ 
\bullet \ e_{ii}^t = 1 - k_i \quad , \ i=1, \dots, n 
\end{cases} $$

The *incidence matrix* of system (2) is then the matrix $E=(e_{ij}^t)$ such that $e_{ij}^t = \sup_{t=1,...,T} (e_{ij}^t)$,

The incidence matrix E is independent of the period t, and gives every interdependence information between variables for each period t=1,...,T.

**Definition 3** : A variable $x_t^j$ is said to be a *loop variable* of system (2) if and only if :

$$ \prod_{i=1}^{j} (1 - e_{ij}^t) = 0 $$
It will be noticed that $k_j = 1 - e_{jj}$. Therefore, $x_t^j$ is a loop variable either if the equation j is not analytically solvable with respect to $x_t^j$, or if the column j in the incidence matrix is a spike column (see Hellerman and Rarick [10]), that is to say, contains non zero elements above main diagonal.

**Remark** : As a matter of fact, the incidence matrix and the set of loop variables of system (2) are not canonically related to that system of equations : every ordering of equations(or variables) will corresponds to probably distinct incidence matrices and loop variables sets. So, we shall only speak of the incidence matrix or loop variables set of a system of equations with reference to a given ordering of those equations.

**Définition 4** : System (2) is said to be recursive if and only if its set of loop variables is empty.

Clearly, a recursive model can be solved at once by a recursive substitution of each variable within the following equation after a possible renumbering of variables and equations such that the incidence matrix becomes low triangular.

Then, indeed, the model will be written :

$$ (2') \qquad x_t^i = f_t^i(x_t^{i-1}, \dots, x_t^{i-1}), \quad i=1, \dots, n $$

One could notice that the Jacobian J of this system is low triangular too. In fact, this relationship between structures of the incidence matrix on one hand and the Jacobi matrix on the other hand, is a general result :

**Proposition 1.**

Every non diagonal entry of the Jacobian of system (2) is a zero element as soon as the corresponding entry in the incidence matrix is a zero element too.

**Proof** : Let us write $e_{ij}=0$ with $i \neq j$. This means :

$$ \frac{\partial}{\partial x_t^j} (f_t^i(x_t)) = 0, \quad t=1, \dots, T $$

or

$$ J_{ij} = 0 $$

if $J=(J_{ij})$ is the jacobian of system (2) at period t. (q.e.d.)
The main result we can draw from this proposition is the following :

**Corollary** : Any line and column permutation in the incidence matrix of system (2) (i.e any renumbering of variables and equations) leading to a structure like the one in figure 2, will allow us to write the system in a quasitriangular form (3).

Finally a simple examination of the incidence matrix of system (3) implies the following proposition :

**Proposition 2.**

The loop variables of system (3) are exactly the variables the indices of which correspond to the ones of loop equations.

### 2.3. *An adapted method for solving a quasi-triangular model.*

For this section, more details are given in [2] and [4].

**Proposition 3.**

For given values of loop variables in system (2), the corresponding reduced (*) system is recursive.

This proposition is a direct consequence of the definition of loop variables. Nevertheless, it becomes obvious, through proposition 2, when the system is written in form (3) (the Jacobian of the reduced system is then a low triangular matrix).

Proposition 3 shows that the set of loop variables is exactly the set of variables we have to initialize (also called iterative set by Kevorkian and Snoek in [14]), so that solving the model equation per equation is reduced to the use of successive formulas (we have already seen that variables $x_t^i$, for which $k_i=0$, are necessary such variables).

Let us call $y_t = (x_t^{m-s+1}, \dots, x_t^m)$ the vector in $\mathbb{R}^s$ of loop variables of system (3).

For a given value of $y_t$, the recursive solution of the reduced system from (3), leads to a vector $x_t(y_t)$ in $\mathbb{R}^n$, such that :

$$x_t(y_t) = (x_t^1, \dots, x_t^{m-s}, y_t, x_t^{m+1}, \dots, x_t^n)$$

We are able to evaluate the corresponding errors upon the loop equations which we shall denote by :

<footnote>(*) after elimination of the corresponding equations.</footnote>
$$ \varphi_t^i(\bar{y}_t) = k_i x_t^i - f_t^i(x_t(\bar{y}_t)), \quad i=m-s+1, \dots, m $$

Then, a solution $\bar{x}_t = x_t(\bar{y}_t)$ of (3) will be given with $\bar{y}_t$ satisfying :

$$ (4) \qquad \varphi_t^i(\bar{y}_t) = 0, \quad i=m-s+1, \dots, m $$

We proposed in [2] a Newton-like method to solve system (4).

So the problem is reduced to the solution of a system with s equations and s unknowns (the loop variables). This method is very interesting especially when n is large and s is much smaller than n. Allocation problems, convergence rate problems are rather scattered. For example, with a Newton method the use of our method (resp. of a global method) leads to the solution of a linear system of (low) dimension s (resp.n) (*).

We give in figure 4 the values of that number s for each mentioned model :

<table><caption>- Figure 4 -</caption><thead><tr><th></th><th>n</th><th>s</th></tr></thead><tbody><tr><td>ANDOMINI (5)</td><td>4</td><td>1</td></tr><tr><td>PIMPON (6)</td><td>14</td><td>2</td></tr><tr><td>FAIR (7)</td><td>83</td><td>6</td></tr><tr><td>STAR (8)</td><td>139</td><td>3</td></tr><tr><td>MPS</td><td>400</td><td>20</td></tr></tbody></table>

Those values have been obtained through our algorithm (see §4.) except for model STAR for which this number s has been given by the constructors.

<footnote>(*) However, there are many good algorithms in order to take into account the sparsity of the matrix for determining its inverse, see Drud [9], Hellerman and Rarick [10], [11], Markowitz [12], Lasdon [13], for instance.</footnote>
### 3. A SIMILAR PROBLEM IN GRAPH THEORY.

In this chapter, we analyse relations between the problem we set out in our paper and some classical problems the solution of which can be found in the literature.

The problem of reduction of a matrix to a quasi-triangular form has already been studied by Kevorkian and Snoek [14], who defined this structure as the LNF-form (lower N-form), and by Cheung and Kuh[15], as the BTF-form (bordered triangular form).

However, D.V. Steward was the first who described the splitting of large systems of non linear equations in [16], [17]. He proposed to reduce the system to several subsystems of minimal size, which could easily be solved one after the other. This method is equivalent to the reduction of the Jacobian J of the system to a block triangular form :

<figure><img src="image_4.png" /></figure>

The non zero pattern of matrix J is hachured.

<caption>- Figure 5 -</caption>

When the dimension of subsystems is still too large, D.V. Steward proposes to get on with that reduction by tearing, that is to say collapsing one or several entries in matrix J (i.e getting fixed some variables in some equations).

Obviously, after a renumbering of variables and equations in order to put those torn variables at the end, the matrix J takes the form :
<figure><img src="image_5.png" /></figure>

The non zero pattern of matrix J is hachured.

<caption>- Figure 6 -</caption>

This structure, called the quasi-block-triangular structure, is the same as the quasi-triangular one in the case when the block J' is low triangular.

Therefore, our problem is equivalent to the following : determine a subset of torn variables with minimal size such that the torn system (i.e. obtained from tearing) should be recursive.

A lot of algorithms were proposed to obtain the different forms corresponding to figures 1, 5, 6.

The *block triangular structure* has been studied by Dulmage and Mendelsohn [18], Harary [19] for systems of linear equations. In [18], it is shown that this structure is unique as soon as each diagonal block is irreducible. Algorithms are proposed by Drud [9], Hellerman and Rarick [11], Steward [17], Kevorkian and Snoek [20], Tarjan [21].

The *quasi-block triangular structure* is rather uneasy to study because of the difficulty to define a criterion of irreducibility. This structure is discussed by Tewarson [22] and Dickson [23]. D.V. Steward studied the problem in terms of tearing [24]. Nathan and Even [25] indirectly tried to solve the problem.

Finally, the *quasi-triangular structure* has been studied by Kevorkian and Snoek [14], Cheung and Kuh [15], Kevorkian [26], [28], Ledet and Himmelblau [27]. It is always defined in relation with the information flow introduced in [16], [17]. This flow diagram sums up interactions between variables by classical means in graph theory. It is also defined in [30] as the associated digraph of the system to be considered.
Let us have a look back to the terminology of Harary [29], [30].

An oriented graph $G$, or digraph, is a set of points, or vertices, indexed by integers $i=1,...,n$ ($n$ is the dimension of $G$), some couples of points $(i,j)$ being joined by oriented paths, or arcs : $i \to j$

**Definition 5** : The adjacency matrix of a digraph $G$ of dimension $n$, is the $n \times n$ boolean matrix $E = E(G)$, such that :

$$ E = (e_{ij}) \quad \text{and} \quad \left\{ \begin{array}{l} e_{ij} = 1 \text{ if and only if } j \to i \text{ in } G, \\ e_{ij} = 0 \text{ otherwise.} \end{array} \right. $$

The associated digraph of a $n \times n$ matrix $E=(e_{ij})$, is the digraph $G = G(E)$ of dimension $n$, such that : $i \to j$ if and only if $e_{ji} \neq 0$.

If $E$ denotes a matrix, $E(G(E))$ is the associated boolean matrix ($^*$), while for every digraph $G$, $G(E(G)) = G$.

Now, let us consider again the general system (2) :

**Definition 6** : The associated digraph of the system (2) is the digraph $G(E)$ where $E$ is the incidence matrix of (2).

In other words, and coming back to definitions 2 and 4, this graph is defined as the digraph the vertices of which are numbered by the indices of equations (or variables) of the system, arcs being defined by :
for $i \neq j$ $i \to j$ if and only if $J_{ji} = \frac{\partial}{\partial x_t^i} (f_t^j(x_t))$ is not identically equal to zero, for at least one period $t$ of simulation, and $i \to i$ (self cycle) if and only if $k_i = 0$.

As an example, we illustrate the different notions we have just introduced on the small model C.S ($^{**}$)(conjectural - structural) that shall be used in the next chapter.

<footnote>($^*$) the boolean matrix whose zero (resp. non zero) entries correspond to zero (resp. non zero) entries in $E$.</footnote>

<footnote>($^{**}$) "A macro model for the Dutch economy", Working Group Conference of European Planners", La Haye, 24-25, April 1969.</footnote>
This model involves 7 equations :

$$1. \quad S = 0.45 * A * (J+L) - 11.9$$

$$2. \quad A = D + 26.69$$

$$3. \quad L = 0.22 * P + 0.395 * H + 0.087$$

$$4. \quad P = 0.25 * (L-H) + 0.04$$

$$5. \quad H = (37.25/B) - 1$$

$$6. \quad B = D + 36.02$$

$$7. \quad D = 0.297 - 0.1036 * L$$

where variables are respectively :

S : Change in disposable income.

A : Paid employment in industry.

L : Income rate per hour.

P : Consumption price.

H : Labor productivity in industry.

B : Total employment in industry.

D : Change in employment in industry.

This model is non linear. It is written in form (2) where K is the identity matrix.

The incidence matrix of the model is the following :

<figure><img src="image_6.png" /><figcaption>- Figure 6 -</figcaption></figure>
The associated digraph is the following :

<figure><img src="image_7.png" /><figcaption>- Figure 7 -</figcaption></figure>

It is shown in [30] that the system (2) is recursive if and only if
the associated digraph is acyclic (i.e. contains no cycle).

The model C.S, for example, is not recursive since it contains the
cycle : {3,7,6,5,4,3} .

Consequently :

**Proposition 4.**

The determination of a quasi-triangular structure in system (2)
is equivalent to the determination of a set of vertices in the associated
digraph such that the corresponding reduced (*) digraph should be acyclic.

Such a set of vertices is called a set of _essential vertices_.
A minimal set of loop variables therefore corresponds to a minimal set of
essential vertices [15] (also called minimal feedback vertex set by
Smith and Walford [31]). The _index_ of the associated digraph is defined as
the lower cardinality of sets of essential vertices.

One could notice that the necessity to cut a self cycle i → i, while
determining a set of essential points, does correspond to the fact that
variable of index i is then a loop variable (kᵢ=0): vertex i belongs to any
set of (minimal) essential vertices as well as variable xᵢ^i with kᵢ = 0
necessarily is a loop variable.

<footnote>(*) obtained by removing corresponding vertices and all arcs joined to those vertices.</footnote>
Most of the algorithms determining a minimal set of essential vertices were discussed in relation with network analysis in signal flow problem as in Nathan and Even [25], Divietti [32], Lempel and Cederbaum [33], where one mainly looks for control points (i.e essential vertices) in circuits (i.e digraphs).

Those algorithms generally involve two parts :

in a first part, the digraph is reduced by local topological rules, while preserving the index,

in a second part, one looks for generating the whole or a part of the cycles of the reduced digraph,
in order to find a minimal (or quasi minimal) number of essential points (i.e cut points within that pertinent circuits). See Smith and Walford [31], Divietti and Grasselli [34], Guardabassi [35], [38], Diaz et al. [36], Yau [39].

Most of algorithms already quoted in relation with the determination of a minimal quasi-triangular structure are direct applications of those algorithms.

It should also be stated that the structures presented above (block triangular, quasi-block triangular, quasi-triangular) actually can be combined to obtain other types of structures, for example after partitionning (*) and looking for a minimal quasi-triangular structure within each diagonal block.

Nevertheless, the use of quasi-triangular decomposition becomes inevitable as long as the system involves an irreducible (i.e with no possible further partitionning) block of large dimension. In FAIR [7], for instance the model involves such a block of dimension 74.

As a conclusion for this chapter, we should like to point out the problems of determining effective algorithms in order to find minimal sets of essential vertices in a digraph. Most of the already quoted

<footnote>(*) In terms of graph theory, this is looking for irreducible (i.e strongly connected) components of the associated digraph [29], [30].</footnote>
algorithms get very expensive as soon as the size of the digraph (i.e of the system) becomes large. The number of cycles in the digraph to be determined can get quite great. The difficulty has been characterized (*) in [40]. It seems that for such problems, the best strategy is to determine a good approximate solution.

That will be done in next chapter.

<footnote>(*) This problem is NP - complet.</footnote>
## 4. A RENUMBERING ALGORITHM FOR GENERATING A MINIMAL SET OF LOOP VARIABLES.

Let us consider an econometric model in form (2). E will denote its incidence matrix.

### 4.1. "Rough" and "fine" analysis.

While looking for a minimal number of loop variables, we could think of minimizing the number of zero coefficients $k_i$, for which the variables $x_t^i$ are necessary loop variables. We shall see that matrix K may be quite different if we make a "rough" or a "fine" analysis.

This problem is related to the determination of the so-called output pairs $(x_t^i, f_t^i)$ in system (2) (i.e which variable $x_t^i$ is associated to equation i), when this choice is not natural as for solvable equations (see §2.1). This problem seems to be essential since each (symmetric) renumbering of equations and variables is to be done according to a given set of output pairs.

For example, equation $f(x_1,x_4,x_8) = 0$ can generally be written in form $x_1 = f_1(x_4,x_8)$ or $x_4 = f_4(x_1,x_8)$ or $x_8 = f_8(x_1,x_4)$. From those three expressions, one is generally analytically known; then choosing the coresponding variable as the depicted one through equation $f(x_1,x_4,x_8)= 0$, a new non-zero entry appears in K, after possible renumbering of variables and equations.

We have a similar problem when some variables are the solutions (within a block) of a linear system which is not analytically solved in the given form. In FAIR [7], for instance, as soon as variables $CG_t$, $DDF_t$, $BORR_t$, $SAVB_t$, $LF_t$ and $SECR_t$ are determined, the five equations: 45, 61, 62, 64, 70 are found to involve a linear system with five unknowns: $BR_t$, $A_t$, $DDB_t$, $LBVBB_t$ and $DDH_t$.

Indexing those equations from 1 to 5, the five last variables called $x_1, \dots, x_5$ may be defined as the solution of the following set of linear equations:

$$1. \quad x_1 = a x_3$$

$$2. \quad x_2 = -x_5 + b$$

$$3. \quad x_3 = x_5 + c$$

$$4. \quad x_4 = -x_1 + x_3 + d$$

$$5. \quad 0 = -x_4 - x_2 + e$$
In a "rough" analysis, on one hand, equation 5 (the equation with number 70 in the model) , that is the only one in FAIR which do not define explicitely one of the endogeneous variables, may, therefore and in a natural manner, determine the last variable : RBILL<sub>t</sub>. This means that a zero coefficient k<sub>70</sub> = 0 and RBILL<sub>t</sub> is necessarily a loop variable.

On the other hand, solution of the described system above leads
to :

$$
1'. \quad ax_3 = b + c + d - e
$$

$$
2'. \quad x_1 = ax_3
$$

$$
3'. \quad x_2 = -x_3 + b + c
$$

$$
4'. \quad x_4 = (1 - a)x_3 + d
$$

$$
5'. \quad x_5 = x_3 - c
$$

So, equation 5, now 5', has been used to determine the variable
DDH<sub>t</sub> and coefficient k<sub>70</sub> is now equal to 1. Similarly, equation 8 in FAIR
which, before, determined the variable DDH<sub>t</sub> in a natural way and which is
written :

$$
\log (\mathrm{DDH}_t) = - 0.0255 * \log (\mathrm{RBILL}_t) + 0.432 * \log (\mathrm{YH}_t) + f \\
\text{may also be written :}
$$

$$
0.0255 * \log(\text{RBILL}_t) = - \log(\text{DDH}_t) + 0.432 * \log(\text{YH}_t) + h \\
\text{and so it may be used to determine the variable } \text{RBILL}_t, \text{ coefficient } k_8 \\
\text{keeping the value } 1.
$$

Anyway, afterwards, our algorithm will allow us to know whether
RBILL<sub>t</sub> is a "true" loop variable or not.

It is clear that this kind of manipulation keeps a very special
character depending on the model to be considered. In an automatic
approach of the problem, one allows only a "rough" analysis. So, for each
equation i, one will look whether or not it is given in form :

$$
x_t^i = \text{some function of the other variables.}
$$

It is easy to come down to this form if the variable is defined
through its growth rate $\frac{x_t^i - x_{t-1}^i}{x_t^i}$, or its increase $x_t^i - x_{t-1}^i$,
or its logarithm $\log (x_t^i)$.
According to the acquired experience, this seems to be rather
satisfying in a first approach.

### 4.2. *The algorithm.*

The first part of the algorithm (step 1 to 5) determines the prologue and the epilogue of the model. It can be drawn out from studies by Drud [9] and Hellerman and Rarick [11]. The second part (step 6 to 10), which determines a smaller set of loop variables is entirely new.

The algorithm is performed on the incidence matrix E of the model :

*Step 1:* Determine RC and CC which are respectively the numbers of non-zero entries in each line and each column of the incidence matrix.
Set $l = 0$.

*Step 2:* If none of the RC(i)'s is zero, go to step 3. Otherwise, set $l = l + 1$ and remove line i and column i. Then add equation i to the set of equations (which is perhaps empty) which already forms the prologue, as element of order $l$ in this set. Revise RC and CC (we recall that line i and column i have been removed) and return to step 2.

*Step 3:* If $l = n$, the incidence matrix can be put in a lower triangular form through the reordering of equations found in the step 2 (prologue). The model is recursive and stop.

Otherwise, set $m = n$.

*Step 4:* If none of the CC(j)'s is zero, go to step 5.

Otherwise, remove line j and column j. Then add equation j to the set of equations (which is perhaps empty) with numbers n to 1, which already forms the epilogue as element of order m in this set. Revise RC and CC, set $m = m-1$ and return to step 4.

*Step 5:* Renumber the equations, starting with those of the prologue (of dimension $l$) and ending with those of the epilogue (with dimension $n-m$), respecting orders that have been found, within those two sets.

We shall now only consider the reduced matrix (i.e corresponding to the reduced system without prologue and epilogue).

We are going to illustrate that first part of the algorithm with model C.S (see §3.)
In figure 8, we show the different values of RC and CC for
step 1 to 4 :

<figure><img src="image_8.png" /><figcaption>- Figure 8 -</figcaption></figure>

The result at end of step 5 is drawn out in Figure 9. The prologue
is empty and epilogue contains equations {2,1}

<figure><img src="image_9.png" /><figcaption>- Figure 9 -</figcaption></figure>
On the reduced system which contains equations {3,4,5,6,7} we now have to take into evidence some set of loop variables as reduced as possible. The algorithm is performed for step 6 to 10 within the reduced matrix as shown in figure 9.

*Step 6:* Determine a first set of s loop variables from the reduced system, the incidence matrix of which has been outlined from step 5. Renumber equations in order to have those equations, which define the s loop variables, at the end.

In general, because of the sparsity of the incidence matrix, this number s is already small compared to n. Nevertheless, in model C.S. we have s = 4; loop variables are: {4,5,6,7} and this result seems still not to be satisfying. Therefore, a finer analysis is necessary.

We have to test each of the loop variables (with numbers from 1 to s) successively.

**Définition 7** : We shall say that variable $x_t^i$ is a predecessor of variable $x_t^{j'}$ if $i < j$ and if there are some variables $x_t^{h_1},...,x_t^{h_p}$ with

$$i < h_1 < h_2 < \dots < h_p < j$$

such that entries $e_{h_1 i}, e_{h_2 h_1}, \dots, e_{h_p j}$ in the incidence matrix are non-zero entries.

<figure><img src="image_10.png" /><figcaption>- Figure 11 -</figcaption></figure>

**Definition 8** : We shall say that variable $x_t^j$ is a successor of variable $x_t^{j'}$ if $i < j$ and if entry $e_{ij}$ in the incidence matrix is a non-zero entry.
<figure><img src="image_11.png" /><figcaption>- Figure 12 -</figcaption></figure>

*Step 7:* Every variable $x_t^i$ for which $e_{ii} = 1$ (i.e. $k_i = 0$) is a "true" loop variable.

If all variables are so, no further reduction is possible ; then stop.

Otherwise, remove lines and columns corresponding to indices of
those variables and go on with the reduced matrix. The other loop
variables are indexed with an integer k. Set k=1.

Step 8 : Determine the sets of predecessors and successors of loop variable $x_t^i$ with index k.

If the intersection of those two sets is empty, go to step 9.

Otherwise, go to step 10.

Step 9 : The loop variable $x_t^i$ is a "false" loop variable. We can remove it by renumbering the equations in such a way that the system involves successively the set of equations defining the set of predecessors of $x_t^i$ (with previous order), equation i and then the set of equations defining the set of successors of $x_t^i$ (with previous order too).

Let us consider the following simple example, with four variables
with numbers 1 to 4 :

<figure><img src="image_12.png" /></figure>
The predecessors of variable 4 are : {1,3} ; there is only one
successor of variable 4, that is variable {2} . The two sets are
distinct and writing equations in order : {1,3} , {4} , and {2} , we
obtain the matrix in low triangular form :

<figure><img src="image_13.png" /></figure>

The main feature is that variable x_t^i, as considered, becomes a
predecessor of its previous successors which means that there are no
longer non zero entries in the incidence matrix above the main
diagonal in column corresponding to that variable. Furthermore, order
is kept within each subset to be considered and no non zero entry is
added above the main diagonal.

Set s=s-1, k=k+1 and return to step 8 except when all loop
variables have been tested, then stop.

Step 10 : The loop variable x_t^i is a "true" loop variable. That means
that the corresponding vertex in the associated digraph belongs to a cycle.
Removing that vertex will cut the cycle. So, remove corresponding line
and column in the incidence matrix.

Set k=k+1 and return to step 8, except when all loop variables have been tested, then stop.
END.

So, in output, the set of loop variables for the new order of
equations is exactly the set with the s "true" loop variables.

Let us come back to our example. Model C.S has been reduced to
equations : {3,4,5,6,7} . In step 8 to 10 we test successively variables
{4,5,6,7} as "false" or "true" loop variables.
So, the only predecessor of variable 4 is variable 3 which is also a successor of variable 4. Variable 4 is a "true" loop variable and matrix is now reduced to :

<figure><img src="image_14.png" /></figure>

Variable 5 has no predecessor. It is surely a "false" loop variable. The set of successors of variable 5 involves only the variable 3. Let us write the set of equations in order : 5, then 3, then 6, 7 :

<figure><img src="image_15.png" /></figure>

There is no longer a 1 above the main diagonal in column corresponding to variable 5.

Variable 6 is also a "false" loop variable. We write equations in order 6, then 5, then 3, 7 :

<figure><img src="image_16.png" /></figure>
The set of predecessors of variable 7 is then : 6,5,3. The only successor is variable 6. The two sets have variable 6 in common and variable 7 is a "true" loop variable. All loop variables have been tested and this is the end of the study.

We then have to reorder equations such that we should successively find :

- prologue : empty

- heart : equations 6, 5, 3

- loop equations : equations 7, 4

- epilogue : equations 2, 1

and with the new ordering, the incidence matrix becomes :

<figure><img src="image_17.png" /><figcaption>- **Figure 13** -</figcaption></figure>

That structure will be compared to that of figure 3.

Programs have been written in Fortran H extended and tested on an IBM 370-168 computer with large models : FAIR (83 equations) and MPS (400 equations).

Results and computing times are given in figure 14. Sparse matrices technics have been used in order to save memory allocation. Indeed, only 2,09 % (resp.: 0,54 %) of the entries of the incidence matrix of the FAIR (resp.: MPS) model are non zero ones. Technics for generating that matrix are discussed in next chapter.
In figures 15 and 17, we give the incidence matrices of FAIR and MPS in output of our algorithm. Incidence matrix of model STAR is presented in figure 16 such as the model has been given by the constructors.

### 4.3. Results.

<table><thead><tr><th></th><th>FAIR</th><th>STAR</th><th>MPS</th></tr></thead><tbody><tr><td>Number of equations</td><td>83</td><td>139</td><td>400</td></tr><tr><td>Prologue</td><td>2</td><td>58</td><td>184</td></tr><tr><td>Epilogue</td><td>7</td><td>5</td><td>34</td></tr><tr><td>Computing time<sup>(*)</sup> steps 1 to 6</td><td>0.02</td><td>-</td><td>1.57</td></tr><tr><td>Initial number of loop var. s<sub>o</sub></td><td>30</td><td>-</td><td>87</td></tr><tr><td>Number of "true" loop var. s</td><td>6</td><td>3</td><td>20</td></tr><tr><td>s / s<sub>o</sub></td><td>0.2</td><td>-</td><td>0.23</td></tr><tr><td>Heart</td><td>68</td><td>73</td><td>162</td></tr><tr><td>Computing time<sup>(*)</sup> steps 7 to 10</td><td>1.12</td><td>-</td><td>45.63</td></tr></tbody></table>

(*) in CPU seconds

- Figure 14 -
<figure><img src="image_18.png" /><figcaption>- Figure 15 -</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>- Figure 16 -</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>- Figure 17 -</figcaption></figure>
### 4.4. A minimality theorem.

The algorithm we have just described allows to obtain a subset of "true" loop variables from the first set of loop variables related to the model in the given form, and that subset is exactly the set of loop variables of the reordered model.

That set is generally much more reduced than the first one. Anyway, there is no reason for the number of those "true" loop variables to be the minimum number of such variables (i.e. equal to the index of the associated digraph).

In model C.S, permutation of variables and equations in order : {7,6,5,4,3,2,1} does clearly draw out only variable 3 as a "true" loop variable in the new incidence matrix (though we found two such variables with our algorithm) :

<figure><img src="image_21.png" /><figcaption>- Figure 18 -</figcaption></figure>

This may easily be explained ; indeed the first set of loop variables we have determined, essentially depends on the order such that equations of the model are given. As there is no reason for this set to contain a set of "true" loop variables with minimum cardinality, there is no reason for the algorithm to build such a set of variables.

That is why we have tested our algorithm, initializing it with some given prescribed random permutations of equations and variables. The number of "true" loop variables thus oscillated from 6 to 9 for FAIR and from 20 to 35 for MPS, the ratio $s_o/s$ keeping approximatively the
same value. The great stability of those results allows to expect a good approximation of the minimum by our algorithm in most of cases.

Nevertheless, if the set of "true" loop variables is not necessarily the smallest set of loop variables that we could find, it is however minimal in the sense of Harary [41] :

**Definition 9** : A set S is said to be **minimal** (resp. **minimum**) with respect to a property "P" if S satisfies "P" and no subset, strictly contained in S, satisfies "P" (resp. no other set T, containing strictly less elements than S, is minimal with respect to "P").

**Theorem.**

The set of "true" loop variables is a minimal set of loop variables.

**Proof** : As a matter of fact, a "true" loop variable is determined by our algorithm (step 10), as soon as it is contained in a cycle within the associated digraph, which contains no other "true" loop variables.

Then, if we could determine a smaller set of "true" loop variable inside the set S we have found, it would exist a variable $x_t^i$ belonging to S and not to S'.

As an element of S, the vertex of index i of the associated digraph is contained in a cycle B, which itself contains no vertex corresponding to indices of elements in S.

But, after removing the vertices corresponding to the variables in S' in the associated digraph, the reduced digraph still contains the cycle B.

It thus cannot be acyclic, which provides a contradiction with the definition of S'. (q.e.d)
## 5. AN AUTOMATIC APPROACH OF A MODEL.

As we have proposed to do, a new ordering of equations and variables can be determined before any simulation, in order to get a model in a "better" form.

From an automatic point of view, that means we have to generate the new code of the renumbered equations of the model, using a preprocessor for example.

So, input could be the Fortran code (for instance) of the given equations of the model.

The preprocessor, then performs a syntaxic analysis equation per equation. In the course of that analysis, a lot of usefull parameters may be determined : the file of variables mnemonics, with the associated types (exogeneous, endogeneous,...) for each variable, the matrix K, the output pairs ($x_t^i, f_t^i$) (see §2.) and so on.

Then, the incidence matrix is generated, analysing whether an endogeneous variable appears in an equation or not.

Furthermore, other formal manipulations may be performed such as : formal derivation in order to compute the Jacobian of the system, formal analysis in order to determine the constant part of each equation (predetermined part), or of the Jacobian (linear part of each equation).

Actually, the Fortran code of the renumbered equations is generated. (see Figure 9).

Of course, generating the incidence matrix seems to be quite easy as soon as the Jacobian of the system is known. Non zero non diagonal entries do correspond in both matrices. Examining whether each equation is solvable or not will complete the information about the main diagonal.

We have proposed formal derivation in order to compute the Jacobian. Anyway, the well known finite difference method gives a good automatic method too, to compute that matrix.
<figure><img src="image_22.png" /><figcaption>- Figure 19 -</figcaption></figure>
REFERENCES

============

(1) D. FOUQUET ; J.M. CHARPIN ; H. GUILLAUME ; P.A. MUET ; D. VALLET :
"DMS, modèle de prévision à moyen terme", Economie et Statistique
n°79, 1976.

(2) P. NEPOMIASTCHY ; A. RAVELLI : "Adapted Methods for solving and optimizing quasi triangular econometric models", Annals of Economics and Social Measurement, Vol. 6, n°6, 1978.

(3) J. DUFOUR ; G. GILLES : "Approche informationnelle de l'analyse structurelle des systèmes complexes", to appear.

(4) D. GABAY ; P. NEPOMIASTCHY ; M. RACHDI ; A. RAVELLI : "Etude, Résolution et Optimisation de Modèles Macroéconomiques", Rapport IRIA/LABORIA, 1978.

(5) K. WALL : NBER Cambridge Memo., 1975.

(6) A. COUTIERE : "Un modèle du système monétaire français", Statistiques et Etudes financières, n°17, INSEE, 1975.

(7) R.C. FAIR : "A model of Macroeconomic Activity", Vol. II, The Empirical Model, Ballinger Publishing Company, 1976.

(8) J. MISTRAL : "STAR, Modèle de prévision à court moyen-terme", Séminaire sur les modèles macroéconomiques, (C.N.R.S. / N.S.F.), Gif/Yvette, 1976.

(9) A. DRUD : "Methods for control of Complex Dynamic systems", IMSOR, n°27, Lingby, 1976.

(10) E. HELLERMAN ; D.C. RARICK : "Reinversion with the Preassigned Pivot Procedure", Mathematical Programming, vol. 1, 1971.

(11) E. HELLERMAN ; D.C. RARICK : "The Partioned Preassigned Pivot Procedure (P⁴)", in Sparce Matrices and Their Applications, D.J. Rose and R.A. Willoughby (eds), Plenum Press, 1972.

(12) H.M. MARKOWITZ : "The Elimination Form of the Inverse and its Application to linear Programming", Management Science, vol.3, 1957.

(13) L.S. LASDON : "Optimization Theory for Large Systems", MacMillan, 1970.

(14) A.K. KEVORKIAN ; J. SNOEK : "Decomposition in Large Scale Systems : Theory and Applications in Solving Large Sets of Non linear Simultaneous Equations", in Decomposition of Large Scale Problems, D.M. Himmelblan (ed;), North Holland, 1973, pp. 491-515.

(15) L.K. CHEUNG ; E.S. KUH : "The Bordered Triangular Matrix and Minimum Essential Sets of a Digraph", IEEE Trans. Circuits and Systems, CAS 21, n°5, 1974.
(16) D.V. STEWARD: "On an approach to techniques for the analysis of the structure of large systems of equations", *SIAM Rev.*, n°4, 1962.

(17) D.V. STEWARD: "Partitionning and Tearing systems of equations", *J. SIAM Numer. Anal.*, vol.2, n°2, 1965.

(18) A.L. DULMAGE; N.S. MENDELSHON: "On the inversion of sparse matrices", *Math. Comput.*, n°16, 1962.

(19) F. HARARY: "A graph Theoretic approach to matrix inversion by partitioning", *Numerische Math.*, n°4, 1962.

(20) A.K. KEVORKIAN; J. SNOEK: "Decomposition in Large scale systems: Theory and Applications of structural analysis in Partitioning, Disjointing and Constructing hierarchical Systems", same as ref. (14), pp. 467-489.

(21) R. TARJAN: "Depth first Search and Linear graph algorithms", *SIAM J. on Computers*, vol.1, 1972.

(22) R.P. TEWARSON: "Computations with Sparse Matrices", *SIAM Rev.*, vol.12, 1972.

(23) J.C. DICKSON: "Finding permutation operations to produce a large triangular submatrix", presented at the *Nat. Meeting Operations Research Soc. Amer.*, Houston, Texas, 1965.

(24) D.V. STEWARD: "Tearing analysis of the structure of disordering sparse matrices in Proc. Sparse Matrix Yorktown Conf.-, R.A. Willoughby (ed.)", 1968, pp. 65-74.

(25) A. NATHAN; R.K. EVEN: "The inversion of sparse matrices by a strategy derived from their graphs", *Comput. J.*, vol. 10, 1967.

(26) A.K. KEVORKIAN: "An optimum Decomposition Method for the Iterative of Nonlinear Algebraic Equations", *J. of Optimization Theory and Applications*, 1976.

(27) W.P. LEDET; D.M. HIMMELBLAU: "Decomposition Procedures for the solving of large-scale systems", *Advan. Chem. Eng.* n°8, 1970.

(28) A.K. KEVORKIAN: "On the bordered Triangular or Lower N-Forms of an Irreducible Matrix", *IEEE Trans. Circuits and Systems*, oct. 1976, Letters to the Editor.

(29) F. HARARY: "Graph Theory", Reading Mass: Addison Wesley, 1969.

(30) F. HARARY; R.Z. NORMAN; D. CARTWRIGHT: "Structural Models", John Wiley, 1965.

(31) G.W. SMITH; R.B. WALFORD: "The Identification of a minimal Feedback Vertex Set of a Directed Graph", *IEEE Trans. Circuits and Systems*, vol. CAS - 22, n°1, 1975.
(32) L. DIVIETTI : "A method for the determination of the paths and the index of a signal flow graph" in Proc. 7th Intern. Conv. Automation and Insturmentation, Milano, 1964.

(33) A. LEMPEL; I. CEDERBAUM : "Minimum feedback arc and vertex sets of a directed graph", IEEE Trans. Circuits Theory, vol.CT;13, 1966.

(34) L. DIVIETTI; A. GRASSELLI : "On the determination of minimum feedback arc and vertex sets", IEEE Trans. Circuit Theory, vol. CT-IS, 1968.

(35) G. GUARDABASSI : "A note on minimal essential sets", IEEE Trans. Circuit Theory, (corresp.) vol. CT-18, 1971.

(36) M. DIAZ; J.P. RICHARD; M. COURVOISIER : "A note on Minimal and Quasi-minimal essential set in complex directed graphs", IEEE Trans. Circuit Theory, (corresp.), vol. CT-19, 1972.

(37) F. HARARY : "A minimal Feedback Vertex Set of a Digraph", IEEE Trans. on Circuits and Systems", oct. 1975, Letters to the Editor.

(38) G. GUARDABASSI : "An indirect Method for Minimal Essential Sets", IEEE Trans. Circuits and Systems, Vol. CAS, 21, n°1, 1974.

(39) S.S. YAU : "Generation of All Hamiltonian Circuits, and center of graph, and related problems", IEEE Trans. Circuit Theory (Corresp.) vol. CT-14 CT-14, 1967.

(40) R.M. KARP : "On the Computational Complexity of Combinatorial Problems", Networks, n°5, 1975.

(41) F. HARARY : "On the minimal Feedback Vertex Sets of a digraph", IEEE Trans. on Circuit and Systems, Letters to the Editor, oct. 1975.